Skip to main content

quantwave_core/indicators/
recursive_median.rs

1use crate::indicators::metadata::{IndicatorMetadata, ParamDef};
2use crate::traits::Next;
3use crate::utils::RingBuffer as VecDeque;
4
5/// Recursive Median Filter
6///
7/// Based on John Ehlers' "Recursive Median Filters" (TASC March 2018).
8/// It applies an EMA smoothing to a 5-bar median filter.
9#[derive(Debug, Clone)]
10pub struct RecursiveMedian {
11    _lp_period: usize,
12    alpha1: f64,
13    window: VecDeque<f64>,
14    prev_rm: f64,
15}
16
17impl RecursiveMedian {
18    pub fn new(lp_period: usize) -> Self {
19        let p = lp_period as f64;
20        let deg_to_rad = std::f64::consts::PI / 180.0;
21        let alpha1 = ((360.0 / p * deg_to_rad).cos() + (360.0 / p * deg_to_rad).sin() - 1.0)
22            / (360.0 / p * deg_to_rad).cos();
23
24        Self {
25            _lp_period: lp_period,
26            alpha1,
27            window: VecDeque::with_capacity(5),
28            prev_rm: 0.0,
29        }
30    }
31}
32
33impl Default for RecursiveMedian {
34    fn default() -> Self {
35        Self::new(12)
36    }
37}
38
39impl Next<f64> for RecursiveMedian {
40    type Output = f64;
41
42    fn next(&mut self, input: f64) -> Self::Output {
43        self.window.push_front(input);
44        if self.window.len() > 5 {
45            self.window.pop_back();
46        }
47
48        if self.window.len() < 5 {
49            self.prev_rm = input;
50            return input;
51        }
52
53        let mut sorted: Vec<f64> = self.window.iter().copied().collect();
54        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
55        let median = sorted[2]; // Median of 5 elements
56
57        let rm = self.alpha1 * median + (1.0 - self.alpha1) * self.prev_rm;
58        self.prev_rm = rm;
59        rm
60    }
61}
62
63/// Recursive Median Oscillator
64///
65/// Based on John Ehlers' "Recursive Median Filters" (TASC March 2018).
66/// It applies a 2nd-order Highpass filter to the Recursive Median Filter output.
67#[derive(Debug, Clone)]
68pub struct RecursiveMedianOscillator {
69    _lp_period: usize,
70    _hp_period: usize,
71    rm_filter: RecursiveMedian,
72    alpha2: f64,
73    prev_rm: [f64; 2],
74    prev_rmo: [f64; 2],
75}
76
77impl RecursiveMedianOscillator {
78    pub fn new(lp_period: usize, hp_period: usize) -> Self {
79        let p = hp_period as f64;
80        let deg_to_rad = std::f64::consts::PI / 180.0;
81        let angle = 0.707 * 360.0 / p * deg_to_rad;
82        let alpha2 = (angle.cos() + angle.sin() - 1.0) / angle.cos();
83
84        Self {
85            _lp_period: lp_period,
86            _hp_period: hp_period,
87            rm_filter: RecursiveMedian::new(lp_period),
88            alpha2,
89            prev_rm: [0.0; 2],
90            prev_rmo: [0.0; 2],
91        }
92    }
93}
94
95impl Default for RecursiveMedianOscillator {
96    fn default() -> Self {
97        Self::new(12, 30)
98    }
99}
100
101impl Next<f64> for RecursiveMedianOscillator {
102    type Output = f64;
103
104    fn next(&mut self, input: f64) -> Self::Output {
105        let rm = self.rm_filter.next(input);
106
107        // Highpass filter calculation:
108        // RMO = (1-alpha2/2)*(1-alpha2/2)*(RM-2*RM[1]+RM[2])+2*(1-alpha2)*RMO[1]-(1-alpha2)*(1-alpha2)*RMO[2]
109        let c1 = (1.0 - self.alpha2 / 2.0) * (1.0 - self.alpha2 / 2.0);
110        let c2 = 2.0 * (1.0 - self.alpha2);
111        let c3 = (1.0 - self.alpha2) * (1.0 - self.alpha2);
112
113        let rmo = c1 * (rm - 2.0 * self.prev_rm[0] + self.prev_rm[1]) + c2 * self.prev_rmo[0]
114            - c3 * self.prev_rmo[1];
115
116        self.prev_rm[1] = self.prev_rm[0];
117        self.prev_rm[0] = rm;
118
119        self.prev_rmo[1] = self.prev_rmo[0];
120        self.prev_rmo[0] = rmo;
121
122        rmo
123    }
124}
125
126pub const RECURSIVE_MEDIAN_METADATA: IndicatorMetadata = IndicatorMetadata {
127    name: "RecursiveMedian",
128    description: "EMA of a 5-bar median filter for smooth tracking with minimal jitter.",
129    usage: "Use to filter out extreme outliers and noise while maintaining trend sensitivity. Excellent as a baseline for other oscillators.",
130    keywords: &["filter", "ehlers", "dsp", "median", "robust", "smoothing"],
131    ehlers_summary: "Standard filters like SMA or EMA are distorted by price spikes. The recursive median filter uses the median to reject outliers and an EMA to provide smoothness, offering a cleaner trend representation than standard moving averages.",
132    params: &[ParamDef {
133        name: "lp_period",
134        default: "12",
135        description: "Low-pass smoothing period",
136    }],
137    formula_source: "https://www.traders.com/Documentation/FEEDbk_docs/2018/03/TradersTips.html",
138    formula_latex: r#"
139\[
140\alpha = \frac{\cos(360/P) + \sin(360/P) - 1}{\cos(360/P)}
141\]
142\[
143RM_t = \alpha \cdot \text{Median}(Price, 5)_t + (1 - \alpha) \cdot RM_{t-1}
144\]
145"#,
146    gold_standard_file: "recursive_median.json",
147    category: "Ehlers DSP",
148};
149
150pub const RECURSIVE_MEDIAN_OSCILLATOR_METADATA: IndicatorMetadata = IndicatorMetadata {
151    name: "RecursiveMedianOscillator",
152    description: "Oscillator derived from the Recursive Median filter using a 2nd-order Highpass filter.",
153    usage: "Identify cyclic turning points with reduced lag and noise. The high-pass component removes the trend, leaving the cycle.",
154    keywords: &["oscillator", "ehlers", "dsp", "median", "cycle", "highpass"],
155    ehlers_summary: "By applying a 2nd-order Highpass filter to the Recursive Median output, we create an oscillator that is specifically tuned to the dominant cycle while remaining immune to the outlier spikes that would otherwise create false signals.",
156    params: &[
157        ParamDef {
158            name: "lp_period",
159            default: "12",
160            description: "Low-pass smoothing period",
161        },
162        ParamDef {
163            name: "hp_period",
164            default: "30",
165            description: "High-pass cutoff period",
166        },
167    ],
168    formula_source: "https://www.traders.com/Documentation/FEEDbk_docs/2018/03/TradersTips.html",
169    formula_latex: r#"
170\[
171\alpha_2 = \frac{\cos(0.707 \cdot 360/HP) + \sin(0.707 \cdot 360/HP) - 1}{\cos(0.707 \cdot 360/HP)}
172\]
173\[
174RMO_t = (1-\alpha_2/2)^2(RM_t - 2RM_{t-1} + RM_{t-2}) + 2(1-\alpha_2)RMO_{t-1} - (1-\alpha_2)^2RMO_{t-2}
175\]
176"#,
177    gold_standard_file: "recursive_median_oscillator.json",
178    category: "Ehlers DSP",
179};
180
181#[cfg(test)]
182mod tests {
183    use super::*;
184    use crate::traits::Next;
185    use proptest::prelude::*;
186
187    #[test]
188    fn test_recursive_median_basic() {
189        let mut rm = RecursiveMedian::new(12);
190        for _ in 0..50 {
191            let val = rm.next(100.0);
192            approx::assert_relative_eq!(val, 100.0, epsilon = 1e-10);
193        }
194    }
195
196    proptest! {
197        #[test]
198        fn test_recursive_median_parity(
199            inputs in prop::collection::vec(1.0..100.0, 50..100),
200        ) {
201            let lp_period = 12;
202            let mut rm = RecursiveMedian::new(lp_period);
203            let streaming_results: Vec<f64> = inputs.iter().map(|&x| rm.next(x)).collect();
204
205            // Batch implementation
206            let mut batch_results = Vec::with_capacity(inputs.len());
207            let p = lp_period as f64;
208            let deg_to_rad = std::f64::consts::PI / 180.0;
209            let alpha1 = ((360.0 / p * deg_to_rad).cos() + (360.0 / p * deg_to_rad).sin() - 1.0)
210                / (360.0 / p * deg_to_rad).cos();
211
212            let mut prev_rm = 0.0;
213            for i in 0..inputs.len() {
214                if i < 4 {
215                    prev_rm = inputs[i];
216                    batch_results.push(inputs[i]);
217                    continue;
218                }
219                let mut window = vec![
220                    inputs[i], inputs[i-1], inputs[i-2], inputs[i-3], inputs[i-4]
221                ];
222                window.sort_by(|a, b| a.partial_cmp(b).unwrap());
223                let median = window[2];
224                let val = alpha1 * median + (1.0 - alpha1) * prev_rm;
225                batch_results.push(val);
226                prev_rm = val;
227            }
228
229            for (s, b) in streaming_results.iter().zip(batch_results.iter()) {
230                approx::assert_relative_eq!(s, b, epsilon = 1e-10);
231            }
232        }
233    }
234}