foxhunt/crates/ml-features/src/bouchaud_features.rs

//! V10 Block BB — Bouchaud-line trade-flow features (light-weight proxies).
//!
//! Two scalars from the trade tape that capture self-excitation / persistence
//! without requiring full Hawkes MLE:
//!
//! ## 1. Trade-sign autocorrelation at lag 1 (Lillo-Mike-Farmer 2005)
//!
//! `ρ(1) = corr(q_t, q_{t-1})` where `q_t ∈ {-1, +1}` is the sign of trade t
//! (buyer-initiated = +1, seller-initiated = -1). Empirically `ρ(1)` runs
//! ~0.3-0.6 in equity markets: trade-sign sequences are far from i.i.d.,
//! exhibiting strong short-term persistence. A bar where `ρ(1)` is unusually
//! high (vs the trailing average) indicates **directional momentum regime**;
//! near-zero values indicate **mean-reverting / dispersed flow**.
//!
//! This is the *fast* version of the Bouchaud γ exponent — a single scalar
//! at one lag rather than a fitted power-law decay.
//!
//! ## 2. Branching ratio proxy via Fano factor (Filimonov-Sornette 2012)
//!
//! For trade arrival times, compute the **Fano factor** `F(T) = Var(N(T)) /
//! E(N(T))` over windows of size T (where N(T) = count of trades in T).
//!
//! - Poisson independent arrivals: `F(T) = 1`
//! - Self-exciting Hawkes: `F(T) > 1`, growing with T
//! - In the limit, `F(T) → 1 / (1 − n)²` where `n` is the branching ratio
//!   (Filimonov-Sornette 2012, *PNAS*)
//!
//! We emit `branching_ratio_proxy = F(T) − 1`, which is:
//! - 0 for independent arrivals
//! - Strictly positive and growing for self-exciting regimes
//! - Bounded for finite-window estimates
//!
//! This is the *fast* version of the Hawkes branching ratio n (full MLE in
//! Block R provides the actual `n = α/β` estimate).
//!
//! ## References
//!
//! - Lillo, Mike & Farmer (2005), "Theory for long memory in supply and demand"
//! - Filimonov & Sornette (2012), "Quantifying reflexivity in financial markets",
//!   *Phys. Rev. E* (also PNAS 2015)
//! - Bouchaud, Bonart, Donier & Gould (2018), *Trades, Quotes and Prices*, Ch. 10

use std::collections::VecDeque;

/// V10 Block BB — streaming Bouchaud trade-flow proxies.
///
/// Maintains:
/// - A rolling buffer of recent trade signs (for autocorrelation)
/// - A rolling buffer of recent trade timestamps (for Fano factor)
#[derive(Debug, Clone)]
pub struct BouchaudFeatures {
    /// Window for autocorrelation estimation (trade signs).
    sign_window: usize,
    /// Window for Fano factor estimation (in trades).
    fano_window: usize,
    /// Fano factor sub-window size T (typically 10-50 trades).
    fano_sub_window: usize,

    signs: VecDeque<f64>,
    arrival_ts_ns: VecDeque<u64>,
}

impl BouchaudFeatures {
    pub fn new(sign_window: usize, fano_window: usize, fano_sub_window: usize) -> Self {
        Self {
            sign_window,
            fano_window,
            fano_sub_window: fano_sub_window.max(2),
            signs: VecDeque::with_capacity(sign_window),
            arrival_ts_ns: VecDeque::with_capacity(fano_window),
        }
    }

    /// Defaults: 200-trade autocorr window, 500-trade Fano buffer, 20-trade
    /// sub-windows. Calibrated for HFT futures (high trade rate).
    pub fn with_defaults() -> Self {
        Self::new(200, 500, 20)
    }

    /// Update with a new trade event.
    pub fn update(&mut self, is_buy: bool, timestamp_ns: u64) {
        let sign: f64 = if is_buy { 1.0 } else { -1.0 };
        self.signs.push_back(sign);
        if self.signs.len() > self.sign_window {
            self.signs.pop_front();
        }
        self.arrival_ts_ns.push_back(timestamp_ns);
        if self.arrival_ts_ns.len() > self.fano_window {
            self.arrival_ts_ns.pop_front();
        }
    }

    /// Returns `[trade_sign_autocorr_lag1, branching_ratio_proxy]`.
    /// Both default to 0 when insufficient data.
    pub fn features(&self) -> [f64; 2] {
        [self.trade_sign_autocorr_lag1(), self.branching_ratio_proxy()]
    }

    /// Lag-1 autocorrelation of trade signs in the rolling window.
    /// Returns 0 for fewer than 10 observations.
    pub fn trade_sign_autocorr_lag1(&self) -> f64 {
        if self.signs.len() < 10 {
            return 0.0;
        }
        let n = self.signs.len() as f64;
        let mean: f64 = self.signs.iter().sum::<f64>() / n;
        let mut numer = 0.0_f64;
        let mut denom = 0.0_f64;
        for i in 1..self.signs.len() {
            let s_t = self.signs[i] - mean;
            let s_tm1 = self.signs[i - 1] - mean;
            numer += s_t * s_tm1;
            denom += s_t * s_t;
        }
        if denom < 1e-12 {
            return 0.0;
        }
        numer / denom
    }

