The Heston Model: Stochastic Volatility for Options

The Heston model makes instantaneous variance a random, mean-reverting factor correlated with the underlying return. That one addition repairs the constant-volatility assumption of Black-Scholes enough to generate skew, term structure, and a volatility-of-volatility risk factor while retaining semi-analytic European option pricing.

Risk-neutral dynamics

Under the pricing measure, the standard Heston specification is:

dS_t = (r - q) S_t dt + sqrt(v_t) S_t dW_S,t
dv_t = kappa (theta - v_t) dt + xi sqrt(v_t) dW_v,t
dW_S dW_v = rho dt

v_t is instantaneous variance, not implied volatility squared. Parameter roles are:

Parameter	Function	Surface consequence
v0	current variance	short-dated ATM level
theta	long-run variance	long-end level
kappa	mean-reversion speed	term-structure transition
xi	vol of variance	smile curvature
rho	spot/variance correlation	skew direction and slope

For equity indices, rho < 0 is typical: falling prices coincide with rising variance, making downside puts rich. A fit with positive equity rho may converge numerically but should be treated as a data or convention failure.

Positivity and the Feller condition

Variance follows a CIR process. A sufficient condition for it to stay strictly positive is:

2 * kappa * theta >= xi^2

This is the Feller condition. Violating it does not automatically invalidate a pricing calibration—zero can be attainable but reflecting—but it increases simulation sensitivity and can produce a parameter set that is unstable under small quote changes. Do not use it as a hard market constraint without testing fit quality; use it as a diagnostic and apply a soft penalty if your simulator cannot robustly handle the boundary.

Why characteristic functions matter

Heston is affine: the conditional characteristic function of log(S_T) has exponential-affine form:

phi(u) = E_Q[exp(i*u*log(S_T)) | F_t]
       = exp(C(u,T) + D(u,T) v_t + i*u*log(S_t))

The Riccati functions C and D are closed form. Fourier inversion then prices a call:

C(K,T) = S0 exp(-qT) P1 - K exp(-rT) P2
Pj = 1/2 + (1/pi) integral_0^infinity Re[ exp(-i*u*log(K)) fj(u)/(i*u) ] du

This avoids Monte Carlo noise during calibration. Implementations must use a stable “little Heston trap” formulation or an equivalent branch-continuous characteristic function; naïvely taking complex square roots can create discontinuities in strikes or maturities.

Calibration is an inverse problem

Fit liquid implied-vol quotes, not option prices with equal weights. Price errors overweight high-premium options, whereas desk risk is normally quoted in volatility points. Use bid-ask-aware weights:

min_p sum_i ((IV_model(p; K_i, T_i) - IV_mid_i) / max(half_spread_i, floor))^2

from scipy.optimize import least_squares

def residuals(raw_params, quotes, heston_iv):
    # transform enforces positivity; rho is bounded by tanh
    kappa, theta, xi = np.exp(raw_params[:3])
    rho = np.tanh(raw_params[3])
    v0 = np.exp(raw_params[4])
    pars = (kappa, theta, xi, rho, v0)
    return [(heston_iv(q.K, q.T, pars) - q.iv) / max(q.half_spread, 0.002)
            for q in quotes]

fit = least_squares(residuals, x0, args=(quotes, heston_iv), loss="soft_l1")

Constrain economic ranges, use multiple starting points, and report parameter uncertainty. Many combinations of kappa, theta, and xi fit vanilla surfaces similarly, especially when only short maturities are liquid.

Simulation for exotics

Plain Euler can make variance negative. “Full truncation” is a practical baseline:

v_plus = np.maximum(v, 0.0)
z1 = rng.standard_normal(n)
z2 = rng.standard_normal(n)
zv = rho * z1 + np.sqrt(1 - rho**2) * z2
v = v + kappa * (theta - v_plus) * dt + xi * np.sqrt(v_plus * dt) * zv
S *= np.exp((r - q - 0.5 * v_plus) * dt + np.sqrt(v_plus * dt) * z1)

For barriers, short maturities, or calibration-quality prices, use Andersen’s QE scheme or exact-CIR sampling. Also add Brownian-bridge barrier correction; daily monitoring in a simulation is not continuous monitoring.

Risk and model governance

Heston captures a single stochastic variance factor, not jumps, dividends, borrow uncertainty, or the full moving surface. Its spot-vol correlation can hedge skew moves locally, but calibrated parameters are not stable trading signals. Track P&L explain by delta, vega, volga, vanna, and parameter moves; a model re-fit can create material “model P&L” even when market quotes barely move.

Compare it to implied-volatility-surface observations daily. A low RMS error can hide static-arbitrage violations between interpolated points, so test calendar monotonicity and convexity in strike after constructing the complete surface.

Key takeaways

• Heston makes variance stochastic and correlated with spot, generating realistic skew.
• rho controls skew, xi curvature, and v0 the front-end variance level.
• Fourier pricing is fast, but complex-branch stability is a production requirement.
• Calibrate to IV using bid-ask weights and diagnose parameter non-identifiability.
• Use positivity-aware simulation and separate model-fit quality from hedge performance.