Heston Model: Calibration and Simulation

Setup

Motivation

Black-Scholes assumes volatility is constant. The empirical reality is that volatility is stochastic: it clusters, mean-reverts, and correlates with the spot. The simplest tractable model that captures these effects is the Heston (1993) model, in which the instantaneous variance follows a mean-reverting square-root (CIR) process correlated with the spot.

The Heston model is analytically semi-tractable: European option prices are available via Fourier inversion of an explicit characteristic function. This makes calibration computationally feasible.

Assumptions and Parameters

Let $(\Omega, \mathcal{F}, \mathbb{Q})$ be the risk-neutral filtered probability space. The Heston dynamics are:

$dS_t = r S_t \, dt + \sqrt{v_t}\, S_t \, dW_t^{(1)}, \qquad S_0 > 0,$

$dv_t = \kappa(\theta - v_t)\, dt + \xi \sqrt{v_t}\, dW_t^{(2)}, \qquad v_0 > 0,$

$d\langle W^{(1)}, W^{(2)} \rangle_t = \rho \, dt.$

Parameters — all must be stated before use:

• $r$: risk-free rate (continuously compounded, constant).
• $\kappa > 0$: mean-reversion speed of variance.
• $\theta > 0$: long-run variance (mean level to which $v_t$ reverts).
• $\xi > 0$: volatility of variance (vol-of-vol).
• $\rho \in (-1, 1)$: correlation between spot and variance Brownian motions.
• $v_0 > 0$: initial instantaneous variance.

The instantaneous volatility is $\sqrt{v_t}$ and the implied vol at ATM at time 0 is approximately $\sqrt{v_0}$ for short maturities.

The Feller Condition

The variance process $v_t$ is a Cox-Ingersoll-Ross (CIR) process. Its boundary behaviour at zero depends on the ratio $2\kappa\theta / \xi^2$:

$\text{Feller condition}: \quad 2\kappa\theta > \xi^2.$

If Feller holds, the boundary $v_t = 0$ is inaccessible: the variance process remains strictly positive almost surely. This ensures $\sqrt{v_t}$ is well-defined at all times.

If Feller fails ($2\kappa\theta \leq \xi^2$), the variance process can reach zero with positive probability and reflect. For $2\kappa\theta < \xi^2$, the boundary is reflecting (the process can touch zero and immediately leave); for $2\kappa\theta = \xi^2$, the boundary is regular.

In practice, calibrated Heston parameters frequently violate the Feller condition — market-implied vol-of-vol $\xi$ is often large relative to mean reversion. This is not a catastrophic failure: the SDE remains well-posed (with reflecting boundary), but numerical simulation requires care. The full-truncation Euler scheme (described below) handles this robustly.

Characteristic Function and the Lewis Formula

Characteristic Function

The key to semi-analytic pricing is that the log-price $x_T = \ln(S_T/S_0)$ has a known conditional characteristic function under $\mathbb{Q}$:

$\varphi_T(u) = \mathbb{E}^{\mathbb{Q}}\!\left[e^{iu x_T}\right] = \exp\!\left(iurT + A(u,T) + B(u,T)\, v_0\right),$

where the functions $A$ and $B$ are given by the Albrecher et al. (2007) stable formulation:

$d(u) = \sqrt{(\rho\xi iu - \kappa)^2 + \xi^2(iu + u^2)},$

$g(u) = \frac{\kappa - \rho\xi iu - d(u)}{\kappa - \rho\xi iu + d(u)},$

$B(u, T) = \frac{\kappa - \rho\xi iu - d(u)}{\xi^2} \cdot \frac{1 - e^{-d(u)T}}{1 - g(u) e^{-d(u)T}},$

$A(u, T) = \frac{\kappa\theta}{\xi^2}\left[(\kappa - \rho\xi iu - d(u))\, T - 2\ln\!\left(\frac{1 - g(u) e^{-d(u)T}}{1 - g(u)}\right)\right].$

Why the Albrecher Formulation?

The original Heston (1993) characteristic function uses a different branch of the square root for $d(u)$. For long maturities or large $|u|$, the original formulation can produce branch-cut discontinuities that cause the characteristic function to be evaluated on the wrong Riemann sheet, leading to mispriced options. The Albrecher formulation avoids this by choosing a consistent branch throughout. Always use the Albrecher version in production.

Lewis Formula for European Options

Lewis (2001) derived a Fourier pricing formula that avoids the dampening-parameter sensitivity of Carr-Madan. For a European call with log-strike $k = \ln(K/F)$ where $F = S_0 e^{rT}$:

$C = S_0 - \frac{\sqrt{FK} \cdot e^{-rT}}{\pi} \int_0^\infty \mathrm{Re}\!\left[\varphi_T\!\left(u - \frac{i}{2}\right) \frac{e^{-iu k}}{u^2 + \frac{1}{4}}\right] du.$

The integrand is real-valued (after taking the real part) and rapidly decaying. Numerical integration via adaptive Gauss-Legendre or Gauss-Kronrod is standard. This integral does not require a dampening parameter $\alpha$ and is numerically stable for $u \in (0, \infty)$.

Derivation sketch. Express the call price as $C = e^{-rT}\mathbb{E}^{\mathbb{Q}}[(S_T - K)^+]$. Write the payoff in terms of the log-price, apply Parseval's identity in Fourier space, and use the characteristic function of the log-price. The $1/(u^2 + 1/4)$ denominator arises from the Fourier transform of the call payoff along a contour $\mathrm{Im}(u) = -1/2$.