    /// Fano-factor-based branching-ratio proxy: `F(T) − 1` where `T` is the
    /// sub-window size in trades. Bins the arrivals into non-overlapping
    /// sub-windows of `fano_sub_window` trades, computes `Var(N) / E(N)`.
    ///
    /// Returns 0 if there are fewer than `3 × fano_sub_window` trades observed.
    pub fn branching_ratio_proxy(&self) -> f64 {
        let n_arrivals = self.arrival_ts_ns.len();
        if n_arrivals < 3 * self.fano_sub_window {
            return 0.0;
        }

        // Determine total observation window (last ts − first ts in nanoseconds)
        // We partition this duration into K equal sub-windows and count arrivals
        // per sub-window. K = floor(n_arrivals / fano_sub_window).
        let first_ts = *self.arrival_ts_ns.front().expect("non-empty");
        let last_ts = *self.arrival_ts_ns.back().expect("non-empty");
        if last_ts <= first_ts {
            return 0.0;
        }
        let total_span_ns = last_ts - first_ts;
        let n_sub_windows = (n_arrivals / self.fano_sub_window).max(2);
        let sub_window_ns = total_span_ns / n_sub_windows as u64;
        if sub_window_ns == 0 {
            return 0.0;
        }

        // Count arrivals per sub-window
        let mut counts: Vec<u32> = vec![0; n_sub_windows];
        for &ts in self.arrival_ts_ns.iter() {
            let bucket = ((ts - first_ts) / sub_window_ns) as usize;
            let bucket = bucket.min(n_sub_windows - 1);
            counts[bucket] += 1;
        }

        // Compute mean + variance of count distribution
        let k = counts.len() as f64;
        let mean: f64 = counts.iter().map(|&c| c as f64).sum::<f64>() / k;
        if mean < 1e-9 {
            return 0.0;
        }
        let var: f64 =
            counts.iter().map(|&c| (c as f64 - mean).powi(2)).sum::<f64>() / k;
        // Fano − 1 ≈ 0 for Poisson; > 0 for self-exciting
        (var / mean - 1.0).max(-1.0).min(50.0) // sanity-bound for numerical stability
    }

    pub fn reset(&mut self) {
        self.signs.clear();
        self.arrival_ts_ns.clear();
    }
}

impl Default for BouchaudFeatures {
    fn default() -> Self {
        Self::with_defaults()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_cold_start_features_zero() {
        let bf = BouchaudFeatures::with_defaults();
        assert_eq!(bf.features(), [0.0, 0.0]);
    }

    #[test]
    fn test_autocorr_lag1_positive_for_persistent_signs() {
        // All same sign → perfect lag-1 autocorrelation
        let mut bf = BouchaudFeatures::new(100, 500, 20);
        for i in 0..50 {
            bf.update(true, (i as u64) * 1_000_000);
        }
        let [ac, _] = bf.features();
        // Variance is zero (all +1) → numer/denom = 0/0 guarded → 0
        // To get nonzero ac we need variance. Let me check the implementation.
        // With all +1: mean=1, deviations all = 0, denom = 0 → returns 0.
        assert_eq!(ac, 0.0, "all-same signs have zero variance → ac=0 by guard");
    }

    #[test]
    fn test_autocorr_lag1_positive_for_persistent_with_variance() {
        // 10× +1 followed by 10× -1 followed by 10× +1, etc. → high persistence
        let mut bf = BouchaudFeatures::new(200, 500, 20);
        let mut t = 0_u64;
        for outer in 0..10 {
            for _ in 0..10 {
                bf.update(outer % 2 == 0, t);
                t += 1_000_000;
            }
        }
        let [ac, _] = bf.features();
        // Same-sign runs → lag-1 ac should be near 1 except at run boundaries
        // 100 trades, 10 sign changes → 90/99 = ~0.91 ac
        assert!(ac > 0.5, "persistent runs → high lag-1 ac, got {ac}");
    }

    #[test]
    fn test_autocorr_lag1_negative_for_alternating_signs() {
        // +1, -1, +1, -1, ... → lag-1 ac → -1
        let mut bf = BouchaudFeatures::new(200, 500, 20);
        for i in 0..100 {
            bf.update(i % 2 == 0, (i as u64) * 1_000_000);
        }
        let [ac, _] = bf.features();
        assert!(ac < -0.9, "alternating signs → ac ≈ -1, got {ac}");
    }

    #[test]
    fn test_branching_ratio_proxy_zero_for_regular_arrivals() {
        // Perfectly regular arrivals: each trade exactly Δt apart → variance
        // of count per sub-window is 0 → Fano = 0 → proxy = -1
        let mut bf = BouchaudFeatures::new(200, 500, 20);
        for i in 0..100 {
            bf.update(i % 2 == 0, (i as u64) * 1_000_000);
        }
        let [_, br] = bf.features();
        // Regular arrivals: each sub-window of 20 trades takes exactly 20Δt
        // → variance ~0 → proxy ~ -1 (less than Poisson)
        assert!(br < 0.0, "regular arrivals are sub-Poisson, got {br}");
    }

    #[test]
    fn test_branching_ratio_proxy_positive_for_clustered_arrivals() {
        // Heavily clustered: 80 trades in first 10% of time, 20 in remaining 90%
        let mut bf = BouchaudFeatures::new(200, 500, 20);
        // Cluster 1: 80 trades quickly
        for i in 0..80 {
            bf.update(i % 2 == 0, (i as u64) * 1_000_000);
        }
        // Sparse: 20 trades spread out
        for i in 0..20 {
            let ts = 80_000_000 + (i as u64) * 100_000_000;
            bf.update(i % 2 == 0, ts);
        }
        let [_, br] = bf.features();
        assert!(br > 0.0, "clustered arrivals are super-Poisson, got {br}");
    }

    #[test]
    fn test_reset_clears_features() {
        let mut bf = BouchaudFeatures::with_defaults();
        for i in 0..50 {
            bf.update(i % 2 == 0, (i as u64) * 1_000_000);
        }
        assert_ne!(bf.features()[0], 0.0);
        bf.reset();
        assert_eq!(bf.features(), [0.0, 0.0]);
    }
}