From 6c2ce3d238ced9fc719892fa88400568b9a78844 Mon Sep 17 00:00:00 2001 From: ThotDjehuty Date: Thu, 19 Feb 2026 19:35:04 +0100 Subject: [PATCH] refactor: remove finance-specific code from point_processes module - Removed order_flow.rs (UnifiedTheoryParams, OrderFlowAnalyzer, MarketImpact) - Removed analyze_order_flow, unified_theory_params, market_impact Python bindings - Updated mod.rs documentation to be generic (no trading references) - Kept general-purpose: Hawkes, fBM, mfBM, Mittag-Leffler, Hurst estimation Finance-specific code moved to rust-hft-arbitrage-lab-internal --- LINKEDIN_POST.md | 68 --- .../algorithms/optimal_control.md.backup | 94 ++++ docs/source/api/hmm.md.backup | 16 + docs/source/api/optimal_control.md.backup | 117 +++++ src/point_processes/mod.rs | 68 +-- src/point_processes/order_flow.rs | 452 ------------------ src/point_processes/python_bindings.rs | 108 +---- 7 files changed, 274 insertions(+), 649 deletions(-) delete mode 100644 LINKEDIN_POST.md create mode 100644 docs/source/algorithms/optimal_control.md.backup create mode 100644 docs/source/api/hmm.md.backup create mode 100644 docs/source/api/optimal_control.md.backup delete mode 100644 src/point_processes/order_flow.rs diff --git a/LINKEDIN_POST.md b/LINKEDIN_POST.md deleted file mode 100644 index 3641dbe..0000000 --- a/LINKEDIN_POST.md +++ /dev/null @@ -1,68 +0,0 @@ -🚀 OptimizR v1.0.0 is LIVE! 🚀 - -I'm thrilled to announce the first stable release of OptimizR - a high-performance Rust library for optimization and statistical inference! 🎉 - -📦 What is OptimizR? - -OptimizR brings blazing-fast implementations of advanced algorithms to data science: -• Differential Evolution (5 strategies with adaptive jDE) -• Hidden Markov Models (Baum-Welch, Viterbi) -• MCMC Sampling (Metropolis-Hastings) -• Mean Field Games (PDE solvers) -• Information Theory tools - -⚡ The Performance Story: -• 50-100× faster than pure Python -• 95% memory reduction vs NumPy/SciPy -• Production-ready with comprehensive testing - -🔥 Real-World Benchmarks: -• HMM Training: 0.03s (Rust) vs 2.4s (Python) = 80× speedup -• Differential Evolution: 0.12s (Rust) vs 8.9s (SciPy) = 74× speedup -• Mean Field Games: 0.4s (Rust) vs 45s (Python) = 112× speedup - -📚 Learn More: -→ Complete documentation: https://optimiz-r.readthedocs.io -→ 6+ working tutorial notebooks -→ API reference with examples -→ Mathematical foundations - -🛠️ Get Started: - -Python: -```bash -pip install optimiz-rs -``` - -Rust: -```bash -cargo add optimiz-rs -``` - -🌟 Why OptimizR? - -Built for researchers and practitioners who need: -✅ Production-grade performance -✅ Stable, well-documented API -✅ Easy Python integration -✅ Robust numerical methods -✅ Active development roadmap - -🤝 Join the Community! - -This project is driven by the belief that high-performance optimization tools should be accessible to everyone. Whether you're: -• A data scientist optimizing portfolios -• A researcher working on Bayesian inference -• An engineer building ML pipelines -• Or just curious about Rust + Python! - -I'd love your feedback, contributions, and star on GitHub! ⭐ - -📖 Explore the docs, try the tutorials, and let me know what you build with it! - -#Rust #Python #DataScience #MachineLearning #Optimization #OpenSource #MCMC #HiddenMarkovModels #NumericalMethods #HighPerformanceComputing - -🔗 GitHub: https://github.com/ThotDjehuty/optimiz-r -📚 Docs: https://optimiz-r.readthedocs.io -📦 PyPI: https://pypi.org/project/optimiz-rs/ -📦 crates.io: https://crates.io/crates/optimiz-rs diff --git a/docs/source/algorithms/optimal_control.md.backup b/docs/source/algorithms/optimal_control.md.backup new file mode 100644 index 0000000..634d016 --- /dev/null +++ b/docs/source/algorithms/optimal_control.md.backup @@ -0,0 +1,94 @@ +# Optimal Control + +Hamilton–Jacobi–Bellman (HJB) solvers, regime-switching thresholds, OU parameter estimation, and Kalman filtering utilities backed by Rust. + +## HJB switching boundaries (OU process) + +```python +from optimizr import solve_hjb_py, solve_hjb_full_py + +lower, upper, residual, iters = solve_hjb_py( + kappa=3.0, + theta=0.0, + sigma=0.2, + rho=0.04, + transaction_cost=0.001, + n_points=400, + max_iter=4000, + tolerance=1e-7, + n_std=5.0, +) +print(f"bounds=({lower:.3f}, {upper:.3f}), residual={residual:.2e}, iters={iters}") +``` + +- Model: $dX_t = \kappa(\theta - X_t)\,dt + \sigma\,dW_t$ with quadratic transaction costs. +- Output: optimal buy/sell thresholds; `solve_hjb_full_py` also returns $V, V_x, V_{xx}$ for diagnostics. +- Diagnostics: plot $V_x$ for smoothness near thresholds; monitor `residual` and increase `max_iter` if not converged. + +## Backtesting optimal switching + +```python +from optimizr import backtest_optimal_switching_py + +metrics = backtest_optimal_switching_py( + spread=spread, + lower_bound=lower, + upper_bound=upper, + transaction_cost=0.001, +) +( + total_return, + sharpe, + max_dd, + n_trades, + win_rate, + pnl_path, +) = metrics +``` + +Inspect `win_rate` vs `max_dd` to tune aggressiveness; combine with HMM regimes for state-aware controls. + +## OU parameter estimation + +```python +import numpy as np +from optimizr import estimate_ou_params_py + +spread = np.random.randn(10_000) +kappa, theta, sigma, half_life = estimate_ou_params_py(spread, dt=1/252) +``` + +Method-of-moments / MLE fit returns $(\kappa, \theta, \sigma, \text{half-life})$. Use a few thousand samples for stability; winsorize heavy tails if needed. + +## Kalman filtering (linear, EKF, UKF) + +```python +import numpy as np +from optimizr import LinearKalmanFilter + +F = [[1.0, 1.0], [0.0, 1.0]] +H = [[1.0, 0.0]] +Q = [[1e-4, 0.0], [0.0, 1e-4]] +R = [[1e-2]] + +kf = LinearKalmanFilter( + f_matrix=F, + h_matrix=H, + q_matrix=Q, + r_matrix=R, + initial_state=[0.0, 0.0], + initial_covariance=[[1.0, 0.0], [0.0, 1.0]], +) + +kf.predict(control=[0.0, 0.0]) +kf.update(observation=[1.2]) +state = kf.get_state() +``` + +- Interfaces: `LinearKalmanFilter`, `UnscentedKalmanFilter`, and `KalmanState` for batch `filter` and smoothing. +- Concept: prediction (dynamics prior) + correction (measurement residual); RTS smoother refines past states. + +## Practical notes +- Rust backend (`optimizr._core`) must be present for control utilities; install from source if wheels are unavailable. +- Grids: for HJB, `n_points≈400` is stable; widen `n_std` for volatile spreads. +- Combine with Mean Field Games: see `mean_field_games.md` for population dynamics; use Kalman estimates as control inputs if needed. diff --git a/docs/source/api/hmm.md.backup b/docs/source/api/hmm.md.backup new file mode 100644 index 0000000..c45f3a1 --- /dev/null +++ b/docs/source/api/hmm.md.backup @@ -0,0 +1,16 @@ +# API: HMM + +```python +from optimizr import HMM + +model = HMM(n_states=2) +model.fit(X, n_iterations=100, tolerance=1e-6) +states = model.predict(X) +logp = model.score(X) +``` + +Parameters +- `n_states`: number of hidden regimes +- `fit(X, n_iterations=100, tolerance=1e-6)`: train with Baum-Welch +- `predict(X)`: Viterbi decoding → `np.ndarray` +- `score(X)`: log-likelihood diff --git a/docs/source/api/optimal_control.md.backup b/docs/source/api/optimal_control.md.backup new file mode 100644 index 0000000..cb1e772 --- /dev/null +++ b/docs/source/api/optimal_control.md.backup @@ -0,0 +1,117 @@ +# API: Optimal Control / Kalman + +High-level bindings exposed by the `optimizr` Python package. All functions require the Rust extension (`optimizr._core`). + +**When to use this module** +- Threshold trading / switching problems solved via HJB (with and without frictions) +- State estimation and smoothing (Kalman, EKF, UKF) +- Parameter inference for mean-reverting spreads (OU) feeding into control logic + +## Hamilton–Jacobi–Bellman (HJB) solvers + +```python +from optimizr import solve_hjb_py, solve_hjb_full_py + +# Switching boundaries for a mean-reverting spread (OU process) +lower, upper, residual, iters = solve_hjb_py( + kappa=3.0, theta=0.0, sigma=0.2, rho=0.04, + transaction_cost=0.001, n_points=400, max_iter=4000, + tolerance=1e-7, n_std=5.0, +) + +# Full state (grid + derivatives) for research/visualization +(lower, upper, residual, iters, x_grid, value, grad, hess) = solve_hjb_full_py( + kappa=3.0, theta=0.0, sigma=0.2, rho=0.04, + transaction_cost=0.001, n_points=400, +) +``` + +### Model + +We assume an Ornstein–Uhlenbeck process $dX_t = \kappa(\theta - X_t)\,dt + \sigma\,dW_t$ with quadratic transaction costs. The HJB on grid $x \in [-n_{std}\,\sigma/\sqrt{\kappa},\; n_{std}\,\sigma/\sqrt{\kappa}]$ solves +$$ +\rho V(x) = \min\Big\{ \tfrac12 \sigma^2 V_{xx}(x) + \kappa(\theta - x) V_x(x),\; V(x) + c_{\text{buy}},\; V(x) + c_{\text{sell}} \Big\}. +$$ +`solve_hjb_py` returns optimal buy/sell thresholds; `solve_hjb_full_py` also returns $V$, $V_x$, and $V_{xx}$ for diagnostics. + +**Diagnostic tips:** +- Plot $V_x$ to verify smoothness near the boundaries; kinks often signal insufficient grid resolution. +- Track `residual` and `iterations` to spot non-convergence; loosen `tolerance` or increase `max_iter` if needed. + +## OU parameter estimation + +```python +import numpy as np +from optimizr import estimate_ou_params_py + +spread = np.random.randn(10_000) +kappa, theta, sigma, half_life = estimate_ou_params_py(spread, dt=1/252) +``` + +Method-of-moments / MLE estimation for +$$ +X_{t+1} = X_t e^{-\kappa \Delta t} + \theta(1-e^{-\kappa \Delta t}) + \eta_t, \quad \eta_t \sim \mathcal{N}\Big(0,\; \tfrac{\sigma^2}{2\kappa}(1-e^{-2\kappa \Delta t})\Big). +$$ +Returns $(\kappa, \theta, \sigma, \text{half\_life})$. + +**Practical guidance:** Use at least a few thousand samples for stable estimates; heavy-tailed series benefit from pre-whitening or winsorizing before fitting. + +## Backtesting optimal switching + +```python +from optimizr import backtest_optimal_switching_py + +metrics = backtest_optimal_switching_py( + spread=spread, + lower_bound=lower, + upper_bound=upper, + transaction_cost=0.001, +) +(total_return, sharpe, max_dd, n_trades, win_rate, pnl_path) = metrics +``` + +Applies HJB thresholds to historical spreads and reports return, Sharpe ratio, drawdown, trade count, win rate, and PnL path. + +**What to inspect:** +- `win_rate` alongside `max_drawdown` to balance aggressiveness +- PnL path for regime shifts; combine with HMM states if you need regime-aware controls + +## Kalman filtering (linear, EKF, UKF) + +```python +import numpy as np +from optimizr import LinearKalmanFilter, KalmanState + +F = [[1.0, 1.0], [0.0, 1.0]] # constant-velocity model +H = [[1.0, 0.0]] # observe position only +Q = [[1e-4, 0.0], [0.0, 1e-4]] +R = [[1e-2]] + +kf = LinearKalmanFilter( + f_matrix=F, + h_matrix=H, + q_matrix=Q, + r_matrix=R, + initial_state=[0.0, 0.0], + initial_covariance=[[1.0, 0.0], [0.0, 1.0]], +) + +kf.predict(control=[0.0, 0.0]) # optional control input via B matrix +kf.update(observation=[1.2]) +state = kf.get_state() # KalmanState with getters for mean/cov + +# Batch filtering +result = kf.filter(observations=[[1.0], [1.4], [1.9]], controls=None) +states = result.get_states() +log_likelihoods = result.get_log_likelihoods() +``` + +### Notes +- `LinearKalmanFilter` implements `predict`, `update`, and batch `filter`. +- `KalmanState` exposes `get_state()`, `get_covariance()`, and `get_log_likelihood()`. +- Extended/Unscented Kalman filters share the same interface (see `UnscentedKalmanFilter` in the Rust module) and are exported through the same bindings. +- For smoothing, use the Rauch–Tung–Striebel smoother (`RTSSmoother`) available in the bindings. + +**Conceptual picture:** Kalman filtering = prediction (dynamics prior) + correction (measurement residual). EKF linearizes $f, h$; UKF propagates sigma points for better nonlinear fidelity. RTS smoothing runs backward in time to refine all past states. + +See [`examples/notebooks/03_optimal_control_tutorial.ipynb`](https://github.com/ThotDjehuty/optimiz-r/blob/main/examples/notebooks/03_optimal_control_tutorial.ipynb) for end-to-end usage combining HJB thresholds, OU estimation, and filtering. diff --git a/src/point_processes/mod.rs b/src/point_processes/mod.rs index 0a645f5..b697909 100644 --- a/src/point_processes/mod.rs +++ b/src/point_processes/mod.rs @@ -1,52 +1,67 @@ -//! Point Processes Module for Order Flow Modeling +//! Point Processes Module //! -//! This module implements the mathematical framework from "A Unified Theory of -//! Order Flow, Market Impact, and Volatility" (Muhle-Karbe et al., 2026). +//! This module provides general-purpose implementations of point processes, +//! fractional Brownian motion, and related mathematical tools. //! //! # Overview //! -//! The module provides tools for modeling order flow in financial markets using: -//! -//! - **Hawkes Processes**: Self-exciting point processes for core and reaction flows +//! - **Hawkes Processes**: Self-exciting point processes with flexible kernels //! - **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 +//! - **Fractional Brownian Motion (fBM)**: Long-memory Gaussian processes +//! - **Mixed Fractional Brownian Motion (mfBM)**: BM + fBM combinations +//! - **Mittag-Leffler Functions**: Special functions for scaling limit analysis +//! - **Hurst Estimation**: R/S analysis and scale-dependent methods //! -//! # Key Relationships (Unified Theory) +//! # Mathematical Background //! -//! All quantities are determined by a single parameter H₀ ≈ 3/4: +//! ## Hawkes Processes //! -//! - 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) +//! A Hawkes process is a self-exciting point process with intensity: +//! +//! $$\lambda(t) = \nu + \sum_{t_i < t} \phi(t - t_i)$$ +//! +//! where $\nu$ is the baseline intensity and $\phi$ is the excitation kernel. +//! +//! Supported kernels: +//! - **Exponential**: $\phi(t) = \alpha e^{-\beta t}$ (short memory) +//! - **Power-law**: $\phi(t) = K_0 (1+t)^{-(1+\alpha_0)}$ (long memory) +//! +//! ## Fractional Brownian Motion +//! +//! fBM $B^H_t$ with Hurst parameter $H \in (0,1)$ has covariance: +//! +//! $$Cov(B^H_s, B^H_t) = \frac{1}{2}(|t|^{2H} + |s|^{2H} - |t-s|^{2H})$$ +//! +//! - $H < 0.5$: Anti-persistent (mean-reverting) +//! - $H = 0.5$: Standard Brownian motion +//! - $H > 0.5$: Persistent (trending) //! //! # Example //! //! ```rust,ignore //! use optimizr::point_processes::{ -//! HawkesProcess, PowerLawKernel, OrderFlowAnalyzer +//! HawkesProcess, PowerLawKernel, FractionalBM //! }; //! -//! // Create a Hawkes process for core order flow +//! // Create a Hawkes process with power-law kernel //! let kernel = PowerLawKernel::new(0.5, 1.0); // α₀ = 0.5, K₀ = 1.0 -//! let hawkes = HawkesProcess::new(0.1, kernel); // ν = 0.1 +//! let hawkes = HawkesProcess::new(0.1, kernel); // baseline ν = 0.1 //! -//! // Simulate order arrivals +//! // Simulate point process //! 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); +//! // Simulate fBM with H = 0.75 +//! let mut fbm = FractionalBM::new(0.75); +//! let path = fbm.simulate_hosking(1000, 1.0, Some(42)); +//! +//! // Estimate Hurst exponent +//! let h_est = FractionalBM::estimate_hurst(&path); //! ``` mod kernels; mod hawkes; mod mittag_leffler; mod mixed_fbm; -mod order_flow; #[cfg(feature = "python-bindings")] pub mod python_bindings; @@ -56,8 +71,3 @@ 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 deleted file mode 100644 index 49c53f7..0000000 --- a/src/point_processes/order_flow.rs +++ /dev/null @@ -1,452 +0,0 @@ -//! 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 index 43980ea..df6667b 100644 --- a/src/point_processes/python_bindings.rs +++ b/src/point_processes/python_bindings.rs @@ -1,14 +1,20 @@ //! Python bindings for point processes module +//! +//! Provides Python access to: +//! - Hawkes process simulation (univariate and bivariate) +//! - Fractional Brownian motion simulation +//! - Mixed fractional Brownian motion +//! - Hurst exponent estimation +//! - Mittag-Leffler special functions 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::hawkes::{HawkesProcess, 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 /// @@ -184,99 +190,6 @@ pub fn scale_dependent_hurst<'py>( 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 @@ -293,10 +206,5 @@ pub fn register_python_functions(m: &Bound<'_, PyModule>) -> PyResult<()> { 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(()) }