Heston Stochastic Volatility Model

Finance

Black-Scholes assumes volatility is CONSTANT — but as the volatility-surfaces page makes clear, real markets price options as if volatility itself were random. The Heston model (Steven Heston, 1993) is the standard parametric specification for stochastic volatility: instantaneous variance follows its own mean-reverting square-root process, correlated with the underlying. The model has SEMI-ANALYTIC vanilla option pricing via the characteristic function — fast calibration to a vol surface, and the natural baseline for exotic-option pricing under stochastic vol. It is the most cited derivative-pricing model after Black-Scholes itself.

The model

Risk-neutral dynamics for the underlying $S_{t}$ and its instantaneous variance $v_{t}$ :

d S_{t} = r S_{t} d t + v_{t} S_{t} d W_{t}^{S},

d v_{t} = κ (θ - v_{t}) d t + ξ v_{t} d W_{t}^{v}, Corr (d W^{S}, d W^{v}) = ρ .

Five parameters with clean interpretations:

$κ > 0$ — MEAN REVERSION speed. Higher κ → variance moves back to its long-run mean faster.
$θ > 0$ — LONG-RUN MEAN of variance. The level $v$ reverts toward; equals long-run unconditional variance.
$ξ > 0$ — VOL-OF-VOL. How volatile the variance process itself is. Drives the CONVEXITY (smile curvature) of the vol surface.
$ρ \in [- 1, 1]$ — CORRELATION between the asset and variance shocks. Negative $ρ$ in equities (vol spikes when stocks fall) drives the SKEW.
$v_{0} > 0$ — initial instantaneous variance. Pins down today's ATM vol.

Three features that make Heston the standard. (1) The variance process is the CIR / Feller process — known to stay positive (under the Feller condition $2 κ θ \geq ξ^{2}$ ), tractable in closed form, mean-reverting in a realistic way. (2) The correlation $ρ$ produces ASYMMETRIC smiles directly — the negative skew of equity options is built in from $ρ < 0$ . (3) The characteristic function of $ln S_{T}$ is KNOWN IN CLOSED FORM, which means vanilla options can be priced by single-integral Fourier inversion — fast enough for daily calibration to market surfaces.

The characteristic function

Heston (1993) solved the PDE for the conditional joint distribution of $(S_{T}, v_{T})$ and obtained the characteristic function of $x = ln S_{T}$ :

ϕ (u) = E [e^{i ux}] = exp (A (u, T) + B (u, T) v_{0} + i u (ln S_{0} + r T))

where (in the Albrecher et al. 2007 "trap" formulation, which avoids the branch-cut issue in the original parameterization):

d = (κ - i ρ ξ u)^{2} + ξ^{2} (u^{2} + i u),

g = \frac{κ - i ρ ξ u - d}{κ - i ρ ξ u + d},

A (u, T) = \frac{κ θ}{ξ ^{2}} [(κ - i ρ ξ u - d) T - 2 ln \frac{1 - g e ^{- d T}}{1 - g}],

B (u, T) = \frac{κ - i ρ ξ u - d}{ξ ^{2}} \frac{1 - e ^{- d T}}{1 - g e ^{- d T}} .

Closed-form (complex-valued, but computable). The original Heston paper used a different equivalent expression involving the complex square root that had numerical issues at long maturities; the "trap" form is robust for all practical parameter ranges.

Vanilla call pricing via inversion

Heston gives the European call price as:

C = S_{0} P_{1} - K e^{- r T} P_{2},

where $P_{1}, P_{2}$ are risk-neutral exercise probabilities under two different measures, each expressed as an inversion integral:

P_{j} = \frac{1}{2} + \frac{1}{π} \int_{0}^{\infty} Re [\frac{e ^{- i u l n K} ϕ _{j} ( u )}{i u}] d u, j = 1, 2,

with appropriate $ϕ_{j}$ built from the same characteristic function. The integrals are typically evaluated by Gauss-Legendre or adaptive quadrature out to $u \sim 200$ (the integrand decays exponentially in $u$ for non-trivial maturities, so high-frequency truncation is benign). A SINGLE-INTEGRAL version exists (Carr-Madan FFT formulation) and is even faster — practically used for production pricing of strike grids.

Code: Heston pricing and the smile

# Heston (1993) vanilla call pricing via the characteristic-function
# inversion formula (Heston's original two-integral approach).
#
#   dS = r S dt + sqrt(v) S dW^S
#   dv = kappa (theta - v) dt + xi sqrt(v) dW^v
#   corr(dW^S, dW^v) = rho

