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
This commit is contained in:
ThotDjehuty
@@ -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
|
||||
@@ -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.
|
||||
@@ -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
|
||||
@@ -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.
|
||||
+39
-29
@@ -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
|
||||
};
|
||||
|
||||
@@ -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<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);
|
||||
}
|
||||
}
|
||||
@@ -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<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
|
||||
@@ -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(())
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user