market_model/process/

hawkes.rs

//! Hawkes self-exciting point process.
//!
//! A stochastic intensity process where each arrival increases the
//! probability of further arrivals, capturing clustering behavior
//! observed in market orders, cancellations, and volatility spikes.

use crate::process::StochasticProcess;
use crate::simulatable::Simulatable;
use rand::Rng;

/// Trait for Hawkes kernel functions.
pub trait Kernel: Clone + Send + Sync {
    fn call(&self, t: f64) -> f64;

    /// Whether this kernel admits an O(1) recursive intensity update.
    /// Only exponential kernels (phi(t) = alpha * exp(-beta * t)) support this.
    fn supports_recursive_intensity(&self) -> bool {
        false
    }

    /// Decay rate beta such that intensity decays as exp(-beta * dt) toward mu.
    /// Only meaningful when `supports_recursive_intensity()` is true.
    fn decay_rate(&self) -> f64 {
        0.0
    }

    /// Excitation jump size added to intensity on each arrival.
    /// Only meaningful when `supports_recursive_intensity()` is true.
    fn excitation_size(&self) -> f64 {
        0.0
    }
}

/// Exponential kernel: phi(t) = alpha * exp(-beta * t).
#[derive(Clone)]
pub struct ExponentialKernel {
    pub alpha: f64,
    pub beta: f64,
}

impl ExponentialKernel {
    pub fn new(alpha: f64, beta: f64) -> Self {
        Self { alpha, beta }
    }
}

impl Kernel for ExponentialKernel {
    fn call(&self, t: f64) -> f64 {
        if t < 0.0 {
            0.0
        } else {
            self.alpha * (-self.beta * t).exp()
        }
    }

    fn supports_recursive_intensity(&self) -> bool {
        true
    }

    fn decay_rate(&self) -> f64 {
        self.beta
    }

    fn excitation_size(&self) -> f64 {
        self.alpha
    }
}

/// Power-law kernel: phi(t) = alpha / (c + t)^p.
#[derive(Clone)]
pub struct PowerLawKernel {
    pub alpha: f64,
    pub c: f64,
    pub p: f64,
}

impl PowerLawKernel {
    pub fn new(alpha: f64, c: f64, p: f64) -> Self {
        Self { alpha, c, p }
    }
}

impl Kernel for PowerLawKernel {
    fn call(&self, t: f64) -> f64 {
        if t < 0.0 {
            0.0
        } else {
            self.alpha / (self.c + t).powf(self.p)
        }
    }
}

/// Hawkes process with arbitrary kernel.
///
/// Simulated via Ogata's thinning algorithm. For exponential kernels,
/// intensity is maintained O(1) via the recursion
/// I(t+dt) = mu + (I(t) - mu) * exp(-beta*dt) + alpha * new_jumps.
/// For non-exponential kernels, intensity is computed by summing over all
/// historical jumps (O(n)).
#[derive(Clone)]
pub struct HawkesProcess<K: Kernel> {
    pub mu: f64,
    pub kernel: K,
    jumps: Vec<f64>,
    current_time: f64,
    state: f64,
    /// Pre-computed decay factor per unit time; used only by recursive path.
    /// Re-computed when `current_time` advances.
    current_intensity: f64,
    /// Time of the oldest jump retained (for pruning).
    prune_threshold: f64,
}

impl<K: Kernel> HawkesProcess<K> {
    pub fn new(mu: f64, kernel: K) -> Self {
        Self {
            mu,
            kernel,
            jumps: Vec::new(),
            current_time: 0.0,
            state: mu,
            current_intensity: mu,
            prune_threshold: f64::NEG_INFINITY,
        }
    }

    pub fn value(&self) -> f64 {
        self.state
    }

