Heston: simulate, price and calibrate

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

with $p = \tfrac12$ for the Heston model. HestonPow::ThreeHalves sets $p = \tfrac32$, a 3/2-type diffusion under the same linear drift.

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

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.

// 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::<f64, _>::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::<f64>() / 2_000.0;
  assert!((mean_st - 100.0 * 0.03_f64.exp()).abs() < 2.0);
}

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),

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.

// 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);
}

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)).

// 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::<Vec<_>>();
  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);
}

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.

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: a whole path per kernel thread, sample_par(m) as one launch. In Rust, re-type the process with .on::<Metal>() (f32 only) or .on::<Cuda>() (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::<f32, _>); 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 — the same model under a calibrated leverage function, and its leverage calibration
• Heston calibration to a market vol surface and the ADI pricer on the quant page
• Heston on the processes page, next to its double, multifactor and rough relatives
• Calibrating with a neural surrogate
• Fit an 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