# Heston: simulate, price and calibrate URL: https://stochastic.rust-dd.com/docs/tutorials/heston > Simulate, price and calibrate the Heston stochastic volatility model in Rust and Python — Euler and QE paths, Fourier and ADI pricers, calibration. The Heston model is a stochastic volatility model in which the variance of an asset follows a mean-reverting square-root (CIR) process driven by a Brownian motion correlated with the asset's own. stochastic-rs simulates it under an Euler or Andersen's quadratic-exponential scheme on the CPU or a GPU, prices European options in semi-closed form, by Fourier inversion and by ADI finite differences, calibrates it to option prices and estimates it from time series, from Rust and from Python. ## Model $$ \begin{aligned} dS_t &= \mu S_t\,dt + \sqrt{v_t}\,S_t\,dW^S_t \\ dv_t &= \kappa(\theta - v_t)\,dt + \sigma\, v_t^{\,p}\,dW^v_t,\qquad d\langle W^S, W^v\rangle_t = \rho\,dt \end{aligned} $$ with $p = \tfrac12$ for the Heston model. `HestonPow::ThreeHalves` sets $p = \tfrac32$, a 3/2-type diffusion under the same linear drift. | Symbol | Argument | Meaning | | -------- | --------- | --------------------------------------------------------------------------------------------------------------- | | $S_0$ | `s0` | Initial price. Pass it: `None` starts an Euler path at 0 | | $v_0$ | `v0` | Initial variance (a variance, not a volatility), $\ge 0$; `None` is 0 | | $\kappa$ | `kappa` | Speed of mean reversion, $\ge 0$ ($> 0$ for the QE scheme) | | $\theta$ | `theta` | Long-run variance, $\ge 0$ | | $\sigma$ | `sigma` | Volatility of the variance, $\ge 0$ | | $\rho$ | `rho` | Correlation of $W^S$ and $W^v$, in $[-1, 1]$ | | $\mu$ | `mu` | Drift of the price; use $r - q$ for risk-neutral paths | | $p$ | `pow` | `HestonPow::Sqrt` ($p = \tfrac12$) or `HestonPow::ThreeHalves` ($p = \tfrac32$); Python `pow="sqrt"` or `"3/2"` | | | `n`, `t` | Grid points including $t = 0$, and the horizon (default 1); the step is $t/(n-1)$ | | | `use_sym` | `Some(true)` reflects the variance, $\lvert v \rvert$, instead of truncating it at 0 | | | `seed` | `Unseeded` or `Deterministic::new(seed)`; Python `seed=` | The constructor asserts the signs above and the range of `rho`. It does not check the Feller condition $2\kappa\theta \ge \sigma^2$: both variance schemes keep the simulated variance non-negative whether it holds or not. ## What the library has for it | Task | Rust | Python | | ------------------------------------------- | --------------------------------------------------------- | ------------------------------------ | | Simulate price and variance paths | `stochastic::volatility::heston::Heston` | `PyHeston` | | Simulate on the log-price | `stochastic::volatility::heston_log::HestonLog` | `PyHestonLog` | | Price European options, with Greeks | `quant::pricing::heston::HestonPricer` | `HestonPricer` | | Price by Fourier inversion | `quant::pricing::fourier::HestonFourier` | `HestonFourier` | | Price vanilla and down-and-out calls by ADI | `quant::pricing::heston_adi::HestonAdiPricer` | `HestonAdiPricer` | | Monte Carlo Malliavin Greeks | `quant::pricing::malliavin_greeks::HestonMalliavinGreeks` | — | | Calibrate to option prices | `quant::calibration::heston::HestonCalibrator` | `HestonCalibrator` | | Estimate from price and variance series | `stats::heston_mle::nmle_heston`, `pmle_heston` | `HestonMLE` | | Estimate from prices alone | `stats::heston_nml_cekf::nmle_cekf_heston` | `HestonNMLECEKF` | | Neural surrogate of the implied-vol surface | `ai::volatility::heston::HestonNn` (feature `ai`) | `HestonNn` (source builds with `ai`) | ## Simulate paths `sample()` returns `[s, v]`, two arrays of `n` points; `sample_par(m)` returns `m` of them. The default variance step is Euler with truncation at zero; `.qe()` switches the same model to Andersen's quadratic-exponential scheme. ```rust title="tests/doctest_tutorials_heston_simulate.rs" // docs: tutorials/heston //! Simulates seeded Heston paths under the default Euler scheme and under Andersen's QE scheme. use stochastic_rs::simd_rng::Deterministic; use stochastic_rs::stochastic::volatility::HestonPow; use stochastic_rs::stochastic::volatility::heston::Heston; use stochastic_rs::traits::ProcessExt; #[test] fn heston_paths_under_euler_and_qe() { let euler = Heston::::new( Some(100.0), // s0 Some(0.04), // v0, a variance 2.0, // kappa 0.04, // theta 0.3, // sigma, the volatility of variance -0.7, // rho 0.03, // mu 252, // n grid points, t = 0 included Some(1.0), // t HestonPow::Sqrt, Some(false), // truncate the variance at zero instead of reflecting it Deterministic::new(42), ); let [s, v] = euler.sample(); assert_eq!((s.len(), v.len()), (252, 252)); assert!(v.iter().all(|&x| x >= 0.0)); // Same parameters, Andersen (2008) quadratic-exponential variance step. let paths = euler.qe().sample_par(2_000); let mean_st = paths.iter().map(|[s, _]| s[251]).sum::() / 2_000.0; assert!((mean_st - 100.0 * 0.03_f64.exp()).abs() < 2.0); } ``` ```python import stochastic_rs as srs p = srs.PyHeston(kappa=2.0, theta=0.04, sigma=0.3, rho=-0.7, mu=0.03, n=252, s0=100.0, v0=0.04, t=1.0, seed=42) s, v = p.sample() # two arrays of shape (252,) S, V = p.sample_par(10_000) # two arrays of shape (10000, 252) print(round(S[:, -1].mean(), 2)) # 103.11 # Log-price variant: the price is exp of an Euler step on ln S, so it stays positive. log_p = srs.PyHestonLog(mu=0.03, kappa=2.0, theta=0.04, xi=0.3, rho=-0.7, n=252, s0=100.0, v0=0.04, t=1.0, seed=42) s_log, v_log = log_p.sample() ``` ## Price options `HestonPricer` evaluates the semi-closed form of Heston (1993), $$ C = S e^{-q\tau} P_1 - K e^{-r\tau} P_2,\qquad P_j = \frac12 + \frac1\pi \int_0^\infty \operatorname{Re}\!\left[\frac{e^{-i\phi \ln K} f_j(\phi)}{i\phi}\right] d\phi, $$ with the characteristic functions $f_j$ in the "little Heston trap" form of Albrecher et al. (2007), which keeps the complex logarithm on its principal branch at long maturities. The struct holds the model only, so one instance prices a whole strike and maturity grid; `greeks(s, k, r, q, tau, option_type)` returns all nine first- and second-order Greeks by central differences, with vega, vanna, volga and veta built on the analytic derivative in $v_0$ and reported per unit of $\sqrt{v_0}$. `HestonFourier` prices the same model through the crate's Fourier pricers (Gil-Pelaez by default, also Carr-Madan, Lewis and COS), and `HestonAdiPricer` solves the Heston PDE on the sinh-stretched meshes of in 't Hout and Foulon, with an optional down-and-out barrier. ```rust title="tests/doctest_tutorials_heston_price.rs" // docs: tutorials/heston //! Prices a Heston European call in semi-closed form and by Gil-Pelaez Fourier inversion. use stochastic_rs::quant::OptionType; use stochastic_rs::quant::pricing::fourier::HestonFourier; use stochastic_rs::quant::pricing::heston::HestonPricer; use stochastic_rs::traits::ModelPricer; #[test] fn heston_call_two_ways() { let (s, k, r, q, tau) = (100.0, 100.0, 0.03, 0.01, 1.0); // Model state only: v0, rho, kappa, theta, sigma, lambda (market price of volatility risk). let model = HestonPricer::new(0.04, -0.7, 2.0, 0.04, 0.3, None); let (call, put) = model.call_put(s, k, r, q, tau); let parity = s * (-q * tau).exp() - k * (-r * tau).exp(); assert!((call - put - parity).abs() < 1e-8); // The Fourier model carries r and q for its characteristic function: pass the same values. let fourier = HestonFourier { v0: 0.04, kappa: 2.0, theta: 0.04, sigma: 0.3, rho: -0.7, r, q, }; assert!((fourier.price_call(s, k, r, q, tau) - call).abs() < 1e-4); let greeks = model.greeks(s, k, r, q, tau, OptionType::Call); assert!(greeks.delta > 0.0 && greeks.delta < 1.0 && greeks.vega > 0.0); } ``` ```python import stochastic_rs as srs # The Python pricers take the query (s, k, r, tau, q) in the constructor. pricer = srs.HestonPricer(s=100.0, v0=0.04, k=100.0, r=0.03, kappa=2.0, theta=0.04, sigma=0.3, rho=-0.7, tau=1.0, q=0.01) call, put = pricer.call_put() print(f"{call:.4f} {put:.4f}") # 8.5996 6.6391 fourier = srs.HestonFourier(v0=0.04, kappa=2.0, theta=0.04, sigma=0.3, rho=-0.7, r=0.03, q=0.01) print(f"{fourier.price_call(100.0, 100.0, 0.03, 0.01, 1.0):.4f}") # 8.5996 adi = srs.HestonAdiPricer(s=100.0, v0=0.04, k=100.0, r=0.03, kappa=2.0, theta=0.04, sigma=0.3, rho=-0.7, tau=1.0, q=0.01, barrier=90.0) print(f"{adi.price():.4f}") # 6.9153, down-and-out at 90 ``` ## Calibrate to market data `HestonCalibrator` fits $(v_0, \kappa, \theta, \sigma, \rho)$ to option prices by Levenberg-Marquardt on the price residuals, jointly over every maturity slice, with the analytic Jacobian of Cui et al. (2017) by default (`HestonJacobianMethod::NumericFiniteDiff` switches to finite differences). The optimiser moves in logistic coordinates, which keeps the parameters inside fixed boxes: $v_0 \in (0.005, 0.25)$, $\kappa \in (0.1, 20)$, $\theta \in (0.001, 0.4)$, $\sigma \in (0.01, 3)$ and $\lvert\rho\rvert \lt 0.9999$. Without an initial guess it starts from an NMLE fit when price and variance series are attached (`HestonCalibrator::new`'s `mle_s`, `mle_v`, `mle_r`), and from fixed defaults otherwise. A Tikhonov pull toward an anchor is available through `with_regularization` (Python `regularization=(anchor, weights)`). ```rust title="tests/doctest_tutorials_heston_calibrate.rs" // docs: tutorials/heston //! Calibrates Heston to call prices generated by known parameters and recovers them. use stochastic_rs::quant::OptionType; use stochastic_rs::quant::calibration::MarketSlice; use stochastic_rs::quant::calibration::heston::HestonCalibrator; use stochastic_rs::quant::pricing::heston::HestonPricer; use stochastic_rs::traits::CalibrationResult; use stochastic_rs::traits::Calibrator; use stochastic_rs::traits::ModelPricer; #[test] fn heston_calibrator_recovers_generating_parameters() { let (s, r) = (100.0, 0.02); let truth = HestonPricer::new(0.05, -0.6, 1.8, 0.06, 0.5, None); let strikes = (0..9).map(|i| 80.0 + 5.0 * i as f64).collect::>(); let slices = [0.25, 0.5, 1.0].map(|tau| MarketSlice { prices: strikes .iter() .map(|&k| truth.price_call(s, k, r, 0.0, tau)) .collect(), is_call: vec![true; strikes.len()], strikes: strikes.clone(), tau, }); // No initial guess: the calibrator starts from its fallback parameters. let calibrator = HestonCalibrator::from_slices(None, &slices, s, r, None, OptionType::Call, false); let result = calibrator.calibrate(None).unwrap(); let p = result.params(); assert!(result.converged() && result.rmse() < 1e-3); assert!((p.v0 - 0.05).abs() < 1e-3 && (p.rho + 0.6).abs() < 1e-2); assert!((p.kappa - 1.8).abs() < 0.05 && (p.sigma - 0.5).abs() < 1e-2); } ``` ```python import stochastic_rs as srs s, r = 100.0, 0.02 strikes = [80.0 + 5.0 * i for i in range(9)] def quote(k, tau): # stands in for market prices: Heston at known parameters return srs.HestonPricer(s=s, v0=0.05, k=k, r=r, kappa=1.8, theta=0.06, sigma=0.5, rho=-0.6, tau=tau).price() # One MarketSlice(strikes, prices, is_call, tau) per expiry. slices = [srs.MarketSlice(strikes, [quote(k, tau) for k in strikes], [True] * 9, tau) for tau in (0.25, 0.5, 1.0)] v0, kappa, theta, sigma, rho, converged, rmse = srs.HestonCalibrator(slices, s=s, r=r).calibrate() print(f"v0={v0:.4f} kappa={kappa:.3f} theta={theta:.4f} sigma={sigma:.3f} rho={rho:.3f}") # v0=0.0500 kappa=1.800 theta=0.0600 sigma=0.500 rho=-0.600 ``` Every quote is priced as the calibrator's `option_type` (`"call"` by default); the slices' `is_call` flags are not read. The result converts to a `HestonFourier` pricer through `to_model(r, q)`. ## Estimate from data `HestonMLE.nmle(s, v, r)` and `HestonMLE.pmle(s, v, r)` estimate the parameters from observed price and variance series with the closed-form estimators of Wang et al. (2018). They take the sampling step as $1/(\text{len} - 1)$, a series spanning one unit of time; in Rust, `nmle_heston_with_delta` and `pmle_heston_with_delta` take the step explicitly. `HestonNMLECEKF(s, r, delta, max_iters)` works from prices alone, alternating a consistent extended Kalman filter for the variance with the NMLE update. ```python import stochastic_rs as srs p = srs.PyHeston(kappa=2.0, theta=0.04, sigma=0.3, rho=-0.7, mu=0.03, n=10_001, s0=100.0, v0=0.04, t=1.0, seed=7) s, v = p.sample() # one year in 10,000 steps est = srs.HestonMLE.nmle(s, v, 0.03) print(round(est.sigma, 3), round(est.rho, 3)) # 0.304 -0.708 ``` A year of data pins down $\sigma$ and $\rho$; the mean reversion is identified by the length of the sample rather than its frequency, and the same run returns $\kappa \approx 5.5$ against a true 2. ## GPU sampling `Heston`, under both variance schemes, and `HestonLog` run on the [Euler engine](/docs/concepts/gpu-support): a whole path per kernel thread, `sample_par(m)` as one launch. In Rust, re-type the process with `.on::()` (`f32` only) or `.on::()` (`f32` or `f64`), calling `.qe()` first for the QE scheme; in Python, pass `device="metal"` together with `dtype="f32"`, or `device="cuda"`. The device draws its own random stream, so its paths agree with the CPU's in distribution, not bit for bit. The published wheels are CPU-only; the device backends need a source build with the `metal` or `cuda` feature. ## Notes * **Euler scheme.** The price takes an Euler step on $S$ itself, $S_{i} = S_{i-1}(1 + \mu\Delta t + \sqrt{v^+_{i-1}}\,\Delta W^S)$, and the variance steps from $v^+ = \max(v, 0)$ and is truncated at zero (or reflected with `use_sym`). `HestonLog` steps $\ln S$ instead, so its price stays positive; its drift comes from the pair `r`, `r_f` (as $r - r_f$), else `b` (cost of carry), else `mu`, one of which is required, and its vol-of-vol argument is named `xi`. * **QE scheme.** `.qe()` is Andersen's (2008) quadratic-exponential step with $\psi_c = 1.5$ and the central discretisation $\gamma_1 = \gamma_2 = \tfrac12$ of the integrated variance, on the log-price. It needs `HestonPow::Sqrt`, $\kappa > 0$ and $S_0 > 0$, and it is Rust-only: `PyHeston` uses the Euler scheme. * **Argument order.** The Rust constructors differ: `HestonPricer::new(v0, rho, kappa, theta, sigma, lambda)`, `HestonAdiPricer::new(v0, kappa, theta, sigma, rho)` and `Heston::new(s0, v0, kappa, theta, sigma, rho, mu, …)`. `lambda` is the market price of volatility risk; `None` is 0. * **`HestonFourier` holds `r` and `q`** for the drift of its characteristic function, while `price_call(s, k, r, q, tau)` discounts with its arguments: pass the same rates to both. * **Deep out of the money.** `HestonPricer`'s price is not floored at zero: at short maturities far out of the money its quadrature can land a few $10^{-5}$ below zero. The Gil-Pelaez pricer behind `HestonFourier::price_call` floors at zero. * **Precision.** The processes are generic over `f32` and `f64` (`Heston::`); the pricers, calibrators and estimators work in `f64`. ## FAQ ### How do I simulate the Heston model in Python? Build `srs.PyHeston(kappa, theta, sigma, rho, mu, n, s0=..., v0=..., t=..., seed=...)` and call `sample()` for one `(s, v)` pair of arrays or `sample_par(m)` for two `(m, n)` arrays. `PyHestonLog` is the log-price variant. ### Does stochastic-rs enforce the Feller condition? No. Neither the process, the pricers nor `HestonCalibrator` imposes $2\kappa\theta \ge \sigma^2$, since market fits often violate it. `HestonParams::satisfies_feller_condition()` reports it and `projected_with_feller_condition()` returns a parameter set that satisfies it. ### How do I calibrate the Heston model to option prices? Collect one `MarketSlice(strikes, prices, is_call, tau)` per expiry and pass the list to `HestonCalibrator(slices, s, r)` in Python, or to `HestonCalibrator::from_slices` in Rust. `calibrate()` returns `(v0, kappa, theta, sigma, rho, converged, rmse)` in Python and a `HestonCalibrationResult` in Rust. ### Which discretisation should I use for Monte Carlo? Euler is the default and the cheapest. The QE scheme (`.qe()` in Rust) has a markedly lower bias when the vol-of-vol is large or the Feller condition is violated, at about the same cost per step. ## Where to go next * [Heston stochastic-local volatility](/docs/processes#heston-stochastic-local-volatility) — the same model under a calibrated leverage function, and its [leverage calibration](/docs/quant#heston-slv-leverage-calibration) * [Heston calibration to a market vol surface](/docs/quant#heston-calibration-to-a-market-vol-surface) and [the ADI pricer](/docs/quant#heston-pde-by-adi-finite-differences) on the quant page * [Heston on the processes page](/docs/processes#heston-stochastic-volatility), next to its double, multifactor and rough relatives * [Calibrating with a neural surrogate](/docs/ai#calibrating-with-a-surrogate) * [Fit an SVI volatility surface](/docs/tutorials/svi-volatility-surface) to the smile a calibrated model produces ## References - Heston, S. L. (1993). A closed-form solution for options with stochastic volatility with applications to bond and currency options. doi:10.1093/rfs/6.2.327 - Andersen, L. (2008). Simple and efficient simulation of the Heston stochastic volatility model. doi:10.21314/JCF.2008.189 - Albrecher, H., Mayer, P., Schoutens, W. and Tistaert, J. (2007). The little Heston trap. - Cui, Y., del Baño Rollin, S. and Germano, G. (2017). Full and fast calibration of the Heston stochastic volatility model. doi:10.1016/j.ejor.2017.05.018 - in 't Hout, K. J. and Foulon, S. (2010). ADI finite difference schemes for option pricing in the Heston model with correlation. arXiv:0811.3427 - Wang, X., He, X., Bao, Y. and Zhao, Y. (2018). Parameter estimates of Heston stochastic volatility model with MLE and consistent EKF algorithm. doi:10.1007/s11432-017-9215-8