Skip to main content

quantwave_core/indicators/
griffiths_dominant_cycle.rs

1use crate::indicators::metadata::{IndicatorMetadata, ParamDef};
2use crate::traits::Next;
3use crate::indicators::high_pass::HighPass;
4use crate::indicators::super_smoother::SuperSmoother;
5use std::collections::VecDeque;
6use std::f64::consts::PI;
7
8/// Griffiths Dominant Cycle
9///
10/// Based on John Ehlers' "Linear Predictive Filters And Instantaneous Frequency" (TASC January 2025).
11/// Uses the Griffiths spectral estimation to find the dominant cycle in the data.
12#[derive(Debug, Clone)]
13pub struct GriffithsDominantCycle {
14    lb: usize,
15    ub: usize,
16    length: usize,
17    mu: f64,
18    hp: HighPass,
19    ss: SuperSmoother,
20    peak: f64,
21    signal_window: VecDeque<f64>,
22    coef: Vec<f64>,
23    prev_cycle: f64,
24}
25
26impl GriffithsDominantCycle {
27    pub fn new(lower_bound: usize, upper_bound: usize, length: usize) -> Self {
28        Self {
29            lb: lower_bound,
30            ub: upper_bound,
31            length,
32            mu: 1.0 / (length as f64),
33            hp: HighPass::new(upper_bound),
34            ss: SuperSmoother::new(lower_bound),
35            peak: 0.1,
36            signal_window: VecDeque::with_capacity(length + 1),
37            coef: vec![0.0; length + 1],
38            prev_cycle: (lower_bound + upper_bound) as f64 / 2.0,
39        }
40    }
41}
42
43impl Default for GriffithsDominantCycle {
44    fn default() -> Self {
45        Self::new(18, 40, 40)
46    }
47}
48
49impl Next<f64> for GriffithsDominantCycle {
50    type Output = f64;
51
52    fn next(&mut self, input: f64) -> Self::Output {
53        let hp_val = self.hp.next(input);
54        let lp_val = self.ss.next(hp_val);
55
56        self.peak *= 0.991;
57        if lp_val.abs() > self.peak {
58            self.peak = lp_val.abs();
59        }
60
61        let signal = if self.peak != 0.0 {
62            lp_val / self.peak
63        } else {
64            0.0
65        };
66
67        self.signal_window.push_front(signal);
68        if self.signal_window.len() > self.length {
69            self.signal_window.pop_back();
70        }
71
72        if self.signal_window.len() < self.length {
73            return self.prev_cycle;
74        }
75
76        let mut xx = vec![0.0; self.length + 1];
77        for (i, val) in xx.iter_mut().enumerate().skip(1).take(self.length) {
78            *val = self.signal_window[self.length - i];
79        }
80
81        let mut x_bar = 0.0;
82        for count in 1..=self.length {
83            x_bar += xx[self.length - count] * self.coef[count];
84        }
85
86        for count in 1..=self.length {
87            self.coef[count] += self.mu * (xx[self.length] - x_bar) * xx[self.length - count];
88        }
89
90        // Spectral scan
91        let mut max_pwr = 0.0;
92        let mut cycle = self.prev_cycle;
93
94        for period_idx in self.lb..=self.ub {
95            let period = period_idx as f64;
96            let mut real = 0.0;
97            let mut imag = 0.0;
98
99            for count in 1..=self.length {
100                let angle = 2.0 * PI * (count as f64) / period;
101                real += self.coef[count] * angle.cos();
102                imag += self.coef[count] * angle.sin();
103            }
104
105            let denom = (1.0 - real).powi(2) + imag.powi(2);
106            let pwr = 0.1 / denom;
107
108            if pwr > max_pwr {
109                max_pwr = pwr;
110                cycle = period;
111            }
112        }
113
114        // Slew rate limiter
115        if cycle > self.prev_cycle + 2.0 {
116            cycle = self.prev_cycle + 2.0;
117        } else if cycle < self.prev_cycle - 2.0 {
118            cycle = self.prev_cycle - 2.0;
119        }
120
121        self.prev_cycle = cycle;
122        cycle
123    }
124}
125
126pub const GRIFFITHS_DOMINANT_CYCLE_METADATA: IndicatorMetadata = IndicatorMetadata {
127    name: "GriffithsDominantCycle",
128    description: "Dominant cycle estimation using Griffiths adaptive spectral analysis.",
129    usage: "Use as a robust dominant cycle estimator less sensitive to amplitude changes than DFT-based methods, making it reliable across different market volatility regimes.",
130    keywords: &["cycle", "dominant-cycle", "ehlers", "dsp", "spectral"],
131    ehlers_summary: "The Griffiths method computes the dominant cycle by solving the real-roots of an autocorrelation polynomial. Adapted by Ehlers in Cycle Analytics for Traders, it remains stable even when market amplitude changes rapidly, unlike power-spectrum methods that can shift with volatility.",
132    params: &[
133        ParamDef {
134            name: "lower_bound",
135            default: "18",
136            description: "Lower period bound",
137        },
138        ParamDef {
139            name: "upper_bound",
140            default: "40",
141            description: "Upper period bound",
142        },
143        ParamDef {
144            name: "length",
145            default: "40",
146            description: "LMS filter length",
147        },
148    ],
149    formula_source: "https://github.com/lavs9/quantwave/blob/main/references/traderstipsreference/TRADERS’%20TIPS%20-%20JANUARY%202025.html",
150    formula_latex: r#"
151\[
152Pwr(Period) = \frac{0.1}{(1-Real)^2 + Imag^2}
153\]
154\[
155Real = \sum coef_i \cos(2\pi i / Period)
156\]
157\[
158Imag = \sum coef_i \sin(2\pi i / Period)
159\]
160"#,
161    gold_standard_file: "griffiths_dominant_cycle.json",
162    category: "Ehlers DSP",
163};
164
165#[cfg(test)]
166mod tests {
167    use super::*;
168    use crate::traits::Next;
169    use proptest::prelude::*;
170
171    #[test]
172    fn test_griffiths_dc_basic() {
173        let mut gdc = GriffithsDominantCycle::new(18, 40, 40);
174        for i in 0..200 {
175            // Sine wave with period 30
176            let val = gdc.next((2.0 * PI * i as f64 / 30.0).sin());
177            if i > 150 {
178                // Should converge towards 30
179                assert!(val > 25.0 && val < 35.0);
180            }
181        }
182    }
183
184    proptest! {
185        #[test]
186        fn test_griffiths_dc_parity(
187            inputs in prop::collection::vec(1.0..100.0, 100..200),
188        ) {
189            let lb = 18;
190            let ub = 40;
191            let length = 40;
192            let mut gdc = GriffithsDominantCycle::new(lb, ub, length);
193            let streaming_results: Vec<f64> = inputs.iter().map(|&x| gdc.next(x)).collect();
194
195            // Batch implementation
196            let mut batch_results = Vec::with_capacity(inputs.len());
197            let mut hp = HighPass::new(ub);
198            let mut ss = SuperSmoother::new(lb);
199            let lp_vals: Vec<f64> = inputs.iter().map(|&x| ss.next(hp.next(x))).collect();
200
201            let mut peak = 0.1;
202            let mut signals = Vec::new();
203            let mut coef = vec![0.0; length + 1];
204            let mu = 1.0 / length as f64;
205            let mut prev_cycle = (lb + ub) as f64 / 2.0;
206
207            for (i, &lp_val) in lp_vals.iter().enumerate() {
208                peak *= 0.991;
209                if lp_val.abs() > peak {
210                    peak = lp_val.abs();
211                }
212                let signal = if peak != 0.0 { lp_val / peak } else { 0.0 };
213                signals.push(signal);
214
215                if signals.len() < length {
216                    batch_results.push(prev_cycle);
217                    continue;
218                }
219
220                let mut xx = vec![0.0; length + 1];
221                for j in 1..=length {
222                    xx[j] = signals[i - (length - j)];
223                }
224
225                let mut x_bar = 0.0;
226                for count in 1..=length {
227                    x_bar += xx[length - count] * coef[count];
228                }
229
230                for count in 1..=length {
231                    coef[count] += mu * (xx[length] - x_bar) * xx[length - count];
232                }
233
234                let mut max_pwr = 0.0;
235                let mut cycle = prev_cycle;
236
237                for period_idx in lb..=ub {
238                    let period = period_idx as f64;
239                    let mut real = 0.0;
240                    let mut imag = 0.0;
241                    for count in 1..=length {
242                        let angle = 2.0 * PI * (count as f64) / period;
243                        real += coef[count] * angle.cos();
244                        imag += coef[count] * angle.sin();
245                    }
246                    let denom = (1.0 - real).powi(2) + imag.powi(2);
247                    let pwr = 0.1 / denom;
248                    if pwr > max_pwr {
249                        max_pwr = pwr;
250                        cycle = period;
251                    }
252                }
253
254                if cycle > prev_cycle + 2.0 {
255                    cycle = prev_cycle + 2.0;
256                } else if cycle < prev_cycle - 2.0 {
257                    cycle = prev_cycle - 2.0;
258                }
259                
260                prev_cycle = cycle;
261                batch_results.push(cycle);
262            }
263
264            for (s, b) in streaming_results.iter().zip(batch_results.iter()) {
265                approx::assert_relative_eq!(s, b, epsilon = 1e-10);
266            }
267        }
268    }
269}