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

// Monte Carlo pricing with full-truncation Euler scheme
// Returns {price, standard_error}
struct MCResult {
    double price;
    double std_error;
};

MCResult heston_mc_call(double S0, double K, double r, double T,
                        const HestonParams& p,
                        int n_paths, int n_steps,
                        unsigned int seed = 42);

HestonModel.cpp

#include "HestonModel.h"
#include <cmath>
#include <complex>
#include <random>
#include <vector>
#include <stdexcept>

using Complex = std::complex<double>;

// ── Characteristic function (Albrecher 2007 stable branch) ──────────────────

static Complex heston_char_fn(const Complex& u, double T,
                               double r, const HestonParams& p) {
    // d = sqrt((rho*xi*iu - kappa)^2 + xi^2*(iu + u^2))
    Complex iu   = Complex(0.0, 1.0) * u;
    Complex d    = std::sqrt(std::pow(p.rho * p.xi * iu - p.kappa, 2)
                           + p.xi * p.xi * (iu + u * u));
    // g = (kappa - rho*xi*iu - d) / (kappa - rho*xi*iu + d)
    Complex kmrd = p.kappa - p.rho * p.xi * iu;
    Complex g    = (kmrd - d) / (kmrd + d);
    Complex edT  = std::exp(-d * T);

    // B(u, T)
    Complex B = (kmrd - d) / (p.xi * p.xi) * (1.0 - edT) / (1.0 - g * edT);

    // A(u, T)  — Albrecher formulation avoids log branch cut
    Complex A = p.kappa * p.theta / (p.xi * p.xi)
                * ((kmrd - d) * T - 2.0 * std::log((1.0 - g * edT) / (1.0 - g)));

    return std::exp(iu * (r * T) + A + B * p.V0);
}

// ── Lewis (2001) formula ─────────────────────────────────────────────────────

double heston_call_price(double S0, double K, double r, double T,
                          const HestonParams& p) {
    if (T <= 0.0)
        throw std::invalid_argument("heston_call_price: T must be positive");

    // Gauss-Legendre quadrature on [0, U_max]
    const int    n     = 200;
    const double U_max = 200.0;

    // Gauss-Legendre nodes and weights on [-1, 1], mapped to [0, U_max]
    // Use midpoint rule for simplicity; production code uses GL nodes.
    double h = U_max / n;
    double integral = 0.0;
    double log_ratio = std::log(S0 / K);

    for (int j = 0; j < n; ++j) {
        double u  = (j + 0.5) * h;
        // phi(u - i/2; T)
        Complex phi = heston_char_fn(Complex(u, -0.5), T, r, p);
        double re   = (std::exp(Complex(0.0, u) * log_ratio) * phi).real();
        integral += re / (u * u + 0.25);
    }
    integral *= h;

    // Lewis formula: C = S0 - sqrt(S0*K) * exp(-r*T/2) / pi * integral
    return S0 - std::sqrt(S0 * K) * std::exp(-0.5 * r * T) / M_PI * integral;
}

double heston_put_price(double S0, double K, double r, double T,
                         const HestonParams& p) {
    // Put-call parity: P = C - S0 + K*exp(-rT)
    return heston_call_price(S0, K, r, T, p) - S0 + K * std::exp(-r * T);
}

// ── Full-Truncation Euler Monte Carlo ────────────────────────────────────────

MCResult heston_mc_call(double S0, double K, double r, double T,
                         const HestonParams& p,
                         int n_paths, int n_steps,
                         unsigned int seed) {
    const double dt   = T / n_steps;
    const double sqdt = std::sqrt(dt);

    std::mt19937_64 rng(seed);
    std::normal_distribution<double> normal(0.0, 1.0);

    double sum  = 0.0;
    double sum2 = 0.0;

    for (int path = 0; path < n_paths; ++path) {
        double S = S0;
        double V = p.V0;

        for (int step = 0; step < n_steps; ++step) {
            double Z1 = normal(rng);
            double Z3 = normal(rng);
            double Z2 = p.rho * Z1 + std::sqrt(1.0 - p.rho * p.rho) * Z3;

            double V_pos = std::max(V, 0.0);  // full truncation

            // Variance step (Euler on V)
            V = V + p.kappa * (p.theta - V_pos) * dt
                  + p.xi * std::sqrt(V_pos) * sqdt * Z1;

            // Asset step (log-Euler — exact for constant vol)
            S = S * std::exp((r - 0.5 * V_pos) * dt
                             + std::sqrt(V_pos) * sqdt * Z2);
        }

        double payoff = std::max(S - K, 0.0) * std::exp(-r * T);
        sum  += payoff;
        sum2 += payoff * payoff;
    }

    double mean = sum / n_paths;
    double var  = (sum2 / n_paths - mean * mean) / (n_paths - 1);

    return {mean, std::sqrt(var)};
}

Validation main.cpp

#include "HestonModel.h"
#include <iostream>
#include <iomanip>

