feat: add the Heston SLV process, its leverage calibration and Python bindings

[dancixx]

Summary

The Heston stochastic-local volatility model was only half present: a Monte Carlo pricer plus a free-function leverage calibration, no ProcessExt process (the processes page already listed a HestonSlv that did not exist), no Calibrator, no Python surface, and a particle method that deviated from its reference (Silverman's rule on the log-spot sample, one Euler step per local-vol maturity, the middle strike used as the spot). This PR completes the model and brings the calibration onto the literature.

Model

• HestonSlv process in stochastic-rs-stochastic::volatility::heston_slv: log-spot Euler with the crate's absorption scheme for the variance, the leverage an Fn2D of (t, S). Written as an Expr it runs on a device through the new Euler family HestonSlv (code 120, program pv); a closure, a Python callable or a tabulated grid keeps it on the host through device_fallback. With eta = 1 and L ≡ 1 it reproduces HestonLog bit for bit at the same seed, which a test pins.
• Grid2D and Fn2D::Grid in stochastic-rs-distributions: a (t, x) lookup table, bilinear inside and flat outside, so a calibrated LeverageSurface converts into a process coefficient (From<LeverageSurface> for Fn2D).

Calibration

• calibrate_leverage now follows Guyon & Henry-Labordère exactly, with the constants as restated by Cozma, Mariapragassam & Reisinger (arXiv:1701.06001 §3.4): Gaussian kernel in spot, bandwidth 1.5 · S₀ · σ_LV(S₀, t) · √max(t, 0.25) · N^(−1/5), the sorted-cloud window ΔK, L(s, 0) = σ_LV(s, 0) / √v₀, a leverage row at every Euler step, and the cloud's spots at every requested maturity so the caller can reprice the grid without a second simulation.
• calibrate_leverage_fokker_planck and heston_slv_density: the finite-volume solution of the forward Kolmogorov equation after Wyns & Du Toit (arXiv:1611.02961) — vertex-centred fluxes with zero boundary flux so the mass is conserved, the four-corner mixed flux with ghost values and the first-order fix above v = 0, the cell-average Dirac start on meshes clustered at the spot and at v₀, Hundsdorfer–Verwer at θ = ½ + √3/6 with Rannacher half-steps, and the inner iteration on the non-linearity (Q = 2).
• HestonSlvCalibrator (Calibrator, CalibrationResult, ToModel): a call-price grid in, a Heston fit (or pinned parameters), a Dupire local volatility (or a supplied grid), the leverage by either LeverageMethod, and the fit repriced from the calibration's own density. to_model builds a HestonSlvPricer anchored to the calibration rates and refuses others; the mixing fraction eta is an input, as in the literature.
• HestonSlvPricer reads the leverage at each step's start rather than its end.

Python

PyHestonSlv, LeverageSurface, HestonSlvPricer, HestonSlvCalibrator (with method="particle" or "fokker_planck") and HestonSlvCalibrationResult, with smoke tests.

Docs

New sections on the quant and processes pages backed by doctests, the GPU support matrix and counts, CLAUDE.md counts, the CUDA notebook cell, and a design spec kept locally under docs/ as the earlier ones are.

Breaking changes (pre-release API)

• calibrate_leverage(params, s0, r, q, &Grid2D, &maturities, &ParticleMethod) -> Result<ParticleCalibration> replaces the old grid-and-seed signature and returns the cloud snapshots beside the surface.
• calibrate_from_dupire is removed; HestonSlvCalibrator is the market-data entry point and takes the spot explicitly.
• HestonSlvCalibrator::method is a LeverageMethod.
• Surfaces calibrated by the previous particle code do not reproduce the new numbers.

Verification

• Gyöngy check: a Heston vanilla surface under its own parameters, through Dupire, calibrated at eta = 1, gives a median |L − 1| below 0.05 over the central strikes and reprices the surface with an RMSE below 0.25 on a spot of 100, by both the particle and the finite-volume route; eta = 0.5 reprices as well.
• Finite-volume solver under L ≡ 1: worst call-price error 0.003 against the closed-form Heston calls at the default mesh, mass conserved to 1e-11, the forward to 1e-6.
• eta = 0 at the long-run variance: the particle leverage equals the closed form σ_LV / √v₀ to 1e-9.
• HestonSlvPricer under L ≡ 1 matches the closed-form Heston call within three standard errors plus the Euler bias.
• Reproducibility guard (132 processes), Metal family parity, the device-law case (sabotage-verified: a term scaled by three fails it), clippy on host / metal / cuda / python, rustdoc, --no-default-features, the docs audit, the full workspace battery (3143 passed, 0 failed) before the Fokker–Planck commit and the quant crate suite after it, and the Python test suite.

Not done, left as follow-ups: Cozma et al.'s control variates for the particle method, a device path for a tabulated leverage (a 2-D table block in the engine), and a Python Expr builder so the process can reach a device from Python.

Validation round 2

• Independent cross-check against QuantLib 1.43 (HestonSLVFDMModel on its own NoExceptLocalVolSurface of the analytic Heston Black surface, the same local volatility handed to our routes so only the leverage calibration is compared): the finite-volume leverage agrees with QuantLib's to a mean of 0.001–0.006 and a worst of 0.003–0.010 over a 4 × 5 probe grid, the particle leverage to a mean of 0.003–0.005 and a worst of 0.009–0.012, at η = 1 and η = 0.5. QuantLib's own FDM and MC models agree with each other to a mean of 0.005–0.010 and a worst of 0.010–0.024 on the same probes. Vanilla prices under the calibrated models: QuantLib's FD engine within 0.035 of the market, our Monte Carlo pricer within 0.04 (200k paths). Kept as stochastic-rs-py/tests/test_slv_quantlib.py, skipped unless QuantLib is installed.
• The finite-difference Dupire read of a call grid is the dominant error: against QuantLib's local volatility it carries 1–3 %, and the in-sample RMSE is 0.036 from that read against 0.004 from the smooth surface. Documented on with_local_vol, which is the recommended input where a smooth surface exists.
• Reference parameter sets of Wyns & Du Toit: set C (Feller satisfied, vol-of-vol 0.9) and set D (Feller violated) reprice the closed-form Heston calls to 0.007 and 0.011 at the default mesh with the mass conserved; a test pins twice those errors.
• η = 0 on the PDE route now takes the closed form L = σ_LV / √V_t exactly instead of estimating a deterministic variance from the density.
• SSVI market route: an arbitrage-free SSVI surface's calls and its own local volatility, under a Heston model that does not generate it, reprice through both routes (RMSE below 0.2, the two leverage surfaces within 0.1 at the money).
• Process versus pricer: a vanilla priced from HestonSlv paths under a calibrated surface agrees with HestonSlvPricer within the two Monte Carlo errors.
• Not validated: the CUDA path (compiled only; no GPU here) and real quoted market data.

[vercel]

The latest updates on your projects. Learn more about Vercel for GitHub.

Project	Deployment	Actions	Updated
stochastic-rs	Ready	Preview	Sep 22, 2026 3:32pm UTC