quant-indicators 0.7.0

//! Hurst Exponent via Rescaled Range (R/S) analysis.
//!
//! Classifies instruments as trending (H > 0.5), mean-reverting (H < 0.5),
//! or random walk (H ≈ 0.5).
//!
//! # Algorithm
//!
//! 1. Compute log-returns from closing prices
//! 2. For each subseries size n in [16, 32, 64], compute the mean R/S statistic
//! 3. Estimate H via log2 slope: H = log2(RS_last / RS_first) / log2(n_last / n_first)
//!
//! # References
//!
//! Hurst (1951) "Long-term storage capacity of reservoirs"

use quant_primitives::Candle;
use rust_decimal::Decimal;

use crate::error::IndicatorError;

/// Hurst Exponent regime classifier via R/S analysis.
///
/// Uses a fixed set of subseries sizes [16, 32, 64] and estimates the
/// Hurst exponent from the log-log slope of average R/S vs subseries size.
#[derive(Debug, Clone)]
pub struct HurstExponent {
    window: usize,
}

/// Subseries sizes used for R/S regression.
const SUBSERIES_SIZES: [usize; 3] = [16, 32, 64];

/// Minimum number of returns required for valid computation.
const MIN_RETURNS: usize = 20;

impl HurstExponent {
    /// Create a new `HurstExponent` with the specified lookback window.
    ///
    /// The window determines how many recent candles are used for computation.
    ///
    /// # Errors
    ///
    /// Returns `InvalidParameter` if window is less than 64 (minimum for
    /// meaningful R/S analysis with subseries sizes up to 64).
    pub fn new(window: usize) -> Result<Self, IndicatorError> {
        if window < 64 {
            return Err(IndicatorError::InvalidParameter {
                message: format!("HurstExponent window must be >= 64, got {window}"),
            });
        }
        Ok(Self { window })
    }

    /// Compute the Hurst exponent from candle closing prices.
    ///
    /// Returns `None` if insufficient data or degenerate price series.
    pub fn compute_from_candles(&self, candles: &[Candle]) -> Option<Decimal> {
        if candles.len() < self.window {
            return None;
        }
        let closes: Vec<Decimal> = candles[candles.len() - self.window..]
            .iter()
            .map(|c| c.close())
            .collect();
        let returns = compute_log_returns(&closes);
        if returns.len() < MIN_RETURNS {
            return None;
        }
        estimate_hurst_from_returns(&returns)
    }
}

/// Compute log-returns as simple percentage changes: (P[i] - P[i-1]) / P[i-1].
///
/// Skips any pair where the denominator is zero.
pub fn compute_log_returns(closes: &[Decimal]) -> Vec<Decimal> {
    closes
        .windows(2)
        .filter_map(|w| {
            if w[0].is_zero() {
                None
            } else {
                Some((w[1] - w[0]) / w[0])
            }
        })
        .collect()
}

/// Compute the average Rescaled Range statistic for a given subseries size `n`.
///
/// Splits the return series into non-overlapping subseries of size `n`,
/// computes R/S for each, and returns the average. Returns `None` if no
/// valid subseries produced a result.
pub fn rescaled_range_for_n(returns: &[Decimal], n: usize) -> Option<Decimal> {
    if n > returns.len() || n == 0 {
        return None;
    }
    let num_subseries = returns.len() / n;
    if num_subseries == 0 {
        return None;
    }

    let mut rs_sum = Decimal::ZERO;
    let mut valid_count = 0u32;

    for k in 0..num_subseries {
        let sub = &returns[k * n..(k + 1) * n];
        let n_dec = Decimal::from(sub.len());
        let mean: Decimal = sub.iter().copied().sum::<Decimal>() / n_dec;
        let deviations: Vec<Decimal> = sub.iter().map(|r| r - mean).collect();

        // Cumulative sum of deviations
        let mut cum = Vec::with_capacity(sub.len());
        let mut running = Decimal::ZERO;
        for d in &deviations {
            running += d;
            cum.push(running);
        }

        let r = cum.iter().copied().max().unwrap_or(Decimal::ZERO)
            - cum.iter().copied().min().unwrap_or(Decimal::ZERO);

        // Standard deviation
        let var: Decimal = deviations.iter().map(|d| d * d).sum::<Decimal>() / n_dec;
        if var.is_zero() {
            continue;
        }
        let s = newton_sqrt(var);
        if s.is_zero() {
            continue;
        }

        rs_sum += r / s;
        valid_count += 1;
    }

    if valid_count == 0 {
        return None;
    }
    Some(rs_sum / Decimal::from(valid_count))
}

/// Estimate the Hurst exponent H from the log-log slope of R/S vs subseries size.
///
/// Uses the first and last valid (n, avg_rs) points to compute:
/// H = log2(RS_last / RS_first) / log2(n_last / n_first)
///
/// Returns `None` if fewer than 2 valid data points.
pub fn estimate_slope(log_ns: &[usize], log_rs: &[Decimal]) -> Option<Decimal> {
    if log_ns.len() < 2 || log_rs.len() < 2 {
        return None;
    }

    let rs1 = log_rs[0];
    let rs2 = log_rs[log_rs.len() - 1];
    let n1 = Decimal::from(log_ns[0]);
    let n2 = Decimal::from(log_ns[log_ns.len() - 1]);

    if rs1.is_zero() || n1.is_zero() {
        return None;
    }

    let rs_ratio = rs2 / rs1;
    let n_ratio = n2 / n1;
    let h = approx_log2(rs_ratio) / approx_log2(n_ratio);
    Some(h.max(Decimal::ZERO).min(Decimal::ONE))
}

/// Full Hurst estimation pipeline from a return series.
fn estimate_hurst_from_returns(returns: &[Decimal]) -> Option<Decimal> {
    let mut log_ns = Vec::new();
    let mut log_rs = Vec::new();

    for &n in &SUBSERIES_SIZES {
        if let Some(avg_rs) = rescaled_range_for_n(returns, n) {
            log_ns.push(n);
            log_rs.push(avg_rs);
        }
    }

    estimate_slope(&log_ns, &log_rs)
}

/// Newton-Raphson square root for Decimal.
fn newton_sqrt(val: Decimal) -> Decimal {
    if val.is_zero() || val < Decimal::ZERO {
        return Decimal::ZERO;
    }
    let mut guess = val / Decimal::TWO;
    for _ in 0..10 {
        if guess.is_zero() {
            return Decimal::ZERO;
        }
        guess = (guess + val / guess) / Decimal::TWO;
    }
    guess
}

/// Approximate log2(x) for positive x using series expansion.
fn approx_log2(x: Decimal) -> Decimal {
    if x <= Decimal::ZERO {
        return Decimal::ZERO;
    }
    let ln2 = Decimal::new(6931, 4); // 0.6931

    let mut reduced = x;
    let mut shifts = 0i64;
    while reduced >= Decimal::TWO {
        reduced /= Decimal::TWO;
        shifts += 1;
    }
    while reduced < Decimal::ONE && !reduced.is_zero() {
        reduced *= Decimal::TWO;
        shifts -= 1;
    }

    let u = (reduced - Decimal::ONE) / (reduced + Decimal::ONE);
    let u2 = u * u;
    let mut term = u;
    let mut sum = u;
    for k in 1..10 {
        term *= u2;
        sum += term / Decimal::from(2 * k + 1);
    }
    let ln_reduced = sum * Decimal::TWO;
    let ln_x = ln_reduced + Decimal::from(shifts) * ln2;
    ln_x / ln2
}

#[cfg(test)]
#[path = "hurst_tests.rs"]
mod tests;