import numpy as np
from scipy.integrate import quad

def heston_char_fn(u, S0, T, r, kappa, theta, xi, rho, v0):
    """Heston characteristic function of x = ln(S_T), 'trap' formulation
    (Albrecher et al. 2007 — avoids branch-cut issues in the log)."""
    d = np.sqrt((kappa - 1j*rho*xi*u)**2 + xi**2 * (u**2 + 1j*u))
    g = (kappa - 1j*rho*xi*u - d) / (kappa - 1j*rho*xi*u + d)
    A = 1j*u*(np.log(S0) + r*T) + kappa*theta/xi**2 * (
        (kappa - 1j*rho*xi*u - d)*T
        - 2*np.log((1 - g*np.exp(-d*T))/(1 - g))
    )
    B = (kappa - 1j*rho*xi*u - d)/xi**2 * (1 - np.exp(-d*T))/(1 - g*np.exp(-d*T))
    return np.exp(A + B*v0)

def heston_call(S0, K, T, r, kappa, theta, xi, rho, v0):
    """Vanilla call via the Heston two-integral formula:
        C = S0 * P1 - K * exp(-rT) * P2
    where P_j are risk-neutral exercise probabilities under two measures."""
    def integrand_P1(u):
        phi_im1 = heston_char_fn(u - 1j, S0, T, r, kappa, theta, xi, rho, v0)
        phi_0   = heston_char_fn(-1j,    S0, T, r, kappa, theta, xi, rho, v0)
        return np.real(np.exp(-1j*u*np.log(K)) * phi_im1 / (1j*u*phi_0))
    def integrand_P2(u):
        phi = heston_char_fn(u, S0, T, r, kappa, theta, xi, rho, v0)
        return np.real(np.exp(-1j*u*np.log(K)) * phi / (1j*u))
    P1 = 0.5 + (1/np.pi)*quad(integrand_P1, 1e-8, 200, limit=200)[0]
    P2 = 0.5 + (1/np.pi)*quad(integrand_P2, 1e-8, 200, limit=200)[0]
    return S0*P1 - K*np.exp(-r*T)*P2

# Typical equity-index parameters
params = dict(kappa=2.0, theta=0.04, xi=0.3, rho=-0.7, v0=0.04)
S0, T, r = 100.0, 0.5, 0.03
print(f"Heston parameters: {params}")
print(f"ATM 6-month call (S=K=100, r=3%): {heston_call(S0, 100, T, r, **params):.4f}")
print(f"BS reference with flat sigma=20%: 6.3710")
print()
print(f"Smile across strikes (vs flat-vol BS):")
from scipy.stats import norm
def bs_call(S, K, T, r, sigma):
    d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    return S*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)

print(f"  {'Strike':>8s}  {'Heston':>10s}  {'BS (20%)':>10s}  {'diff':>10s}")
for K in [85, 90, 95, 100, 105, 110, 115]:
    h = heston_call(S0, K, T, r, **params)
    b = bs_call(S0, K, T, r, 0.20)
    print(f"  {K:>8d}  {h:>10.4f}  {b:>10.4f}  {h - b:>10.4f}")
print("(rho < 0 → low-strike puts cost MORE than flat-vol BS predicts:")
print(" the equity-index 'volatility skew'.)")

Output:

Heston parameters: {'kappa': 2.0, 'theta': 0.04, 'xi': 0.3, 'rho': -0.7, 'v0': 0.04}
ATM 6-month call (S=K=100, r=3%): 6.2647
BS reference with flat sigma=20%: 6.3710

Smile across strikes (vs flat-vol BS):
    Strike      Heston    BS (20%)        diff
        85     17.2482     16.9119      0.3363
        90     13.1036     12.7993      0.3043
        95      9.3988      9.2511      0.1477
       100      6.2647      6.3710     -0.1063
       105      3.8061      4.1783     -0.3722
       110      2.0631      2.6119     -0.5488
       115      0.9789      1.5593     -0.5804
(rho < 0 → low-strike puts cost MORE than flat-vol BS predicts:
 the equity-index 'volatility skew'.)

