Skip to main content

quantwave_core/indicators/
correlation_cycle.rs

1use crate::indicators::metadata::{IndicatorMetadata, ParamDef};
2use crate::traits::Next;
3use std::collections::VecDeque;
4use std::f64::consts::PI;
5
6/// Correlation Cycle Indicator
7///
8/// Based on John Ehlers' "Correlation as a Cycle Indicator".
9/// Computes the correlation of price against a fixed-period Cosine and negative Sine wave
10/// to derive a phase angle. This angle is used to identify market cycles and trends.
11#[derive(Debug, Clone)]
12pub struct CorrelationCycle {
13    period: usize,
14    price_window: VecDeque<f64>,
15    cosine_wave: Vec<f64>,
16    sine_wave: Vec<f64>,
17    prev_angle: f64,
18    count: usize,
19}
20
21impl CorrelationCycle {
22    pub fn new(period: usize) -> Self {
23        let mut cosine_wave = Vec::with_capacity(period);
24        let mut sine_wave = Vec::with_capacity(period);
25        for n in 0..period {
26            let angle = 2.0 * PI * n as f64 / period as f64;
27            cosine_wave.push(angle.cos());
28            sine_wave.push(-angle.sin());
29        }
30
31        Self {
32            period,
33            price_window: VecDeque::with_capacity(period),
34            cosine_wave,
35            sine_wave,
36            prev_angle: 0.0,
37            count: 0,
38        }
39    }
40
41    fn pearson_correlation(n: usize, x: &std::collections::VecDeque<f64>, y: &[f64]) -> f64 {
42        let mut sx = 0.0;
43        let mut sy = 0.0;
44        let mut sxx = 0.0;
45        let mut syy = 0.0;
46        let mut sxy = 0.0;
47
48        for i in 0..n {
49            let xi = x[i];
50            let yi = y[i];
51            sx += xi;
52            sy += yi;
53            sxx += xi * xi;
54            syy += yi * yi;
55            sxy += xi * yi;
56        }
57
58        let nf = n as f64;
59        let num = nf * sxy - sx * sy;
60        let den = ((nf * sxx - sx * sx) * (nf * syy - sy * sy)).sqrt();
61
62        if den > 0.0 {
63            num / den
64        } else {
65            0.0
66        }
67    }
68}
69
70impl Default for CorrelationCycle {
71    fn default() -> Self {
72        Self::new(20)
73    }
74}
75
76impl Next<f64> for CorrelationCycle {
77    type Output = (f64, f64, f64); // (Real, Imag, Angle)
78
79    fn next(&mut self, input: f64) -> Self::Output {
80        self.count += 1;
81        self.price_window.push_front(input);
82        if self.price_window.len() > self.period {
83            self.price_window.pop_back();
84        }
85
86        if self.price_window.len() < self.period {
87            return (0.0, 0.0, 0.0);
88        }
89
90        let real = Self::pearson_correlation(self.period, &self.price_window, &self.cosine_wave);
91        let imag = Self::pearson_correlation(self.period, &self.price_window, &self.sine_wave);
92
93        // Angle in degrees: 90 + atan(Real / Imag)
94        // Ehlers resolve ambiguity: if Imag > 0 then Angle = Angle - 180
95        let mut angle = if imag != 0.0 {
96            (real / imag).atan().to_degrees() + 90.0
97        } else {
98            90.0
99        };
100
101        if imag > 0.0 {
102            angle -= 180.0;
103        }
104
105        // Do not allow rate change of angle to go negative
106        // If Angle[1] - Angle < 270 and Angle < Angle[1] Then Angle = Angle[1];
107        if self.count > self.period + 1 && self.prev_angle - angle < 270.0 && angle < self.prev_angle {
108            angle = self.prev_angle;
109        }
110
111        self.prev_angle = angle;
112
113        (real, imag, angle)
114    }
115}
116
117pub const CORRELATION_CYCLE_METADATA: IndicatorMetadata = IndicatorMetadata {
118    name: "CorrelationCycle",
119    description: "Determines cycle phase angle by correlating price with orthogonal sinusoids.",
120    usage: "Use to measure the dominant cycle period via autocorrelation in an amplitude-independent way. Prefer over DFT methods when price amplitude varies significantly across the measurement window.",
121    keywords: &["cycle", "dominant-cycle", "ehlers", "dsp", "spectral"],
122    ehlers_summary: "Ehlers introduces Correlation Cycle measurement in Cycle Analytics for Traders (2013) as an improvement on DFT. By normalizing autocorrelation coefficients to unity variance, the resulting periodogram is independent of price amplitude variations, producing more consistent cycle period estimates.",
123    params: &[
124        ParamDef {
125            name: "period",
126            default: "20",
127            description: "Correlation wavelength",
128        },
129    ],
130    formula_source: "https://github.com/lavs9/quantwave/blob/main/references/Ehlers%20Papers/CORRELATION%20AS%20A%20CYCLE%20INDICATOR.pdf",
131    formula_latex: r#"
132\[
133R = \text{Corr}(Price, \cos(2\pi n/P)), I = \text{Corr}(Price, -\sin(2\pi n/P))
134\]
135\[
136\text{Angle} = 90 + \arctan(R/I) \text{ (with quadrant resolution)}
137\]
138"#,
139    gold_standard_file: "correlation_cycle.json",
140    category: "Ehlers DSP",
141};
142
143#[cfg(test)]
144mod tests {
145    use super::*;
146    use crate::traits::Next;
147    use proptest::prelude::*;
148
149    #[test]
150    fn test_correlation_cycle_basic() {
151        let mut cc = CorrelationCycle::new(20);
152        for i in 0..100 {
153            let (r, im, a) = cc.next((2.0 * PI * i as f64 / 20.0).sin());
154            assert!(!r.is_nan());
155            assert!(!im.is_nan());
156            assert!(!a.is_nan());
157        }
158    }
159
160    proptest! {
161        #[test]
162        fn test_correlation_cycle_parity(
163            inputs in prop::collection::vec(1.0..100.0, 100..150),
164        ) {
165            let period = 20;
166            let mut cc = CorrelationCycle::new(period);
167            let streaming_results: Vec<(f64, f64, f64)> = inputs.iter().map(|&x| cc.next(x)).collect();
168
169            // Batch implementation
170            let mut batch_results = Vec::with_capacity(inputs.len());
171            let mut cos_w = Vec::new();
172            let mut sin_w = Vec::new();
173            for n in 0..period {
174                let ang = 2.0 * PI * n as f64 / period as f64;
175                cos_w.push(ang.cos());
176                sin_w.push(-ang.sin());
177            }
178
179            let mut prev_a = 0.0;
180            for i in 0..inputs.len() {
181                if i < period - 1 {
182                    batch_results.push((0.0, 0.0, 0.0));
183                    continue;
184                }
185                
186                let mut sx = 0.0;
187                let mut sy_c = 0.0;
188                let mut sy_s = 0.0;
189                let mut sxx = 0.0;
190                let mut syy_c = 0.0;
191                let mut syy_s = 0.0;
192                let mut sxy_c = 0.0;
193                let mut sxy_s = 0.0;
194
195                for j in 0..period {
196                    let xi = inputs[i - j];
197                    let yc = cos_w[j];
198                    let ys = sin_w[j];
199                    sx += xi;
200                    sy_c += yc;
201                    sy_s += ys;
202                    sxx += xi * xi;
203                    syy_c += yc * yc;
204                    syy_s += ys * ys;
205                    sxy_c += xi * yc;
206                    sxy_s += xi * ys;
207                }
208
209                let nf = period as f64;
210                let den_c = ((nf * sxx - sx * sx) * (nf * syy_c - sy_c * sy_c)).sqrt();
211                let real = if den_c > 0.0 { (nf * sxy_c - sx * sy_c) / den_c } else { 0.0 };
212                
213                let den_s = ((nf * sxx - sx * sx) * (nf * syy_s - sy_s * sy_s)).sqrt();
214                let imag = if den_s > 0.0 { (nf * sxy_s - sx * sy_s) / den_s } else { 0.0 };
215
216                let mut angle = if imag != 0.0 {
217                    (real / imag).atan().to_degrees() + 90.0
218                } else {
219                    90.0
220                };
221                if imag > 0.0 { angle -= 180.0; }
222
223                if i > period {
224                    if prev_a - angle < 270.0 && angle < prev_a {
225                        angle = prev_a;
226                    }
227                }
228                
229                prev_a = angle;
230                batch_results.push((real, imag, angle));
231            }
232
233            for (s, b) in streaming_results.iter().zip(batch_results.iter()) {
234                approx::assert_relative_eq!(s.0, b.0, epsilon = 1e-10);
235                approx::assert_relative_eq!(s.1, b.1, epsilon = 1e-10);
236                approx::assert_relative_eq!(s.2, b.2, epsilon = 1e-10);
237            }
238        }
239    }
240}