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 self.rx.len() < p {
410            return None;
411        }
412        let n = p as f64;
413        let denom = n * self.sxx - self.sx * self.sx;
414        Some(if denom > 0.0 {
415            (n * self.sxy - self.sx * self.sy) / denom
416        } else {
417            0.0
418        })
419    }
420}
421
422impl Next<(f64, f64)> for TaBETA {
423    type Output = f64;
424
425    fn next(&mut self, (v0, v1): (f64, f64)) -> Self::Output {
426        let (Some(p0), Some(p1)) = (self.prev0, self.prev1) else {
427            self.prev0 = Some(v0);
428            self.prev1 = Some(v1);
429            return f64::NAN;
430        };
431        self.prev0 = Some(v0);
432        self.prev1 = Some(v1);
433        if p0 == 0.0 || p1 == 0.0 {
434            return f64::NAN;
435        }
436        let rx = (v0 - p0) / p0;
437        let ry = (v1 - p1) / p1;
438        self.push_return(rx, ry).unwrap_or(f64::NAN)
439    }
440}
441
442#[cfg(test)]
443mod tests {
444    use super::*;
445    use crate::test_utils::{
446        assert_bitwise_batch_streaming_parity, log_price_series, reference_stddev,
447    };
448    use proptest::prelude::*;
449
450    // ---- quantwave-qkft: numerical stability of the stddev accumulator ------
451
452    /// The failing case from quantwave-qkft: on log-price the one-pass
453    /// `E[X²] - E[X]²` form disagreed with a two-pass reference by ~1e-6
454    /// relative, blowing a strict 1e-7 parity gate. TaSTDDEV must now clear it
455    /// by orders of magnitude.
456    #[test]
457    fn ta_stddev_matches_two_pass_reference_on_log_price() {
458        let period = 20;
459        let data = log_price_series(600);
460        let mut ind = TaSTDDEV::new(period, 1.0);
461        let mut worst: f64 = 0.0;
462        for (i, &x) in data.iter().enumerate() {
463            let got = ind.next(x);
464            if got.is_nan() {
465                continue;
466            }
467            let reference = reference_stddev(&data[i + 1 - period..=i]);
468            assert!(reference > 0.0);
469            worst = worst.max(((got - reference) / reference).abs());
470        }
471        assert!(
472            worst < 1e-7,
473            "TaSTDDEV vs two-pass reference on log-price: {worst:e} (gate 1e-7)"
474        );
475        // Not just inside the gate — at reference precision.
476        assert!(worst < 1e-12, "expected ~1e-15 precision, got {worst:e}");
477    }
478
479    #[test]
480    fn ta_var_matches_two_pass_reference_on_log_price() {
481        let period = 20;
482        let data = log_price_series(600);
483        let mut ind = TaVAR::new(period, 1.0);
484        let mut worst: f64 = 0.0;
485        for (i, &x) in data.iter().enumerate() {
486            let got = ind.next(x);
487            if got.is_nan() {
488                continue;
489            }
490            let s = reference_stddev(&data[i + 1 - period..=i]);
491            let reference = s * s;
492            assert!(reference > 0.0);
493            worst = worst.max(((got - reference) / reference).abs());
494        }
495        assert!(worst < 1e-12, "TaVAR vs reference on log-price: {worst:e}");
496    }
497
498    /// Batch and streaming must stay bit-identical: the accumulator's periodic
499    /// exact refresh is driven by bar index only, never by data or by how the
500    /// series was chunked.
501    #[test]
502    fn ta_stddev_batch_streaming_bitwise_parity() {
503        let data = log_price_series(500);
504        assert_bitwise_batch_streaming_parity(&data, || TaSTDDEV::new(20, 2.0));
505        assert_bitwise_batch_streaming_parity(&data, || TaVAR::new(14, 1.0));
506        let prices: Vec<f64> = (0..500)
507            .map(|i| 2980.0 + (i as f64 * 0.31).sin() * 40.0)
508            .collect();
509        assert_bitwise_batch_streaming_parity(&prices, || TaSTDDEV::new(20, 2.0));
510        assert_bitwise_batch_streaming_parity(&prices, || TaVAR::new(14, 1.0));
511    }
512
513    /// The LINEARREG family shares the same rolling window; make sure swapping
514    /// the accumulator underneath did not disturb it.
515    #[test]
516    fn linreg_family_batch_streaming_bitwise_parity() {
517        let data = log_price_series(400);
518        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG::new(14));
519        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG_SLOPE::new(14));
520        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG_INTERCEPT::new(14));
521        assert_bitwise_batch_streaming_parity(&data, || TaLINEARREG_ANGLE::new(14));
522        assert_bitwise_batch_streaming_parity(&data, || TaTSF::new(14));
523    }
524
525    proptest! {
526        #[test]
527        fn test_ta_beta_parity(
528            x in prop::collection::vec(1.0..100.0, 15..100),
529            y in prop::collection::vec(1.0..100.0, 15..100)
530        ) {
531            let len = x.len().min(y.len());
532            let period = 14;
533            let mut ind = TaBETA::new(period);
534            let streaming: Vec<f64> = (0..len).map(|i| ind.next((x[i], y[i]))).collect();
535            let batch = talib_rs::statistic::beta(&x[..len], &y[..len], period)
536                .unwrap_or_else(|_| vec![f64::NAN; len]);
537            for (s, b) in streaming.iter().zip(batch.iter()) {
538                if s.is_nan() { assert!(b.is_nan()); }
539                else if !b.is_nan() { approx::assert_relative_eq!(s, b, epsilon = 1e-5); }
540            }
541        }
542        #[test]
543        fn test_ta_stddev_parity(input in prop::collection::vec(0.1..100.0, 1..100)) {
544            let period = 10;
545            let nbdev = 1.0;
546            let mut s = TaSTDDEV::new(period, nbdev);
547            let streaming: Vec<f64> = input.iter().map(|&x| s.next(x)).collect();
548            let batch = talib_rs::statistic::stddev(&input, period, nbdev)
549                .unwrap_or_else(|_| vec![f64::NAN; input.len()]);
550            for (a, b) in streaming.iter().zip(batch.iter()) {
551                if a.is_nan() { assert!(b.is_nan()); }
552                else { approx::assert_relative_eq!(a, b, epsilon = 1e-6); }
553            }
554        }
555
556        #[test]
557        fn test_ta_linearreg_parity(input in prop::collection::vec(0.1..100.0, 1..100)) {
558            let period = 10;
559            let mut s = TaLINEARREG::new(period);
560            let streaming: Vec<f64> = input.iter().map(|&x| s.next(x)).collect();
561            let batch = talib_rs::statistic::linearreg(&input, period)
562                .unwrap_or_else(|_| vec![f64::NAN; input.len()]);
563            for (a, b) in streaming.iter().zip(batch.iter()) {
564                if a.is_nan() { assert!(b.is_nan()); }
565                else { approx::assert_relative_eq!(a, b, epsilon = 1e-6); }
566            }
567        }
568
569        #[test]
570        fn test_ta_correl_parity(
571            x in prop::collection::vec(0.1..100.0, 1..100),
572            y in prop::collection::vec(0.1..100.0, 1..100)
573        ) {
574            let len = x.len().min(y.len());
575            let period = 10;
576            let mut s = TaCORREL::new(period);
577            let streaming: Vec<f64> = (0..len).map(|i| s.next((x[i], y[i]))).collect();
578            let batch = talib_rs::statistic::correl(&x[..len], &y[..len], period)
579                .unwrap_or_else(|_| vec![f64::NAN; len]);
580            for (a, b) in streaming.iter().zip(batch.iter()) {
581                if a.is_nan() { assert!(b.is_nan()); }
582                else { approx::assert_relative_eq!(a, b, epsilon = 1e-6); }
583            }
584        }
585    }
586}
quantwave_core/indicators/incremental/statistics_ta.rs

quantwave_core/indicators/incremental/
statistics_ta.rs