    fn compute_intensity(&self, t: f64) -> f64 {
        if self.kernel.supports_recursive_intensity() {
            // O(1) for exponential kernel: compute directly from jump times.
            // The contribution of each jump at tj is alpha * exp(-beta * (t - tj)).
            // Since we only ever need the O(n) scan once per step() call (at the
            // very end), this is simpler and numerically identical to the
            // recursive O(1) path inside the Ogata loop.
            let alpha = self.kernel.excitation_size();
            let beta = self.kernel.decay_rate();
            let mut intensity = self.mu;
            for &tj in &self.jumps {
                if tj < t {
                    intensity += alpha * (-beta * (t - tj)).exp();
                } else {
                    break;
                }
            }
            intensity
        } else {
            // O(n): sum over all historical jumps using the kernel call.
            let mut sum = 0.0;
            for &tj in &self.jumps {
                if tj < t {
                    sum += self.kernel.call(t - tj);
                } else {
                    break;
                }
            }
            self.mu + sum
        }
    }

    /// Drop jumps older than ~10 decay half-lives.
    ///
    /// A jump at time t_j contributes alpha * exp(-beta * (t - t_j)) to
    /// the intensity at time t. After 10 half-lives, exp(-10) ~ 4.5e-5.
    /// Removing such jumps has negligible impact on the intensity but bounds
    /// the Vec growth from O(steps) to O(1).
    fn prune_jumps(&mut self) {
        let cutoff = self.current_time - 10.0 / self.kernel.decay_rate();
        if cutoff <= self.prune_threshold {
            return;
        }
        self.prune_threshold = cutoff;
        let first = self.jumps.partition_point(|tj| *tj < cutoff);
        if first > 0 {
            self.jumps.drain(..first);
        }
    }
}

impl<K: Kernel + 'static> Simulatable for HawkesProcess<K> {
    type Output = f64;

    fn dim(&self) -> usize {
        0 // No Brownian driver; all randomness from Ogata thinning via rng.
    }

    fn step(&mut self, dt: f64, _dw: &[f64], rng: &mut impl Rng) -> f64 {
        let end_time = self.current_time + dt;

        if self.kernel.supports_recursive_intensity() {
            // O(1) path for exponential kernel.
            //
            // Two intensity values are maintained:
            //   lambda_t: true intensity I(t), including all accepted jumps.
            //   lambda_p: proposal upper bound for Ogata, which after a jump
            //     at t does NOT include the jump (matches compute_intensity(t)
            //     strict < semantics in the O(n) path).
            //
            // After a jump: lambda_t += alpha, lambda_p stays at pre-jump.
            // Both decay at rate beta toward mu between events.
            let beta = self.kernel.decay_rate();
            let alpha = self.kernel.excitation_size();
            let mut lambda_t = self.current_intensity;
            let mut lambda_p = self.current_intensity;
            let mut t_now = self.current_time;

            while t_now < end_time {
                if lambda_p <= 1e-12 {
                    break;
                }

                let u: f64 = rng.random();
                let tau = -u.ln() / lambda_p;
                let candidate = t_now + tau;

                if candidate > end_time {
                    // Decay both to end_time.
                    let dtau = end_time - t_now;
                    let decay = (-beta * dtau).exp();
                    lambda_t = self.mu + (lambda_t - self.mu) * decay;
                    t_now = end_time;
                    break;
                }

                let dtau = candidate - t_now;
                let decay = (-beta * dtau).exp();
                let lambda_t_candidate = self.mu + (lambda_t - self.mu) * decay;
                let lambda_p_candidate = self.mu + (lambda_p - self.mu) * decay;

                let u_accept: f64 = rng.random();
                // Ogata acceptance uses lambda_p (the proposal upper bound)
                // against lambda_t_candidate (the true intensity at candidate).
                if u_accept * lambda_p <= lambda_t_candidate {
                    self.jumps.push(candidate);
                    lambda_t = lambda_t_candidate + alpha;
                    // lambda_p for next interval: intensity at candidate
                    // WITHOUT the just-accepted jump (matches O(n)).
                    lambda_p = lambda_t_candidate;
                } else {
                    lambda_t = lambda_t_candidate;
                    lambda_p = lambda_p_candidate;
                }
                t_now = candidate;
            }

            self.current_time = t_now;
            self.state = lambda_t;
            self.current_intensity = lambda_t;
            self.prune_jumps();
        } else {
            // O(n) path: Ogata thinning with full intensity scan.
            while self.current_time < end_time {
                let lambda_star = self.compute_intensity(self.current_time);

                if lambda_star <= 1e-12 {
                    self.current_time = end_time;
                    break;
                }

                let u: f64 = rng.random();
                let tau = -u.ln() / lambda_star;
                let candidate = self.current_time + tau;

                if candidate > end_time {
                    self.current_time = end_time;
                    break;
                }

                self.current_time = candidate;
                let lambda_s = self.compute_intensity(self.current_time);
                let u_accept: f64 = rng.random();

                if u_accept * lambda_star <= lambda_s {
                    self.jumps.push(self.current_time);
                }
            }

            self.state = self.compute_intensity(self.current_time);
            self.current_intensity = self.state;
        }

        self.state
    }

    fn current(&self) -> f64 {
        self.state
    }
}

