quant_indicators/
hurst.rs

1//! Hurst Exponent via Rescaled Range (R/S) analysis.
2//!
3//! Classifies instruments as trending (H > 0.5), mean-reverting (H < 0.5),
4//! or random walk (H ≈ 0.5).
5//!
6//! # Algorithm
7//!
8//! 1. Compute log-returns from closing prices
9//! 2. For each subseries size n in [16, 32, 64], compute the mean R/S statistic
10//! 3. Estimate H via log2 slope: H = log2(RS_last / RS_first) / log2(n_last / n_first)
11//!
12//! # References
13//!
14//! Hurst (1951) "Long-term storage capacity of reservoirs"
15
16use quant_primitives::Candle;
17use rust_decimal::Decimal;
18
19use crate::error::IndicatorError;
20
21/// Hurst Exponent regime classifier via R/S analysis.
22///
23/// Uses a fixed set of subseries sizes [16, 32, 64] and estimates the
24/// Hurst exponent from the log-log slope of average R/S vs subseries size.
25#[derive(Debug, Clone)]
26pub struct HurstExponent {
27    window: usize,
28}
29
30/// Subseries sizes used for R/S regression.
31const SUBSERIES_SIZES: [usize; 3] = [16, 32, 64];
32
33/// Minimum number of returns required for valid computation.
34const MIN_RETURNS: usize = 20;
35
36impl HurstExponent {
37    /// Create a new `HurstExponent` with the specified lookback window.
38    ///
39    /// The window determines how many recent candles are used for computation.
40    ///
41    /// # Errors
42    ///
43    /// Returns `InvalidParameter` if window is less than 64 (minimum for
44    /// meaningful R/S analysis with subseries sizes up to 64).
45    pub fn new(window: usize) -> Result<Self, IndicatorError> {
46        if window < 64 {
47            return Err(IndicatorError::InvalidParameter {
48                message: format!("HurstExponent window must be >= 64, got {window}"),
49            });
50        }
51        Ok(Self { window })
52    }
53
54    /// Compute the Hurst exponent from candle closing prices.
55    ///
56    /// Returns `None` if insufficient data or degenerate price series.
57    pub fn compute_from_candles(&self, candles: &[Candle]) -> Option<Decimal> {
58        if candles.len() < self.window {
59            return None;
60        }
61        let closes: Vec<Decimal> = candles[candles.len() - self.window..]
62            .iter()
63            .map(|c| c.close())
64            .collect();
65        let returns = compute_log_returns(&closes);
66        if returns.len() < MIN_RETURNS {
67            return None;
68        }
69        estimate_hurst_from_returns(&returns)
70    }
71}
72
73/// Compute log-returns as simple percentage changes: (P[i] - P[i-1]) / P[i-1].
74///
75/// Skips any pair where the denominator is zero.
76pub fn compute_log_returns(closes: &[Decimal]) -> Vec<Decimal> {
77    closes
78        .windows(2)
79        .filter_map(|w| {
80            if w[0].is_zero() {
81                None
82            } else {
83                Some((w[1] - w[0]) / w[0])
84            }
85        })
86        .collect()
87}
88
89/// Compute the average Rescaled Range statistic for a given subseries size `n`.
90///
91/// Splits the return series into non-overlapping subseries of size `n`,
92/// computes R/S for each, and returns the average. Returns `None` if no
93/// valid subseries produced a result.
94pub fn rescaled_range_for_n(returns: &[Decimal], n: usize) -> Option<Decimal> {
95    if n > returns.len() || n == 0 {
96        return None;
97    }
98    let num_subseries = returns.len() / n;
99    if num_subseries == 0 {
100        return None;
101    }
102
103    let mut rs_sum = Decimal::ZERO;
104    let mut valid_count = 0u32;
105
106    for k in 0..num_subseries {
107        let sub = &returns[k * n..(k + 1) * n];
108        let n_dec = Decimal::from(sub.len());
109        let mean: Decimal = sub.iter().copied().sum::<Decimal>() / n_dec;
110        let deviations: Vec<Decimal> = sub.iter().map(|r| r - mean).collect();
111
112        // Cumulative sum of deviations
113        let mut cum = Vec::with_capacity(sub.len());
114        let mut running = Decimal::ZERO;
115        for d in &deviations {
116            running += d;
117            cum.push(running);
118        }
119
120        let r = cum.iter().copied().max().unwrap_or(Decimal::ZERO)
121            - cum.iter().copied().min().unwrap_or(Decimal::ZERO);
122
123        // Standard deviation
124        let var: Decimal = deviations.iter().map(|d| d * d).sum::<Decimal>() / n_dec;
125        if var.is_zero() {
126            continue;
127        }
128        let s = newton_sqrt(var);
129        if s.is_zero() {
130            continue;
131        }
132
133        rs_sum += r / s;
134        valid_count += 1;
135    }
136
137    if valid_count == 0 {
138        return None;
139    }
140    Some(rs_sum / Decimal::from(valid_count))
141}
142
143/// Estimate the Hurst exponent H from the log-log slope of R/S vs subseries size.
144///
145/// Uses the first and last valid (n, avg_rs) points to compute:
146/// H = log2(RS_last / RS_first) / log2(n_last / n_first)
147///
148/// Returns `None` if fewer than 2 valid data points.
149pub fn estimate_slope(log_ns: &[usize], log_rs: &[Decimal]) -> Option<Decimal> {
150    if log_ns.len() < 2 || log_rs.len() < 2 {
151        return None;
152    }
153
154    let rs1 = log_rs[0];
155    let rs2 = log_rs[log_rs.len() - 1];
156    let n1 = Decimal::from(log_ns[0]);
157    let n2 = Decimal::from(log_ns[log_ns.len() - 1]);
158
159    if rs1.is_zero() || n1.is_zero() {
160        return None;
161    }
162
163    let rs_ratio = rs2 / rs1;
164    let n_ratio = n2 / n1;
165    let h = approx_log2(rs_ratio) / approx_log2(n_ratio);
166    Some(h.max(Decimal::ZERO).min(Decimal::ONE))
167}
168
169/// Full Hurst estimation pipeline from a return series.
170fn estimate_hurst_from_returns(returns: &[Decimal]) -> Option<Decimal> {
171    let mut log_ns = Vec::new();
172    let mut log_rs = Vec::new();
173
174    for &n in &SUBSERIES_SIZES {
175        if let Some(avg_rs) = rescaled_range_for_n(returns, n) {
176            log_ns.push(n);
177            log_rs.push(avg_rs);
178        }
179    }
180
181    estimate_slope(&log_ns, &log_rs)
182}
183
184/// Newton-Raphson square root for Decimal.
185fn newton_sqrt(val: Decimal) -> Decimal {
186    if val.is_zero() || val < Decimal::ZERO {
187        return Decimal::ZERO;
188    }
189    let mut guess = val / Decimal::TWO;
190    for _ in 0..10 {
191        if guess.is_zero() {
192            return Decimal::ZERO;
193        }
194        guess = (guess + val / guess) / Decimal::TWO;
195    }
196    guess
197}
198
199/// Approximate log2(x) for positive x using series expansion.
200fn approx_log2(x: Decimal) -> Decimal {
201    if x <= Decimal::ZERO {
202        return Decimal::ZERO;
203    }
204    let ln2 = Decimal::new(6931, 4); // 0.6931
205
206    let mut reduced = x;
207    let mut shifts = 0i64;
208    while reduced >= Decimal::TWO {
209        reduced /= Decimal::TWO;
210        shifts += 1;
211    }
212    while reduced < Decimal::ONE && !reduced.is_zero() {
213        reduced *= Decimal::TWO;
214        shifts -= 1;
215    }
216
217    let u = (reduced - Decimal::ONE) / (reduced + Decimal::ONE);
218    let u2 = u * u;
219    let mut term = u;
220    let mut sum = u;
221    for k in 1..10 {
222        term *= u2;
223        sum += term / Decimal::from(2 * k + 1);
224    }
225    let ln_reduced = sum * Decimal::TWO;
226    let ln_x = ln_reduced + Decimal::from(shifts) * ln2;
227    ln_x / ln2
228}
229
230#[cfg(test)]
231#[path = "hurst_tests.rs"]
232mod tests;
quant_indicators/hurst.rs

quant_indicators/
hurst.rs
