quantwave_core/indicators/incremental/
statistics_ta.rs

1//! Native O(1) TA-Lib statistics: STDDEV, VAR, LINEARREG family, BETA, CORREL, TSF.
2
3use crate::indicators::incremental::rolling_variance::RollingVariance;
4use crate::traits::Next;
5use crate::utils::RingBuffer;
6
7#[inline]
8fn linreg_coeffs(ys: &[f64]) -> (f64, f64) {
9    let n = ys.len() as f64;
10    let mut sum_x = 0.0;
11    let mut sum_x2 = 0.0;
12    let mut sum_y = 0.0;
13    let mut sum_xy = 0.0;
14    for (i, &y) in ys.iter().enumerate() {
15        let x = i as f64;
16        sum_x += x;
17        sum_x2 += x * x;
18        sum_y += y;
19        sum_xy += x * y;
20    }
21    let denom = n.mul_add(sum_x2, -sum_x * sum_x);
22    if denom == 0.0 {
23        return (f64::NAN, f64::NAN);
24    }
25    let slope = n.mul_add(sum_xy, -sum_x * sum_y) / denom;
26    let intercept = (sum_y - slope * sum_x) / n;
27    (slope, intercept)
28}
29
30/// Fixed-length sliding window with a numerically stable variance accumulator
31/// (see [`RollingVariance`] — quantwave-qkft).
32#[derive(Debug, Clone)]
33struct RollingWindow {
34    acc: RollingVariance,
35}
36
37impl RollingWindow {
38    fn new(period: usize) -> Self {
39        Self {
40            acc: RollingVariance::new(period),
41        }
42    }
43
44    fn push(&mut self, v: f64) -> bool {
45        self.acc.push(v)
46    }
47
48    fn values(&self) -> Vec<f64> {
49        self.acc.values()
50    }
51}
52
53#[derive(Debug, Clone)]
54#[allow(non_camel_case_types)]
55pub struct TaSTDDEV {
56    pub timeperiod: usize,
57    pub nbdev: f64,
58    inner: RollingWindow,
59}
60
61impl TaSTDDEV {
62    pub fn new(timeperiod: usize, nbdev: f64) -> Self {
63        Self {
64            timeperiod,
65            nbdev,
66            inner: RollingWindow::new(timeperiod),
67        }
68    }
69}
70
71impl Next<f64> for TaSTDDEV {
72    type Output = f64;
73
74    fn next(&mut self, input: f64) -> Self::Output {
75        if !self.inner.push(input) {
76            return f64::NAN;
77        }
78        self.inner.acc.stddev() * self.nbdev
79    }
80}
81
82#[derive(Debug, Clone)]
83#[allow(non_camel_case_types)]
84pub struct TaVAR {
85    pub timeperiod: usize,
86    pub nbdev: f64,
87    inner: RollingWindow,
88}
89
90impl TaVAR {
91    pub fn new(timeperiod: usize, nbdev: f64) -> Self {
92        Self {
93            timeperiod,
94            nbdev,
95            inner: RollingWindow::new(timeperiod),
96        }
97    }
98}
99
100impl Next<f64> for TaVAR {
101    type Output = f64;
102
103    fn next(&mut self, input: f64) -> Self::Output {
104        if !self.inner.push(input) {
105            return f64::NAN;
106        }
107        // NOTE: `nbdev` is deliberately unused here — TA-Lib's VAR ignores it
108        // too, and the batch oracle agrees. Behaviour preserved from before the
109        // quantwave-qkft accumulator swap.
110        self.inner.acc.variance()
111    }
112}
113
114#[derive(Debug, Clone)]
115struct LinRegCore {
116    _period: usize,
117    inner: RollingWindow,
118}
119
120impl LinRegCore {
121    fn new(period: usize) -> Self {
122        Self {
123            _period: period,
124            inner: RollingWindow::new(period),
125        }
126    }
127
128    fn push(&mut self, v: f64) -> Option<(f64, f64)> {
129        if !self.inner.push(v) {
130            return None;
131        }
132        Some(linreg_coeffs(&self.inner.values()))
133    }
134}
135
136#[derive(Debug, Clone)]
137#[allow(non_camel_case_types)]
138pub struct TaLINEARREG {
139    pub timeperiod: usize,
140    core: LinRegCore,
141}
142
143impl TaLINEARREG {
144    pub fn new(timeperiod: usize) -> Self {
145        Self {
146            timeperiod,
147            core: LinRegCore::new(timeperiod),
148        }
149    }
150}
151
152impl Next<f64> for TaLINEARREG {
153    type Output = f64;
154
155    fn next(&mut self, input: f64) -> Self::Output {
156        let Some((slope, intercept)) = self.core.push(input) else {
157            return f64::NAN;
158        };
159        let x = (self.timeperiod - 1) as f64;
160        intercept + slope * x
161    }
162}
163
164#[derive(Debug, Clone)]
165#[allow(non_camel_case_types)]
166pub struct TaLINEARREG_SLOPE {
167    pub timeperiod: usize,
168    core: LinRegCore,
169}
170
171impl TaLINEARREG_SLOPE {
172    pub fn new(timeperiod: usize) -> Self {
173        Self {
174            timeperiod,
175            core: LinRegCore::new(timeperiod),
176        }
177    }
178}
179
180impl Next<f64> for TaLINEARREG_SLOPE {
181    type Output = f64;
182
183    fn next(&mut self, input: f64) -> Self::Output {
184        match self.core.push(input) {
185            Some((slope, _)) => slope,
186            None => f64::NAN,
187        }
188    }
189}
190
191#[derive(Debug, Clone)]
192#[allow(non_camel_case_types)]
193pub struct TaLINEARREG_INTERCEPT {
194    pub timeperiod: usize,
195    core: LinRegCore,
196}
197
198impl TaLINEARREG_INTERCEPT {
199    pub fn new(timeperiod: usize) -> Self {
200        Self {
201            timeperiod,
202            core: LinRegCore::new(timeperiod),
203        }
204    }
205}
206
207impl Next<f64> for TaLINEARREG_INTERCEPT {
208    type Output = f64;
209
210    fn next(&mut self, input: f64) -> Self::Output {
211        match self.core.push(input) {
212            Some((_, intercept)) => intercept,
213            None => f64::NAN,
214        }
215    }
216}
217
218#[derive(Debug, Clone)]
219#[allow(non_camel_case_types)]
220pub struct TaLINEARREG_ANGLE {
221    pub timeperiod: usize,
222    core: LinRegCore,
223}
224
225impl TaLINEARREG_ANGLE {
226    pub fn new(timeperiod: usize) -> Self {
227        Self {
228            timeperiod,
229            core: LinRegCore::new(timeperiod),
230        }
231    }
232}
233
234impl Next<f64> for TaLINEARREG_ANGLE {
235    type Output = f64;
236
237    fn next(&mut self, input: f64) -> Self::Output {
238        match self.core.push(input) {
239            Some((slope, _)) => slope.atan().to_degrees(),
240            None => f64::NAN,
241        }
242    }
243}
244
245#[derive(Debug, Clone)]
246#[allow(non_camel_case_types)]
247pub struct TaTSF {
248    pub timeperiod: usize,
249    core: LinRegCore,
250}
251
252impl TaTSF {
253    pub fn new(timeperiod: usize) -> Self {
254        Self {
255            timeperiod,
256            core: LinRegCore::new(timeperiod),
257        }
258    }
259}
260
261impl Next<f64> for TaTSF {
262    type Output = f64;
263
264    fn next(&mut self, input: f64) -> Self::Output {
265        let Some((slope, intercept)) = self.core.push(input) else {
266            return f64::NAN;
267        };
268        intercept + slope * self.timeperiod as f64
269    }
270}
271
272#[derive(Debug, Clone)]
273struct DualWindow {
274    period: usize,
275    xs: RingBuffer<f64>,
276    ys: RingBuffer<f64>,
277}
278
279impl DualWindow {
280    fn new(period: usize) -> Self {
281        Self {
282            period,
283            xs: RingBuffer::with_capacity(period),
284            ys: RingBuffer::with_capacity(period),
285        }
286    }
287
288    fn push(&mut self, x: f64, y: f64) -> bool {
289        if self.xs.len() >= self.period {
290            let _ = self.xs.pop_front();
291            let _ = self.ys.pop_front();
292        }
293        self.xs.push_back(x);
294        self.ys.push_back(y);
295        self.xs.len() >= self.period
296    }
297
298    fn xs_values(&self) -> Vec<f64> {
299        self.xs.iter().cloned().collect()
300    }
301
302    fn ys_values(&self) -> Vec<f64> {
303        self.ys.iter().cloned().collect()
304    }
305}
306
307#[inline]
308fn correl_beta(xs: &[f64], ys: &[f64]) -> (f64, f64) {
309    let n = xs.len() as f64;
310    let mut sum_x = 0.0;
311    let mut sum_y = 0.0;
312    let mut sum_x2 = 0.0;
313    let mut sum_y2 = 0.0;
314    let mut sum_xy = 0.0;
315    for i in 0..xs.len() {
316        let x = xs[i];
317        let y = ys[i];
318        sum_x += x;
319        sum_y += y;
320        sum_x2 += x * x;
321        sum_y2 += y * y;
322        sum_xy += x * y;
323    }
324    let cov = n.mul_add(sum_xy, -sum_x * sum_y);
325    let var_x = n.mul_add(sum_x2, -sum_x * sum_x);
326    let var_y = n.mul_add(sum_y2, -sum_y * sum_y);
327    if var_x == 0.0 || var_y == 0.0 {
328        return (f64::NAN, f64::NAN);
329    }
330    let correl = cov / (var_x.sqrt() * var_y.sqrt());
331    // Linear regression Y = alpha + beta * X  =>  beta = cov(X,Y) / var(X).
332    let beta = cov / var_x;
333    (correl, beta)
334}
335
336#[derive(Debug, Clone)]
337#[allow(non_camel_case_types)]
338pub struct TaCORREL {
339    pub timeperiod: usize,
340    inner: DualWindow,
341}
342
343impl TaCORREL {
344    pub fn new(timeperiod: usize) -> Self {
345        Self {
346            timeperiod,
347            inner: DualWindow::new(timeperiod),
348        }
349    }
350}
351
352impl Next<(f64, f64)> for TaCORREL {
353    type Output = f64;
354
355    fn next(&mut self, (x, y): (f64, f64)) -> Self::Output {
356        if !self.inner.push(x, y) {
357            return f64::NAN;
358        }
359        correl_beta(&self.inner.xs_values(), &self.inner.ys_values()).0
360    }
361}
362
363/// Beta on percentage returns of two series (TA-Lib `statistic::beta`).
364#[derive(Debug, Clone)]
365#[allow(non_camel_case_types)]
366pub struct TaBETA {
367    pub timeperiod: usize,
368    prev0: Option<f64>,
369    prev1: Option<f64>,
370    rx: RingBuffer<f64>,
371    ry: RingBuffer<f64>,
372    sx: f64,
373    sy: f64,
374    sxx: f64,
375    sxy: f64,
376}
377
378impl TaBETA {
379    pub fn new(timeperiod: usize) -> Self {
380        Self {
381            timeperiod,
382            prev0: None,
383            prev1: None,
384            rx: RingBuffer::with_capacity(timeperiod),
385            ry: RingBuffer::with_capacity(timeperiod),
386            sx: 0.0,
387            sy: 0.0,
388            sxx: 0.0,
389            sxy: 0.0,
390        }
391    }
392
393    fn push_return(&mut self, x: f64, y: f64) -> Option<f64> {
394        let p = self.timeperiod;
395        if self.rx.len() >= p
396            && let (Some(ox), Some(oy)) = (self.rx.pop_front(), self.ry.pop_front())
397        {
398            self.sx -= ox;
399            self.sy -= oy;
400            self.sxx -= ox * ox;
401            self.sxy -= ox * oy;
402        }
403        self.rx.push_back(x);
404        self.ry.push_back(y);
405        self.sx += x;
406        self.sy += y;
407        self.sxx += x * x;
408        self.sxy += x * y;
409        // If accumulators are contaminated by NaN, attempt self-heal from the ring buffer
410        // so that once the NaN ages out of the rolling window, beta recovers.
411        if self.sx.is_nan() || self.sy.is_nan() || self.sxx.is_nan() || self.sxy.is_nan() {
412            self.sx = self.rx.iter().copied().sum();
413            self.sy = self.ry.iter().copied().sum();
414            self.sxx = self.rx.iter().map(|&v| v * v).sum();
415            self.sxy = self
416                .rx
417                .iter()
418                .zip(self.ry.iter())
419                .map(|(&a, &b)| a * b)
420                .sum();
421        }
422        if self.rx.len() < p {
423            return None;
424        }
425        let n = p as f64;
426        let denom = n * self.sxx - self.sx * self.sx;
427        Some(if denom > 0.0 {
428            (n * self.sxy - self.sx * self.sy) / denom
429        } else {
430            0.0
431        })
432    }
433}
434
435impl Next<(f64, f64)> for TaBETA {
436    type Output = f64;
437
438    fn next(&mut self, (v0, v1): (f64, f64)) -> Self::Output {
439        let (Some(p0), Some(p1)) = (self.prev0, self.prev1) else {
440            self.prev0 = Some(v0);
441            self.prev1 = Some(v1);
442            return f64::NAN;
443        };
444        self.prev0 = Some(v0);
445        self.prev1 = Some(v1);
446        if p0 == 0.0 || p1 == 0.0 || v0.is_nan() || v1.is_nan() || p0.is_nan() || p1.is_nan() {
447            self.push_return(f64::NAN, f64::NAN);
448            return f64::NAN;
449        }
450        let rx = (v0 - p0) / p0;
451        let ry = (v1 - p1) / p1;
452        self.push_return(rx, ry).unwrap_or(f64::NAN)
453    }
454}
455
456#[cfg(test)]
457mod tests {
458    use super::*;
459    use crate::test_utils::{
460        assert_bitwise_batch_streaming_parity, log_price_series, reference_stddev,
461    };
462    use proptest::prelude::*;
463
464    // ---- quantwave-qkft: numerical stability of the stddev accumulator ------
465
466    /// The failing case from quantwave-qkft: on log-price the one-pass
467    /// `E[X²] - E[X]²` form disagreed with a two-pass reference by ~1e-6
468    /// relative, blowing a strict 1e-7 parity gate. TaSTDDEV must now clear it
469    /// by orders of magnitude.
470    #[test]
471    fn ta_stddev_matches_two_pass_reference_on_log_price() {
472        let period = 20;
473        let data = log_price_series(600);
474        let mut ind = TaSTDDEV::new(period, 1.0);
475        let mut worst: f64 = 0.0;
476        for (i, &x) in data.iter().enumerate() {
477            let got = ind.next(x);
478            if got.is_nan() {
479                continue;
480            }
481            let reference = reference_stddev(&data[i + 1 - period..=i]);
482            assert!(reference > 0.0);
483            worst = worst.max(((got - reference) / reference).abs());
484        }
485        assert!(
486            worst < 1e-7,
487            "TaSTDDEV vs two-pass reference on log-price: {worst:e} (gate 1e-7)"
488        );
489        // Not just inside the gate — at reference precision.
490        assert!(worst < 1e-12, "expected ~1e-15 precision, got {worst:e}");
491    }
492
493    #[test]
494    fn ta_var_matches_two_pass_reference_on_log_price() {
495        let period = 20;
496        let data = log_price_series(600);
497        let mut ind = TaVAR::new(period, 1.0);
498        let mut worst: f64 = 0.0;
499        for (i, &x) in data.iter().enumerate() {
500            let got = ind.next(x);
501            if got.is_nan() {
502                continue;
503            }
504            let s = reference_stddev(&data[i + 1 - period..=i]);
505            let reference = s * s;
506            assert!(reference > 0.0);
507            worst = worst.max(((got - reference) / reference).abs());
508        }
509        assert!(worst < 1e-12, "TaVAR vs reference on log-price: {worst:e}");
510    }
511
512    /// Batch and streaming must stay bit-identical: the accumulator's periodic
513    /// exact refresh is driven by bar index only, never by data or by how the
514    /// series was chunked.
515    #[test]
516    fn ta_stddev_batch_streaming_bitwise_parity() {
517        let data = log_price_series(500);
518        assert_bitwise_batch_streaming_parity(&data, || TaSTDDEV::new(20, 2.0));
519        assert_bitwise_batch_streaming_parity(&data, || TaVAR::new(14, 1.0));
520        let prices: Vec<f64> = (0..500)
521            .map(|i| 2980.0 + (i as f64 * 0.31).sin() * 40.0)
522            .collect();
523        assert_bitwise_batch_streaming_parity(&prices, || TaSTDDEV::new(20, 2.0));
524        assert_bitwise_batch_streaming_parity(&prices, || TaVAR::new(14, 1.0));
525    }
526
527    /// The LINEARREG family shares the same rolling window; make sure swapping
528    /// the accumulator underneath did not disturb it.
529    #[test]
530    fn linreg_family_batch_streaming_bitwise_parity() {
531        let data = log_price_series(400);
532        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG::new(14));
533        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG_SLOPE::new(14));
534        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG_INTERCEPT::new(14));
535        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG_ANGLE::new(14));
536        assert_bitwise_batch_streaming_parity(&data, || TaTSF::new(14));
537    }
538
539    #[test]
540    fn test_ta_beta_recovers_after_nan() {
541        let period = 5;
542        let mut beta = TaBETA::new(period);
543        // Feed valid series
544        for i in 1..=10 {
545            beta.next((i as f64 * 10.0, i as f64 * 15.0));
546        }
547        // Poison with a single NaN
548        let nan_out = beta.next((110.0, f64::NAN));
549        assert!(nan_out.is_nan());
550
551        // For the next `period` bars, output should be NaN
552        for i in 12..12 + period {
553            beta.next((i as f64 * 10.0, i as f64 * 15.0));
554        }
555
556        // Once the NaN has aged out of the rolling window, it must recover to finite
557        let recovered = beta.next((200.0, 300.0));
558        assert!(
559            !recovered.is_nan(),
560            "Beta must recover once NaN ages out of rolling window, got {recovered}"
561        );
562        approx::assert_relative_eq!(recovered, 1.0, epsilon = 1e-4);
563    }
564
565    proptest! {
566        #[test]
567        fn test_ta_beta_parity(
568            x in prop::collection::vec(1.0..100.0, 15..100),
569            y in prop::collection::vec(1.0..100.0, 15..100)
570        ) {
571            let len = x.len().min(y.len());
572            let period = 14;
573            let mut ind = TaBETA::new(period);
574            let streaming: Vec<f64> = (0..len).map(|i| ind.next((x[i], y[i]))).collect();
575            let batch = talib_rs::statistic::beta(&x[..len], &y[..len], period)
576                .unwrap_or_else(|_| vec![f64::NAN; len]);
577            for (s, b) in streaming.iter().zip(batch.iter()) {
578                if s.is_nan() { assert!(b.is_nan()); }
579                else if !b.is_nan() { approx::assert_relative_eq!(s, b, epsilon = 1e-5); }
580            }
581        }
582        #[test]
583        fn test_ta_stddev_parity(input in prop::collection::vec(0.1..100.0, 1..100)) {
584            let period = 10;
585            let nbdev = 1.0;
586            let mut s = TaSTDDEV::new(period, nbdev);
587            let streaming: Vec<f64> = input.iter().map(|&x| s.next(x)).collect();
588            let batch = talib_rs::statistic::stddev(&input, period, nbdev)
589                .unwrap_or_else(|_| vec![f64::NAN; input.len()]);
590            for (a, b) in streaming.iter().zip(batch.iter()) {
591                if a.is_nan() { assert!(b.is_nan()); }
592                else { approx::assert_relative_eq!(a, b, epsilon = 1e-6); }
593            }
594        }
595
596        #[test]
597        fn test_ta_linearreg_parity(input in prop::collection::vec(0.1..100.0, 1..100)) {
598            let period = 10;
599            let mut s = TaLINEARREG::new(period);
600            let streaming: Vec<f64> = input.iter().map(|&x| s.next(x)).collect();
601            let batch = talib_rs::statistic::linearreg(&input, period)
602                .unwrap_or_else(|_| vec![f64::NAN; input.len()]);
603            for (a, b) in streaming.iter().zip(batch.iter()) {
604                if a.is_nan() { assert!(b.is_nan()); }
605                else { approx::assert_relative_eq!(a, b, epsilon = 1e-6); }
606            }
607        }
608
609        #[test]
610        fn test_ta_correl_parity(
611            x in prop::collection::vec(0.1..100.0, 1..100),
612            y in prop::collection::vec(0.1..100.0, 1..100)
613        ) {
614            let len = x.len().min(y.len());
615            let period = 10;
616            let mut s = TaCORREL::new(period);
617            let streaming: Vec<f64> = (0..len).map(|i| s.next((x[i], y[i]))).collect();
618            let batch = talib_rs::statistic::correl(&x[..len], &y[..len], period)
619                .unwrap_or_else(|_| vec![f64::NAN; len]);
620            for (a, b) in streaming.iter().zip(batch.iter()) {
621                if a.is_nan() { assert!(b.is_nan()); }
622                else { approx::assert_relative_eq!(a, b, epsilon = 1e-6); }
623            }
624        }
625    }
626}
quantwave_core/indicators/incremental/statistics_ta.rs

quantwave_core/indicators/incremental/
statistics_ta.rs