/// Simulate the Hawkes process and return the intensity path on a grid.
impl<K: Kernel + 'static> HawkesProcess<K> {
    pub fn simulate_intensity(&self, total_time: f64, n_steps: usize) -> Vec<f64> {
        let mut sim = self.clone();
        let dt = total_time / n_steps as f64;
        let mut path = Vec::with_capacity(n_steps + 1);
        path.push(sim.current());
        for _ in 0..n_steps {
            let v = Simulatable::step(&mut sim, dt, &[], &mut rand::rng());
            path.push(v);
        }
        path
    }
}

impl<K: Kernel + 'static> StochasticProcess for HawkesProcess<K> {
    type State = f64;
    type Diffusion = f64;
    type BrownianIncrement = f64;

    fn drift(&self, x: &f64, _t: f64) -> f64 {
        self.mu - x
    }

    fn diffusion(&self, _x: &f64, _t: f64) -> f64 {
        0.0
    }

    fn jump_intensity(&self, _x: &f64, _t: f64) -> f64 {
        self.mu
    }

    fn jump_size<R: Rng + ?Sized>(&self, _x: &f64, t: f64, _rng: &mut R) -> f64 {
        self.kernel.call(t)
    }

    fn expected_value(&self, x0: &f64, _t: f64) -> f64 {
        *x0
    }

    fn variance(&self, _x0: &f64, _t: f64) -> f64 {
        self.mu
    }

    fn step<R: Rng + ?Sized>(&self, x: &f64, _t: f64, dt: f64, _dw: &f64, rng: &mut R) -> f64 {
        let lambda = self.mu;
        if lambda <= 1e-12 {
            return *x;
        }
        let u: f64 = rng.random();
        let tau = -u.ln() / lambda;
        if tau < dt { *x + 1.0 } else { *x }
    }

    fn generate_increment<R: Rng + ?Sized>(&self, _rng: &mut R, _dt: f64) -> f64 {
        0.0
    }

    fn simulate<R: Rng + ?Sized>(
        &self,
        x0: f64,
        t0: f64,
        t_end: f64,
        dt: f64,
        rng: &mut R,
    ) -> Vec<(f64, f64)> {
        let n_steps = ((t_end - t0) / dt).ceil() as usize;
        let mut path = Vec::with_capacity(n_steps + 1);
        let mut t = t0;
        let mut x = x0;
        path.push((t, x));

        for _ in 0..n_steps {
            x = self.step(&x, t, dt, &0.0, rng);
            t += dt;
            path.push((t, x));
        }
        path
    }
}

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

    #[test]
    fn test_hawkes_no_brownian_dimension() {
        let kernel = ExponentialKernel::new(0.5, 2.0);
        let h = HawkesProcess::new(0.1, kernel);
        assert_eq!(h.dim(), 0);
    }

    #[test]
    fn test_hawkes_intensity_positive() {
        let kernel = ExponentialKernel::new(0.5, 2.0);
        let h = HawkesProcess::new(0.5, kernel);
        let path = h.simulate_intensity(1.0, 100);
        for v in path {
            assert!(v >= 0.0);
        }
    }

    #[test]
    fn test_exponential_kernel() {
        let k = ExponentialKernel::new(1.0, 2.0);
        assert_eq!(k.call(0.0), 1.0);
        assert!(k.call(1.0) < 1.0);
        assert_eq!(k.call(-0.5), 0.0);
    }

    #[test]
    fn test_power_law_kernel() {
        let k = PowerLawKernel::new(1.0, 1.0, 2.0);
        assert_eq!(k.call(0.0), 1.0);
        assert!(k.call(1.0) < 0.5);
        assert_eq!(k.call(-0.5), 0.0);
    }

    #[test]
    fn test_recursive_intensity_same_seed_match() {
        // Same seed trace: the recursive path must consume RNG identically.
        // We run a long single-path trace with StdRng.  If the two paths ever
        // diverge, the RNGs desync and all subsequent steps differ.
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        let kernel = ExponentialKernel::new(2.0, 3.0);
        let mu = 1.0;
        let dt = 0.01;
        let n_steps = 5000;
        let seed = 42u64;

        #[derive(Clone)]
        struct NonRecursiveKernel(ExponentialKernel);
        impl Kernel for NonRecursiveKernel {
            fn call(&self, t: f64) -> f64 {
                self.0.call(t)
            }
        }

        let mut ha = HawkesProcess::new(mu, kernel.clone());
        let mut rng_a = StdRng::seed_from_u64(seed);
        let mut hb = HawkesProcess::new(mu, NonRecursiveKernel(kernel.clone()));
        let mut rng_b = StdRng::seed_from_u64(seed);

        let mut first_diverge = None;
        for i in 0..n_steps {
            let a = Simulatable::step(&mut ha, dt, &[], &mut rng_a);
            let b = Simulatable::step(&mut hb, dt, &[], &mut rng_b);
            if first_diverge.is_none() && (a - b).abs() > 1e-12 {
                first_diverge = Some((i, a, b));
            }
        }

        if let Some((i, a, b)) = first_diverge {
            panic!(
                "recursive and exact diverged at step {}: rec={:.15e}, exact={:.15e}",
                i, a, b
            );
        }
    }

    #[test]
    fn test_recursive_intensity_o_n_self_consistency() {
        // Two independent runs of the O(n) path: grand means should agree
        // within Monte Carlo error.  This establishes the baseline noise floor
        // for any cross-algorithm comparison.
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        let kernel = ExponentialKernel::new(2.0, 3.0);
        let mu = 1.0;
        let dt = 0.01;
        let n_steps = 2000;
        let n_paths = 1000;

        #[derive(Clone)]
        struct NonRecursiveKernel(ExponentialKernel);
        impl Kernel for NonRecursiveKernel {
            fn call(&self, t: f64) -> f64 {
                self.0.call(t)
            }
        }

        let mut sum_a = 0.0f64;
        let mut sum_b = 0.0f64;
        let mut count = 0u64;

        for path_i in 0u64..n_paths {
            let seed_a = path_i;
            let seed_b = path_i + 1_000_000;

            let mut ha = HawkesProcess::new(mu, NonRecursiveKernel(kernel.clone()));
            let mut rng_a = StdRng::seed_from_u64(seed_a);
            for _ in 0..n_steps {
                sum_a += Simulatable::step(&mut ha, dt, &[], &mut rng_a);
                count += 1;
            }

            let mut hb = HawkesProcess::new(mu, NonRecursiveKernel(kernel.clone()));
            let mut rng_b = StdRng::seed_from_u64(seed_b);
            for _ in 0..n_steps {
                sum_b += Simulatable::step(&mut hb, dt, &[], &mut rng_b);
            }
        }

        let mean_a = sum_a / count as f64;
        let mean_b = sum_b / count as f64;
        let rel_diff = (mean_a - mean_b).abs() / mean_b.max(1e-12);

        eprintln!(
            "O(n) self-consistency: mean_a={:.6}, mean_b={:.6}, rel_diff={:.4}%",
            mean_a,
            mean_b,
            rel_diff * 100.0
        );

        // The exponential Hawkes with branching ratio 2/3 has long-range
        // dependence (bursty).  1000 paths x 2000 steps may not be enough
        // for the grand mean to converge tighter than ~1%.  This is a
        // property of the process, not of the algorithm.
        assert!(
            rel_diff < 0.01,
            "O(n) self-consistency: means differ by {:.4}%",
            rel_diff * 100.0
        );
    }
}