diff --git a/src/lib.rs b/src/lib.rs index b7b71fc..7f1112e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -45,6 +45,7 @@ pub mod optimal_control; pub mod risk_metrics; pub mod sparse_optimization; pub mod mean_field; // Mean Field Games and Mean Field Type Control +pub mod point_processes; // Point processes for order flow modeling (Hawkes, fBM) // Python bindings for legacy compatibility #[cfg(feature = "python-bindings")] @@ -115,6 +116,9 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> { // Mean Field Games functions mean_field::python_bindings::register_python_functions(m)?; + // Point Processes functions (Hawkes, fBM, order flow) + point_processes::python_bindings::register_python_functions(m)?; + // Optimal Control functions (includes Kalman Filter) optimal_control::py_bindings::register_py_module(m)?; diff --git a/src/point_processes/hawkes.rs b/src/point_processes/hawkes.rs new file mode 100644 index 0000000..6c028fb --- /dev/null +++ b/src/point_processes/hawkes.rs @@ -0,0 +1,462 @@ +//! Hawkes Process Implementation +//! +//! Self-exciting point processes where past events increase the probability +//! of future events. Used to model order flow clustering in financial markets. + +use super::kernels::ExcitationKernel; +use rand::prelude::*; +use rand_distr::{Exp, Uniform}; +use std::collections::VecDeque; + +/// Configuration for a Hawkes process +#[derive(Clone, Debug)] +pub struct HawkesProcessConfig { + /// Baseline intensity ν > 0 + pub baseline_intensity: f64, + /// Excitation kernel + pub kernel: K, + /// Maximum history to track (for efficiency) + pub max_history: usize, + /// Numerical tolerance + pub tolerance: f64, +} + +impl HawkesProcessConfig { + pub fn new(baseline_intensity: f64, kernel: K) -> Self { + assert!(baseline_intensity > 0.0); + Self { + baseline_intensity, + kernel, + max_history: 10000, + tolerance: 1e-10, + } + } + + pub fn with_max_history(mut self, max_history: usize) -> Self { + self.max_history = max_history; + self + } +} + +/// Univariate Hawkes process +/// +/// Models self-exciting point processes where the intensity is: +/// λ(t) = ν + ∫₀ᵗ⁻ φ(t - s) dN(s) +/// +/// where ν is the baseline intensity and φ is the excitation kernel. +#[derive(Clone, Debug)] +pub struct HawkesProcess { + config: HawkesProcessConfig, + /// Event times + events: Vec, + /// Current time + current_time: f64, +} + +impl HawkesProcess { + /// Create a new Hawkes process + pub fn new(baseline_intensity: f64, kernel: K) -> Self { + let config = HawkesProcessConfig::new(baseline_intensity, kernel); + Self { + config, + events: Vec::new(), + current_time: 0.0, + } + } + + /// Create from configuration + pub fn from_config(config: HawkesProcessConfig) -> Self { + Self { + config, + events: Vec::new(), + current_time: 0.0, + } + } + + /// Get the intensity at time t + pub fn intensity(&self, t: f64) -> f64 { + let mut lambda = self.config.baseline_intensity; + for &event_time in &self.events { + if event_time < t { + lambda += self.config.kernel.evaluate(t - event_time); + } + } + lambda + } + + /// Simulate the process up to time T + pub fn simulate(&mut self, t_max: f64, seed: Option) -> Vec { + let mut rng = match seed { + Some(s) => StdRng::seed_from_u64(s), + None => StdRng::from_entropy(), + }; + + self.events.clear(); + self.current_time = 0.0; + + // Use Ogata's thinning algorithm + let lambda_max_factor = 1.5; // Safety factor for upper bound + + while self.current_time < t_max { + // Upper bound for intensity + let lambda_max = self.intensity(self.current_time) * lambda_max_factor + + self.config.baseline_intensity; + + // Generate candidate inter-arrival time + let exp_dist = Exp::new(lambda_max).unwrap(); + let tau: f64 = rng.sample(exp_dist); + let candidate_time = self.current_time + tau; + + if candidate_time >= t_max { + break; + } + + // Accept/reject + let u: f64 = rng.gen(); + let lambda_at_candidate = self.intensity(candidate_time); + + if u <= lambda_at_candidate / lambda_max { + // Accept the event + self.events.push(candidate_time); + self.current_time = candidate_time; + + // Prune old events if needed + if self.events.len() > self.config.max_history { + let recent_events: Vec = self.events + .iter() + .rev() + .take(self.config.max_history) + .cloned() + .collect(); + self.events = recent_events.into_iter().rev().collect(); + } + } else { + self.current_time = candidate_time; + } + } + + self.events.clone() + } + + /// Estimate kernel parameters from event times using MLE + pub fn fit(&mut self, events: &[f64]) -> Result { + if events.is_empty() { + return Err("No events provided"); + } + + self.events = events.to_vec(); + let t_max = *events.last().unwrap(); + let n = events.len(); + + // Compute log-likelihood + let log_likelihood = self.log_likelihood(t_max); + + // Estimate baseline intensity + let lambda_avg = n as f64 / t_max; + let branching_ratio = self.config.kernel.l1_norm(); + let estimated_baseline = lambda_avg * (1.0 - branching_ratio); + + Ok(FitResult { + log_likelihood, + n_events: n, + duration: t_max, + estimated_baseline, + estimated_branching_ratio: branching_ratio, + }) + } + + /// Compute log-likelihood of the process + pub fn log_likelihood(&self, t_max: f64) -> f64 { + if self.events.is_empty() { + return 0.0; + } + + // Log-likelihood = Σ log(λ(tᵢ)) - ∫₀ᵀ λ(t) dt + let mut ll = 0.0; + + // Sum of log intensities at event times + for &t in &self.events { + let lambda_t = self.intensity(t); + if lambda_t > self.config.tolerance { + ll += lambda_t.ln(); + } + } + + // Compensator (integrated intensity) + // ∫₀ᵀ λ(t) dt = νT + Σᵢ ∫₀^{T-tᵢ} φ(s) ds + let compensator = self.config.baseline_intensity * t_max + + self.events.iter() + .map(|&ti| self.config.kernel.integrate(t_max - ti)) + .sum::(); + + ll - compensator + } + + /// Get all event times + pub fn events(&self) -> &[f64] { + &self.events + } + + /// Get number of events + pub fn n_events(&self) -> usize { + self.events.len() + } +} + +/// Result of fitting a Hawkes process +#[derive(Clone, Debug)] +pub struct FitResult { + pub log_likelihood: f64, + pub n_events: usize, + pub duration: f64, + pub estimated_baseline: f64, + pub estimated_branching_ratio: f64, +} + +/// Bivariate Hawkes process for reaction order flow +/// +/// Models the interplay between buy and sell reaction orders. +/// N = (N⁺, N⁻) with intensity: +/// +/// λ⁺(t) = μ⁺(t) + ∫ [φ₁(t-s)dN⁺(s) + φ₂(t-s)dN⁻(s)] +/// λ⁻(t) = μ⁻(t) + ∫ [φ₂(t-s)dN⁺(s) + φ₁(t-s)dN⁻(s)] +/// +/// where μ(t) is driven by the core order flow. +#[derive(Clone, Debug)] +pub struct BivariateHawkes { + /// Same-side kernel φ₁ (buy->buy, sell->sell) + pub phi_1: K, + /// Cross-side kernel φ₂ (buy->sell, sell->buy) + pub phi_2: K, + /// Buy events N⁺ + pub buy_events: Vec, + /// Sell events N⁻ + pub sell_events: Vec, + /// Current time + current_time: f64, +} + +impl BivariateHawkes { + pub fn new(phi_1: K, phi_2: K) -> Self { + Self { + phi_1, + phi_2, + buy_events: Vec::new(), + sell_events: Vec::new(), + current_time: 0.0, + } + } + + /// Get buy intensity at time t given external baseline μ⁺(t) + pub fn buy_intensity(&self, t: f64, baseline: f64) -> f64 { + let mut lambda = baseline; + + // Self-excitation from buy events + for &ti in &self.buy_events { + if ti < t { + lambda += self.phi_1.evaluate(t - ti); + } + } + + // Cross-excitation from sell events + for &ti in &self.sell_events { + if ti < t { + lambda += self.phi_2.evaluate(t - ti); + } + } + + lambda + } + + /// Get sell intensity at time t given external baseline μ⁻(t) + pub fn sell_intensity(&self, t: f64, baseline: f64) -> f64 { + let mut lambda = baseline; + + // Cross-excitation from buy events + for &ti in &self.buy_events { + if ti < t { + lambda += self.phi_2.evaluate(t - ti); + } + } + + // Self-excitation from sell events + for &ti in &self.sell_events { + if ti < t { + lambda += self.phi_1.evaluate(t - ti); + } + } + + lambda + } + + /// Simulate given core order flow as driver + /// + /// # Arguments + /// * `core_buy_events` - Core buy order arrival times F⁺ + /// * `core_sell_events` - Core sell order arrival times F⁻ + /// * `t_max` - Maximum simulation time + /// * `seed` - Optional random seed + pub fn simulate_driven( + &mut self, + core_buy_events: &[f64], + core_sell_events: &[f64], + t_max: f64, + seed: Option, + ) -> (Vec, Vec) { + let mut rng = match seed { + Some(s) => StdRng::seed_from_u64(s), + None => StdRng::from_entropy(), + }; + + self.buy_events.clear(); + self.sell_events.clear(); + self.current_time = 0.0; + + // Compute baseline intensity from core flow reaction + let baseline_buy = |t: f64| -> f64 { + let mut mu = 0.0; + for &ti in core_buy_events { + if ti < t { mu += self.phi_1.evaluate(t - ti); } + } + for &ti in core_sell_events { + if ti < t { mu += self.phi_2.evaluate(t - ti); } + } + mu.max(0.001) // Ensure positive + }; + + let baseline_sell = |t: f64| -> f64 { + let mut mu = 0.0; + for &ti in core_buy_events { + if ti < t { mu += self.phi_2.evaluate(t - ti); } + } + for &ti in core_sell_events { + if ti < t { mu += self.phi_1.evaluate(t - ti); } + } + mu.max(0.001) + }; + + // Use thinning algorithm for bivariate process + let safety_factor = 2.0; + + while self.current_time < t_max { + let lambda_buy = self.buy_intensity(self.current_time, baseline_buy(self.current_time)); + let lambda_sell = self.sell_intensity(self.current_time, baseline_sell(self.current_time)); + let lambda_max = (lambda_buy + lambda_sell) * safety_factor + 0.1; + + let exp_dist = Exp::new(lambda_max).unwrap(); + let tau: f64 = rng.sample(exp_dist); + let candidate_time = self.current_time + tau; + + if candidate_time >= t_max { + break; + } + + // Accept/reject and determine type + let u: f64 = rng.gen(); + let lambda_buy_cand = self.buy_intensity(candidate_time, baseline_buy(candidate_time)); + let lambda_sell_cand = self.sell_intensity(candidate_time, baseline_sell(candidate_time)); + let total_intensity = lambda_buy_cand + lambda_sell_cand; + + if u <= total_intensity / lambda_max { + // Accept - determine buy or sell + let p_buy = lambda_buy_cand / total_intensity; + if rng.gen::() < p_buy { + self.buy_events.push(candidate_time); + } else { + self.sell_events.push(candidate_time); + } + } + + self.current_time = candidate_time; + } + + (self.buy_events.clone(), self.sell_events.clone()) + } + + /// Spectral radius of kernel matrix (stability condition) + pub fn spectral_radius(&self) -> f64 { + // For symmetric 2x2 [[φ₁, φ₂], [φ₂, φ₁]], eigenvalues are φ₁+φ₂ and φ₁-φ₂ + // Spectral radius = max(|φ₁+φ₂|, |φ₁-φ₂|) in L¹ norm + let l1_phi1 = self.phi_1.l1_norm(); + let l1_phi2 = self.phi_2.l1_norm(); + (l1_phi1 + l1_phi2).max((l1_phi1 - l1_phi2).abs()) + } + + /// Check if the process is stable (spectral radius < 1) + pub fn is_stable(&self) -> bool { + self.spectral_radius() < 1.0 + } + + /// Get signed reaction flow N⁺ - N⁻ + pub fn signed_flow(&self) -> Vec<(f64, i32)> { + let mut events: Vec<(f64, i32)> = Vec::new(); + events.extend(self.buy_events.iter().map(|&t| (t, 1))); + events.extend(self.sell_events.iter().map(|&t| (t, -1))); + events.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + events + } + + /// Get unsigned reaction volume N⁺ + N⁻ + pub fn unsigned_volume(&self) -> usize { + self.buy_events.len() + self.sell_events.len() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::point_processes::kernels::{ExponentialKernel, PowerLawKernel}; + + #[test] + fn test_hawkes_simulation() { + let kernel = ExponentialKernel::new(0.3, 1.0); + let mut process = HawkesProcess::new(1.0, kernel); + + let events = process.simulate(100.0, Some(42)); + assert!(!events.is_empty()); + + // Check events are sorted + for i in 1..events.len() { + assert!(events[i] > events[i-1]); + } + } + + #[test] + fn test_hawkes_intensity() { + let kernel = ExponentialKernel::new(0.5, 1.0); + let mut process = HawkesProcess::new(1.0, kernel); + + // Initial intensity should be baseline + assert!((process.intensity(0.0) - 1.0).abs() < 1e-10); + + // Add an event and check intensity increases + process.events.push(1.0); + let intensity_after = process.intensity(1.001); + assert!(intensity_after > 1.0); + } + + #[test] + fn test_power_law_hawkes() { + let kernel = PowerLawKernel::nearly_critical(0.375, 0.05); // H₀ ≈ 0.75 + let mut process = HawkesProcess::new(0.1, kernel); + + let events = process.simulate(100.0, Some(123)); + println!("Power-law Hawkes: {} events in [0, 100]", events.len()); + assert!(!events.is_empty()); + } + + #[test] + fn test_bivariate_stability() { + let phi_1 = ExponentialKernel::new(0.3, 1.0); // L¹ = 0.3 + let phi_2 = ExponentialKernel::new(0.2, 1.0); // L¹ = 0.2 + + let process = BivariateHawkes::new(phi_1, phi_2); + assert!(process.is_stable()); // 0.3 + 0.2 = 0.5 < 1 + + // Unstable case + let phi_1_unstable = ExponentialKernel::new(0.8, 1.0); + let phi_2_unstable = ExponentialKernel::new(0.5, 1.0); + let process_unstable = BivariateHawkes::new(phi_1_unstable, phi_2_unstable); + assert!(!process_unstable.is_stable()); // 0.8 + 0.5 = 1.3 > 1 + } +} diff --git a/src/point_processes/kernels.rs b/src/point_processes/kernels.rs new file mode 100644 index 0000000..b0ae1a3 --- /dev/null +++ b/src/point_processes/kernels.rs @@ -0,0 +1,307 @@ +//! Excitation Kernels for Hawkes Processes +//! +//! This module provides different kernel functions that describe how past events +//! influence the intensity of future arrivals in a Hawkes process. +//! +//! # Kernels +//! +//! - **ExponentialKernel**: φ(t) = α * exp(-β * t), for short-range dependence +//! - **PowerLawKernel**: φ(t) = K₀ * (1 + t)^(-1-α₀), for long-range dependence + +use std::f64::consts::PI; + +/// Trait for excitation kernels in Hawkes processes +pub trait ExcitationKernel: Clone + Send + Sync { + /// Evaluate the kernel at time t (t >= 0) + fn evaluate(&self, t: f64) -> f64; + + /// Integrate the kernel from 0 to t + fn integrate(&self, t: f64) -> f64; + + /// L¹ norm of the kernel (total mass) + fn l1_norm(&self) -> f64; + + /// Check if the process is stable (L¹ norm < 1 for stability) + fn is_stable(&self) -> bool { + self.l1_norm() < 1.0 + } + + /// Tail exponent α₀ (for asymptotic analysis) + fn tail_exponent(&self) -> Option; +} + +/// Exponential kernel: φ(t) = α * exp(-β * t) +/// +/// Suitable for processes with short-range temporal dependence. +/// The L¹ norm is α/β. +#[derive(Clone, Debug)] +pub struct ExponentialKernel { + /// Peak intensity α > 0 + pub alpha: f64, + /// Decay rate β > 0 + pub beta: f64, +} + +impl ExponentialKernel { + pub fn new(alpha: f64, beta: f64) -> Self { + assert!(alpha > 0.0, "alpha must be positive"); + assert!(beta > 0.0, "beta must be positive"); + Self { alpha, beta } + } + + /// Create a stable kernel with given L¹ norm (< 1) + pub fn with_branching_ratio(branching_ratio: f64, beta: f64) -> Self { + assert!(branching_ratio > 0.0 && branching_ratio < 1.0); + assert!(beta > 0.0); + Self { + alpha: branching_ratio * beta, + beta, + } + } +} + +impl ExcitationKernel for ExponentialKernel { + fn evaluate(&self, t: f64) -> f64 { + if t < 0.0 { + return 0.0; + } + self.alpha * (-self.beta * t).exp() + } + + fn integrate(&self, t: f64) -> f64 { + if t <= 0.0 { + return 0.0; + } + (self.alpha / self.beta) * (1.0 - (-self.beta * t).exp()) + } + + fn l1_norm(&self) -> f64 { + self.alpha / self.beta + } + + fn tail_exponent(&self) -> Option { + // Exponential kernel decays faster than any power law + None + } +} + +/// Power-law kernel: φ(t) = K₀ * (1 + t)^(-1-α₀) +/// +/// Suitable for processes with long-range temporal dependence (memory). +/// Used in the unified theory paper for core order flow modeling. +/// +/// Parameters: +/// - α₀ ∈ (0, 1): tail exponent controlling memory persistence +/// - K₀ > 0: scaling constant +/// +/// Smaller α₀ means stronger persistence (longer memory). +#[derive(Clone, Debug)] +pub struct PowerLawKernel { + /// Tail exponent α₀ ∈ (0, 1) + pub alpha_0: f64, + /// Scaling constant K₀ > 0 + pub k_0: f64, + /// Normalization factor to achieve unit L¹ norm + norm_factor: f64, +} + +impl PowerLawKernel { + /// Create a power-law kernel with given parameters + /// + /// # Arguments + /// * `alpha_0` - Tail exponent in (0, 1) + /// * `k_0` - Scaling constant > 0 + pub fn new(alpha_0: f64, k_0: f64) -> Self { + assert!(alpha_0 > 0.0 && alpha_0 < 1.0, "alpha_0 must be in (0, 1)"); + assert!(k_0 > 0.0, "k_0 must be positive"); + + // L¹ norm = K₀ / α₀ (integral of (1+t)^{-1-α₀} from 0 to ∞) + let norm_factor = alpha_0 / k_0; + + Self { alpha_0, k_0, norm_factor } + } + + /// Create a kernel with unit L¹ norm (critical regime) + pub fn unit_norm(alpha_0: f64) -> Self { + assert!(alpha_0 > 0.0 && alpha_0 < 1.0); + Self { + alpha_0, + k_0: alpha_0, // This gives L¹ norm = 1 + norm_factor: 1.0, + } + } + + /// Create a nearly-critical kernel (L¹ norm = 1 - ε) + pub fn nearly_critical(alpha_0: f64, epsilon: f64) -> Self { + assert!(alpha_0 > 0.0 && alpha_0 < 1.0); + assert!(epsilon > 0.0 && epsilon < 1.0); + let k_0 = alpha_0 * (1.0 - epsilon); + Self::new(alpha_0, k_0) + } + + /// Get the corresponding Hurst exponent H₀ = 2 * α₀ + pub fn hurst_exponent(&self) -> f64 { + 2.0 * self.alpha_0 + } +} + +impl ExcitationKernel for PowerLawKernel { + fn evaluate(&self, t: f64) -> f64 { + if t < 0.0 { + return 0.0; + } + self.k_0 * (1.0 + t).powf(-1.0 - self.alpha_0) + } + + fn integrate(&self, t: f64) -> f64 { + if t <= 0.0 { + return 0.0; + } + // ∫₀ᵗ K₀(1+s)^{-1-α₀} ds = (K₀/α₀) * [1 - (1+t)^{-α₀}] + (self.k_0 / self.alpha_0) * (1.0 - (1.0 + t).powf(-self.alpha_0)) + } + + fn l1_norm(&self) -> f64 { + self.k_0 / self.alpha_0 + } + + fn tail_exponent(&self) -> Option { + Some(self.alpha_0) + } +} + +/// Completely monotone power-law kernel (as in Assumption A of the paper) +/// +/// This satisfies the complete monotonicity requirement for the scaling limit +/// theorems. φ(t) = K₀ * t^{-α₀} * E_{1-α₀}(-λ * t^{1-α₀}) +/// where E is the Mittag-Leffler function. +#[derive(Clone, Debug)] +pub struct CompletelyMonotoneKernel { + pub alpha_0: f64, + pub k_0: f64, + pub lambda: f64, +} + +impl CompletelyMonotoneKernel { + pub fn new(alpha_0: f64, k_0: f64, lambda: f64) -> Self { + assert!(alpha_0 > 0.0 && alpha_0 < 1.0); + assert!(k_0 > 0.0); + assert!(lambda > 0.0); + Self { alpha_0, k_0, lambda } + } +} + +impl ExcitationKernel for CompletelyMonotoneKernel { + fn evaluate(&self, t: f64) -> f64 { + if t <= 0.0 { + return 0.0; + } + // Use the Mittag-Leffler approximation for now + // Full implementation requires mittag_leffler function + let ml_arg = -self.lambda * t.powf(1.0 - self.alpha_0); + + // Approximate E_{1-α₀}(z) ≈ exp(z^{1/(1-α₀)}) / (1-α₀) for large |z| + // For small arguments, E_α(z) ≈ 1 + z/Γ(1+α) + let ml_approx = if ml_arg.abs() < 0.1 { + 1.0 + ml_arg / gamma_fn(2.0 - self.alpha_0) + } else { + // Asymptotic expansion + (-ml_arg).powf(-1.0) / gamma_fn(self.alpha_0) + }; + + self.k_0 * t.powf(-self.alpha_0) * ml_approx + } + + fn integrate(&self, t: f64) -> f64 { + // Numerical integration fallback + let n = 1000; + let dt = t / n as f64; + let mut sum = 0.0; + for i in 0..n { + let ti = (i as f64 + 0.5) * dt; + sum += self.evaluate(ti); + } + sum * dt + } + + fn l1_norm(&self) -> f64 { + // For completely monotone kernels, need numerical approximation + // or analytical form based on Mittag-Leffler integral + 1.0 // Placeholder for unit norm case + } + + fn tail_exponent(&self) -> Option { + Some(self.alpha_0) + } +} + +/// Gamma function approximation (Lanczos approximation) +fn gamma_fn(z: f64) -> f64 { + // Use Lanczos approximation for Γ(z) + if z < 0.5 { + PI / ((PI * z).sin() * gamma_fn(1.0 - z)) + } else { + let z = z - 1.0; + let g = 7; + let c = [ + 0.99999999999980993, + 676.5203681218851, + -1259.1392167224028, + 771.32342877765313, + -176.61502916214059, + 12.507343278686905, + -0.13857109526572012, + 9.9843695780195716e-6, + 1.5056327351493116e-7, + ]; + + let mut x = c[0]; + for i in 1..=(g + 1) { + x += c[i] / (z + i as f64); + } + + let t = z + g as f64 + 0.5; + (2.0 * PI).sqrt() * t.powf(z + 0.5) * (-t).exp() * x + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_exponential_kernel() { + let kernel = ExponentialKernel::new(0.5, 1.0); + assert!((kernel.l1_norm() - 0.5).abs() < 1e-10); + assert!(kernel.is_stable()); + + // Check decay + let v0 = kernel.evaluate(0.0); + let v1 = kernel.evaluate(1.0); + assert!(v1 < v0); + assert!((v0 - 0.5).abs() < 1e-10); + } + + #[test] + fn test_power_law_kernel() { + let kernel = PowerLawKernel::new(0.375, 0.375); // H₀ = 0.75 + assert!((kernel.l1_norm() - 1.0).abs() < 1e-10); + assert!((kernel.hurst_exponent() - 0.75).abs() < 1e-10); + + // Check power-law decay + let v1 = kernel.evaluate(1.0); + let v10 = kernel.evaluate(10.0); + let v100 = kernel.evaluate(100.0); + + // φ(t) ~ t^{-1-α₀}, so φ(10)/φ(1) ≈ 10^{-1-α₀} + let expected_ratio = 10.0_f64.powf(-1.375); + assert!((v10 / v1 - expected_ratio).abs() < 0.1); + } + + #[test] + fn test_nearly_critical() { + let kernel = PowerLawKernel::nearly_critical(0.5, 0.01); + assert!((kernel.l1_norm() - 0.99).abs() < 1e-10); + assert!(kernel.is_stable()); + } +} diff --git a/src/point_processes/mittag_leffler.rs b/src/point_processes/mittag_leffler.rs new file mode 100644 index 0000000..b8c404d --- /dev/null +++ b/src/point_processes/mittag_leffler.rs @@ -0,0 +1,304 @@ +//! Mittag-Leffler Functions +//! +//! The Mittag-Leffler function E_{α,β}(z) is a generalization of the +//! exponential function that plays a crucial role in fractional calculus +//! and the scaling limits of Hawkes processes. +//! +//! # Definition +//! +//! E_{α,β}(z) = Σ_{k=0}^∞ z^k / Γ(αk + β) +//! +//! Special cases: +//! - E_{1,1}(z) = exp(z) +//! - E_{2,1}(z²) = cosh(z) + +use std::f64::consts::PI; + +/// Compute the Mittag-Leffler function E_{α,β}(z) +/// +/// Uses series expansion for |z| < R and asymptotic expansion for |z| > R. +/// +/// # Arguments +/// * `alpha` - Parameter α > 0 +/// * `beta` - Parameter β > 0 +/// * `z` - Complex argument (real part only) +/// +/// # Returns +/// The value E_{α,β}(z) +pub fn mittag_leffler(alpha: f64, beta: f64, z: f64) -> f64 { + assert!(alpha > 0.0, "alpha must be positive"); + assert!(beta > 0.0, "beta must be positive"); + + if z.abs() < 1e-15 { + return 1.0 / gamma(beta); + } + + // Use different methods based on |z| + if z.abs() < 10.0 { + mittag_leffler_series(alpha, beta, z, 100) + } else { + mittag_leffler_asymptotic(alpha, beta, z) + } +} + +/// Series expansion for Mittag-Leffler function +fn mittag_leffler_series(alpha: f64, beta: f64, z: f64, max_terms: usize) -> f64 { + let mut sum: f64 = 0.0; + let mut z_pow: f64 = 1.0; // z^k + + for k in 0..max_terms { + let term = z_pow / gamma(alpha * k as f64 + beta); + + // Check for convergence + if term.abs() < 1e-15 * sum.abs() && k > 10 { + break; + } + + sum += term; + z_pow *= z; + + // Prevent overflow + if z_pow.abs() > 1e100 { + break; + } + } + + sum +} + +/// Asymptotic expansion for large |z| +fn mittag_leffler_asymptotic(alpha: f64, beta: f64, z: f64) -> f64 { + if alpha < 1.0 { + // For 0 < α < 1, different behavior for positive/negative z + if z > 0.0 { + // Leading term: (1/α) * z^{(1-β)/α} * exp(z^{1/α}) + let z_pow = z.powf(1.0 / alpha); + (1.0 / alpha) * z.powf((1.0 - beta) / alpha) * z_pow.exp() + } else { + // For z < 0, algebraic decay + // E_{α,β}(-x) ~ -Σ_{k=1}^∞ (-x)^{-k} / Γ(β - αk) + let mut sum = 0.0; + let z_abs = z.abs(); + for k in 1..=10 { + let term = (-1.0_f64).powi(k as i32 + 1) + * z_abs.powi(-(k as i32)) + / gamma(beta - alpha * k as f64); + sum += term; + } + sum + } + } else if alpha == 1.0 { + // E_{1,β}(z) ~ z^{1-β} * exp(z) for large z > 0 + if z > 0.0 { + z.powf(1.0 - beta) * z.exp() + } else { + // E_{1,β}(-x) ~ 0 for large x + 0.0 + } + } else { + // For α > 1, more complex behavior + // Use series with damping + mittag_leffler_series(alpha, beta, z, 200) + } +} + +/// Derivative of Mittag-Leffler function: d/dz E_{α,β}(z) +pub fn mittag_leffler_derivative(alpha: f64, beta: f64, z: f64) -> f64 { + // d/dz E_{α,β}(z) = (1/α) * (E_{α,β-1}(z) - (β-1) * E_{α,β}(z)) + // More stable: use series directly + + if z.abs() < 1e-15 { + return 1.0 / gamma(alpha + beta); + } + + let mut sum = 0.0; + let mut z_pow = 1.0; + + for k in 1..100 { + let term = k as f64 * z_pow / gamma(alpha * k as f64 + beta); + sum += term; + z_pow *= z; + + if z_pow.abs() > 1e100 || (term.abs() < 1e-15 * sum.abs() && k > 10) { + break; + } + } + + sum / z // Factor out one z from derivative +} + +/// The function f_{α₀,λ₀}(x) from Theorem 3.1 of the paper +/// +/// f_{α₀,λ₀}(x) = λ₀ * x^{α₀-1} * E_{α₀,α₀}(-λ₀ * x^{α₀}) +/// +/// This function controls the scaling limit of the Hawkes process. +pub fn f_alpha_lambda(alpha_0: f64, lambda_0: f64, x: f64) -> f64 { + assert!(alpha_0 > 0.0 && alpha_0 < 1.0); + assert!(lambda_0 > 0.0); + + if x <= 0.0 { + return 0.0; + } + + let x_alpha = x.powf(alpha_0); + let ml_arg = -lambda_0 * x_alpha; + + lambda_0 * x.powf(alpha_0 - 1.0) * mittag_leffler(alpha_0, alpha_0, ml_arg) +} + +/// Integral of f_{α₀,λ₀} from 0 to t +/// +/// ∫₀ᵗ f_{α₀,λ₀}(s) ds = t^{α₀} * E_{α₀,α₀+1}(-λ₀ * t^{α₀}) +pub fn f_alpha_lambda_integral(alpha_0: f64, lambda_0: f64, t: f64) -> f64 { + if t <= 0.0 { + return 0.0; + } + + let t_alpha = t.powf(alpha_0); + let ml_arg = -lambda_0 * t_alpha; + + t_alpha * mittag_leffler(alpha_0, alpha_0 + 1.0, ml_arg) +} + +/// Gamma function using Lanczos approximation +pub fn gamma(z: f64) -> f64 { + if z < 0.5 { + // Reflection formula: Γ(z) * Γ(1-z) = π / sin(πz) + PI / ((PI * z).sin() * gamma(1.0 - z)) + } else { + let z = z - 1.0; + let g = 7; + let c = [ + 0.99999999999980993, + 676.5203681218851, + -1259.1392167224028, + 771.32342877765313, + -176.61502916214059, + 12.507343278686905, + -0.13857109526572012, + 9.9843695780195716e-6, + 1.5056327351493116e-7, + ]; + + let mut x = c[0]; + for i in 1..=(g + 1) { + x += c[i] / (z + i as f64); + } + + let t = z + g as f64 + 0.5; + (2.0 * PI).sqrt() * t.powf(z + 0.5) * (-t).exp() * x + } +} + +/// Log-gamma function for numerical stability +pub fn lgamma(z: f64) -> f64 { + gamma(z).abs().ln() +} + +/// Incomplete gamma function γ(s, x) = ∫₀ˣ t^{s-1} e^{-t} dt +/// Used for various probability computations +pub fn incomplete_gamma_lower(s: f64, x: f64) -> f64 { + if x < 0.0 || s <= 0.0 { + return 0.0; + } + + if x < s + 1.0 { + // Series expansion + let mut sum = 0.0; + let mut term = 1.0 / s; + sum += term; + + for n in 1..100 { + term *= x / (s + n as f64); + sum += term; + if term.abs() < 1e-12 * sum.abs() { + break; + } + } + + sum * x.powf(s) * (-x).exp() + } else { + // Continued fraction for large x + gamma(s) - incomplete_gamma_upper(s, x) + } +} + +/// Upper incomplete gamma Γ(s, x) = ∫ₓ^∞ t^{s-1} e^{-t} dt +pub fn incomplete_gamma_upper(s: f64, x: f64) -> f64 { + if x < 0.0 { + return gamma(s); + } + + // Continued fraction representation (Lentz's algorithm) + let mut f = 1.0 + x - s; + if f.abs() < 1e-30 { + f = 1e-30; + } + + let mut c = f; + let mut d = 0.0; + + for i in 1..100 { + let an = (i as f64) * (s - i as f64); + let bn = 2.0 * i as f64 + 1.0 + x - s; + + d = bn + an * d; + if d.abs() < 1e-30 { d = 1e-30; } + + c = bn + an / c; + if c.abs() < 1e-30 { c = 1e-30; } + + d = 1.0 / d; + let delta = c * d; + f *= delta; + + if (delta - 1.0).abs() < 1e-10 { + break; + } + } + + x.powf(s) * (-x).exp() / f +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_gamma() { + // Γ(1) = 1 + assert!((gamma(1.0) - 1.0).abs() < 1e-10); + + // Γ(2) = 1! = 1 + assert!((gamma(2.0) - 1.0).abs() < 1e-10); + + // Γ(3) = 2! = 2 + assert!((gamma(3.0) - 2.0).abs() < 1e-10); + + // Γ(1/2) = √π + assert!((gamma(0.5) - PI.sqrt()).abs() < 1e-10); + } + + #[test] + fn test_mittag_leffler_special_cases() { + // E_{1,1}(z) = exp(z) + let z = 1.5; + let ml = mittag_leffler(1.0, 1.0, z); + assert!((ml - z.exp()).abs() < 1e-8); + + // E_{α,β}(0) = 1/Γ(β) + let ml_zero = mittag_leffler(0.5, 1.0, 0.0); + assert!((ml_zero - 1.0).abs() < 1e-10); + } + + #[test] + fn test_f_alpha_lambda() { + // Test that f is positive for positive arguments + let val = f_alpha_lambda(0.375, 1.0, 1.0); + assert!(val > 0.0); + + // Test integral property + let integral = f_alpha_lambda_integral(0.5, 1.0, 2.0); + assert!(integral > 0.0); + } +} diff --git a/src/point_processes/mixed_fbm.rs b/src/point_processes/mixed_fbm.rs new file mode 100644 index 0000000..41a053b --- /dev/null +++ b/src/point_processes/mixed_fbm.rs @@ -0,0 +1,444 @@ +//! Mixed Fractional Brownian Motion +//! +//! The mixed fractional Brownian motion (mfBM) is the sum of a standard +//! Brownian motion and an independent fractional Brownian motion. +//! +//! MH(t) = a * B(t) + b * BH(t) +//! +//! where: +//! - B(t) is standard Brownian motion (Hurst H = 1/2) +//! - BH(t) is fractional Brownian motion with Hurst index H +//! - a, b are mixing coefficients +//! +//! This is the limiting process for aggregate order flow in the unified theory. + +use rand::prelude::*; +use rand_distr::Normal; +use std::f64::consts::PI; + +/// Fractional Brownian Motion with Hurst parameter H +#[derive(Clone, Debug)] +pub struct FractionalBM { + /// Hurst parameter H ∈ (0, 1) + pub hurst: f64, + /// Time grid + pub times: Vec, + /// Sample path values + pub values: Vec, +} + +impl FractionalBM { + /// Create a new fBM specification + pub fn new(hurst: f64) -> Self { + assert!(hurst > 0.0 && hurst < 1.0, "Hurst parameter must be in (0, 1)"); + Self { + hurst, + times: Vec::new(), + values: Vec::new(), + } + } + + /// Covariance function of fBM: E[BH(s) * BH(t)] + pub fn covariance(&self, s: f64, t: f64) -> f64 { + let h = self.hurst; + let two_h = 2.0 * h; + 0.5 * (s.powf(two_h) + t.powf(two_h) - (s - t).abs().powf(two_h)) + } + + /// Variance of increment: Var[BH(t) - BH(s)] + pub fn increment_variance(&self, s: f64, t: f64) -> f64 { + let h = self.hurst; + (t - s).abs().powf(2.0 * h) + } + + /// Simulate a sample path using Cholesky method + /// + /// # Arguments + /// * `times` - Time points + /// * `seed` - Optional random seed + pub fn simulate(&mut self, times: &[f64], seed: Option) -> Vec { + let n = times.len(); + if n == 0 { + return Vec::new(); + } + + let mut rng = match seed { + Some(s) => StdRng::seed_from_u64(s), + None => StdRng::from_entropy(), + }; + + // Build covariance matrix + let mut cov = vec![vec![0.0; n]; n]; + for i in 0..n { + for j in 0..n { + cov[i][j] = self.covariance(times[i], times[j]); + } + } + + // Cholesky decomposition + let chol = cholesky(&cov).expect("Cholesky decomposition failed"); + + // Generate standard normal samples + let normal = Normal::new(0.0, 1.0).unwrap(); + let z: Vec = (0..n).map(|_| rng.sample(normal)).collect(); + + // Apply Cholesky factor + let mut values = vec![0.0; n]; + for i in 0..n { + for j in 0..=i { + values[i] += chol[i][j] * z[j]; + } + } + + self.times = times.to_vec(); + self.values = values.clone(); + values + } + + /// Simulate using Hosking's method (more efficient for regular grids) + pub fn simulate_hosking(&mut self, n: usize, dt: f64, seed: Option) -> Vec { + let mut rng = match seed { + Some(s) => StdRng::seed_from_u64(s), + None => StdRng::from_entropy(), + }; + + let h = self.hurst; + let normal = Normal::new(0.0, 1.0).unwrap(); + + // Hosking's algorithm for fractional Gaussian noise + // Then cumsum to get fBM + + // Autocovariance of fGn: γ(k) = 0.5 * (|k-1|^{2H} - 2|k|^{2H} + |k+1|^{2H}) + let acf = |k: i64| -> f64 { + let k_abs = k.abs() as f64; + let h2 = 2.0 * h; + 0.5 * ((k_abs - 1.0).abs().powf(h2) - 2.0 * k_abs.powf(h2) + (k_abs + 1.0).powf(h2)) + }; + + // Durbin-Levinson algorithm for prediction coefficients + let mut phi = vec![vec![0.0; n]; n]; + let mut v = vec![0.0; n]; + + v[0] = acf(0); + phi[0][0] = 0.0; + + let mut fgn = vec![0.0; n]; // Fractional Gaussian noise + fgn[0] = v[0].sqrt() * rng.sample(normal); + + for i in 1..n { + // Compute phi_ii + let num: f64 = acf(i as i64) - (0..i).map(|j| phi[i-1][j] * acf((i - 1 - j) as i64)).sum::(); + phi[i][i] = num / v[i-1]; + + // Update other coefficients + for j in 0..i { + phi[i][j] = phi[i-1][j] - phi[i][i] * phi[i-1][i-1-j]; + } + + // Update variance + v[i] = v[i-1] * (1.0 - phi[i][i].powi(2)); + + // Generate fGn sample + let pred: f64 = (0..i).map(|j| phi[i][j] * fgn[i-1-j]).sum(); + fgn[i] = pred + v[i].sqrt() * rng.sample(normal); + } + + // Cumulative sum to get fBM + let scale = dt.powf(h); + let mut values = vec![0.0; n + 1]; + for i in 0..n { + values[i + 1] = values[i] + scale * fgn[i]; + } + + self.times = (0..=n).map(|i| i as f64 * dt).collect(); + self.values = values.clone(); + values + } + + /// Estimate Hurst exponent from data using rescaled range (R/S) analysis + pub fn estimate_hurst(data: &[f64]) -> f64 { + let n = data.len(); + if n < 20 { + return 0.5; // Not enough data + } + + let mut log_ns = Vec::new(); + let mut log_rs = Vec::new(); + + // Try different subseries lengths + let min_n = 10; + let mut sub_n = min_n; + while sub_n <= n / 2 { + let k = n / sub_n; // Number of subseries + let mut rs_values = Vec::new(); + + for i in 0..k { + let start = i * sub_n; + let end = start + sub_n; + let subseries = &data[start..end]; + + // Mean and standard deviation + let mean: f64 = subseries.iter().sum::() / sub_n as f64; + let var: f64 = subseries.iter().map(|x| (x - mean).powi(2)).sum::() / sub_n as f64; + let std_dev = var.sqrt(); + + if std_dev > 1e-10 { + // Cumulative deviations + let mut cumdev = vec![0.0; sub_n]; + cumdev[0] = subseries[0] - mean; + for j in 1..sub_n { + cumdev[j] = cumdev[j-1] + (subseries[j] - mean); + } + + // Range + let max_dev = cumdev.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + let min_dev = cumdev.iter().cloned().fold(f64::INFINITY, f64::min); + let r = max_dev - min_dev; + + rs_values.push(r / std_dev); + } + } + + if !rs_values.is_empty() { + let avg_rs: f64 = rs_values.iter().sum::() / rs_values.len() as f64; + log_ns.push((sub_n as f64).ln()); + log_rs.push(avg_rs.ln()); + } + + sub_n += sub_n / 4 + 1; + } + + // Linear regression to estimate H + if log_ns.len() < 3 { + return 0.5; + } + + let n_points = log_ns.len() as f64; + let mean_x: f64 = log_ns.iter().sum::() / n_points; + let mean_y: f64 = log_rs.iter().sum::() / n_points; + + let mut num = 0.0; + let mut den = 0.0; + for (x, y) in log_ns.iter().zip(log_rs.iter()) { + num += (x - mean_x) * (y - mean_y); + den += (x - mean_x).powi(2); + } + + let h = num / den; + h.clamp(0.01, 0.99) + } +} + +/// Mixed Fractional Brownian Motion +/// +/// MH(t) = a * B(t) + b * BH(t) +/// +/// Combines a standard BM (diffusive component) with an fBM (persistent component). +/// This is the limiting process for aggregate order flow. +#[derive(Clone, Debug)] +pub struct MixedFractionalBM { + /// Coefficient for standard BM component + pub a: f64, + /// Coefficient for fBM component + pub b: f64, + /// Hurst index H of the fBM component + pub hurst: f64, + /// Time grid + pub times: Vec, + /// Sample path values + pub values: Vec, +} + +impl MixedFractionalBM { + /// Create a new mixed fBM + pub fn new(a: f64, b: f64, hurst: f64) -> Self { + assert!(hurst > 0.0 && hurst < 1.0); + Self { + a, + b, + hurst, + times: Vec::new(), + values: Vec::new(), + } + } + + /// Create with Hurst index H₀ from unified theory (typically ~3/4) + /// Using a = 1, b = 1 (unit mixing) + pub fn from_h0(h0: f64) -> Self { + assert!(h0 > 0.5 && h0 < 1.0, "H0 should be in (0.5, 1) for persistent flow"); + Self::new(1.0, 1.0, h0) + } + + /// Covariance function + pub fn covariance(&self, s: f64, t: f64) -> f64 { + // Cov(MH(s), MH(t)) = a² * min(s,t) + b² * ρH(s,t) + // where ρH is the fBM covariance + let min_st = s.min(t); + let fbm_cov = { + let h = self.hurst; + let two_h = 2.0 * h; + 0.5 * (s.powf(two_h) + t.powf(two_h) - (s - t).abs().powf(two_h)) + }; + self.a.powi(2) * min_st + self.b.powi(2) * fbm_cov + } + + /// Simulate the mixed fBM + pub fn simulate(&mut self, n: usize, dt: f64, seed: Option) -> Vec { + let mut rng = match seed { + Some(s) => StdRng::seed_from_u64(s), + None => StdRng::from_entropy(), + }; + + // Simulate standard BM component + let normal = Normal::new(0.0, 1.0).unwrap(); + let sqrt_dt = dt.sqrt(); + let mut bm = vec![0.0; n + 1]; + for i in 0..n { + bm[i + 1] = bm[i] + sqrt_dt * rng.sample(normal); + } + + // Simulate fBM component + let mut fbm = FractionalBM::new(self.hurst); + let fbm_values = fbm.simulate_hosking(n, dt, seed.map(|s| s + 12345)); + + // Combine + let mut values = vec![0.0; n + 1]; + for i in 0..=n { + values[i] = self.a * bm[i] + self.b * fbm_values[i.min(fbm_values.len() - 1)]; + } + + self.times = (0..=n).map(|i| i as f64 * dt).collect(); + self.values = values.clone(); + values + } + + /// Check if process is semi-martingale (H > 3/4 for mfBM) + pub fn is_semimartingale(&self) -> bool { + self.hurst > 0.75 + } + + /// Estimate effective Hurst at different time scales + pub fn scale_dependent_hurst(data: &[f64], scales: &[usize]) -> Vec<(usize, f64)> { + let mut results = Vec::new(); + + for &scale in scales { + if scale >= data.len() / 4 { + continue; + } + + // Compute increments at this scale + let increments: Vec = (scale..data.len()) + .map(|i| data[i] - data[i - scale]) + .collect(); + + if increments.is_empty() { + continue; + } + + // Estimate local Hurst using variance ratio + let var_1: f64 = increments.iter().map(|x| x.powi(2)).sum::() / increments.len() as f64; + + // Compare with increments at double scale + let double_scale = 2 * scale; + if double_scale >= data.len() { + continue; + } + + let increments_2: Vec = (double_scale..data.len()) + .map(|i| data[i] - data[i - double_scale]) + .collect(); + + if increments_2.is_empty() { + continue; + } + + let var_2: f64 = increments_2.iter().map(|x| x.powi(2)).sum::() / increments_2.len() as f64; + + // E[|X_t - X_s|^2] ~ |t-s|^{2H} + // var_2 / var_1 ~ 2^{2H} + let h = (var_2 / var_1).ln() / (2.0_f64.ln() * 2.0); + results.push((scale, h.clamp(0.01, 0.99))); + } + + results + } +} + +/// Cholesky decomposition of a positive semi-definite matrix +fn cholesky(a: &[Vec]) -> Option>> { + let n = a.len(); + let mut l = vec![vec![0.0; n]; n]; + + for i in 0..n { + for j in 0..=i { + let mut sum = 0.0; + for k in 0..j { + sum += l[i][k] * l[j][k]; + } + + if i == j { + let diag = a[i][i] - sum; + if diag < -1e-10 { + return None; // Not positive semi-definite + } + l[i][j] = diag.max(0.0).sqrt(); + } else { + if l[j][j].abs() < 1e-10 { + l[i][j] = 0.0; + } else { + l[i][j] = (a[i][j] - sum) / l[j][j]; + } + } + } + } + + Some(l) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_fbm_covariance() { + let fbm = FractionalBM::new(0.75); + + // Variance at t=1 should be 1 + assert!((fbm.covariance(1.0, 1.0) - 1.0).abs() < 1e-10); + + // Cov should be non-negative + assert!(fbm.covariance(1.0, 2.0) >= 0.0); + } + + #[test] + fn test_fbm_simulation() { + let mut fbm = FractionalBM::new(0.7); + let values = fbm.simulate_hosking(100, 0.01, Some(42)); + + assert_eq!(values.len(), 101); + assert!((values[0]).abs() < 1e-10); // Starts at 0 + } + + #[test] + fn test_mixed_fbm() { + let mut mfbm = MixedFractionalBM::new(1.0, 0.5, 0.8); + let values = mfbm.simulate(100, 0.01, Some(42)); + + assert_eq!(values.len(), 101); + assert!((values[0]).abs() < 1e-10); + } + + #[test] + fn test_hurst_estimation() { + // Generate fBM with known H + let mut fbm = FractionalBM::new(0.7); + let values = fbm.simulate_hosking(1000, 1.0, Some(42)); + + // Estimate H + let h_est = FractionalBM::estimate_hurst(&values); + + // R/S method is noisy, allow wider tolerance (within 0.3) + // The estimation is approximate and depends on sample size + assert!(h_est > 0.3 && h_est < 1.0, "Estimated H={:.3} out of reasonable range", h_est); + } +} diff --git a/src/point_processes/mod.rs b/src/point_processes/mod.rs new file mode 100644 index 0000000..0a645f5 --- /dev/null +++ b/src/point_processes/mod.rs @@ -0,0 +1,63 @@ +//! Point Processes Module for Order Flow Modeling +//! +//! This module implements the mathematical framework from "A Unified Theory of +//! Order Flow, Market Impact, and Volatility" (Muhle-Karbe et al., 2026). +//! +//! # Overview +//! +//! The module provides tools for modeling order flow in financial markets using: +//! +//! - **Hawkes Processes**: Self-exciting point processes for core and reaction flows +//! - **Excitation Kernels**: Power-law and exponential kernels for temporal dependence +//! - **Mixed Fractional Brownian Motion**: Scaling limits of aggregate order flow +//! - **Mittag-Leffler Functions**: Key functions for scaling limit analysis +//! - **Order Flow Analysis**: Signed/unsigned flow, Hurst estimation, market impact +//! +//! # Key Relationships (Unified Theory) +//! +//! All quantities are determined by a single parameter H₀ ≈ 3/4: +//! +//! - Signed order flow: Hurst index H₀ +//! - Unsigned volume: Hurst index H₀ - 1/2 (rough, ~0.25) +//! - Volatility: Hurst index 2H₀ - 3/2 (~0 for H₀=3/4) +//! - Market impact: power law exponent 2 - 2H₀ (~0.5, square-root law) +//! +//! # Example +//! +//! ```rust,ignore +//! use optimizr::point_processes::{ +//! HawkesProcess, PowerLawKernel, OrderFlowAnalyzer +//! }; +//! +//! // Create a Hawkes process for core order flow +//! let kernel = PowerLawKernel::new(0.5, 1.0); // α₀ = 0.5, K₀ = 1.0 +//! let hawkes = HawkesProcess::new(0.1, kernel); // ν = 0.1 +//! +//! // Simulate order arrivals +//! let arrivals = hawkes.simulate(1000.0, Some(42)); +//! +//! // Analyze order flow +//! let analyzer = OrderFlowAnalyzer::new(); +//! let h0 = analyzer.estimate_h0(&arrivals); +//! println!("Estimated H₀: {:.4}", h0); +//! ``` + +mod kernels; +mod hawkes; +mod mittag_leffler; +mod mixed_fbm; +mod order_flow; + +#[cfg(feature = "python-bindings")] +pub mod python_bindings; + +// Re-export public API +pub use kernels::{ExcitationKernel, PowerLawKernel, ExponentialKernel}; +pub use hawkes::{HawkesProcess, HawkesProcessConfig, BivariateHawkes}; +pub use mittag_leffler::{mittag_leffler, mittag_leffler_derivative, f_alpha_lambda}; +pub use mixed_fbm::{MixedFractionalBM, FractionalBM}; +pub use order_flow::{ + OrderFlowAnalyzer, OrderFlowMetrics, UnifiedTheoryParams, + signed_order_flow, unsigned_volume, market_impact_exponent, + volatility_hurst, volume_hurst +}; diff --git a/src/point_processes/order_flow.rs b/src/point_processes/order_flow.rs new file mode 100644 index 0000000..49c53f7 --- /dev/null +++ b/src/point_processes/order_flow.rs @@ -0,0 +1,452 @@ +//! Order Flow Analysis Module +//! +//! Implements the unified theory framework for analyzing signed and unsigned +//! order flow, estimating H₀, and deriving market impact and volatility. +//! +//! # Key Relationships (from unified theory) +//! +//! Given H₀ ≈ 3/4 (persistence of core flow): +//! - Signed order flow: Hurst index H₀ +//! - Unsigned volume: Hurst index H₁ = H₀ - 1/2 ≈ 0.25 (rough) +//! - Volatility: Hurst index H_vol = 2H₀ - 3/2 ≈ 0 (very rough) +//! - Market impact: power law exponent δ = 2 - 2H₀ ≈ 0.5 (square root) + +use crate::point_processes::mixed_fbm::{FractionalBM, MixedFractionalBM}; +use std::f64::consts::PI; + +/// Parameters from the unified theory +#[derive(Clone, Debug)] +pub struct UnifiedTheoryParams { + /// H₀: Hurst index of signed order flow / core flow persistence + /// Typically H₀ ≈ 0.75 empirically + pub h0: f64, + + /// μ₀: Baseline intensity scaling constant + pub mu0: f64, + + /// λ₀: Decay rate parameter + pub lambda0: f64, +} + +impl UnifiedTheoryParams { + pub fn new(h0: f64) -> Self { + assert!(h0 > 0.5 && h0 < 1.0, "H0 must be in (0.5, 1)"); + Self { + h0, + mu0: 1.0, + lambda0: 1.0, + } + } + + /// Create with typical empirical value H₀ ≈ 0.75 + pub fn empirical() -> Self { + Self::new(0.75) + } + + /// Hurst index of unsigned volume: H₁ = H₀ - 1/2 + pub fn volume_hurst(&self) -> f64 { + self.h0 - 0.5 + } + + /// Hurst index of volatility: H_vol = 2H₀ - 3/2 + pub fn volatility_hurst(&self) -> f64 { + 2.0 * self.h0 - 1.5 + } + + /// Market impact exponent: δ = 2 - 2H₀ + pub fn impact_exponent(&self) -> f64 { + 2.0 - 2.0 * self.h0 + } + + /// Check if mfBM is semimartingale (H₀ > 3/4) + pub fn is_semimartingale(&self) -> bool { + self.h0 > 0.75 + } + + /// α₀ (tail exponent): H₀ = 2α₀, so α₀ = H₀/2 + pub fn alpha0(&self) -> f64 { + self.h0 / 2.0 + } +} + +/// Calculate signed order flow Hurst index (= H₀) +pub fn signed_order_flow(h0: f64) -> f64 { + h0 +} + +/// Calculate unsigned volume Hurst index: H₁ = H₀ - 1/2 +pub fn unsigned_volume(h0: f64) -> f64 { + h0 - 0.5 +} + +/// Alias for unsigned_volume +pub fn volume_hurst(h0: f64) -> f64 { + unsigned_volume(h0) +} + +/// Calculate volatility Hurst index: H_vol = 2H₀ - 3/2 +pub fn volatility_hurst(h0: f64) -> f64 { + 2.0 * h0 - 1.5 +} + +/// Calculate market impact exponent: δ = 2 - 2H₀ +/// Impact ~ Q^δ where Q is order size +pub fn market_impact_exponent(h0: f64) -> f64 { + 2.0 - 2.0 * h0 +} + +/// Order flow metrics computed from data +#[derive(Clone, Debug)] +pub struct OrderFlowMetrics { + /// Estimated H₀ (core flow persistence) + pub h0: f64, + + /// Estimated Hurst of signed flow (under fBM assumption) + pub h_signed_fbm: f64, + + /// Estimated Hurst of signed flow (under mfBM assumption) + pub h_signed_mfbm: f64, + + /// Estimated Hurst of unsigned volume + pub h_unsigned: f64, + + /// Total signed flow + pub total_signed: f64, + + /// Total unsigned volume + pub total_unsigned: f64, + + /// Implied volatility Hurst + pub h_volatility: f64, + + /// Implied impact exponent + pub impact_exponent: f64, + + /// Autocorrelation at lag 1 + pub acf_1: f64, + + /// Scale-dependent Hurst estimates + pub scale_hurst: Vec<(usize, f64)>, +} + +/// Analyzer for order flow data +#[derive(Clone, Debug, Default)] +pub struct OrderFlowAnalyzer { + /// Whether to compute scale-dependent statistics + pub compute_scales: bool, + + /// Scales for multi-scale analysis + pub scales: Vec, +} + +impl OrderFlowAnalyzer { + pub fn new() -> Self { + Self { + compute_scales: true, + scales: vec![10, 50, 100, 500, 1000, 2000, 5000], + } + } + + /// Analyze signed order flow + /// + /// # Arguments + /// * `flow` - Signed order flow data (positive = buy, negative = sell) + pub fn analyze_signed_flow(&self, flow: &[f64]) -> OrderFlowMetrics { + let n = flow.len(); + if n < 100 { + return self.default_metrics(); + } + + // Cumulative signed flow + let mut cum_flow = vec![0.0; n + 1]; + for i in 0..n { + cum_flow[i + 1] = cum_flow[i] + flow[i]; + } + + // Estimate H under fBM assumption (R/S analysis) + let h_fbm = FractionalBM::estimate_hurst(&cum_flow); + + // Estimate H under mfBM assumption (scale-dependent) + let scale_hurst = MixedFractionalBM::scale_dependent_hurst( + &cum_flow, + &self.scales, + ); + + // Average of high-frequency estimates (mfBM martingale component dominates) + let h_mfbm_hf: f64 = scale_hurst + .iter() + .filter(|(s, _)| *s < 100) + .map(|(_, h)| *h) + .sum::() / scale_hurst.iter().filter(|(s, _)| *s < 100).count().max(1) as f64; + + // Average of low-frequency estimates (fBM component dominates) + let h_mfbm_lf: f64 = scale_hurst + .iter() + .filter(|(s, _)| *s >= 500) + .map(|(_, h)| *h) + .sum::() / scale_hurst.iter().filter(|(s, _)| *s >= 500).count().max(1) as f64; + + // H₀ estimate: use low-frequency Hurst (persistent component) + let h0 = if h_mfbm_lf > 0.5 { h_mfbm_lf } else { h_fbm }; + + // Unsigned volume (absolute values) + let unsigned: Vec = flow.iter().map(|x| x.abs()).collect(); + let cum_unsigned: Vec = unsigned.iter() + .scan(0.0, |acc, &x| { *acc += x; Some(*acc) }) + .collect(); + let h_unsigned = estimate_hurst_variance_ratio(&cum_unsigned); + + // Autocorrelation + let mean: f64 = flow.iter().sum::() / n as f64; + let var: f64 = flow.iter().map(|x| (x - mean).powi(2)).sum::() / n as f64; + let acf_1 = if var > 1e-10 { + let cov: f64 = flow[..n-1].iter().zip(flow[1..].iter()) + .map(|(x, y)| (x - mean) * (y - mean)) + .sum::() / (n - 1) as f64; + cov / var + } else { + 0.0 + }; + + OrderFlowMetrics { + h0, + h_signed_fbm: h_fbm, + h_signed_mfbm: h_mfbm_lf, + h_unsigned, + total_signed: cum_flow.last().copied().unwrap_or(0.0), + total_unsigned: cum_unsigned.last().copied().unwrap_or(0.0), + h_volatility: volatility_hurst(h0), + impact_exponent: market_impact_exponent(h0), + acf_1, + scale_hurst, + } + } + + /// Analyze order arrival times + /// + /// # Arguments + /// * `buy_times` - Buy order arrival times + /// * `sell_times` - Sell order arrival times + /// * `t_max` - Maximum time + /// * `bin_size` - Time bin for aggregation + pub fn analyze_arrivals( + &self, + buy_times: &[f64], + sell_times: &[f64], + t_max: f64, + bin_size: f64, + ) -> OrderFlowMetrics { + // Bin the arrivals + let n_bins = (t_max / bin_size).ceil() as usize; + let mut signed_flow = vec![0.0; n_bins]; + + for &t in buy_times { + let bin = ((t / bin_size).floor() as usize).min(n_bins - 1); + signed_flow[bin] += 1.0; + } + for &t in sell_times { + let bin = ((t / bin_size).floor() as usize).min(n_bins - 1); + signed_flow[bin] -= 1.0; + } + + self.analyze_signed_flow(&signed_flow) + } + + /// Estimate H₀ from signed order flow data + pub fn estimate_h0(&self, flow: &[f64]) -> f64 { + self.analyze_signed_flow(flow).h0 + } + + fn default_metrics(&self) -> OrderFlowMetrics { + OrderFlowMetrics { + h0: 0.75, + h_signed_fbm: 0.5, + h_signed_mfbm: 0.75, + h_unsigned: 0.25, + total_signed: 0.0, + total_unsigned: 0.0, + h_volatility: 0.0, + impact_exponent: 0.5, + acf_1: 0.0, + scale_hurst: Vec::new(), + } + } +} + +/// Estimate Hurst exponent using variance ratio method +fn estimate_hurst_variance_ratio(data: &[f64]) -> f64 { + let n = data.len(); + if n < 50 { + return 0.5; + } + + let scales = [5, 10, 20, 40, 80, 160]; + let mut log_scales = Vec::new(); + let mut log_vars = Vec::new(); + + for &s in &scales { + if s >= n / 4 { + break; + } + + // Compute increments at scale s + let increments: Vec = (s..n).map(|i| data[i] - data[i - s]).collect(); + if increments.is_empty() { + continue; + } + + let var: f64 = increments.iter().map(|x| x.powi(2)).sum::() / increments.len() as f64; + if var > 1e-15 { + log_scales.push((s as f64).ln()); + log_vars.push(var.ln()); + } + } + + if log_scales.len() < 3 { + return 0.5; + } + + // Linear regression: log(var) = 2H * log(scale) + const + let n_pts = log_scales.len() as f64; + let mean_x: f64 = log_scales.iter().sum::() / n_pts; + let mean_y: f64 = log_vars.iter().sum::() / n_pts; + + let mut num = 0.0; + let mut den = 0.0; + for (x, y) in log_scales.iter().zip(log_vars.iter()) { + num += (x - mean_x) * (y - mean_y); + den += (x - mean_x).powi(2); + } + + let slope = num / den; + (slope / 2.0).clamp(0.01, 0.99) +} + +/// Market impact function: Impact(Q) ~ Q^δ where δ = 2 - 2H₀ +#[derive(Clone, Debug)] +pub struct MarketImpact { + /// Impact exponent δ + pub delta: f64, + /// Scaling constant + pub scale: f64, +} + +impl MarketImpact { + /// Create from unified theory parameter H₀ + pub fn from_h0(h0: f64, scale: f64) -> Self { + Self { + delta: 2.0 - 2.0 * h0, + scale, + } + } + + /// Create square-root impact (H₀ = 0.75) + pub fn square_root(scale: f64) -> Self { + Self { + delta: 0.5, + scale, + } + } + + /// Calculate market impact for order size Q + pub fn impact(&self, q: f64) -> f64 { + self.scale * q.abs().powf(self.delta) * q.signum() + } + + /// Inverse: order size needed for target impact + pub fn order_size(&self, target_impact: f64) -> f64 { + (target_impact.abs() / self.scale).powf(1.0 / self.delta) * target_impact.signum() + } +} + +/// Temporary price impact with decay +#[derive(Clone, Debug)] +pub struct TransientImpact { + /// Impact function + pub impact: MarketImpact, + /// Decay kernel exponent (controls how impact dissipates) + pub decay_exponent: f64, +} + +impl TransientImpact { + pub fn from_h0(h0: f64, scale: f64) -> Self { + Self { + impact: MarketImpact::from_h0(h0, scale), + decay_exponent: 2.0 * h0 - 1.0, // Derived from no-arbitrage + } + } + + /// Decay kernel G(t) ~ t^{-(2H₀-1)} + pub fn decay(&self, t: f64) -> f64 { + if t <= 0.0 { + return 1.0; + } + t.powf(-self.decay_exponent) + } + + /// Total impact at time t from orders (times, sizes) + pub fn total_impact(&self, t: f64, orders: &[(f64, f64)]) -> f64 { + let mut total = 0.0; + for (order_time, size) in orders { + if *order_time < t { + let elapsed = t - order_time; + total += self.impact.impact(*size) * self.decay(elapsed); + } + } + total + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_unified_theory_params() { + let params = UnifiedTheoryParams::new(0.75); + + assert!((params.volume_hurst() - 0.25).abs() < 1e-10); + assert!((params.volatility_hurst() - 0.0).abs() < 1e-10); + assert!((params.impact_exponent() - 0.5).abs() < 1e-10); + assert!(params.is_semimartingale() == false); // H₀ = 0.75 is boundary + + let params_high = UnifiedTheoryParams::new(0.8); + assert!(params_high.is_semimartingale()); + } + + #[test] + fn test_market_impact() { + let impact = MarketImpact::square_root(1.0); + + // Impact should be proportional to sqrt(Q) + let i1 = impact.impact(100.0); + let i4 = impact.impact(400.0); + + // i4 / i1 ≈ 2 (since sqrt(400) / sqrt(100) = 2) + assert!((i4 / i1 - 2.0).abs() < 0.01); + } + + #[test] + fn test_order_flow_analysis() { + // Generate some synthetic order flow + use rand::prelude::*; + use rand_distr::Normal; + + let mut rng = StdRng::seed_from_u64(42); + let normal = Normal::new(0.0, 1.0).unwrap(); + + // Add some persistence + let mut flow = vec![0.0; 1000]; + flow[0] = rng.sample(normal); + for i in 1..1000 { + flow[i] = 0.3 * flow[i-1] + rng.sample(normal); + } + + let analyzer = OrderFlowAnalyzer::new(); + let metrics = analyzer.analyze_signed_flow(&flow); + + // Hurst should be between 0 and 1 + assert!(metrics.h0 > 0.0 && metrics.h0 < 1.0); + assert!(metrics.h_signed_fbm > 0.0 && metrics.h_signed_fbm < 1.0); + } +} diff --git a/src/point_processes/python_bindings.rs b/src/point_processes/python_bindings.rs new file mode 100644 index 0000000..43980ea --- /dev/null +++ b/src/point_processes/python_bindings.rs @@ -0,0 +1,302 @@ +//! Python bindings for point processes module + +use pyo3::prelude::*; +use pyo3::types::PyDict; +use numpy::{PyArray1, PyReadonlyArray1}; + +use super::kernels::{ExcitationKernel, PowerLawKernel, ExponentialKernel}; +use super::hawkes::{HawkesProcess, HawkesProcessConfig, BivariateHawkes}; +use super::mittag_leffler::{mittag_leffler, f_alpha_lambda}; +use super::mixed_fbm::{FractionalBM, MixedFractionalBM}; +use super::order_flow::{OrderFlowAnalyzer, UnifiedTheoryParams, MarketImpact}; + +/// Simulate a univariate Hawkes process +/// +/// # Arguments +/// * `baseline` - Baseline intensity ν +/// * `alpha` - Kernel amplitude (for exponential) or K₀ (for power-law) +/// * `beta` - Decay rate (for exponential) or α₀ (for power-law) +/// * `t_max` - Maximum simulation time +/// * `kernel_type` - "exponential" or "power_law" +/// * `seed` - Optional random seed +/// +/// # Returns +/// Array of event times +#[pyfunction] +#[pyo3(signature = (baseline, alpha, beta, t_max, kernel_type="exponential", seed=None))] +pub fn simulate_hawkes<'py>( + py: Python<'py>, + baseline: f64, + alpha: f64, + beta: f64, + t_max: f64, + kernel_type: &str, + seed: Option, +) -> PyResult>> { + let events = match kernel_type { + "exponential" => { + let kernel = ExponentialKernel::new(alpha, beta); + let mut process = HawkesProcess::new(baseline, kernel); + process.simulate(t_max, seed) + } + "power_law" => { + let kernel = PowerLawKernel::new(beta, alpha); // beta = α₀, alpha = K₀ + let mut process = HawkesProcess::new(baseline, kernel); + process.simulate(t_max, seed) + } + _ => return Err(PyErr::new::( + format!("Unknown kernel type: {}. Use 'exponential' or 'power_law'", kernel_type) + )), + }; + + Ok(PyArray1::from_vec_bound(py, events)) +} + +/// Simulate a bivariate Hawkes process for buy/sell reaction flow +/// +/// # Arguments +/// * `core_buy_times` - Core buy order times (driver process) +/// * `core_sell_times` - Core sell order times (driver process) +/// * `phi1_alpha` - Same-side kernel amplitude +/// * `phi1_beta` - Same-side kernel decay +/// * `phi2_alpha` - Cross-side kernel amplitude +/// * `phi2_beta` - Cross-side kernel decay +/// * `t_max` - Maximum simulation time +/// * `seed` - Optional random seed +/// +/// # Returns +/// Tuple of (buy_times, sell_times) +#[pyfunction] +#[pyo3(signature = (core_buy_times, core_sell_times, phi1_alpha, phi1_beta, phi2_alpha, phi2_beta, t_max, seed=None))] +pub fn simulate_bivariate_hawkes<'py>( + py: Python<'py>, + core_buy_times: PyReadonlyArray1, + core_sell_times: PyReadonlyArray1, + phi1_alpha: f64, + phi1_beta: f64, + phi2_alpha: f64, + phi2_beta: f64, + t_max: f64, + seed: Option, +) -> PyResult<(Bound<'py, PyArray1>, Bound<'py, PyArray1>)> { + let core_buys = core_buy_times.as_slice()?; + let core_sells = core_sell_times.as_slice()?; + + let phi_1 = ExponentialKernel::new(phi1_alpha, phi1_beta); + let phi_2 = ExponentialKernel::new(phi2_alpha, phi2_beta); + + let mut process = BivariateHawkes::new(phi_1, phi_2); + let (buys, sells) = process.simulate_driven(core_buys, core_sells, t_max, seed); + + Ok(( + PyArray1::from_vec_bound(py, buys), + PyArray1::from_vec_bound(py, sells), + )) +} + +/// Compute Mittag-Leffler function E_{α,β}(z) +#[pyfunction] +#[pyo3(signature = (alpha, beta, z))] +pub fn mittag_leffler_py(alpha: f64, beta: f64, z: f64) -> PyResult { + Ok(mittag_leffler(alpha, beta, z)) +} + +/// Compute f_{α₀,λ₀}(x) function from unified theory scaling limits +#[pyfunction] +#[pyo3(signature = (alpha0, lambda0, x))] +pub fn f_alpha_lambda_py(alpha0: f64, lambda0: f64, x: f64) -> PyResult { + Ok(f_alpha_lambda(alpha0, lambda0, x)) +} + +/// Simulate fractional Brownian motion +/// +/// # Arguments +/// * `hurst` - Hurst parameter H ∈ (0, 1) +/// * `n` - Number of time steps +/// * `dt` - Time step size +/// * `seed` - Optional random seed +#[pyfunction] +#[pyo3(signature = (hurst, n, dt=1.0, seed=None))] +pub fn simulate_fbm<'py>( + py: Python<'py>, + hurst: f64, + n: usize, + dt: f64, + seed: Option, +) -> PyResult>> { + let mut fbm = FractionalBM::new(hurst); + let values = fbm.simulate_hosking(n, dt, seed); + Ok(PyArray1::from_vec_bound(py, values)) +} + +/// Simulate mixed fractional Brownian motion +/// +/// # Arguments +/// * `a` - Coefficient for standard BM component +/// * `b` - Coefficient for fBM component +/// * `hurst` - Hurst parameter of fBM component +/// * `n` - Number of time steps +/// * `dt` - Time step size +/// * `seed` - Optional random seed +#[pyfunction] +#[pyo3(signature = (a, b, hurst, n, dt=1.0, seed=None))] +pub fn simulate_mixed_fbm<'py>( + py: Python<'py>, + a: f64, + b: f64, + hurst: f64, + n: usize, + dt: f64, + seed: Option, +) -> PyResult>> { + let mut mfbm = MixedFractionalBM::new(a, b, hurst); + let values = mfbm.simulate(n, dt, seed); + Ok(PyArray1::from_vec_bound(py, values)) +} + +/// Estimate Hurst exponent using R/S analysis +#[pyfunction] +pub fn estimate_hurst<'py>( + _py: Python<'py>, + data: PyReadonlyArray1, +) -> PyResult { + let slice = data.as_slice()?; + Ok(FractionalBM::estimate_hurst(slice)) +} + +/// Compute scale-dependent Hurst exponents (for mfBM identification) +#[pyfunction] +#[pyo3(signature = (data, scales=None))] +pub fn scale_dependent_hurst<'py>( + py: Python<'py>, + data: PyReadonlyArray1, + scales: Option>, +) -> PyResult> { + let slice = data.as_slice()?; + let scales = scales.unwrap_or_else(|| vec![10, 50, 100, 500, 1000, 2000, 5000]); + + let results = MixedFractionalBM::scale_dependent_hurst(slice, &scales); + + let dict = PyDict::new_bound(py); + for (scale, h) in results { + dict.set_item(scale, h)?; + } + Ok(dict) +} + +/// Analyze order flow data using unified theory framework +/// +/// # Arguments +/// * `flow` - Signed order flow data (positive = buy, negative = sell) +/// +/// # Returns +/// Dictionary with metrics: +/// - h0: Estimated H₀ (core flow persistence) +/// - h_signed_fbm: Hurst under pure fBM assumption +/// - h_signed_mfbm: Hurst under mfBM assumption +/// - h_unsigned: Hurst of unsigned volume +/// - h_volatility: Implied volatility Hurst (2H₀ - 3/2) +/// - impact_exponent: Implied market impact exponent (2 - 2H₀) +/// - acf_1: First-order autocorrelation +/// - scale_hurst: Scale-dependent Hurst estimates +#[pyfunction] +pub fn analyze_order_flow<'py>( + py: Python<'py>, + flow: PyReadonlyArray1, +) -> PyResult> { + let slice = flow.as_slice()?; + + let analyzer = OrderFlowAnalyzer::new(); + let metrics = analyzer.analyze_signed_flow(slice); + + let dict = PyDict::new_bound(py); + dict.set_item("h0", metrics.h0)?; + dict.set_item("h_signed_fbm", metrics.h_signed_fbm)?; + dict.set_item("h_signed_mfbm", metrics.h_signed_mfbm)?; + dict.set_item("h_unsigned", metrics.h_unsigned)?; + dict.set_item("h_volatility", metrics.h_volatility)?; + dict.set_item("impact_exponent", metrics.impact_exponent)?; + dict.set_item("total_signed", metrics.total_signed)?; + dict.set_item("total_unsigned", metrics.total_unsigned)?; + dict.set_item("acf_1", metrics.acf_1)?; + + // Scale-dependent Hurst as nested dict + let scale_dict = PyDict::new_bound(py); + for (scale, h) in metrics.scale_hurst { + scale_dict.set_item(scale, h)?; + } + dict.set_item("scale_hurst", scale_dict)?; + + Ok(dict) +} + +/// Get unified theory derived quantities from H₀ +/// +/// # Arguments +/// * `h0` - Hurst index of signed order flow (typically ~0.75) +/// +/// # Returns +/// Dictionary with derived parameters: +/// - h0: Input H₀ +/// - alpha0: Tail exponent α₀ = H₀/2 +/// - h_volume: Volume Hurst H₁ = H₀ - 0.5 +/// - h_volatility: Volatility Hurst = 2H₀ - 1.5 +/// - impact_exponent: Market impact exponent δ = 2 - 2H₀ +/// - is_semimartingale: Whether mfBM is a semimartingale (H₀ > 3/4) +#[pyfunction] +pub fn unified_theory_params<'py>( + py: Python<'py>, + h0: f64, +) -> PyResult> { + let params = UnifiedTheoryParams::new(h0); + + let dict = PyDict::new_bound(py); + dict.set_item("h0", params.h0)?; + dict.set_item("alpha0", params.alpha0())?; + dict.set_item("h_volume", params.volume_hurst())?; + dict.set_item("h_volatility", params.volatility_hurst())?; + dict.set_item("impact_exponent", params.impact_exponent())?; + dict.set_item("is_semimartingale", params.is_semimartingale())?; + + Ok(dict) +} + +/// Compute market impact for given order size +/// +/// Impact(Q) = scale * |Q|^δ * sign(Q) +/// where δ = 2 - 2*H₀ +#[pyfunction] +#[pyo3(signature = (q, h0=0.75, scale=1.0))] +pub fn market_impact<'py>( + _py: Python<'py>, + q: f64, + h0: f64, + scale: f64, +) -> PyResult { + let impact = MarketImpact::from_h0(h0, scale); + Ok(impact.impact(q)) +} + +/// Register all point process functions to PyO3 module +pub fn register_python_functions(m: &Bound<'_, PyModule>) -> PyResult<()> { + // Hawkes processes + m.add_function(wrap_pyfunction!(simulate_hawkes, m)?)?; + m.add_function(wrap_pyfunction!(simulate_bivariate_hawkes, m)?)?; + + // Special functions + m.add_function(wrap_pyfunction!(mittag_leffler_py, m)?)?; + m.add_function(wrap_pyfunction!(f_alpha_lambda_py, m)?)?; + + // Fractional Brownian motion + m.add_function(wrap_pyfunction!(simulate_fbm, m)?)?; + m.add_function(wrap_pyfunction!(simulate_mixed_fbm, m)?)?; + m.add_function(wrap_pyfunction!(estimate_hurst, m)?)?; + m.add_function(wrap_pyfunction!(scale_dependent_hurst, m)?)?; + + // Order flow analysis + m.add_function(wrap_pyfunction!(analyze_order_flow, m)?)?; + m.add_function(wrap_pyfunction!(unified_theory_params, m)?)?; + m.add_function(wrap_pyfunction!(market_impact, m)?)?; + + Ok(()) +}