int main() {
    // Parameters from Heston (1993) test case
    HestonParams p{
        .V0    = 0.04,   // initial variance: 20% vol
        .kappa = 2.0,
        .theta = 0.04,   // long-run variance: 20% vol
        .xi    = 0.3,
        .rho   = -0.7
    };

    double S0 = 100.0, K = 100.0, r = 0.05, T = 1.0;

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Feller condition: "
              << (p.feller_satisfied() ? "satisfied" : "VIOLATED") << "\n";

    double call = heston_call_price(S0, K, r, T, p);
    double put  = heston_put_price(S0, K, r, T, p);
    std::cout << "Lewis call: " << call << "\n";
    std::cout << "Lewis put:  " << put  << "\n";

    // Put-call parity check: C - P should equal S - K*exp(-rT)
    double parity = call - put - (S0 - K * std::exp(-r * T));
    std::cout << "PCP error:  " << parity << " (should be ~0)\n";

    // Monte Carlo comparison
    auto mc = heston_mc_call(S0, K, r, T, p, 100000, 252);
    std::cout << "MC call:    " << mc.price
              << " +/- " << mc.std_error << "\n";

    return 0;
}

Validation

Put-Call Parity

For any pricing model, $C - P = S_0 - K e^{-rT}$ (European options, no dividends). This is a model-free no-arbitrage relation — it must hold to machine precision:

$|C_{\text{Lewis}} - P_{\text{Lewis}} - (S_0 - Ke^{-rT})| < 10^{-10}$

Convergence: Lewis vs. Monte Carlo

The standard error decays as $1/\sqrt{N}$. The MC price should be within 3 standard errors of the Lewis price for any seed.

Black-Scholes Limit

When $\xi \to 0$, the Heston model degenerates to Black-Scholes with constant vol $\sigma = \sqrt{V_0}$. The Lewis formula should then reproduce the Black-Scholes formula to machine precision:

HestonParams bs{.V0=0.04, .kappa=1.0, .theta=0.04, .xi=1e-10, .rho=0.0};
// heston_call_price(..., bs) ≈ bs_call(S0, K, r, T, sigma=0.20)

Limitations

Calibration instability. The Heston model has 5 parameters ($V_0, \kappa, \theta, \xi, \rho$) plus initial conditions. Calibration to the implied vol surface is a non-linear optimisation problem that is frequently ill-posed: multiple parameter sets produce very similar surface fits. In practice, $\kappa$ and $\xi$ are highly correlated in the calibration — small perturbations in market vols can cause large parameter jumps. Regularisation (penalising distance from a prior) is standard.

Feller condition violations. Most calibrated parameter sets violate the Feller condition ($2\kappa\theta \leq \xi^2$), meaning the simulated variance process can hit zero. Full-truncation Euler handles this without crashing but introduces discretisation bias that does not vanish as $\Delta t \to 0$. The Andersen (2008) QE (Quadratic Exponential) scheme samples the exact conditional distribution of the CIR process and eliminates this bias at a modest computational cost.

No jumps. Equity index implied vol surfaces exhibit a pronounced short-dated skew driven by jump risk. The Heston model fits the term structure of skew but not the absolute level for short maturities. Bates (1996) augments Heston with Poisson jumps; Carr-Wu (2003) adds a stochastic time change. These extensions retain the characteristic function pricing framework but add parameters.

The Lewis integral needs careful tuning. The upper limit $U_{\max}$ and number of quadrature points depend on the parameters. For small $T$ (short-dated options), the characteristic function decays slowly and $U_{\max} = 500$ or higher may be needed. Adaptive quadrature (Gauss-Kronrod) is more robust than a fixed-grid midpoint rule.

Interview Angle

Junior (L1): What is the Feller condition and why does it matter for Monte Carlo simulation of the Heston model?

The Feller condition $2\kappa\theta > \xi^2$ ensures the CIR variance process stays strictly positive. When it is violated, the process can reach zero, making $\sqrt{V_t}$ undefined. The full-truncation fix (replace $V$ with $\max(V, 0)$ in the diffusion term) prevents a crash but introduces a systematic bias: paths where vol hits zero are "stuck" longer than the true process would allow.

Senior (L2): Why does the original Heston (1993) characteristic function formulation cause pricing errors for long maturities, and what is the Albrecher fix?

The original formulation involves $\ln(1 - g e^{-dT})$ where $g$ is a ratio of complex numbers. For long $T$ or certain parameter combinations, $g e^{-dT}$ rotates around the unit circle in the complex plane — the complex logarithm is multi-valued and a naive implementation picks the wrong branch, causing discontinuous price errors. Albrecher et al. (2007) rewrite the formula using a different algebraic identity for the $e^{dT}$ terms that avoids the branch-cut discontinuity entirely, at the cost of a slightly less intuitive expression.

Researcher (L3): Compare the full-truncation Euler scheme with the Andersen QE scheme for simulating the Heston variance process. When does each scheme produce the largest discretisation error?

Full-truncation Euler is first-order in the discretisation step $\Delta t$. Its error is largest when the Feller condition is strongly violated (high $\xi$, low $\kappa\theta$) or when the time step is large relative to the mean-reversion time $1/\kappa$. The QE scheme (Andersen 2008) samples from the exact conditional distribution of the CIR process given $V_t$, switching between a quadratic approximation (when $V_t$ is large, variance is approximately Gaussian) and an exponential distribution (when $V_t$ is near zero). The QE scheme has effectively zero discretisation bias for the variance process regardless of Feller — the remaining error is purely in the Euler step for the asset. At equal path count, QE reduces MC bias by 10–50x near the Feller boundary.