diff --git a/CHANGELOG.md b/CHANGELOG.md index 20217a39..627c1ce7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 site makes the invariant explicit. ### Changed +- `LinearRegression`, `LinRegSlope` and `LinRegAngle` (via composition over + `LinRegSlope`) now run their rolling ordinary-least-squares fit + **incrementally** in O(1) per update (audit finding R2). Previously every + tick refit the line from scratch in O(period). The OLS denominators (`Σx` + and `Σxx`) depend only on `period`, so they were already precomputed; this + release adds running `Σy` and `Σxy` accumulators and slides them in closed + form via the identity + `new_Σxy = old_Σxy − old_Σy + popped_y₀` (then `Σxy += (n − 1) · new_value` + and `Σy += new_value`). New per-bar equivalence tests compare the O(1) + output against a fresh O(n) refit on noisy ramps, step functions, and + constants — values agree to within 1e-9. - Fuzz suite expanded from 2 indicators to the full catalogue (audit finding R9). The existing `indicator_update` target now exercises every scalar-input indicator (~33 classes including MACD and Bollinger Bands); a new diff --git a/crates/wickra-core/src/indicators/linreg.rs b/crates/wickra-core/src/indicators/linreg.rs index 7fe36eea..2795ba30 100644 --- a/crates/wickra-core/src/indicators/linreg.rs +++ b/crates/wickra-core/src/indicators/linreg.rs @@ -19,8 +19,20 @@ use crate::traits::Indicator; /// /// This is TA-Lib's `LINEARREG`: a smoothed price that lags less than an SMA /// because it extrapolates the *local trend* forward to the current bar -/// instead of averaging it away. The `Σx` terms depend only on `period`, so -/// they are computed once; each `update` is O(period). +/// instead of averaging it away. +/// +/// Each `update` is O(1): the `Σx` and `Σxx` terms depend only on `period` and +/// are precomputed once, while `Σy` and `Σxy` are maintained incrementally as +/// the window slides. The closed-form sliding-window identity for +/// `x = 0, 1, …, period − 1` is +/// +/// ```text +/// new_sum_xy = old_sum_xy − old_sum_y + popped_y0 // index shift by −1 +/// new_sum_y = old_sum_y − popped_y0 +/// // then push the new value at index n−1: +/// sum_xy += (n − 1) · new_value +/// sum_y += new_value +/// ``` /// /// # Example /// @@ -38,8 +50,16 @@ use crate::traits::Indicator; pub struct LinearRegression { period: usize, window: VecDeque, + /// Closed form of `Σx` over `x = 0, 1, …, period − 1` — constant in `period`. sum_x: f64, + /// Closed form of `n · Σxx − (Σx)²` — constant in `period`, the OLS + /// denominator. denom: f64, + /// Running sum of the values currently in the window. + sum_y: f64, + /// Running `Σ(x · y)` where `x` is the position of each value within the + /// trailing window (`0` for the oldest, `period − 1` for the newest). + sum_xy: f64, } impl LinearRegression { @@ -63,6 +83,8 @@ impl LinearRegression { window: VecDeque::with_capacity(period), sum_x, denom: n * sum_xx - sum_x * sum_x, + sum_y: 0.0, + sum_xy: 0.0, }) } @@ -70,20 +92,6 @@ impl LinearRegression { pub const fn period(&self) -> usize { self.period } - - /// Ordinary-least-squares `(slope, endpoint)` over the current full window. - fn fit(&self) -> (f64, f64) { - let n = self.period as f64; - let mut sum_y = 0.0; - let mut sum_xy = 0.0; - for (x, &y) in self.window.iter().enumerate() { - sum_y += y; - sum_xy += x as f64 * y; - } - let slope = (n * sum_xy - self.sum_x * sum_y) / self.denom; - let intercept = (sum_y - slope * self.sum_x) / n; - (slope, intercept + slope * (n - 1.0)) - } } impl Indicator for LinearRegression { @@ -92,17 +100,35 @@ impl Indicator for LinearRegression { fn update(&mut self, value: f64) -> Option { if self.window.len() == self.period { - self.window.pop_front(); + // Sliding phase: pop the oldest, then shift every remaining index + // down by 1 in the running `sum_xy`. The identity + // Σ((i − 1) · y_i for i = 1..n−1) = Σ(i · y_i) − Σ(y_i) + y_0 + // gives the closed-form update below. + let y0 = self.window.pop_front().expect("non-empty"); + self.sum_xy = self.sum_xy - self.sum_y + y0; + self.sum_y -= y0; } + // Append at position `k = current length` before the push. During + // warmup `k` ranges over `0..period − 1`; once the window is full it + // is always `period − 1`. + let k = self.window.len() as f64; self.window.push_back(value); + self.sum_y += value; + self.sum_xy += k * value; + if self.window.len() < self.period { return None; } - Some(self.fit().1) + let n = self.period as f64; + let slope = (n * self.sum_xy - self.sum_x * self.sum_y) / self.denom; + let intercept = (self.sum_y - slope * self.sum_x) / n; + Some(intercept + slope * (n - 1.0)) } fn reset(&mut self) { self.window.clear(); + self.sum_y = 0.0; + self.sum_xy = 0.0; } fn warmup_period(&self) -> usize { @@ -195,4 +221,65 @@ mod tests { prices.iter().map(|x| b.update(*x)).collect::>() ); } + + /// Incremental OLS equivalence: the O(1) implementation must agree to + /// `1e-9` with a fresh-from-scratch O(n) refit on every bar, on inputs + /// chosen to stress every code path: a noisy ramp (sliding phase + /// dominates), a step function (the new value differs sharply from the + /// popped one), and constants (the floating-point accumulators must not + /// drift). + #[test] + fn incremental_matches_naive_fit_bar_by_bar() { + fn naive_endpoint(window: &[f64]) -> f64 { + let n = window.len() as f64; + let mut sum_y = 0.0; + let mut sum_xy = 0.0; + let mut sum_x = 0.0; + let mut sum_xx = 0.0; + for (i, &y) in window.iter().enumerate() { + let x = i as f64; + sum_y += y; + sum_xy += x * y; + sum_x += x; + sum_xx += x * x; + } + let denom = n * sum_xx - sum_x * sum_x; + let slope = (n * sum_xy - sum_x * sum_y) / denom; + let intercept = (sum_y - slope * sum_x) / n; + intercept + slope * (n - 1.0) + } + + fn check(prices: &[f64], period: usize) { + let mut lr = LinearRegression::new(period).unwrap(); + for (t, p) in prices.iter().enumerate() { + let streaming = lr.update(*p); + if t + 1 >= period { + let lo = t + 1 - period; + let expected = naive_endpoint(&prices[lo..=t]); + let got = streaming.expect("warmed up"); + assert!( + (got - expected).abs() < 1e-9, + "endpoint diverges at t={t}, period={period}: got={got}, expected={expected}", + ); + } + } + } + + let noisy_ramp: Vec = (0..120) + .map(|i| 100.0 + f64::from(i) * 0.5 + (f64::from(i) * 0.7).sin() * 3.0) + .collect(); + check(&noisy_ramp, 5); + check(&noisy_ramp, 14); + check(&noisy_ramp, 30); + + let mut step = vec![1.0; 30]; + step.extend(std::iter::repeat_n(100.0, 30)); + step.extend(std::iter::repeat_n(0.001, 30)); + check(&step, 5); + check(&step, 14); + + let constant = vec![42.0; 50]; + check(&constant, 8); + check(&constant, 25); + } } diff --git a/crates/wickra-core/src/indicators/linreg_slope.rs b/crates/wickra-core/src/indicators/linreg_slope.rs index 5b0c4343..aab13211 100644 --- a/crates/wickra-core/src/indicators/linreg_slope.rs +++ b/crates/wickra-core/src/indicators/linreg_slope.rs @@ -17,8 +17,12 @@ use crate::traits::Indicator; /// This is TA-Lib's `LINEARREG_SLOPE`: a momentum-like reading of how steeply /// price is trending over the window — positive while it rises, negative /// while it falls, near zero when it is flat — without the band-pass quirks -/// of a difference-based oscillator. The `Σx` terms depend only on `period`, -/// so they are computed once; each `update` is O(period). +/// of a difference-based oscillator. +/// +/// Each `update` is O(1): the same incremental OLS state as +/// [`LinearRegression`](crate::LinearRegression) is maintained — `Σx` and +/// `Σxx` are precomputed once from `period`, while `Σy` and `Σxy` are slid +/// forward in closed form on every push. /// /// # Example /// @@ -36,8 +40,14 @@ use crate::traits::Indicator; pub struct LinRegSlope { period: usize, window: VecDeque, + /// Closed form of `Σx` over `x = 0, 1, …, period − 1` — constant in `period`. sum_x: f64, + /// Closed form of `n · Σxx − (Σx)²` — constant in `period`. denom: f64, + /// Running sum of the values currently in the window. + sum_y: f64, + /// Running `Σ(x · y)` where `x` is the position within the trailing window. + sum_xy: f64, } impl LinRegSlope { @@ -61,6 +71,8 @@ impl LinRegSlope { window: VecDeque::with_capacity(period), sum_x, denom: n * sum_xx - sum_x * sum_x, + sum_y: 0.0, + sum_xy: 0.0, }) } @@ -76,24 +88,30 @@ impl Indicator for LinRegSlope { fn update(&mut self, value: f64) -> Option { if self.window.len() == self.period { - self.window.pop_front(); + // Sliding-window identity: when the window slides one step forward + // the indices `x` for every kept entry shift down by 1, so + // new_sum_xy = old_sum_xy − old_sum_y + y0 + // (`y0` is the popped front value). + let y0 = self.window.pop_front().expect("non-empty"); + self.sum_xy = self.sum_xy - self.sum_y + y0; + self.sum_y -= y0; } + let k = self.window.len() as f64; self.window.push_back(value); + self.sum_y += value; + self.sum_xy += k * value; + if self.window.len() < self.period { return None; } let n = self.period as f64; - let mut sum_y = 0.0; - let mut sum_xy = 0.0; - for (x, &y) in self.window.iter().enumerate() { - sum_y += y; - sum_xy += x as f64 * y; - } - Some((n * sum_xy - self.sum_x * sum_y) / self.denom) + Some((n * self.sum_xy - self.sum_x * self.sum_y) / self.denom) } fn reset(&mut self) { self.window.clear(); + self.sum_y = 0.0; + self.sum_xy = 0.0; } fn warmup_period(&self) -> usize { @@ -192,4 +210,52 @@ mod tests { prices.iter().map(|x| b.update(*x)).collect::>() ); } + + /// Incremental OLS equivalence for the slope: the O(1) implementation must + /// agree bar-by-bar with a fresh-from-scratch O(n) refit, on a noisy ramp + /// (sliding-phase dominated) and a step function (large pop/push deltas). + #[test] + fn incremental_matches_naive_slope_bar_by_bar() { + fn naive_slope(window: &[f64]) -> f64 { + let n = window.len() as f64; + let mut sum_y = 0.0; + let mut sum_xy = 0.0; + let mut sum_x = 0.0; + let mut sum_xx = 0.0; + for (i, &y) in window.iter().enumerate() { + let x = i as f64; + sum_y += y; + sum_xy += x * y; + sum_x += x; + sum_xx += x * x; + } + (n * sum_xy - sum_x * sum_y) / (n * sum_xx - sum_x * sum_x) + } + + fn check(prices: &[f64], period: usize) { + let mut ls = LinRegSlope::new(period).unwrap(); + for (t, p) in prices.iter().enumerate() { + let streaming = ls.update(*p); + if t + 1 >= period { + let lo = t + 1 - period; + let expected = naive_slope(&prices[lo..=t]); + let got = streaming.expect("warmed up"); + assert!( + (got - expected).abs() < 1e-9, + "slope diverges at t={t}, period={period}: got={got}, expected={expected}", + ); + } + } + } + + let noisy_ramp: Vec = (0..120) + .map(|i| 100.0 + f64::from(i) * 0.5 + (f64::from(i) * 0.7).sin() * 3.0) + .collect(); + check(&noisy_ramp, 5); + check(&noisy_ramp, 14); + + let mut step = vec![1.0; 30]; + step.extend(std::iter::repeat_n(100.0, 30)); + check(&step, 7); + } }