ta_statistics/
rolling_moments.rs

1use num_traits::Float;
2
3use crate::{Kbn, RingBuffer};
4
5/// This module provides functionality for calculating rolling statistical moments over a time series.
6///
7/// Rolling moments are essential statistics that help analyze the characteristics of data over
8/// a moving window. These include measures like mean, variance, skewness, and kurtosis, which
9/// provide insights into the distribution and behavior of financial time series data.
10///
11/// The implementation uses Kahan-Babuska-Neumaier summation algorithm for numerical stability
12/// when computing these statistics over potentially large datasets with floating-point values.
13#[derive(Debug, Clone)]
14pub struct RollingMoments<T> {
15    /// Statistics period
16    period: usize,
17    /// Ring buffer to maintain the window
18    buf: RingBuffer<T>,
19    /// Most recent value pushed into the rolling window.
20    value: Option<T>,
21    /// Most recent value popped out of the rolling window (if full).
22    popped: Option<T>,
23    /// Delta Degrees of Freedom
24    ddof: bool,
25    /// Sum of inputs
26    sum: Kbn<T>,
27    /// Sum of squares
28    sum_sq: Kbn<T>,
29    /// Sum of cubes
30    sum_cube: Kbn<T>,
31    /// Sum of fourth powers
32    sum_quad: Kbn<T>,
33    /// Current mean
34    mean: T,
35    /// Second central moment
36    m2: T,
37    /// Third central moment
38    m3: T,
39    /// Fourth central moment
40    m4: T,
41}
42
43impl<T: Float + Default> RollingMoments<T> {
44    /// Creates a new `RollingMoments` instance with the specified period.
45    ///
46    /// # Arguments
47    ///
48    /// * `period` - The period of the statistics
49    ///
50    /// # Returns
51    ///
52    /// * `Self` - The statistics object
53    pub fn new(period: usize) -> Self {
54        Self {
55            period,
56            buf: RingBuffer::new(period),
57            value: None,
58            popped: None,
59            ddof: false,
60            sum: Kbn::default(),
61            sum_sq: Kbn::default(),
62            sum_cube: Kbn::default(),
63            sum_quad: Kbn::default(),
64            mean: T::zero(),
65            m2: T::zero(),
66            m3: T::zero(),
67            m4: T::zero(),
68        }
69    }
70
71    /// Resets the sums
72    #[inline]
73    fn reset_sums(&mut self) {
74        self.sum = Kbn::default();
75        self.sum_sq = Kbn::default();
76        self.sum_cube = Kbn::default();
77        self.sum_quad = Kbn::default();
78    }
79
80    /// Resets the moments
81    #[inline]
82    fn reset_moments(&mut self) {
83        self.mean = T::zero();
84        self.m2 = T::zero();
85        self.m3 = T::zero();
86        self.m4 = T::zero();
87    }
88
89    /// Updates the central moments
90    ///
91    /// # Returns
92    ///
93    /// * `Option<()>` - `None` if the window is not full, `Some(())` otherwise
94    fn update_central_moments(&mut self) -> Option<()> {
95        let n = T::from(self.buf.len())?;
96        if n == T::zero() {
97            self.reset_moments();
98            return None;
99        }
100
101        self.mean = self.sum.total() / n;
102
103        let m1 = self.mean;
104        let m2_raw = self.sum_sq.total() / n;
105        let m3_raw = self.sum_cube.total() / n;
106        let m4_raw = self.sum_quad.total() / n;
107
108        let m1_sq = m1 * m1;
109        let m1_cb = m1_sq * m1;
110        let m1_4 = m1_cb * m1;
111
112        let _2 = T::from(2.0)?;
113        let _3 = T::from(3.0)?;
114        let _4 = T::from(4.0)?;
115        let _6 = T::from(6.0)?;
116        self.m2 = m2_raw - m1_sq;
117        self.m3 = m3_raw - _3 * m1 * m2_raw + _2 * m1_cb;
118        self.m4 = m4_raw - _4 * m1 * m3_raw + _6 * m1_sq * m2_raw - _3 * m1_4;
119        Some(())
120    }
121
122    /// Returns the Delta Degrees of Freedom
123    ///
124    /// # Returns
125    ///
126    /// * `bool` - The Delta Degrees of Freedom
127    #[inline]
128    pub const fn ddof(&self) -> bool {
129        self.ddof
130    }
131
132    /// Sets the Delta Degrees of Freedom
133    ///
134    /// # Arguments
135    ///
136    /// * `ddof` - The Delta Degrees of Freedom
137    ///
138    /// # Returns
139    ///
140    /// * `&mut Self` - The statistics object
141    #[inline]
142    pub const fn set_ddof(&mut self, ddof: bool) -> &mut Self {
143        self.ddof = ddof;
144        self
145    }
146
147    /// Resets the rolling moments
148    ///
149    /// # Returns
150    ///
151    /// * `&mut Self` - The rolling moments object
152    #[inline]
153    pub fn reset(&mut self) -> &mut Self {
154        self.buf.reset();
155        self.value = None;
156        self.popped = None;
157        self.reset_sums();
158        self.reset_moments();
159        self
160    }
161
162    /// Updates the value to the rolling moments
163    ///
164    /// # Arguments
165    ///
166    /// * `value` - The value to update the rolling moments with
167    ///
168    /// # Returns
169    ///
170    /// * `&mut Self` - The rolling moments object
171    ///
172    #[inline]
173    pub fn next(&mut self, value: T) -> &mut Self {
174        self.value = Some(value);
175        self.popped = self.buf.push(value);
176        if let Some(popped) = self.popped {
177            self.sum -= popped;
178            self.sum_sq -= popped * popped;
179            self.sum_cube -= popped * popped * popped;
180            self.sum_quad -= popped * popped * popped * popped;
181        }
182
183        self.sum += value;
184        self.sum_sq += value * value;
185        self.sum_cube += value * value * value;
186        self.sum_quad += value * value * value * value;
187
188        self.update_central_moments();
189
190        self
191    }
192
193    /// Recomputes the rolling statistics, could be called to avoid
194    /// prolonged compounding of floating rounding errors
195    ///
196    /// # Returns
197    ///
198    /// * `&mut Self` - The rolling moments object
199    #[inline]
200    pub fn recompute(&mut self) {
201        self.reset_sums();
202
203        for &v in self.buf.iter() {
204            self.sum += v;
205            self.sum_sq += v * v;
206            self.sum_cube += v * v * v;
207            self.sum_quad += v * v * v * v;
208        }
209
210        self.update_central_moments();
211    }
212
213    /// Returns the value that was removed from the window
214    ///
215    /// # Returns
216    ///
217    /// * `Option<T>` - The value that was removed from the window
218    pub const fn popped(&self) -> Option<T> {
219        self.popped
220    }
221
222    /// Returns the value that was added to the window
223    ///
224    /// # Returns
225    ///
226    /// * `Option<T>` - The value that was added to the window
227    pub const fn value(&self) -> Option<T> {
228        self.value
229    }
230
231    /// Returns the maximum value in the ring buffer
232    #[inline]
233    pub fn max(&self) -> Option<T> {
234        self.buf.max()
235    }
236
237    /// Returns the minimum value in the ring buffer    
238    #[inline]
239    pub fn min(&self) -> Option<T> {
240        self.buf.min()
241    }
242
243    /// Returns the window period
244    ///
245    /// # Returns
246    ///
247    /// * `usize` - The window period
248    #[inline]
249    pub const fn period(&self) -> usize {
250        self.period
251    }
252
253    /// Returns true of the calculation was ready
254    ///
255    /// # Returns
256    ///
257    /// * `bool` - True if the calculation was ready
258    #[inline]
259    pub const fn is_ready(&self) -> bool {
260        self.buf.is_full()
261    }
262
263    /// Returns the number of elements in the buffer
264    ///
265    /// # Returns
266    ///
267    /// * `usize` - The number of elements in the buffer
268    #[inline]
269    pub const fn count(&self) -> usize {
270        self.buf.len()
271    }
272
273    /// Returns an iterator over the elements in the ring buffer
274    #[inline]
275    pub fn iter(&self) -> impl Iterator<Item = &T> {
276        self.buf.iter()
277    }
278
279    /// Returns a slice of the elements in the buffer
280    #[inline]
281    pub fn as_slice(&self) -> &[T] {
282        self.buf.as_slice()
283    }
284
285    /// Returns the sum of all values in the rolling window
286    ///
287    /// # Returns
288    ///
289    /// * `Option<T>` - The sum of all values if the window is ready, None otherwise
290    #[inline]
291    pub fn sum(&self) -> Option<T> {
292        self.is_ready().then_some(self.sum.total())
293    }
294
295    /// Returns the sum of squares of all values in the rolling window
296    ///
297    /// # Returns
298    ///
299    /// * `Option<T>` - The sum of squares of all values if the window is ready, None otherwise
300    #[inline]
301    pub fn sum_sq(&self) -> Option<T> {
302        self.is_ready().then_some(self.sum_sq.total())
303    }
304
305    /// Returns the mean of all values in the rolling window
306    ///
307    /// # Returns
308    ///
309    /// * `Option<T>` - The mean of all values if the window is ready, None otherwise
310    #[inline]
311    pub fn mean(&self) -> Option<T> {
312        self.is_ready().then_some(self.mean)
313    }
314
315    /// Returns the mean of squared values in the rolling window
316    ///
317    /// # Returns
318    ///
319    /// * `Option<T>` - The mean of squared values if the window is ready, None otherwise
320    #[inline]
321    pub fn mean_sq(&self) -> Option<T> {
322        self.sum_sq()
323            .zip(T::from(self.count()))
324            .map(|(ss, n)| ss / n)
325    }
326
327    /// Returns the variance of values in the rolling window
328    ///
329    /// Variance quantifies the dispersion of data points around the mean,
330    /// providing essential volatility measurement for financial time series:
331    ///
332    /// - Serves as the foundation for volatility-based position sizing
333    /// - Enables normalization of returns for cross-asset comparison
334    /// - Provides critical input for statistical arbitrage models
335    /// - Forms the basis for numerous technical indicators like Bollinger Bands
336    ///
337    /// # Returns
338    ///
339    /// * `Option<T>` - The variance, or `None` if the window is not full
340    #[inline]
341    pub fn variance(&self) -> Option<T> {
342        if !self.is_ready() {
343            return None;
344        }
345        let n = T::from(self.count())?;
346        let denom = if self.ddof { n - T::one() } else { n };
347        if denom > T::zero() {
348            Some(self.m2 * n / denom)
349        } else {
350            None
351        }
352    }
353
354    /// Returns the standard deviation of values in the rolling window
355    ///
356    /// Standard deviation quantifies the dispersion of data points around the mean,
357    /// providing essential volatility measurement for financial time series:
358    ///
359    /// - Serves as the foundation for volatility-based position sizing
360    /// - Enables normalization of returns for cross-asset comparison
361    /// - Provides critical input for statistical arbitrage models
362    /// - Forms the basis for numerous technical indicators like Bollinger Bands
363    ///
364    /// # Returns
365    ///
366    /// * `Option<T>` - The standard deviation, or `None` if the window is not full
367    #[inline]
368    pub fn stddev(&self) -> Option<T> {
369        self.variance().and_then(|var| {
370            if var >= T::zero() {
371                Some(var.sqrt())
372            } else {
373                None
374            }
375        })
376    }
377
378    /// Returns the Zscore of the most recent value
379    ///
380    /// Zscore measures the number of standard deviations a value is from the mean.
381    ///
382    /// # Returns
383    ///
384    /// * `Option<T>` - The Zscore if the window is ready and standard deviation is positive, None otherwise
385    #[inline]
386    pub fn zscore(&self) -> Option<T> {
387        let value = self.value?;
388        let mean = self.mean()?;
389        let stddev = self.stddev()?;
390
391        if stddev > T::zero() {
392            Some((value - mean) / stddev)
393        } else {
394            None
395        }
396    }
397
398    /// Returns the skewness of values in the rolling window
399    ///
400    /// Skewness measures the asymmetry of the distribution of values.
401    /// - Positive values indicate a right-skewed distribution (long tail on the right)
402    /// - Negative values indicate a left-skewed distribution (long tail on the left)
403    /// - A normal distribution has a skewness of 0
404    ///
405    /// When `ddof` is true, applies a bias correction for sample skewness.
406    ///
407    /// # Returns
408    ///
409    /// * `Option<T>` - The skewness if the window is ready and variance is positive, None otherwise
410    #[inline]
411    pub fn skew(&self) -> Option<T> {
412        if !self.is_ready() || self.m2 <= T::zero() {
413            return None;
414        }
415
416        let n = T::from(self.count())?;
417        let m3 = self.m3;
418        let m2 = self.m2;
419
420        let denominator = m2 * m2.sqrt();
421        if denominator <= T::zero() {
422            return None;
423        }
424        let g1 = m3 / denominator;
425
426        if self.ddof {
427            if n <= T::from(2)? {
428                return None;
429            }
430            let correction = (n * (n - T::one())).sqrt() / (n - T::from(2)?);
431            Some(correction * g1)
432        } else {
433            Some(g1)
434        }
435    }
436
437    /// Returns the excess kurtosis of values in the rolling window
438    ///
439    /// Excess kurtosis measures the "tailedness" of a distribution compared to a normal distribution.
440    /// - Positive values indicate a distribution with heavier tails (more outliers)
441    /// - Negative values indicate a distribution with lighter tails (fewer outliers)
442    /// - A normal distribution has an excess kurtosis of 0
443    ///
444    /// When `ddof` is true, applies a bias correction for sample kurtosis.
445    ///
446    /// # Returns
447    ///
448    /// * `Option<T>` - The excess kurtosis if the window is ready and variance is positive, None otherwise
449    #[inline]
450    pub fn kurt(&self) -> Option<T> {
451        if !self.is_ready() || self.m2 <= T::zero() {
452            return None;
453        }
454
455        let n = T::from(self.count())?;
456        if n < T::from(4.0)? {
457            return None;
458        }
459
460        let _1 = T::from(1.0)?;
461        let _2 = T::from(2.0)?;
462        let _3 = T::from(3.0)?;
463
464        if !self.ddof {
465            let kurt_pop = self.m4 / (self.m2 * self.m2);
466            Some(kurt_pop - _3)
467        } else {
468            let sample_var = self.m2 * n / (n - _1);
469            let numerator = n * n * (n + _1);
470            let denominator = (n - _1) * (n - _2) * (n - _3);
471            let correction = (_3 * (n - _1) * (n - _1)) / ((n - _2) * (n - _3));
472
473            let g2 = (numerator / denominator) * (self.m4 / (sample_var * sample_var)) - correction;
474
475            Some(g2)
476        }
477    }
478}
479
480#[cfg(test)]
481mod tests {
482    use assert_approx_eq::assert_approx_eq;
483
484    use super::*;
485
486    #[test]
487    fn sum_works() {
488        let mut stats = RollingMoments::new(3);
489        let inputs = [
490            1_000_000.1,
491            1_000_000.2,
492            1_000_000.3,
493            1_000_000.4,
494            1_000_000.5,
495            1_000_000.6,
496            1_000_000.7,
497        ];
498        let mut results = vec![];
499
500        inputs.iter().for_each(|i| {
501            if let Some(v) = stats.next(*i).sum() {
502                results.push(v)
503            }
504        });
505
506        let expected = [3000000.6, 3000000.9, 3000001.2, 3000001.5, 3000001.8];
507        assert_eq!(&results, &expected);
508    }
509
510    #[test]
511    fn sum_sq_works() {
512        let mut stats = RollingMoments::new(3);
513        let inputs = [
514            100_000.1, 100_000.2, 100_000.3, 100_000.4, 100_000.5, 100_000.6, 100_000.7,
515        ];
516        let mut results = vec![];
517
518        inputs.iter().for_each(|i| {
519            if let Some(v) = stats.next(*i).sum_sq() {
520                results.push(v)
521            }
522        });
523        let expected = [
524            30000120000.14,
525            30000180000.289997,
526            30000240000.5,
527            30000300000.769997,
528            30000360001.1,
529        ];
530        assert_eq!(&results, &expected);
531    }
532
533    #[test]
534    fn mean_works() {
535        let mut stats = RollingMoments::new(3);
536        let inputs = [
537            1_000_000.1,
538            1_000_000.2,
539            1_000_000.3,
540            1_000_000.4,
541            1_000_000.5,
542            1_000_000.6,
543            1_000_000.7,
544        ];
545        let mut results = vec![];
546
547        inputs.iter().for_each(|i| {
548            if let Some(v) = stats.next(*i).mean() {
549                results.push(v)
550            }
551        });
552
553        let expected: [f64; 5] = [1000000.2, 1000000.3, 1000000.4, 1000000.5, 1000000.6];
554        for (i, e) in expected.iter().enumerate() {
555            assert_approx_eq!(e, results[i], 0.0001);
556        }
557    }
558
559    #[test]
560    fn mean_sq_works() {
561        let mut stats = RollingMoments::new(3);
562        let inputs = [
563            1_000_000.1,
564            1_000_000.2,
565            1_000_000.3,
566            1_000_000.4,
567            1_000_000.5,
568            1_000_000.6,
569            1_000_000.7,
570        ];
571        let mut results = vec![];
572
573        inputs.iter().for_each(|i| {
574            if let Some(v) = stats.next(*i).mean_sq() {
575                results.push(v)
576            }
577        });
578
579        let expected: [f64; 5] = [
580            1000000400000.05,
581            1000000600000.1,
582            1000000800000.17,
583            1000001000000.26,
584            1000001200000.37,
585        ];
586        for (i, e) in expected.iter().enumerate() {
587            assert_approx_eq!(e, results[i], 0.01);
588        }
589    }
590
591    #[test]
592    fn variance_works() {
593        let mut stats = RollingMoments::new(3);
594        let mut results = vec![];
595        let inputs = [25.4, 26.2, 26.0, 26.1, 25.8, 25.9, 26.3, 26.2, 26.5];
596        inputs.iter().for_each(|i| {
597            if let Some(v) = stats.next(*i).variance() {
598                results.push(v)
599            }
600        });
601
602        let expected: [f64; 7] = [0.1156, 0.0067, 0.0156, 0.0156, 0.0467, 0.0289, 0.0156];
603        for (i, e) in expected.iter().enumerate() {
604            assert_approx_eq!(e, results[i], 0.0001);
605        }
606
607        stats.reset().set_ddof(true);
608        results = vec![];
609        inputs.iter().for_each(|i| {
610            if let Some(v) = stats.next(*i).variance() {
611                results.push(v)
612            }
613        });
614
615        let expected: [f64; 7] = [0.1733, 0.01, 0.0233, 0.0233, 0.07, 0.0433, 0.0233];
616        for (i, e) in expected.iter().enumerate() {
617            assert_approx_eq!(e, results[i], 0.0001);
618        }
619    }
620
621    #[test]
622    fn stddev_works() {
623        let mut stats = RollingMoments::new(3);
624        let mut results = vec![];
625        let inputs = [25.4, 26.2, 26.0, 26.1, 25.8, 25.9, 26.3, 26.2, 26.5];
626        inputs.iter().for_each(|i| {
627            if let Some(v) = stats.next(*i).stddev() {
628                results.push(v)
629            }
630        });
631
632        let expected: [f64; 7] = [0.3399, 0.0816, 0.1247, 0.1247, 0.216, 0.17, 0.1247];
633        for (i, e) in expected.iter().enumerate() {
634            assert_approx_eq!(e, results[i], 0.0001);
635        }
636
637        stats.reset().set_ddof(true);
638        results = vec![];
639        inputs.iter().for_each(|i| {
640            if let Some(v) = stats.next(*i).stddev() {
641                results.push(v)
642            }
643        });
644
645        let expected: [f64; 7] = [0.4163, 0.1, 0.1528, 0.1528, 0.2646, 0.2082, 0.1528];
646        for (i, e) in expected.iter().enumerate() {
647            assert_approx_eq!(e, results[i], 0.0001);
648        }
649    }
650
651    #[test]
652    fn zscore_works() {
653        let mut stats = RollingMoments::new(3);
654        let mut results = vec![];
655        let inputs = [1.2, -0.7, 3.4, 2.1, -1.5, 0.0, 2.2, -0.3, 1.5, -2.0];
656        inputs.iter().for_each(|i| {
657            if let Some(v) = stats.next(*i).zscore() {
658                results.push(v)
659            }
660        });
661
662        let expected: [f64; 8] = [
663            1.2535, 0.2923, -1.3671, -0.1355, 1.2943, -0.8374, 0.3482, -1.2129,
664        ];
665        for (i, e) in expected.iter().enumerate() {
666            assert_approx_eq!(e, results[i], 0.0001);
667        }
668
669        stats.reset().set_ddof(true);
670        results = vec![];
671        inputs.iter().for_each(|i| {
672            if let Some(v) = stats.next(*i).zscore() {
673                results.push(v)
674            }
675        });
676
677        let expected: [f64; 8] = [
678            1.0235, 0.2386, -1.1162, -0.1106, 1.0568, -0.6837, 0.2843, -0.9903,
679        ];
680        for (i, e) in expected.iter().enumerate() {
681            assert_approx_eq!(e, results[i], 0.0001);
682        }
683    }
684
685    #[test]
686    fn skew_works() {
687        let mut stats = RollingMoments::new(4);
688        let mut results = vec![];
689        let inputs = [25.4, 26.2, 26.0, 26.1, 25.8, 25.9, 26.3, 26.2, 26.5];
690        inputs.iter().for_each(|i| {
691            if let Some(v) = stats.next(*i).skew() {
692                results.push(v)
693            }
694        });
695
696        let expected: [f64; 6] = [-0.9794, -0.4347, 0.0000, 0.2780, 0.0000, -0.3233];
697        for (i, e) in expected.iter().enumerate() {
698            assert_approx_eq!(e, results[i], 0.0001);
699        }
700
701        stats.reset().set_ddof(true);
702        results = vec![];
703        inputs.iter().for_each(|i| {
704            if let Some(v) = stats.next(*i).skew() {
705                results.push(v)
706            }
707        });
708
709        let expected: [f64; 6] = [-1.6964, -0.7528, 0.0000, 0.4816, 0.0000, -0.5600];
710        for (i, e) in expected.iter().enumerate() {
711            assert_approx_eq!(e, results[i], 0.0001);
712        }
713    }
714
715    #[test]
716    fn kurt_works() {
717        let mut stats = RollingMoments::new(4);
718        let mut results = vec![];
719        let inputs = [25.4, 26.2, 26.0, 26.1, 25.8, 25.9, 26.3, 26.2, 26.5];
720        inputs.iter().for_each(|i| {
721            if let Some(v) = stats.next(*i).kurt() {
722                results.push(v)
723            }
724        });
725
726        let expected: [f64; 6] = [-0.7981, -1.1543, -1.3600, -1.4266, -1.7785, -1.0763];
727        for (i, e) in expected.iter().enumerate() {
728            assert_approx_eq!(e, results[i], 0.0001);
729        }
730
731        stats.reset().set_ddof(true);
732        results = vec![];
733        inputs.iter().for_each(|i| {
734            if let Some(v) = stats.next(*i).kurt() {
735                results.push(v);
736            }
737        });
738
739        let expected: [f64; 6] = [3.0144, 0.3429, -1.2, -1.6995, -4.3391, 0.928];
740        for (i, e) in expected.iter().enumerate() {
741            assert_approx_eq!(e, results[i], 0.0001);
742        }
743    }
744}
ta_statistics/rolling_moments.rs

ta_statistics/
rolling_moments.rs
