140 lines of μ, 119 ms, 98,000 rows

The problem: detecting hidden market states

Stock correlations are not constant. In calm markets, large-cap tech equities move together strongly — a long position in AAPL and a short in MSFT earns nothing. During factor rotations they decouple — pairs trades become profitable. In crisis periods correlations break down or invert unpredictably, making naïve hedges dangerous. A strategy that ignores regime will make money in one state and give it all back in another.

The regime is not directly observable. You see prices. From prices you compute correlations, and from correlations you infer which of a small number of hidden states the market is in. That inference problem is a Hidden Markov Model.

What does each state look like? At the correlation-matrix level the difference is structural:

The synthetic dataset contains 16 regime blocks in sequence — alternating between HIGH_CORR, LOW_CORR, and BREAKDOWN states — with known ground-truth durations. This makes correctness easy to verify: the recovered state means should match the design values exactly.

Three scalar features compress the correlation matrix into something an EM classifier can act on:

• mean_abs_corr — mean absolute off-diagonal correlation: 0.75 / 0.15 / 0.00
• dispersion — standard deviation of correlations: 0.08 / 0.25 / 0.35
• eigenratio — ratio of first to second eigenvalue: 5.0 / 1.5 / 0.8

HIGH_CORR is easy to detect: one dominant factor drives everything. BREAKDOWN is also clear: correlations collapse toward zero. LOW_CORR is the hard case — correlations exist but are weak and noisy, sitting between the other two states with no sharp boundary.

Simulated mean_abs_corr signal across 252 trading days — the feature the HMM reads to detect regime. Shading shows ground-truth regime blocks.

The μ code: two implementations

The experiment has two μ files. The first (correlation_hmm.mu, 167 lines) is analytic: the regime schedule is hardcoded and the EM means are computed algebraically from ground truth. It runs in ~7 ms and confirms the language evaluates conditional lambdas and global mutation correctly.

The second (correlation_hmm_v2.mu, ~140 lines) learns from real data. It loads the full 98k-row CSV, computes log returns, extracts rolling window features, classifies states via EM, and runs a trading simulation — all in plain μ.

The pipeline in v2, top to bottom:

# Phase 1: Load and compute log returns
csv_data: (csv "experiments/correlation_hmm/stock_prices.csv")
columns:  csv_data@1

ratios:  {[p] p % 1,-1 _ p}
logrets: {[p] 1 _ log (ratios p)}

rAAPL: (logrets columns@1)
rAMZN: (logrets columns@2)
# ... (8 more stocks)

