gavindiaz/optimiz-rs
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: add point_processes module implementing unified theory framework

Browse Source
This commit is contained in:
ThotDjehuty
2026-02-19 19:01:29 +01:00
parent fc8a8d036e
commit f7e62cbe6d
8 changed files with 2338 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/lib.rs
+4
View File
@@ -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)?;
src/point_processes/hawkes.rs
+462
View File
@@ -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<K: ExcitationKernel> {
/// 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<K: ExcitationKernel> HawkesProcessConfig<K> {
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<K: ExcitationKernel> {
config: HawkesProcessConfig<K>,
/// Event times
events: Vec<f64>,
/// Current time
current_time: f64,
}
impl<K: ExcitationKernel> HawkesProcess<K> {
/// 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<K>) -> 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<u64>) -> Vec<f64> {
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<f64> = 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<FitResult, &'static str> {
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::<f64>();
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<K: ExcitationKernel> {
/// 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<f64>,
/// Sell events N⁻
pub sell_events: Vec<f64>,
/// Current time
current_time: f64,
}
impl<K: ExcitationKernel> BivariateHawkes<K> {
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<u64>,
) -> (Vec<f64>, Vec<f64>) {
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::<f64>() < 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
}
}
src/point_processes/kernels.rs
+307
View File
@@ -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<f64>;
}
/// 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<f64> {
// 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<f64> {
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<f64> {
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());
}
}
src/point_processes/mittag_leffler.rs
+304
View File
@@ -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);
}
}
src/point_processes/mixed_fbm.rs
+444
View File
@@ -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<f64>,
/// Sample path values
pub values: Vec<f64>,
}
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<u64>) -> Vec<f64> {
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<f64> = (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<u64>) -> Vec<f64> {
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::<f64>();
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::<f64>() / sub_n as f64;
let var: f64 = subseries.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / 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::<f64>() / 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::<f64>() / n_points;
let mean_y: f64 = log_rs.iter().sum::<f64>() / 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<f64>,
/// Sample path values
pub values: Vec<f64>,
}
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<u64>) -> Vec<f64> {
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<f64> = (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::<f64>() / 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<f64> = (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::<f64>() / 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<f64>]) -> Option<Vec<Vec<f64>>> {
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);
}
}
src/point_processes/mod.rs
+63
View File
@@ -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
};
src/point_processes/order_flow.rs
+452
View File
@@ -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<usize>,
}
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::<f64>() / 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::<f64>() / 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<f64> = flow.iter().map(|x| x.abs()).collect();
let cum_unsigned: Vec<f64> = 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::<f64>() / n as f64;
let var: f64 = flow.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / 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::<f64>() / (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<f64> = (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::<f64>() / 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::<f64>() / n_pts;
let mean_y: f64 = log_vars.iter().sum::<f64>() / 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);
}
}
src/point_processes/python_bindings.rs
+302
View File
@@ -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<u64>,
) -> PyResult<Bound<'py, PyArray1<f64>>> {
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::<pyo3::exceptions::PyValueError, _>(
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<f64>,
core_sell_times: PyReadonlyArray1<f64>,
phi1_alpha: f64,
phi1_beta: f64,
phi2_alpha: f64,
phi2_beta: f64,
t_max: f64,
seed: Option<u64>,
) -> PyResult<(Bound<'py, PyArray1<f64>>, Bound<'py, PyArray1<f64>>)> {
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<f64> {
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<f64> {
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<u64>,
) -> PyResult<Bound<'py, PyArray1<f64>>> {
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<u64>,
) -> PyResult<Bound<'py, PyArray1<f64>>> {
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<f64>,
) -> PyResult<f64> {
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<f64>,
scales: Option<Vec<usize>>,
) -> PyResult<Bound<'py, PyDict>> {
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<f64>,
) -> PyResult<Bound<'py, PyDict>> {
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<Bound<'py, PyDict>> {
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<f64> {
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(())
}