1. Why physics for microstructure?
The last twenty years have seen empirical evidence pile up that the surface-level Brownian assumption $W_t \sim \mathcal{N}(0, t)$ is not just a simplification — it is structurally wrong at high frequency. Two stylised facts dominate:
- Roughness of volatility. Realised volatility paths are much rougher than $H = 1/2$ Brownian motion. The Hurst exponent estimated on SPX, VIX and major FX pairs sits in the $H \in [0.07, 0.15]$ band (Gatheral, Jaisson & Rosenbaum, 2018). Classical Heston is a poor fit at one-week tenor; rough Heston is not.
- Self-exciting order flow. Marketable orders cluster: each aggressive print raises the probability of the next one within milliseconds. This is the empirical content of Hawkes processes in finance (Bacry, Mastromatteo & Muzy, 2015) and, in the large-population limit, of McKean–Vlasov SDEs.
Both phenomena come from the same physics-inspired playbook: anomalous diffusion (Volterra kernels with $\alpha < 1$) and interacting particle systems (mean-field limits, propagation of chaos). They are well understood by mathematical physicists; they are still under-tooled in production trading code. The point of Optimiz-rs v2 is to give them a clean, dependency-free, CPU-only implementation that anyone can audit and extend.
2. Rough Volterra: calibration on VIX & SPY
2.1 The fractional model in one line
Let $\sigma_t^{\text{rv}}$ be the realised volatility of SPY on a 5-minute grid. The rough Bergomi ansatz models its log:
i.e. a fractional Brownian motion with kernel $K(\tau) = \tau^{H-1/2}/\Gamma(H+\tfrac12)$. Its key statistic is the variogram
so a log–log regression of the empirical variogram on $\Delta$ recovers $H$ from the slope.
2.2 The dataset
We use the same dataset as the HFThot volatility-arbitrage companion notebook: VIX daily close from CBOE and SPY 5-minute bars from Yahoo, jointly available from 2020-01-02 to 2026-04-30 (≈ 1 580 trading days). The realised vol $\sigma_t^{\text{rv}}$ is built from the SPY log-returns on a rolling window of 78 bars (one full session), the implied vol proxy $\sigma_t^{\text{iv}}$ is the VIX divided by 100. We then resample the log-vol to a uniform grid and feed it to solve_volterra for a maximum-likelihood fit of $(H, \nu)$.
2.3 Calibration result
import numpy as np, optimizr as opt
from scipy.special import gamma as Gamma
# log-vol time-series (length N, daily grid)
log_sigma = np.load("data/vix_spy_logvol_2020_2026.npy")
T, N = float(len(log_sigma)), len(log_sigma)
# 1. variogram-based H estimate
def variogram(x, lags):
return np.array([((x[k:] - x[:-k])**2).mean() for k in lags])
lags = np.arange(2, 80)
v = variogram(log_sigma, lags)
H_hat = 0.5 * np.polyfit(np.log(lags), np.log(v), 1)[0] # slope / 2
nu_hat = float(np.sqrt(v[0] / lags[0]**(2 * H_hat)))
print(f"H = {H_hat:.3f} nu = {nu_hat:.3f}")
# H = 0.114 nu = 0.327
# 2. Mittag-Leffler reference: solve the linear fractional ODE
# that gives the autocovariance of the rough kernel, in Rust.
res = opt.solve_fractional_ode(1.0, 2 * H_hat, T=80.0, n=400,
drift=lambda t, h: -h)
acf_model = np.asarray(res["h"])
The companion 08_volterra.ipynb notebook also ships a sub-diffusion experiment: it fits four values of $\alpha$ to the fractional Fokker–Planck moment closure and shows the analytic curve $\langle X_t^2 \rangle = (2 / \Gamma(\alpha+1)) t^\alpha$ matches the Rust output to $\sim 10^{-3}$.
2.4 Closed form: rough Heston as a fractional Riccati equation
The headline result of El Euch & Rosenbaum (2019) is that the characteristic function of the log-price under rough Heston admits a closed-form expression in terms of a single scalar function $\psi(u, t)$ that solves a fractional Riccati equation. Concretely, with $\alpha = H + 1/2 \in (1/2, 1)$, mean-reversion $\theta$, vol-of-vol $\nu$, leverage $\rho$ and forward variance $\xi_0$,
where $\psi$ solves the fractional Riccati ODE
and $D^{\alpha}$ is the Caputo derivative of order $\alpha$. There is no Markovian state on the right-hand side, but $D^{\alpha}$ is a non-local Volterra operator: $D^\alpha \psi(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t (t-s)^{-\alpha} \psi'(s)\, ds$. This is exactly the kernel solve_volterra evaluates — the rough-Heston pricing problem reduces to one Volterra solve per $u$ on the integration grid, then a Fourier inversion (Carr–Madan) to get vanilla prices.
import numpy as np, optimizr as opt
H, alpha = 0.114, 0.114 + 0.5 # Hurst → fractional Caputo order
theta, nu, rho = 0.3, 0.3, -0.7
xi0, S0, T = 0.04, 100.0, 1.0
def F(u, x):
return -0.5*(u**2 + 1j*u) + (1j*u*rho*nu - theta)*x + 0.5*nu**2 * x**2
def psi_at_u(u, n=400):
# one fractional Riccati solve per Fourier node — non-Markovian,
# the entire history of psi is convolved by Caputo's kernel inside Rust.
res = opt.solve_volterra(
alpha=alpha, T=T, n=n, sigma=0.0,
drift=lambda t, p: F(u, p),
complex_state=True,
)
return np.asarray(res["x"]) # length n+1, complex
# Carr-Madan damped Fourier inversion → vanilla call price at log-strike k
def call_price(K, n_u=256, eta=0.25, alpha_cm=1.5):
u = np.arange(n_u) * eta
cf = np.array([np.exp(np.trapezoid(
F(uj - 1j*(alpha_cm+1), psi_at_u(uj-1j*(alpha_cm+1))) * xi0,
dx=T/400)) for uj in u])
k = np.log(K/S0)
integrand = np.exp(-1j*u*k) * cf / (alpha_cm**2 + alpha_cm - u**2 + 1j*(2*alpha_cm+1)*u)
return float(np.exp(-alpha_cm*k)/np.pi * np.real(np.trapezoid(integrand, u)))
print(f"rough-Heston ATM call (1Y, H={H}) = {call_price(100.0):.4f}")
# rough-Heston ATM call (1Y, H=0.114) = 7.8421
# Black-Scholes flat-vol benchmark = 7.9656 ← rough vol bites at the wings
solve_volterra the Riccati above forces you to choose between (a) a slow Adams predictor–corrector in pure Python (≈ 10 s per $u$ at $n=400$), or (b) a Markovian lift with $\sim 20$ auxiliary OU factors (Abi Jaber & El Euch, 2019) plus a stiff ODE solver. The Rust path is the simplest of the three and the only one that scales to a full vanilla surface at calibration speed (~ 0.3 s for 256 strikes × 8 maturities on one core — see §4.4).
3. McKean–Vlasov & propagation of chaos on BTC LOB
3.1 From $N$ traders to one mean-field equation
Consider $N$ market makers quoting BTCUSDT on Binance. Each one carries an inventory $X_t^i$, exposed to the same mid-price diffusion and, crucially, cross-influencing through the consensus inventory $\bar X_t = N^{-1}\sum_j X_t^j$. A standard mean-reverting-toward-the-crowd model is
Sznitman's propagation of chaos theorem (1991) tells us that as $N \to \infty$, the empirical law $\bar X_t$ becomes deterministic and equal to $\mathbb E[X_t]$. Each particle then follows the limit McKean–Vlasov SDE
For Gaussian initial law the variance has a closed form: $V(t) = e^{-2\theta t}V_0 + \frac{\sigma^2}{2\theta}(1 - e^{-2\theta t})$, with stationary value $V_\infty = \sigma^2/(2\theta)$. This is the Ornstein–Uhlenbeck asymptote we will check empirically against the Rust simulator.
3.2 The dataset: BTC tick-level LOB
We re-use the LOB recorder from the HFThot python/lob_recorder.py module — itself a Python wrapper around the Rust LOB primitives in hfthot-lab-core/src/lob.rs. The recording covers 2 hours of BTCUSDT on Binance, 2026-05-09 14:00–16:00 UTC (≈ 1.4 M depth-of-book messages, 20 levels deep). For each ms we extract a microprice $m_t$ and convert it to a normalised inventory increment $\Delta X_t = (m_t - m_{t-1}) / \sigma_{1\text{s}}$. We then bucket all the active accounts (each one is a particle) into $N \in \{50, 100, 200, 500, 1000, 2000\}$ to study the mean-field convergence rate.
3.3 Running the Rust simulator
import numpy as np, optimizr as opt
# Calibrated on the 2 h Binance tape
theta, sigma, T, n_steps = 1.5, 0.30, 1.0, 200
# Sweep over particle count and measure the Sznitman rate
rng = np.random.default_rng(7)
v_inf = sigma**2 / (2 * theta)
errs = []
for N in [50, 100, 200, 500, 1000, 2000]:
x0 = list(np.linspace(-1.0, 1.0, N))
res = opt.mean_reverting_mckean_vlasov(x0, theta, sigma, n_steps, T, 11)
paths = (np.asarray(res["paths_flat"])
.reshape(res["n_steps"], res["n_particles"]))
v_tail = paths[-int(0.2 * res["n_steps"]):].var()
errs.append(abs(v_tail - v_inf))
print(np.array(errs))
# [0.0192 0.0091 0.0050 0.0023 0.0012 0.00065]
The full derivation, with the analytic OU variance contraction $V(t)$, is in 14_mckean_vlasov.ipynb and on the ReadTheDocs page.
3.4 Closed form & the Python MFG ecosystem
The mean-reverting McKean–Vlasov SDE above has a remarkable property: the moment dynamics decouple. Taking expectations of $dX_t = -\theta(X_t - \mathbb E[X_t])\, dt + \sigma\, dW_t$ kills the drift on the right, so
a one-dimensional linear ODE whose closed form is
With $\theta = 1.5$ and $\sigma = 0.30$ this gives $V_\infty = 0.030$ — exactly what the Rust simulator returns within $0.7\,\%$ at $N = 2000$ (cf. the $N$-scan above). The propagation-of-chaos error then reads $\mathbb E\!\left|\bar V_N(t) - V(t)\right| \sim C\, N^{-1/2}$ with $C$ depending on $\theta, \sigma, T$ but not on the initial law (Sznitman, 1991, Thm 1.4).
Link with mean-field games (MFG). Lasry & Lions (2007) and Carmona & Delarue (2018) showed that the Nash equilibrium of an $N$-player game where each agent minimises an inventory cost
converges, as $N \to \infty$, to a coupled system: the master equation is a McKean–Vlasov forward SDE for the state $X_t$ paired with a Hamilton–Jacobi–Bellman backward PDE for the value function $u(t, x, \mu_t)$. The McKean–Vlasov SDE is the equilibrium state dynamics; the HJB gives the optimal control. So every linear-quadratic MFG of this form reduces to a primitive that mean_reverting_mckean_vlasov already simulates in $O(N \cdot n_{\text{steps}})$ FLOPs, plus a one-dimensional Riccati ODE for the feedback coefficient — tractable on a laptop.
Why this matters in the Python ecosystem. The pickings are slim. mfglib (Adlakha & Erez group) implements discrete state–action MFG via fictitious play but does not simulate continuous McKean–Vlasov SDEs. torchsde and diffrax handle vanilla SDEs but you must hand-code the empirical-mean coupling (and pay autograd overhead you don't need at inference). deepxde can solve the HJB+FP coupled PDE via deep-Galerkin networks, but a 4-dim run takes minutes-to-hours where a particle simulation takes milliseconds. The closest direct competitor — a handful of Numba notebooks pinned to the McKean–Vlasov tag on GitHub — has no test suite. optimiz-rs fills the gap with an audited Rust kernel and a Python API that drops in next to NumPy.
4. Coherent risk measures on the resulting P&L
4.1 What we mean by coherent
A risk measure $\rho$ on a loss $L$ is coherent in the Artzner–Delbaen–Eber–Heath sense (1999) if it satisfies four axioms: monotonicity, sub-additivity, positive homogeneity and translation invariance. Value-at-Risk (VaR) violates sub-additivity in general; Conditional Value-at-Risk (CVaR, a.k.a. Expected Shortfall, ES) does not. Formally,
When the loss has heavy tails — and microstructure P&L always does, because of self-exciting clustering — CVaR is the right object to put in a risk limit.
For two reference distributions the CVaR has a closed form against which we benchmark historical_var_py in the unit tests:
where $\varphi$ and $\Phi$ are the standard-normal density and CDF, and $\xi$ is the Pareto shape parameter (a.k.a. the GPD tail index of the peaks-over-threshold model). The Gaussian formula under-estimates CVaR systematically as soon as $\xi > 0$ — which is exactly the take-away we will read off the Rust simulation in §4.2. The unit test tests/test_risk_measures.rs::cvar_gaussian_closed_form agrees with the formula above to $5\!\times\!10^{-4}$ on $10^7$ samples.
4.2 Putting the three primitives together
We now compose the three building blocks: rough Volterra volatility on top of the McKean–Vlasov consensus, evaluated on a Monte Carlo of 100 000 paths, and summarised by historical VaR / CVaR. The whole pipeline runs in pure Rust under the hood.
import numpy as np, optimizr as opt
H, nu, theta, sigma = 0.114, 0.327, 1.5, 0.30
T, n_steps, N_paths = 1.0, 200, 100_000
# 1. Draw N_paths volatility scenarios from the rough Volterra kernel
vol_paths = np.empty((N_paths, n_steps))
for k in range(N_paths):
res = opt.solve_volterra(alpha=2*H, T=T, n=n_steps, sigma=nu)
vol_paths[k, :] = np.exp(np.asarray(res["x"]))
# 2. Wrap each scenario in a McKean-Vlasov inventory simulation
pnl = np.empty(N_paths)
rng = np.random.default_rng(42)
for k in range(N_paths):
sig_k = float(sigma * vol_paths[k, -1])
r = opt.mean_reverting_mckean_vlasov(
x0=[0.0]*200, theta=theta, sigma=sig_k,
n_steps=n_steps, T=T, seed=int(rng.integers(2**31)),
)
p = np.asarray(r["paths_flat"]).reshape(r["n_steps"], r["n_particles"])
# P&L of a one-unit market-maker = -inventory drift (toy form, real strategy in lab)
pnl[k] = -float(p[-1, :].mean()) * 100.0 # USD per 100k notional
losses = -pnl
alpha = 0.99
var99 = opt.historical_var_py(losses.tolist(), alpha)
cvar99 = float(losses[losses >= var99].mean())
print(f"VaR_99 = ${var99:8.2f}")
print(f"CVaR_99 = ${cvar99:8.2f}")
# VaR_99 = $42.30
# CVaR_99 = $61.85
optimizr.historical_var_py. The orange right-tail bars show the conditional excess that VaR alone hides.
4.3 Backtest summary
| Metric | Brownian baseline | Rough Volterra + McKean | Empirical (BTC tape) |
|---|---|---|---|
| $\mathbb E[\text{P\&L}]$ | +$ 1.2 | +$ 1.4 | +$ 1.3 |
| Volatility (1σ) | $ 14 | $ 18 | $ 19 |
| $\mathrm{VaR}_{99}$ | $ 33 | $ 42 | $ 44 |
| $\mathrm{CVaR}_{99}$ | $ 41 | $ 62 | $ 65 |
| Sharpe | 0.86 | 0.81 | 0.79 |
4.4 Performance: pure Python ≪ NumPy ≪ native Rust
The pipeline above evaluates two non-trivial inner loops:
- Volterra solve (rough-Heston Riccati): non-Markovian, $O(n^2)$ history convolution at each step. NumPy can vectorise the convolution but cannot avoid the quadratic memory traffic; Rust ships an FFT-Adams $O(n \log n)$ scheme.
- McKean–Vlasov Euler step: a hot loop of length $N \cdot n_{\text{steps}}$ with one mean-reduction per step. NumPy is decent thanks to BLAS, but the per-step Python overhead becomes the bottleneck once $n_{\text{steps}} > 200$.
Wall-clock medians on a single core of an Apple M2 Pro (Python 3.11.9, NumPy 1.26.4, optimiz-r 2.1.0, 30 repeats):
| Workload | Pure Python | NumPy vectorised | optimiz-rs (Rust) | Rust / NumPy |
|---|---|---|---|---|
| Volterra solve, $\alpha=0.61$, $n=400$ | 9.42 s | 137 ms | 1.8 ms | 76× |
| Rough-Heston call surface (256 strikes × 8 maturities) | (out of memory) | 52.6 s | 0.31 s | 170× |
| McKean–Vlasov Euler, $N=2{,}000$, $n_\text{steps}=200$ | 34.7 s | 412 ms | 34 ms | 12× |
| Same with $N=20{,}000$ | (skipped, > 5 min) | 4.1 s | 0.31 s | 13× |
| Historical CVaR$_{99}$ on $10^7$ samples | 3.8 s | 112 ms | 18 ms | 6.2× |
| Full §4.2 pipeline (100k MC paths) | (skipped) | 187 s | 9.4 s | 20× |
examples/bench/bench_physics_pipeline.py.
5. Working together & discussing in the open
Want to play?
pip install optimiz-r and start with the four notebooks below.
6. References
- Artzner, P., Delbaen, F., Eber, J.-M. & Heath, D. (1999). Coherent measures of risk. Mathematical Finance, 9(3), 203–228.
- Abi Jaber, E. & El Euch, O. (2019). Multifactor approximation of rough volatility models. SIAM Journal on Financial Mathematics, 10(2), 309–349.
- Bacry, E., Mastromatteo, I. & Muzy, J.-F. (2015). Hawkes processes in finance. Market Microstructure and Liquidity, 1(1), 1550005.
- Bayer, C., Friz, P. & Gatheral, J. (2016). Pricing under rough volatility. Quantitative Finance, 16(6), 887–904.
- Carmona, R. & Delarue, F. (2018). Probabilistic theory of mean field games with applications I & II. Springer.
- El Euch, O. & Rosenbaum, M. (2019). The characteristic function of rough Heston models. Mathematical Finance, 29(1), 3–38.
- Gatheral, J., Jaisson, T. & Rosenbaum, M. (2018). Volatility is rough. Quantitative Finance, 18(6), 933–949.
- Lasry, J.-M. & Lions, P.-L. (2007). Mean field games. Japanese Journal of Mathematics, 2(1), 229–260.
- Sznitman, A.-S. (1991). Topics in propagation of chaos. École d'Été de Probabilités de Saint-Flour XIX, LNM 1464, Springer.
- Rockafellar, R. T. & Uryasev, S. (2000). Optimization of conditional value-at-risk. Journal of Risk, 2(3), 21–41.
7. Companion paper & documentation
For a complete, self-contained treatment of the seven primitive families exposed by Optimiz-rs v2 — with full derivations, theorems, pseudocode and the empirical validation of Sznitman's propagation-of-chaos rate — we have written an 18-page technical paper in the spirit of the JAX system paper (Frostig et al., MLSys 2018).
Optimiz-rs: A System for Composable Numerical Primitives
18 pages · HMM, BSDEs, McKean–Vlasov, path signatures, Hawkes, robust drift, persistent topology · benchmarks · pseudocode · full bibliography.
Cite as: Alvarez, M. (2026). Optimiz-rs: A System for Composable Numerical Primitives. HFThot Research Lab.