Skip to main content

quantwave_core/indicators/
recursive_median.rs

1use crate::indicators::metadata::{IndicatorMetadata, ParamDef};
2use crate::traits::Next;
3use std::collections::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])
114            + c2 * self.prev_rmo[0]
115            - c3 * self.prev_rmo[1];
116
117        self.prev_rm[1] = self.prev_rm[0];
118        self.prev_rm[0] = rm;
119
120        self.prev_rmo[1] = self.prev_rmo[0];
121        self.prev_rmo[0] = rmo;
122
123        rmo
124    }
125}
126
127pub const RECURSIVE_MEDIAN_METADATA: IndicatorMetadata = IndicatorMetadata {
128    name: "RecursiveMedian",
129    description: "EMA of a 5-bar median filter for smooth tracking with minimal jitter.",
130    usage: "Use to filter out extreme outliers and noise while maintaining trend sensitivity. Excellent as a baseline for other oscillators.",
131    keywords: &["filter", "ehlers", "dsp", "median", "robust", "smoothing"],
132    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.",
133    params: &[
134        ParamDef {
135            name: "lp_period",
136            default: "12",
137            description: "Low-pass smoothing period",
138        },
139    ],
140    formula_source: "https://www.traders.com/Documentation/FEEDbk_docs/2018/03/TradersTips.html",
141    formula_latex: r#"
142\[
143\alpha = \frac{\cos(360/P) + \sin(360/P) - 1}{\cos(360/P)}
144\]
145\[
146RM_t = \alpha \cdot \text{Median}(Price, 5)_t + (1 - \alpha) \cdot RM_{t-1}
147\]
148"#,
149    gold_standard_file: "recursive_median.json",
150    category: "Ehlers DSP",
151};
152
153pub const RECURSIVE_MEDIAN_OSCILLATOR_METADATA: IndicatorMetadata = IndicatorMetadata {
154    name: "RecursiveMedianOscillator",
155    description: "Oscillator derived from the Recursive Median filter using a 2nd-order Highpass filter.",
156    usage: "Identify cyclic turning points with reduced lag and noise. The high-pass component removes the trend, leaving the cycle.",
157    keywords: &["oscillator", "ehlers", "dsp", "median", "cycle", "highpass"],
158    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.",
159    params: &[
160        ParamDef {
161            name: "lp_period",
162            default: "12",
163            description: "Low-pass smoothing period",
164        },
165        ParamDef {
166            name: "hp_period",
167            default: "30",
168            description: "High-pass cutoff period",
169        },
170    ],
171    formula_source: "https://www.traders.com/Documentation/FEEDbk_docs/2018/03/TradersTips.html",
172    formula_latex: r#"
173\[
174\alpha_2 = \frac{\cos(0.707 \cdot 360/HP) + \sin(0.707 \cdot 360/HP) - 1}{\cos(0.707 \cdot 360/HP)}
175\]
176\[
177RMO_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}
178\]
179"#,
180    gold_standard_file: "recursive_median_oscillator.json",
181    category: "Ehlers DSP",
182};
183
184#[cfg(test)]
185mod tests {
186    use super::*;
187    use crate::traits::Next;
188    use proptest::prelude::*;
189
190    #[test]
191    fn test_recursive_median_basic() {
192        let mut rm = RecursiveMedian::new(12);
193        for _ in 0..50 {
194            let val = rm.next(100.0);
195            approx::assert_relative_eq!(val, 100.0, epsilon = 1e-10);
196        }
197    }
198
199    proptest! {
200        #[test]
201        fn test_recursive_median_parity(
202            inputs in prop::collection::vec(1.0..100.0, 50..100),
203        ) {
204            let lp_period = 12;
205            let mut rm = RecursiveMedian::new(lp_period);
206            let streaming_results: Vec<f64> = inputs.iter().map(|&x| rm.next(x)).collect();
207
208            // Batch implementation
209            let mut batch_results = Vec::with_capacity(inputs.len());
210            let p = lp_period as f64;
211            let deg_to_rad = std::f64::consts::PI / 180.0;
212            let alpha1 = ((360.0 / p * deg_to_rad).cos() + (360.0 / p * deg_to_rad).sin() - 1.0)
213                / (360.0 / p * deg_to_rad).cos();
214
215            let mut prev_rm = 0.0;
216            for i in 0..inputs.len() {
217                if i < 4 {
218                    prev_rm = inputs[i];
219                    batch_results.push(inputs[i]);
220                    continue;
221                }
222                let mut window = vec![
223                    inputs[i], inputs[i-1], inputs[i-2], inputs[i-3], inputs[i-4]
224                ];
225                window.sort_by(|a, b| a.partial_cmp(b).unwrap());
226                let median = window[2];
227                let val = alpha1 * median + (1.0 - alpha1) * prev_rm;
228                batch_results.push(val);
229                prev_rm = val;
230            }
231
232            for (s, b) in streaming_results.iter().zip(batch_results.iter()) {
233                approx::assert_relative_eq!(s, b, epsilon = 1e-10);
234            }
235        }
236    }
237}