# Phase 2: Rolling window features
var_at: {[t] (variance (WINDOW_MINUTES # ((t-WINDOW_MINUTES) _ rAAPL)))}
vars_sample: (var_at' sample_ts)

# Phase 3: EM — classify each sample point
classify:     {[v] (v < threshold); 0; 2}
states_sample: (classify' vars_sample)

# Phase 4: Trading simulation
n_high: (+/ (states_sample = 0))
n_break: (+/ (states_sample = 2))

No libraries. No NumPy. No special matrix verb. Just csv, log, take/drop, each', and arithmetic — the core of the language applied to a production-scale problem.

Analytic EM — verifying correctness

The v1 analytic implementation is worth examining on its own. Given the known regime schedule, it computes weighted EM state means in a dozen lines:

# Indicator arrays: 1.0 where each block is in that state
is0: {[b] (get[RSTATES; b]) = 0.0}' !16

# Weighted minutes per state
wm0: {[b] get[block_mins; b] * get[is0; b]}' !16
.cnt0: +/ wm0

# Weighted feature sums
f1: {[b] corr_of[get[block_state; b]]}' !16
.sum0_f1: +/ {[b] get[wm0; b] * get[f1; b]}' !16

# EM means
mean0_f1: cif[(.cnt0 > 0.0); (.sum0_f1 % .cnt0); 0.0]

The output matches design values exactly — 0.75, 0.15, 0.0 for mean_abs_corr; 0.08, 0.25, 0.35 for dispersion; 5.0, 1.5, 0.8 for eigenratio. This confirms μ evaluates conditional branches, global accumulation inside each', and weighted aggregation correctly.

Speed: 138× faster than Python

Three implementations, same dataset, different algorithmic scope. The table shows what each one actually does and how long it takes:

The 138× number is real but misleading in isolation — it reflects different algorithms, not just different runtimes. The number that matters is 119 ms: that is how long a working regime classifier takes in μ on a full year of tick data. Here is what drives the gaps:

• Python loops over 98,280 timesteps in interpreted Python. Every rolling-window computation is a Python for loop calling NumPy on small slices (a 1,950-minute window). NumPy's per-call overhead dominates when the arrays are small.
• C processes every timestep with a full correlation matrix. O(window × n²) per timestep = ~18M inner-loop iterations × 16s. Thorough but expensive.
• μ samples 18 timesteps and uses variance as a proxy. It does not compute the full correlation matrix. It computes scalar variance via {variance (WINDOW_MINUTES # ((t-W) _ rAAPL))} at 18 spread-out points. csv reads directly into packed float[] columns (no char-list intermediate), and all arithmetic runs on flat contiguous arrays with a single loop per operation.

The structural reason μ is fast here is the same as in the Monte Carlo note: packed column storage eliminates per-element boxing, and the language's functional style means the hot path is a handful of tight C loops with no Python interpreter overhead.

Wall-clock time, log scale. μ v2 processes the same 98k rows in 119 ms vs Python's 136,000 ms.

Simplicity: code size

Fewer lines is not the goal — correctness and clarity are. But code size is a reasonable proxy for cognitive load, and the comparison is striking:

The C implementation is 400 lines because it must manage memory manually: allocate rolling window buffers, free them, track array bounds, accumulate into pre-allocated structs. The Python version is 280 lines because NumPy requires explicit shape management and the Viterbi decoder is non-trivial to vectorise.

The μ v2 implementation has none of that. take # (drop _ data) is the entire rolling window. log (ratios p) is log returns. classify' vars_sample is the EM pass. No allocations to manage, no shapes to track.

Lines of code. μ v2 does more per line — no memory management, no shape tracking.

What was not done

Variance proxy instead of full correlation matrix. The v2 implementation computes variance of AAPL log returns as a scalar feature. The Python and C versions compute the full 10×10 correlation matrix at every timestep (O(window × n²)). Variance is a coarser signal — it detects regime changes but cannot distinguish HIGH_CORR from LOW_CORR as precisely.

Sampled timesteps instead of every timestep. v2 evaluates features at 18 spread-out sample points across the 98k dataset. Python and C evaluate at every one of the ~96k valid timesteps. This is the primary reason for the speed difference — 18 window operations vs 96k.

2-state threshold instead of 3-state Gaussian mixture EM. v2 classifies each sample as HIGH_CORR or BREAKDOWN based on a fixed variance threshold. The full algorithm uses a 3-state Gaussian mixture model with iterative EM and Viterbi decoding — which the Python and C versions implement.

These are practical trade-offs, not oversights. They demonstrate that μ can express and run the algorithmic structure correctly. The path to a full HMM is clear; the language is not the bottleneck.

Variance proxy (y) vs true mean_abs_corr (x) at the 18 sampled timesteps. HIGH_CORR and BREAKDOWN cluster separately — the proxy works, but a threshold classifier cannot recover LOW_CORR precisely.

What comes next

The remaining steps to a full HMM implementation in μ, in rough order of difficulty:

1. Full correlation matrix per window. Currently using scalar variance. The full 10×10 matrix requires O(window × n²) per timestep. μ has eigen built in; the question is whether rolling matrix computation fits the language ergonomically without per-element boxing killing performance.
2. 3-state Gaussian mixture EM. Replace the threshold classifier with a proper E-step (soft assignments) and M-step (mean/variance updates). Iterative, but convergence is fast (<30 iterations for this problem).
3. Viterbi decoding with learned transition matrix. Dynamic programming over the full timestep sequence. Requires mutable state per path — currently the most structurally challenging piece.
4. Per-timestamp trading simulation with mark-to-market. The v1 simulation uses block-level returns. A full simulation requires per-minute P&L tracking tied to regime transitions.

The current 119 ms baseline leaves substantial headroom. If the full correlation matrix computation is O(18 samples × 1950 window × 100 stock pairs), that is ~3.5M multiplications — well within what μ's packed arithmetic can handle without per-element overhead.