Three things to read off. (1) ATM Heston price (6.26) is slightly BELOW the flat-vol BS (6.37) at $σ = θ = 20%$ — a small effect from the convexity of the vol-of-vol term. (2) For low-strike puts (K = 85, 90), Heston prices ABOVE flat BS — fat left tail driven by $ρ < 0$ , i.e., variance spikes when the stock falls, amplifying downside risk. (3) For high-strike calls (K = 110, 115), Heston prices BELOW flat BS — thin right tail, same $ρ < 0$ mechanism. The skew pattern (low-strike pricier than flat BS, high-strike cheaper) is the equity-index volatility skew built into Heston by construction.

The Feller condition and the boundary at zero

The variance process $d v = κ (θ - v) d t + ξ v d W$ is a Cox-Ingersoll-Ross / Feller process. Its boundary behavior at $v = 0$ depends on the FELLER CONDITION:

2 κ θ \geq ξ^{2} .

When satisfied, $v_{t}$ almost surely never reaches zero (the boundary is INACCESSIBLE). When violated, the process can hit zero in finite time (still non-negative, just reflecting). Empirically calibrated equity-index parameters often violate Feller — typical $ξ$ is high enough that $ξ^{2} > 2 κ θ$ . Numerical schemes for simulating Heston in the violation regime need care (the standard Euler scheme can go negative on isolated paths, requiring reflection or truncation; QE / log-normal schemes handle the boundary cleanly).

Calibration in practice

Given the daily vol surface (matrix of implied vols across strikes and maturities), find the Heston parameters $(κ, θ, ξ, ρ, v_{0})$ that best reproduce it. The optimization:

κ, θ, ξ, ρ, v_{0} min (K, T) \sum w_{K T} (σ_{Heston}^{imp} (K, T) - σ_{mkt}^{imp} (K, T))^{2} .

Five-parameter nonlinear least-squares. Standard tools: Levenberg-Marquardt or differential evolution. Each evaluation of the loss requires Heston-pricing every $(K, T)$ on the surface; with Carr-Madan FFT this is fast. Typical calibration time: 1-10 seconds per surface. Done at every market update during trading hours.

KNOWN CALIBRATION FAILURES: Heston cannot fit the EXTREME WINGS of the surface (very far OTM options) — the model's tail behavior is exponential, while the market's is more like power-law. For equity options out to maturities of a year or so, Heston is excellent; for longer-dated or far-wing pricing, more complex models (SABR with stochastic skew, Bates with jumps, full local-stochastic-vol) are needed.

Heston extensions

Bates model: Heston + jumps in the asset price (Merton-style). Captures short-dated skew better; the abrupt smile of weekly options needs jumps to reproduce.
Double Heston: two variance factors (fast and slow). Captures the term structure of vol better — short-dated and long-dated options can have different fitted Heston parameters under single-factor; double Heston resolves the tension.
Heston-Hull-White: Heston for equity vol + Hull-White for stochastic rates. Standard for long-dated equity-linked products.
SABR (Hagan et al. 2002): different SV specification — backbone parameter and a much simpler asymptotic expansion. Industry standard for interest-rate volatility (cap/floor and swaption smile).
Local-stochastic vol (LSV): Heston-style dynamics times a local-vol leverage function that ensures exact calibration to the vanilla surface. Production standard for path-dependent exotics under stochastic vol.

Greeks and risk management

Heston Greeks beyond those of Black-Scholes:

Vega^2 / volga: $\partial^{2} C / \partial v_{0}^{2}$ . Quadratic exposure to current vol; captures convexity in the SHORT-DATED smile.
Vol-of-vol Greek: $\partial C / \partial ξ$ . Sensitivity to the curvature of the smile; key for trading vol-of-vol products (variance dispersion, VIX options).
Skew Greek: $\partial C / \partial ρ$ . Sensitivity to the slope of the smile; key for trading equity skew via OTM puts vs ATM.

Computed by finite differences on the calibrated model. A volatility-trading desk has positions whose Greeks are aggregated across the surface and risk-managed against limits — same operational structure as a vanilla options book, but with the additional vol-of-vol and skew exposures explicit.

Volatility surfaces and calibration — what Heston is calibrated TO; the empirical fact Heston explains.
Black-Scholes model — the constant-vol baseline that Heston extends.
Greeks and delta hedging — the vanilla option Greeks; Heston adds vega-of-vega and skew Greeks.
GARCH and realized volatility — discrete-time analogue: Heston is the continuous-time stochastic-vol model; GARCH is its discrete-time counterpart for time-series modelling.
American options — LSM under Heston dynamics is the standard for American exotics with stochastic vol.