Heston Stochastic Volatility: Theory and C++ Implementation

Setup

Why Stochastic Volatility?

The Black-Scholes model assumes constant volatility. The market implied volatility surface $\sigma^*(T,K)$ is demonstrably not flat — it exhibits the volatility smile (a U-shape) and a downward skew for equity index options. These shapes cannot be reproduced by any constant-volatility model.

Two model families address this:

1. Local volatility (Dupire): $\sigma = \sigma(t, S_t)$. Fits the surface exactly but produces a term structure of skew that flattens too fast — the model forward smile is too flat relative to observed dynamics.

2. Stochastic volatility: the volatility itself is a random process with its own SDE. The correlation between the asset and vol processes generates skew that persists across maturities.

The Heston (1993) model is the industry benchmark stochastic vol model. Under the risk-neutral measure:

$dS_t = r S_t\,dt + \sqrt{V_t}\, S_t\,dW_t^S$ $dV_t = \kappa(\theta - V_t)\,dt + \xi\sqrt{V_t}\,dW_t^V$ $d\langle W^S, W^V \rangle_t = \rho\,dt$

Parameters:

• $V_0$: initial variance (not volatility — note $\sqrt{V_0}$ is the initial vol)
• $\kappa > 0$: mean-reversion speed
• $\theta > 0$: long-run variance
• $\xi > 0$: volatility of variance (vol-of-vol)
• $\rho \in (-1, 1)$: correlation between spot and vol shocks; negative for equity (larger moves down are associated with vol spikes)

Conventions:

• $r$: continuously compounded risk-free rate, annualised.
• $T$: time to expiry in years.
• Volatility parameters in decimal (not percent).

Theory

1. The Feller Condition

The CIR variance process $dV_t = \kappa(\theta - V_t)\,dt + \xi\sqrt{V_t}\,dW_t^V$ can hit zero. Whether it stays at zero or bounces off depends on the Feller condition:

$2\kappa\theta > \xi^2 \quad \Longrightarrow \quad V_t > 0 \text{ almost surely}$

When the Feller condition is violated ($2\kappa\theta \leq \xi^2$), the process can reach zero and the square root is undefined during simulation. This is not merely numerical — it reflects model instability. On a derivatives desk, parameters are routinely calibrated to surfaces that violate Feller; the simulation scheme must handle this gracefully.

2. Characteristic Function Pricing (Lewis Formula)

The Heston model has a known characteristic function for $\ln(S_T/S_0)$:

$\phi(u; T) = \mathbb{E}^{\mathbb{Q}}\!\left[e^{iu\ln(S_T/S_0)}\right] = \exp\!\left(iu(r - q)T + A(u, T) + B(u, T) V_0\right)$

where (using the Albrecher et al. 2007 numerically stable formulation):

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

The Lewis (2001) formula prices a European call using a single real-valued integral:

$C(S_0, K, T) = S_0 - \frac{\sqrt{S_0 K}\, e^{-rT/2}}{\pi} \int_0^\infty \text{Re}\!\left[\frac{e^{iu\ln(S_0/K)} \phi(u - i/2; T)}{u^2 + 1/4}\right] du$

The integrand is smooth and rapidly decaying; standard Gauss-Legendre quadrature on $[0, U_{\max}]$ with $U_{\max} \approx 200$ and $n \approx 200$ points gives machine precision.

The stable branch selection (Albrecher 2007 vs. the original Heston 1993 formulation) matters critically: the original formulation has a branch-cut discontinuity in the complex logarithm that causes pricing errors for long maturities or extreme parameters. Always use the Albrecher formulation.

3. Monte Carlo Simulation

For path-dependent products (barriers, Asians) that the Fourier method cannot price analytically, Monte Carlo is required. The variance process must be discretised carefully.

Full truncation Euler (Lord et al. 2010): $V_{t+\Delta t} = V_t + \kappa(\theta - V_t^+)\Delta t + \xi\sqrt{V_t^+\,\Delta t}\,Z_1$ where $V^+ = \max(V, 0)$ and $Z_1 \sim \mathcal{N}(0,1)$.

The asset step uses the current (potentially zero-floored) variance: $S_{t+\Delta t} = S_t \exp\!\left[\left(r - \tfrac{1}{2}V_t^+\right)\Delta t + \sqrt{V_t^+\,\Delta t}\,Z_2\right]$

where $Z_2 = \rho Z_1 + \sqrt{1-\rho^2}\, Z_3$ with $Z_3 \sim \mathcal{N}(0,1)$ independent of $Z_1$.

The log-Euler scheme for $S$ avoids negative asset prices without any flooring. The correlation between $Z_1$ and $Z_2$ is exactly $\rho$.

Variance reduction — control variate: use the Black-Scholes price with current spot vol $\sqrt{V_0}$ as a control. This reduces standard error by 80-95% near at-the-money.

Implementation

HestonModel.h

#pragma once
#include <complex>

// Heston model parameters
struct HestonParams {
    double V0;       // Initial variance (vol = sqrt(V0))
    double kappa;    // Mean reversion speed
    double theta;    // Long-run variance
    double xi;       // Vol of variance (vol-of-vol)
    double rho;      // Spot-vol correlation

    // Returns true if the Feller condition 2*kappa*theta > xi^2 is satisfied
    bool feller_satisfied() const {
        return 2.0 * kappa * theta > xi * xi;
    }
};

// Computes European call/put price via Lewis (2001) formula
// using the Albrecher et al. (2007) numerically stable characteristic function.
double heston_call_price(double S0, double K, double r, double T,
                         const HestonParams& p);

double heston_put_price(double S0, double K, double r, double T,
                        const HestonParams& p);