diff --git a/crates/ml-alpha/src/multi_horizon_labels.rs b/crates/ml-alpha/src/multi_horizon_labels.rs index 3600208a8..c1801eecf 100644 --- a/crates/ml-alpha/src/multi_horizon_labels.rs +++ b/crates/ml-alpha/src/multi_horizon_labels.rs @@ -165,48 +165,57 @@ pub fn generate_outcome_labels_ab( continue; } - // Pass 1: rolling Welford std of log-returns r[t] = ln(p[t+K] / p[t]). - // Maintain a window of the last SIGMA_WINDOW finite log-returns and - // recompute mean/var on the window. Window is small (1000) and N - // typically <= a few hundred K, so the O(W) per step is acceptable - // and avoids the numerical instability of subtracting popped values - // from a streaming sum-of-squares. + // Pass 1 (SPEED-C, 2026-05-22): rolling sample-std of log-returns + // r[t] = ln(p[t+K] / p[t]). Maintains f64 running sum_x and sum_x2 + // updated incrementally on push/pop — O(1) per step vs the prior + // O(SIGMA_WINDOW) recompute. f64 accumulators contain ~16 decimal + // digits of precision; over a 5M-step file the accumulated rounding + // error stays well below the f32 output precision (`.as f32` + // truncation absorbs it). Every RECOMPUTE_PERIOD pops we still do + // a full window sweep to reset the running totals as defense in + // depth against pathological drift. + const RECOMPUTE_PERIOD: usize = 10_000; let mut window: VecDeque = VecDeque::with_capacity(SIGMA_WINDOW); + let mut sum_x: f64 = 0.0; + let mut sum_x2: f64 = 0.0; + let mut pops_since_reset: usize = 0; for t in 0..n - k { let p_t = prices[t]; let p_kt = prices[t + k]; // Compute σ_K[t] BEFORE pushing r[t] into the window, so σ_K[t] is - // strictly causal (no peek at the future bar's return). Window - // therefore must have warmed up via prior iterations. + // strictly causal (no peek at the future bar's return). if window.len() >= 2 { - // Sample variance via two-pass mean/var on the window (small W, - // accuracy > speed). Welford's online form would be marginal here. - let w_len = window.len() as f32; - let mean: f32 = window.iter().sum::() / w_len; - let var: f32 = window - .iter() - .map(|&x| { - let d = x - mean; - d * d - }) - .sum::() - / (w_len - 1.0); - let std = var.sqrt(); - sigma_k[h][t] = std.max(sigma_floor); + let w_len = window.len() as f64; + let mean = sum_x / w_len; + // Sample variance via centered form: E[x²] - mean² scaled by + // Bessel's correction. Computed in f64; cast to f32 at the end. + let var = (sum_x2 / w_len - mean * mean) * (w_len / (w_len - 1.0)); + // Clamp negatives that can appear from f64 rounding when var≈0. + let var = var.max(0.0); + sigma_k[h][t] = ((var as f32).sqrt()).max(sigma_floor); } else if sigma_floor > 0.0 { - // Even with no window data, expose the structural floor so - // downstream code can normalize meaningfully once cost > 0. - // With cost == 0 (test path) we leave NaN so y_size is masked. sigma_k[h][t] = sigma_floor; } + // Push r[t] into the window for FUTURE steps. if p_t.is_finite() && p_kt.is_finite() && p_t > 0.0 && p_kt > 0.0 { let r = (p_kt / p_t).ln(); if r.is_finite() { if window.len() == SIGMA_WINDOW { - window.pop_front(); + let old = window.pop_front().expect("window non-empty"); + sum_x -= old as f64; + sum_x2 -= (old as f64) * (old as f64); + pops_since_reset += 1; + if pops_since_reset >= RECOMPUTE_PERIOD { + // Periodic full recompute against accumulated rounding. + sum_x = window.iter().map(|&v| v as f64).sum(); + sum_x2 = window.iter().map(|&v| (v as f64) * (v as f64)).sum(); + pops_since_reset = 0; + } } window.push_back(r); + sum_x += r as f64; + sum_x2 += (r as f64) * (r as f64); } } } @@ -571,4 +580,74 @@ mod tests { assert!(labels.y_size_short[0][t].is_nan()); } } + + #[test] + fn online_sigma_matches_naive_full_window_recompute_within_tolerance() { + // Build a synthetic 30k-snapshot price series with stochastic returns. + // Verifies that the online (sum_x, sum_x2) maintenance matches the + // straightforward two-pass mean+var on every step, within f32 tolerance. + use rand::{Rng, SeedableRng}; + let mut rng = rand_chacha::ChaCha8Rng::seed_from_u64(0xDEADBEEF); + let mut prices = vec![5000.0_f32; 30_000]; + for i in 1..prices.len() { + // Geometric Brownian with σ=0.0001 step-noise. + let r: f32 = rng.gen_range(-0.0001..0.0001); + prices[i] = (prices[i - 1] * (1.0 + r)).max(1.0); + } + let horizons = [100_usize, 300, 1000]; + let cost = 0.5_f32; + let out = generate_outcome_labels_ab(&prices, &horizons, cost) + .expect("generate_outcome_labels_ab should succeed"); + + // For each horizon, recompute σ_K[t] independently via the OLD naive + // O(W) algorithm and check the new fast path agrees within f32 tol. + for (h_idx, &k) in horizons.iter().enumerate() { + let mut window: std::collections::VecDeque = + std::collections::VecDeque::with_capacity(SIGMA_WINDOW); + let sigma_floor = cost / 4.0; + for t in 0..prices.len() - k { + let p_t = prices[t]; + let p_kt = prices[t + k]; + let naive_sigma = if window.len() >= 2 { + let w_len = window.len() as f32; + let mean: f32 = window.iter().sum::() / w_len; + let var: f32 = window + .iter() + .map(|&x| { + let d = x - mean; + d * d + }) + .sum::() + / (w_len - 1.0); + var.sqrt().max(sigma_floor) + } else if sigma_floor > 0.0 { + sigma_floor + } else { + f32::NAN + }; + if p_t.is_finite() && p_kt.is_finite() && p_t > 0.0 && p_kt > 0.0 { + let r = (p_kt / p_t).ln(); + if r.is_finite() { + if window.len() == SIGMA_WINDOW { + window.pop_front(); + } + window.push_back(r); + } + } + let fast_sigma = out.sigma_k[h_idx][t]; + if naive_sigma.is_finite() && fast_sigma.is_finite() { + // Relative tolerance 1e-4 — f32 accumulation drift over + // 30k snapshots stays well within this even with naive + // push-pop. The fix maintains f64 accumulators so the + // drift is even smaller in practice. + let rel_err = ((fast_sigma - naive_sigma).abs()) / + naive_sigma.abs().max(1e-12); + assert!( + rel_err < 1e-4, + "h_idx={h_idx} k={k} t={t} fast={fast_sigma} naive={naive_sigma} rel_err={rel_err}", + ); + } + } + } + } }