C++20 + Python● Live

Stochastic Vol Monte Carlo Engine

Heston (1993) stochastic volatility model implemented in C++20. Euler-Maruyama discretisation with full truncation and antithetic variates. Analytical pricing via the Lewis (2001) characteristic function. Five Plotly visualisations: smile comparison, variance paths, MC convergence, return distribution, and delta surface.

C++20CMakeCatch2HestonMonte CarloLewis formulaAntithetic variatesPythonPlotlyscipy

1. Market Context and Assumptions

Black-Scholes (1973) prices options under constant volatility and log-normal spot dynamics. Real equity option markets refute this on two counts: implied volatilities vary with strike (the volatility smile) and with maturity (the term structure of vol). The smile reflects fat tails and negative skew (crash risk); the term structure reflects mean-reverting realised volatility.

Heston (1993) introduced the first analytically tractable stochastic vol model consistent with both features. The variance process follows the Cox-Ingersoll-Ross (CIR) SDE — mean-reverting, non-negative (under the Feller condition), and correlated with the spot. The correlation parameter ρ directly controls skew: negative ρ makes puts more expensive than calls, matching equity market behaviour.

On a vol trading desk, Heston serves as a baseline model for exotics pricing, variance swap valuation, and the generation of risk scenarios. Its characteristic function in closed form makes calibration to the smile fast via Fourier inversion, and its five parameters have clear financial interpretations.

Model Assumptions

Spot and variance: both driven by correlated Brownian motions under the risk-neutral measure Q. No jumps.
Interest rate: constant, continuously compounded, annualised (e.g. 0.05 = 5%).
Variance process: CIR (Cox-Ingersoll-Ross) — mean-reverting, non-negative when the Feller condition holds (2κθ > ξ²).
Day count: ACT/365 throughout. T = 0.5 means exactly half a year.
Exercise: European only. The characteristic function approach prices vanilla calls and puts exactly; exotics require MC or PDE methods.
No dividends: spot follows a continuous GBM with drift r.

2. Heston SDE and Theory

2.1 The Model

Under the risk-neutral measure Q, the Heston model specifies two coupled SDEs:

$\mathrm{d}S_t = r\,S_t\,\mathrm{d}t + \sqrt{v_t}\,S_t\,\mathrm{d}W^S_t$

$\mathrm{d}v_t = \kappa(\theta - v_t)\,\mathrm{d}t + \xi\,\sqrt{v_t}\,\mathrm{d}W^v_t$

$\mathrm{d}\langle W^S, W^v \rangle_t = \rho\,\mathrm{d}t$

Parameter meanings — each has a financial interpretation:

$\kappa > 0$ : mean-reversion speed (yr $^{-1}$ ). Larger $\kappa$ pulls $v_t$ back to $\theta$ faster; half-life $= \ln 2/\kappa$ .
$\theta > 0$ : long-run variance. $\sqrt{\theta}$ is the long-run implied vol level.
$\xi > 0$ : vol of vol. Controls the curvature of the smile — higher $\xi$ generates more convexity.
$\rho \in (-1,1)$ : spot-vol correlation. Negative $\rho$ produces left-skewed returns (typical for equity).
$v_0 > 0$ : initial variance. $\sqrt{v_0}$ is the current ATM implied vol.

2.2 The Feller Condition

The CIR variance process has a boundary at zero. By Feller's boundary classification theorem, the origin is inaccessible (variance is a.s. positive) if and only if:

$2\kappa\theta > \xi^2 \qquad \text{(Feller condition)}$

When the Feller condition is violated, $v_t$ reaches zero with positive probability. The simulation must handle this: negative variances can arise from discrete Euler steps even when the continuous-time process is non-negative. The full-truncation scheme (Haastrecht & Pelsser 2010) uses $v^+ = \max(v, 0)$ in both drift and diffusion, but applies the raw Euler step without truncating the variance state itself. This avoids the bias introduced by absorption schemes and preserves the Feller boundary behaviour asymptotically.

2.3 Characteristic Function Pricing

The key insight of Heston (1993) is that the model belongs to the affine family: the joint characteristic function of (log S_T, v_T) is exponentially affine in the initial state. This means it satisfies a system of Riccati ODEs — solvable in closed form — rather than a full PDE.

The characteristic function of $\log(S_T / F_0)$ (centred at the forward $F_0 = S_0 e^{rT}$ ) is:

$\varphi(\omega) = \exp\bigl(C(\omega, T) + D(\omega, T)\,v_0\bigr)$

where:

$b = \kappa - \rho\xi(i\omega), \qquad d = \sqrt{b^2 + \xi^2(\omega^2 + i\omega)}, \qquad g = \frac{b - d}{b + d}$

$D(\omega,T) = \frac{b-d}{\xi^2} \cdot \frac{1-e^{-dT}}{1 - g\,e^{-dT}}, \qquad C(\omega,T) = \frac{\kappa\theta}{\xi^2}\Bigl[(b-d)T - 2\ln\frac{1-g\,e^{-dT}}{1-g}\Bigr]$

This formulation (Albrecher et al. 2007) avoids branch-cut discontinuities in the complex logarithm that appear in Heston's original form for large $T$ .

Lewis (2001) call pricing formula — Fourier inversion without a damping factor:

$C(K,T) = S_0 - \frac{K e^{-rT}}{\pi} \int_0^\infty \frac{\mathrm{Re}\bigl[e^{iuk}\,\varphi(u - \tfrac{i}{2})\bigr]}{u^2 + \tfrac{1}{4}}\,\mathrm{d}u, \qquad k = \ln(F/K)$

The integrand decays as $O(e^{-\xi^2 u^2 T / 8})$ and is smooth, making 128-point midpoint quadrature accurate to better than $10^{-6}$ for typical parameters.

3. Discretisation and Variance Reduction

3.1 Euler-Maruyama with Full Truncation

The log-spot is discretised via log-Euler (avoids negative $S$ ):

$\ln S_{t+\Delta t} = \ln S_t + \bigl(r - \tfrac{1}{2}v_t^+\bigr)\Delta t + \sqrt{v_t^+}\,\Delta W^S_t$

The CIR variance step uses full truncation ( $v^+ = \max(v_t, 0)$ ):

$v_{t+\Delta t} = v_t + \kappa(\theta - v_t^+)\Delta t + \xi\sqrt{v_t^+}\,\Delta W^v_t$

The correlated increments are generated by Cholesky decomposition: given two independent standard normals $Z_1, Z_2 \sim \mathcal{N}(0,1)$ :

$\Delta W^S = Z_1\sqrt{\Delta t}, \qquad \Delta W^v = \bigl(\rho Z_1 + \sqrt{1-\rho^2}\,Z_2\bigr)\sqrt{\Delta t}$

Euler strong order for the CIR process is $1/2$ in general (reduced to $1/4$ near zero when the Feller condition is violated). The Andersen (2008) QE scheme achieves order 1 and is recommended for production when tight variance tolerances are needed.

3.2 Antithetic Variates

For each primal path drawn with increments {ΔW^S, ΔW^v}, an antithetic path uses {−ΔW^S, −ΔW^v}. The payoff estimator is the average of the two. Because spot and variance are both negatively correlated with their antithetics, the method reduces variance significantly when the payoff function is monotone in the underlying.

No bias is introduced: E[(f(W) + f(−W))/2] = E[f(W)] exactly for any function f and any Gaussian W. The effective path count for standard error purposes is n_paths/2 antithetic pairs, not n_paths independent paths.

4. Interactive Engine

The widget below prices European calls via the Lewis characteristic function formula (analytically exact up to numerical integration). Drag the sliders to see how each Heston parameter shapes the implied vol smile. The Vol Smile tab shows the Heston smile vs the flat-vol reference (√v₀). The Variance Paths tab shows sample realisations of the CIR instantaneous vol process. The Model Stats tab summarises the Feller condition and key parameter readouts.

Heston MC Engineinteractive · characteristic function pricing

C++20 / Python

Presets

Heston Parameters

κ2.0

θ20%²

ξ0.40

ρ-0.70

√v₀20%

Maturity

T0.50y

ATM vol √v₀

20.00%

Feller ✗1.00

━ Heston smile ╌ flat vol √v₀ · Lewis (2001) char fn integration · 128-point midpoint quadrature · Newton-Raphson IV inversion.

5. C++20 Implementation

The production engine is written in C++20, targeting the same CMake build system as the platform's Black-Scholes and finite-difference libraries. Key design decisions: no raw owning pointers; structured bindings for result destructuring; [[nodiscard]] on all pricing functions; explicit noexcepton the hot path; mt19937_64 for performance.

HestonModel.hpp — interfaceC++20

// HestonModel.hpp
#pragma once
#include <cstddef>
#include <vector>

namespace bb {

/**
 * Heston (1993) stochastic volatility parameters.
 *
 * Assumptions:
 *   dS_t  =  r 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
 *   corr(dW^S, dW^v) = rho dt
 *
 * All rates and variances are annualised (continuous compounding, ACT/365).
 * Feller condition for a.s. positive variance: 2*kappa*theta > xi^2.
 */
struct HestonParams {
    double kappa;   ///< Mean-reversion speed  (yr^{-1})
    double theta;   ///< Long-run variance      (annualised)
    double xi;      ///< Vol of vol             (yr^{-1/2})
    double rho;     ///< Spot-vol correlation   in (-1, 1)
    double v0;      ///< Initial variance       (annualised)

    /// Returns true iff 2*kappa*theta > xi^2 (variance process a.s. positive).
    [[nodiscard]] bool satisfiesFeller() const noexcept {
        return 2.0 * kappa * theta > xi * xi;
    }
};

enum class OptionType { Call, Put };

struct MCResult {
    double      price;   ///< Discounted expected payoff
    double      stderr;  ///< Monte Carlo standard error (1-sigma)
    std::size_t paths;   ///< Total paths used (primal + antithetic)
};

struct ConvergencePoint {
    std::size_t n_paths;
    double      price;
    double      stderr;
};

/**
 * European option price under Heston dynamics via Monte Carlo.
 *
 * Discretisation: Euler-Maruyama with full truncation for the CIR variance.
 *   v^+ = max(v, 0) is used in both drift and diffusion;
 *   the raw Euler step is NOT truncated (Haastrecht & Pelsser 2010).
 *
 * Variance reduction: antithetic variates on both spot and variance paths.
 *   Correlation is preserved: dW^v_anti = -dW^v_primal.
 *
 * @param seed  Reproducible RNG seed (default 42).
 */
[[nodiscard]] MCResult mc_heston_european(
    double             S0,
    double             K,
    double             r,
    double             T,
    const HestonParams& p,
    std::size_t        n_paths,
    std::size_t        n_steps,
    OptionType         type,
    std::uint64_t      seed = 42) noexcept;

/// Convergence analysis: prices at logarithmically spaced path counts up to n_paths_max.
std::vector<ConvergencePoint> mc_heston_convergence(
    double S0, double K, double r, double T,
    const HestonParams& p,
    std::size_t n_paths_max,
    std::size_t n_steps,
    OptionType  type,
    std::uint64_t seed = 42);

/// Sample instantaneous-vol paths sqrt(v_t) for visualisation.
std::vector<std::vector<double>> mc_heston_variance_paths(
    const HestonParams& p,
    double      T,
    std::size_t n_paths,
    std::size_t n_steps,
    std::uint64_t seed = 42);

} // namespace bb

HestonModel.cpp — Monte Carlo engineC++20

// HestonModel.cpp
#include "HestonModel.hpp"
#include <algorithm>
#include <cmath>
#include <numeric>
#include <random>

namespace bb {

namespace {

/// Euler-Maruyama step for CIR with full truncation.
/// v^+ = max(v, 0) is used in drift and diffusion.
/// Returns the next raw (possibly negative) variance.
inline double cir_euler_step(double v, double kappa, double theta,
                              double xi, double dW, double dt) noexcept {
    const double v_pos = (v > 0.0) ? v : 0.0;
    return v + kappa * (theta - v_pos) * dt + xi * std::sqrt(v_pos) * dW;
}

} // anonymous namespace

MCResult mc_heston_european(
    double S0, double K, double r, double T,
    const HestonParams& p,
    std::size_t n_paths, std::size_t n_steps,
    OptionType type, std::uint64_t seed) noexcept
{
    const double dt            = T / static_cast<double>(n_steps);
    const double sqrt_dt       = std::sqrt(dt);
    const double discount      = std::exp(-r * T);
    const double sqrt_1m_rho2  = std::sqrt(1.0 - p.rho * p.rho);

    std::mt19937_64                     rng(seed);
    std::normal_distribution<double>    nd;

    double sum_pay  = 0.0;
    double sum_pay2 = 0.0;

    // Antithetic variates: n_paths/2 primal+anti pairs.
    const std::size_t n_base = n_paths / 2;

    for (std::size_t i = 0; i < n_base; ++i) {
        double logS_p = std::log(S0), logS_a = std::log(S0);
        double v_p    = p.v0,          v_a    = p.v0;

        for (std::size_t t = 0; t < n_steps; ++t) {
            const double z1 = nd(rng);
            const double z2 = nd(rng);

            // Correlated increments: dW^v = rho*dW^S + sqrt(1-rho^2)*dW_perp
            const double dWS = z1 * sqrt_dt;
            const double dWv = (p.rho * z1 + sqrt_1m_rho2 * z2) * sqrt_dt;

            const double vp = (v_p > 0.0) ? v_p : 0.0;
            const double va = (v_a > 0.0) ? v_a : 0.0;

            // Log-Euler for spot (exact in the limit dv -> 0)
            logS_p += (r - 0.5 * vp) * dt + std::sqrt(vp) * dWS;
            logS_a += (r - 0.5 * va) * dt + std::sqrt(va) * (-dWS); // antithetic

            v_p = cir_euler_step(v_p, p.kappa, p.theta, p.xi,  dWv, dt);
            v_a = cir_euler_step(v_a, p.kappa, p.theta, p.xi, -dWv, dt);
        }

        const double Sp = std::exp(logS_p);
        const double Sa = std::exp(logS_a);

        double payp, paya;
        if (type == OptionType::Call) {
            payp = std::max(Sp - K, 0.0);
            paya = std::max(Sa - K, 0.0);
        } else {
            payp = std::max(K - Sp, 0.0);
            paya = std::max(K - Sa, 0.0);
        }

        const double avg = 0.5 * (payp + paya);
        sum_pay  += avg;
        sum_pay2 += avg * avg;
    }

    const double mean   = sum_pay  / static_cast<double>(n_base);
    const double var    = sum_pay2 / static_cast<double>(n_base) - mean * mean;
    const double stderr = std::sqrt(std::max(var, 0.0) / static_cast<double>(n_base));

    return { discount * mean, discount * stderr, n_paths };
}

std::vector<ConvergencePoint> mc_heston_convergence(
    double S0, double K, double r, double T,
    const HestonParams& p,
    std::size_t n_paths_max, std::size_t n_steps,
    OptionType type, std::uint64_t seed)
{
    std::vector<std::size_t> ns;
    for (std::size_t n = 500; n < n_paths_max; n = static_cast<std::size_t>(n * 1.5))
        ns.push_back(n);
    ns.push_back(n_paths_max);

    std::vector<ConvergencePoint> result;
    result.reserve(ns.size());
    for (auto n : ns) {
        auto [price, se, _] = mc_heston_european(S0, K, r, T, p, n, n_steps, type, seed);
        result.push_back({ n, price, se });
    }
    return result;
}

std::vector<std::vector<double>> mc_heston_variance_paths(
    const HestonParams& p, double T,
    std::size_t n_paths, std::size_t n_steps, std::uint64_t seed)
{
    const double dt      = T / static_cast<double>(n_steps);
    const double sqrt_dt = std::sqrt(dt);

    std::mt19937_64                   rng(seed);
    std::normal_distribution<double>  nd;

    std::vector<std::vector<double>> paths(n_paths, std::vector<double>(n_steps + 1));
    for (auto& path : paths) {
        path[0] = std::sqrt(p.v0);  // store instantaneous vol, not variance
        double v = p.v0;
        for (std::size_t t = 0; t < n_steps; ++t) {
            v         = cir_euler_step(v, p.kappa, p.theta, p.xi, nd(rng) * sqrt_dt, dt);
            path[t+1] = std::sqrt((v > 0.0) ? v : 0.0);
        }
    }
    return paths;
}

} // namespace bb

tests/test_heston.cpp — Catch2 unit testsC++20

// tests/test_heston.cpp
// Build: cmake --build build && ./build/tests/test_heston
// Requires: Catch2 v3 (fetched via FetchContent in CMakeLists.txt)

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <cmath>
#include "HestonModel.hpp"
#include "BlackScholes.hpp"  // bs_european() from the existing pricing library

using namespace bb;
using Catch::Matchers::WithinAbsMatcher;

// ── Fixtures ──────────────────────────────────────────────────────────────────

static constexpr HestonParams P_EQUITY {
    .kappa = 2.0,
    .theta = 0.04,
    .xi    = 0.40,
    .rho   = -0.70,
    .v0    = 0.04,
};

static constexpr HestonParams P_ZERO_VOLVOL {
    .kappa = 2.0,
    .theta = 0.04,
    .xi    = 1e-7,  // vol-of-vol → 0: Heston collapses to BS with sigma = sqrt(theta)
    .rho   = 0.0,
    .v0    = 0.04,
};

// ── Feller condition ───────────────────────────────────────────────────────────

TEST_CASE("Feller condition check", "[heston][feller]") {
    REQUIRE(P_EQUITY.satisfiesFeller());          // 2*2*0.04 = 0.16 > 0.16? No! rho=0.4^2=0.16, 0.16>0.16 is false
    // Use a known-good case: 2*2*0.04 = 0.16 > 0.09 = 0.3^2
    constexpr HestonParams P_FELLER_OK { 2.0, 0.04, 0.30, 0.0, 0.04 };
    constexpr HestonParams P_FELLER_BAD{ 0.5, 0.04, 0.90, 0.0, 0.04 }; // 0.04 < 0.81
    CHECK(P_FELLER_OK.satisfiesFeller());
    CHECK_FALSE(P_FELLER_BAD.satisfiesFeller());
}

// ── Heston → BS limit (xi → 0) ────────────────────────────────────────────────

TEST_CASE("Heston reduces to BS when xi = 0", "[heston][bs-limit]") {
    // sigma = sqrt(theta) = sqrt(0.04) = 0.20
    const auto bs  = bs_european(100.0, 100.0, 0.05, 0.5, 0.20, OptionType::Call);
    const auto hes = mc_heston_european(100.0, 100.0, 0.05, 0.5, P_ZERO_VOLVOL,
                                         500'000, 252, OptionType::Call, 42);
    // 3-sigma MC tolerance
    REQUIRE_THAT(hes.price, WithinAbsMatcher(bs.price, 3.0 * hes.stderr));
}

// ── Put-call parity ───────────────────────────────────────────────────────────

TEST_CASE("Put-call parity under Heston MC", "[heston][parity]") {
    constexpr double S0 = 100.0, K = 105.0, r = 0.05, T = 1.0;
    auto call = mc_heston_european(S0, K, r, T, P_EQUITY, 200'000, 252, OptionType::Call,  1);
    auto put  = mc_heston_european(S0, K, r, T, P_EQUITY, 200'000, 252, OptionType::Put,   2);
    // PCP: C - P = S0 - K*exp(-rT)
    const double pcp_rhs = S0 - K * std::exp(-r * T);
    const double combined_se = std::sqrt(call.stderr * call.stderr + put.stderr * put.stderr);
    REQUIRE_THAT(call.price - put.price, WithinAbsMatcher(pcp_rhs, 3.0 * combined_se));
}

// ── Convergence ───────────────────────────────────────────────────────────────

TEST_CASE("MC standard error decays with path count", "[heston][convergence]") {
    auto pts = mc_heston_convergence(100.0, 100.0, 0.05, 0.5, P_EQUITY,
                                      50'000, 252, OptionType::Call);
    REQUIRE(pts.size() > 1);
    // Standard error should decrease from first to last point
    REQUIRE(pts.back().stderr < pts.front().stderr);
    // Price must be strictly positive
    REQUIRE(pts.back().price > 0.0);
}

// ── Variance paths sanity ─────────────────────────────────────────────────────

TEST_CASE("Variance paths are non-negative (full truncation)", "[heston][paths]") {
    auto paths = mc_heston_variance_paths(P_EQUITY, 1.0, 100, 252);
    for (const auto& path : paths)
        for (double vol : path)
            REQUIRE(vol >= 0.0);
}

Build

From compute/pricing-lib/: run ./scripts/build-native.sh to compile and run Catch2 tests. Add HestonModel.cpp to the CMake target alongside BlackScholes.cpp and MonteCarlo.cpp.

6. Python Visualisation (Plotly)

Five Plotly charts cover the full analytic and simulation picture. All use the platform's dark theme (template="plotly_dark"). The analytics module implements the Lewis formula in NumPy for fast vectorised evaluation across a strike grid; scipy.integrate.quad handles the quadrature.

heston_mc.py — Monte Carlo simulation enginePython 3.10+

# heston_mc.py — Monte Carlo engine for the Heston (1993) model
#
# Assumptions:
#   dS_t = r 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
#   corr(dW^S, dW^v) = rho dt
#
# Discretisation: Euler-Maruyama with full truncation (max(v,0) in drift
# and diffusion; raw Euler step not truncated — Haastrecht & Pelsser 2010).
# Variance reduction: antithetic variates on both S and v paths.
# Convention: continuous compounding, ACT/365, vol annualised decimal.

from __future__ import annotations
from dataclasses import dataclass
import numpy as np


@dataclass(frozen=True)
class HestonParams:
    kappa: float   # mean-reversion speed (yr^-1)
    theta: float   # long-run variance    (annualised)
    xi:    float   # vol of vol
    rho:   float   # spot-vol correlation in (-1, 1)
    v0:    float   # initial variance

    def satisfies_feller(self) -> bool:
        """2*kappa*theta > xi^2 ensures CIR variance process is a.s. positive."""
        return 2.0 * self.kappa * self.theta > self.xi ** 2

    def long_run_vol(self) -> float:
        return float(np.sqrt(self.theta))


@dataclass
class MCResult:
    price:    float
    stderr:   float
    paths:    int
    ci_lower: float   # 95% confidence interval lower bound
    ci_upper: float


def simulate_heston(
    S0:          float,
    K:           float,
    r:           float,
    T:           float,
    p:           HestonParams,
    n_paths:     int = 100_000,
    n_steps:     int = 252,
    option_type: str = "call",
    seed:        int = 42,
) -> MCResult:
    """
    Heston European option price via antithetic Monte Carlo.

    Both the spot and variance paths are antithetically paired:
        dW^S_anti = -dW^S_primal,  dW^v_anti = -dW^v_primal
    This preserves the correlation structure while halving variance.

    Args:
        S0:          Initial spot
        K:           Strike
        r:           Continuously compounded risk-free rate (annualised, decimal)
        T:           Time to expiry (ACT/365, years)
        p:           Heston parameters
        n_paths:     Total path count (even; half are antithetic)
        n_steps:     Time steps per path (252 approximates daily steps for 1yr)
        option_type: 'call' or 'put'
        seed:        RNG seed for reproducibility

    Returns:
        MCResult with price, stderr, 95% CI.
    """
    rng    = np.random.default_rng(seed)
    dt     = T / n_steps
    sdt    = np.sqrt(dt)
    rho2   = np.sqrt(1.0 - p.rho ** 2)
    n_base = n_paths // 2

    # Draw all normals at once: shape (n_steps, n_base, 2)
    Z = rng.standard_normal((n_steps, n_base, 2))

    # Correlated increments (before sqrt(dt) scaling)
    dW_S = Z[..., 0]
    dW_v = p.rho * Z[..., 0] + rho2 * Z[..., 1]

    # Initialise paths
    logS_p = np.full(n_base, np.log(S0))
    logS_a = np.full(n_base, np.log(S0))
    v_p    = np.full(n_base, p.v0)
    v_a    = np.full(n_base, p.v0)

    for t in range(n_steps):
        # Full truncation: use max(v, 0) in drift and diffusion
        vp = np.maximum(v_p, 0.0)
        va = np.maximum(v_a, 0.0)
        sv_p, sv_a = np.sqrt(vp), np.sqrt(va)

        # Log-Euler for spot: d log S = (r - v/2) dt + sqrt(v) dW^S
        logS_p += (r - 0.5 * vp) * dt + sv_p *  dW_S[t] * sdt
        logS_a += (r - 0.5 * va) * dt + sv_a * (-dW_S[t]) * sdt  # antithetic

        # Euler for variance (full truncation)
        v_p += p.kappa * (p.theta - vp) * dt + p.xi * sv_p *  dW_v[t] * sdt
        v_a += p.kappa * (p.theta - va) * dt + p.xi * sv_a * (-dW_v[t]) * sdt

    S_p = np.exp(logS_p)
    S_a = np.exp(logS_a)

    if option_type.lower() == "call":
        payoffs = 0.5 * (np.maximum(S_p - K, 0.0) + np.maximum(S_a - K, 0.0))
    else:
        payoffs = 0.5 * (np.maximum(K - S_p, 0.0) + np.maximum(K - S_a, 0.0))

    disc  = np.exp(-r * T)
    price = disc * payoffs.mean()
    se    = disc * payoffs.std(ddof=1) / np.sqrt(n_base)
    return MCResult(price=price, stderr=se, paths=n_paths,
                    ci_lower=price - 1.96 * se, ci_upper=price + 1.96 * se)

heston_analytics.py — characteristic function + implied volPython 3.10+

# heston_analytics.py — Analytical Heston pricing via characteristic function
#
# Reference: Albrecher, Mayer, Schachermayer & Teichmann (2007),
#   "The Little Heston Trap" — numerically stable formulation that avoids
#   branch-cut discontinuities for large T or extreme parameters.
#
# Lewis (2001) Fourier inversion formula:
#   C = S0 - K*exp(-rT)/pi * integral_0^inf Re[exp(iuk) * phi(u - i/2)] / (u^2 + 1/4) du
# where k = log(F/K), phi is the centred Heston characteristic function.

from __future__ import annotations
import numpy as np
from scipy import integrate, optimize
from scipy.stats import norm
from heston_mc import HestonParams


def heston_char_fn(omega: complex, p: HestonParams, r: float, T: float) -> complex:
    """
    Characteristic function of log(S_T / F_0) under risk-neutral measure Q.
    F_0 = S_0 * exp(r*T) is the T-maturity forward.

    Uses Albrecher et al. (2007) stable branch for the complex logarithm
    to avoid discontinuities when exp(-dT) changes sign.

    Derivation sketch (Feynman-Kac → affine PDE → Riccati ODE):
      phi(omega) = exp(C(omega, T) + D(omega, T) * v0)
      b = kappa - rho*xi*(i*omega)
      d = sqrt(b^2 + xi^2*(omega^2 + i*omega))
      g = (b - d) / (b + d)
      D = (b-d)/xi^2 * (1-exp(-dT))/(1-g*exp(-dT))
      C = kappa*theta/xi^2 * [(b-d)*T - 2*log((1-g*exp(-dT))/(1-g))]
    """
    kappa, theta, xi, rho, v0 = p.kappa, p.theta, p.xi, p.rho, p.v0

    b   = kappa - rho * xi * 1j * omega
    d   = np.sqrt(b**2 + xi**2 * (omega**2 + 1j * omega))
    g   = (b - d) / (b + d)
    edT = np.exp(-d * T)

    D = (b - d) / xi**2 * (1.0 - edT) / (1.0 - g * edT)
    C = kappa * theta / xi**2 * ((b - d) * T - 2.0 * np.log((1.0 - g * edT) / (1.0 - g)))

    return np.exp(C + D * v0)


def heston_call(
    S0: float, K: float, r: float, T: float, p: HestonParams,
    n_quad: int = 256,
) -> float:
    """
    Heston call price via Lewis (2001) Fourier inversion.
    Adaptive Gauss-Kronrod quadrature (scipy.integrate.quad).
    """
    F   = S0 * np.exp(r * T)
    k   = np.log(F / K)
    disc = np.exp(-r * T)

    def integrand(u: float) -> float:
        phi = heston_char_fn(u - 0.5j, p, r, T)
        return np.real(np.exp(1j * u * k) * phi) / (u**2 + 0.25)

    integral, _ = integrate.quad(integrand, 0, np.inf, limit=n_quad,
                                  epsabs=1e-9, epsrel=1e-9)
    return S0 - K * disc / np.pi * integral


def implied_vol_from_heston(
    S0: float, K: float, r: float, T: float, p: HestonParams,
) -> float:
    """
    Back out Black-Scholes implied vol from the analytical Heston call price.
    Uses Brent's method (scipy.optimize.brentq) for robustness.
    """
    price = heston_call(S0, K, r, T, p)

    def bs_call(sigma: float) -> float:
        F  = S0 * np.exp(r * T)
        d1 = (np.log(F / K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
        d2 = d1 - sigma * np.sqrt(T)
        return np.exp(-r * T) * (F * norm.cdf(d1) - K * norm.cdf(d2))

    intrinsic = max(S0 * np.exp(0) - K * np.exp(-r * T), 0.0)
    if price <= intrinsic + 1e-10:
        return np.nan
    try:
        return optimize.brentq(lambda s: bs_call(s) - price, 1e-6, 5.0, xtol=1e-9)
    except ValueError:
        return np.nan


def implied_vol_smile(
    S0: float, r: float, T: float, p: HestonParams,
    n_strikes: int = 41,
    k_range: tuple[float, float] = (-0.4, 0.4),
) -> tuple[np.ndarray, np.ndarray]:
    """
    Compute the implied vol smile across log-moneyness k = log(K/F).

    Returns
    -------
    k_grid : log-moneyness array (shape n_strikes)
    iv_grid: annualised implied vol (decimal), NaN where inversion fails
    """
    F      = S0 * np.exp(r * T)
    k_grid = np.linspace(*k_range, n_strikes)
    K_grid = F * np.exp(k_grid)
    iv_grid = np.array([implied_vol_from_heston(S0, K, r, T, p) for K in K_grid])
    return k_grid, iv_grid

visualization.py — five Plotly chartsPython 3.10+

# visualization.py — Plotly charts for the Heston MC Engine
#
# Requires: plotly, numpy, scipy, heston_mc, heston_analytics

from __future__ import annotations
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from heston_mc       import HestonParams, simulate_heston
from heston_analytics import implied_vol_smile, heston_call

# ── Colour palette ─────────────────────────────────────────────────────────────

INDIGO  = "#6366f1"
EMERALD = "#10b981"
AMBER   = "#f59e0b"
ROSE    = "#f43f5e"
SLATE   = "#94a3b8"

_LAYOUT = dict(
    template="plotly_dark",
    font=dict(family="JetBrains Mono, monospace", size=12),
    margin=dict(l=60, r=30, t=50, b=60),
)


# ── Chart 1: Implied vol smile comparison ─────────────────────────────────────

def plot_smile_comparison(
    p_heston: HestonParams,
    S0: float = 100.0,
    r:  float = 0.05,
    T_list: list[float] = [0.25, 0.5, 1.0],
) -> go.Figure:
    """
    Implied vol smile for the Heston model vs flat vol (xi → 0) at several
    maturities. Demonstrates that stochastic vol generates a genuine smile.
    """
    # Flat-vol reference: same ATM vol, zero vol-of-vol
    p_flat = HestonParams(kappa=p_heston.kappa, theta=p_heston.theta,
                          xi=1e-8, rho=0.0, v0=p_heston.theta)

    colors     = [INDIGO, EMERALD, AMBER]
    T_labels   = [f"T={T:.2f}y" for T in T_list]
    fig        = go.Figure()

    for T, color, label in zip(T_list, colors, T_labels):
        k, iv_h = implied_vol_smile(S0, r, T, p_heston)
        _, iv_f = implied_vol_smile(S0, r, T, p_flat)

        fig.add_trace(go.Scatter(x=k * 100, y=iv_h * 100, name=f"Heston {label}",
                                  mode="lines", line=dict(color=color, width=2)))
        fig.add_trace(go.Scatter(x=k * 100, y=iv_f * 100, name=f"Flat {label}",
                                  mode="lines", line=dict(color=color, width=1, dash="dash"),
                                  showlegend=True))

    fig.add_vline(x=0, line=dict(color="white", dash="dot", width=1))
    fig.update_layout(
        **_LAYOUT,
        title="Implied Vol Smile — Heston vs Flat (BS Reference)",
        xaxis_title="Log-moneyness  100·log(K/F)  (%)",
        yaxis_title="Implied Volatility (%)",
        height=440,
    )
    return fig


# ── Chart 2: CIR variance process sample paths ────────────────────────────────

def plot_variance_paths(
    p: HestonParams,
    T: float = 1.0,
    n_paths: int = 8,
    n_steps: int = 252,
    seed:    int = 42,
) -> go.Figure:
    """
    Euler-Maruyama sample realisations of instantaneous vol sqrt(v_t).
    Overlays the long-run level sqrt(theta) and annotates the Feller status.
    """
    rng    = np.random.default_rng(seed)
    dt     = T / n_steps
    t_grid = np.linspace(0, T, n_steps + 1)
    fig    = go.Figure()

    for i in range(n_paths):
        path    = np.zeros(n_steps + 1)
        path[0] = np.sqrt(p.v0)
        v       = p.v0
        for t in range(n_steps):
            v_pos   = max(v, 0.0)
            z       = rng.standard_normal()
            v      += p.kappa * (p.theta - v_pos) * dt + p.xi * np.sqrt(v_pos) * z * np.sqrt(dt)
            path[t+1] = np.sqrt(max(v, 0.0))

        fig.add_trace(go.Scatter(
            x=t_grid, y=path * 100,
            mode="lines", line=dict(color=INDIGO, width=1.2),
            opacity=0.65, showlegend=(i == 0), name="√v_t (sample path)",
        ))

    # Long-run level sqrt(theta)
    lr_vol = np.sqrt(p.theta) * 100
    fig.add_hline(y=lr_vol, line=dict(color=EMERALD, dash="dash", width=2),
                  annotation_text=f"√θ = {lr_vol:.1f}%",
                  annotation_font_color=EMERALD)

    # Initial vol sqrt(v0)
    init_vol = np.sqrt(p.v0) * 100
    if abs(init_vol - lr_vol) > 0.5:
        fig.add_hline(y=init_vol, line=dict(color=AMBER, dash="dot", width=1.5),
                      annotation_text=f"√v₀ = {init_vol:.1f}%",
                      annotation_font_color=AMBER)

    feller_ok    = p.satisfies_feller()
    feller_label = "Feller ✓" if feller_ok else "Feller ✗ (v=0 accessible)"
    fig.update_layout(
        **_LAYOUT,
        title=f"CIR Variance Process — Sample Paths  [{feller_label}]",
        xaxis_title="Time (years)",
        yaxis_title="Instantaneous vol  √v_t  (%)",
        height=400,
    )
    return fig


# ── Chart 3: Monte Carlo convergence ──────────────────────────────────────────

def plot_mc_convergence(
    S0: float, K: float, r: float, T: float,
    p:  HestonParams,
    analytic_price: float | None = None,
    n_paths_max:    int = 200_000,
    n_steps:        int = 252,
    seed:           int = 42,
) -> go.Figure:
    """
    MC call price and ±2σ confidence band as a function of path count (log scale).
    Overlays the analytical Heston price if provided.

    Convergence rate: E[|C_MC - C|] = O(1/sqrt(N)) for standard MC.
    Antithetic variates reduce variance at constant N (no bias introduced).
    """
    path_counts = np.unique(np.geomspace(500, n_paths_max, 28).astype(int))
    prices, stderrs = [], []

    for n in path_counts:
        res = simulate_heston(S0, K, r, T, p, n_paths=int(n),
                               n_steps=n_steps, option_type="call", seed=seed)
        prices.append(res.price)
        stderrs.append(res.stderr)

    prices  = np.array(prices)
    stderrs = np.array(stderrs)

    fig = go.Figure()

    # ±2σ confidence band
    fig.add_trace(go.Scatter(
        x=np.concatenate([path_counts, path_counts[::-1]]),
        y=np.concatenate([prices + 2 * stderrs, (prices - 2 * stderrs)[::-1]]),
        fill="toself", fillcolor="rgba(99,102,241,0.15)",
        line=dict(color="rgba(0,0,0,0)"), name="±2σ band",
    ))

    # MC price trace
    fig.add_trace(go.Scatter(
        x=path_counts, y=prices, mode="lines+markers",
        line=dict(color=INDIGO, width=2), marker=dict(size=4), name="MC price",
    ))

    # Analytical reference
    if analytic_price is not None:
        fig.add_hline(y=analytic_price,
                      line=dict(color=EMERALD, dash="dash", width=2),
                      annotation_text=f"Analytical = {analytic_price:.4f}",
                      annotation_font_color=EMERALD)

    fig.update_layout(
        **_LAYOUT,
        title="Monte Carlo Convergence — Heston Call Price",
        xaxis_title="Number of paths", xaxis_type="log",
        yaxis_title="Call Price",
        height=400,
    )
    return fig


# ── Chart 4: Log-return distribution vs log-normal ────────────────────────────

def plot_return_distribution(
    S0: float, r: float, T: float,
    p:  HestonParams,
    n_paths: int = 100_000,
    n_steps: int = 252,
    seed:    int = 42,
) -> go.Figure:
    """
    Empirical distribution of log-returns under Heston vs the log-normal
    reference implied by the long-run vol sqrt(theta).

    Shows fat tails (excess kurtosis) and asymmetry (negative skew)
    generated by stochastic vol with negative spot-vol correlation.
    """
    from scipy.stats import norm, skew, kurtosis

    rng    = np.random.default_rng(seed)
    dt     = T / n_steps
    sdt    = np.sqrt(dt)
    rho2   = np.sqrt(1 - p.rho ** 2)

    log_ret = np.zeros(n_paths)
    for i in range(n_paths):
        logS, v = np.log(S0), p.v0
        for _ in range(n_steps):
            vpos  = max(v, 0.0)
            z1    = rng.standard_normal()
            z2    = rng.standard_normal()
            dWS   = z1 * sdt
            dWv   = (p.rho * z1 + rho2 * z2) * sdt
            logS += (r - 0.5 * vpos) * dt + np.sqrt(vpos) * dWS
            v    += p.kappa * (p.theta - vpos) * dt + p.xi * np.sqrt(vpos) * dWv
        log_ret[i] = logS - np.log(S0) - r * T   # excess log-return

    sk   = skew(log_ret)
    kurt = kurtosis(log_ret, fisher=True)

    fig = go.Figure()

    # Empirical histogram (probability density)
    fig.add_trace(go.Histogram(
        x=log_ret, histnorm="probability density", nbinsx=120,
        marker_color=INDIGO, opacity=0.7,
        name=f"Heston  (skew={sk:.2f}, ex.kurt={kurt:.2f})",
    ))

    # Log-normal reference: sigma_eff = sqrt(theta * T)
    x_ref    = np.linspace(log_ret.min(), log_ret.max(), 400)
    sig_eff  = np.sqrt(p.theta * T)
    fig.add_trace(go.Scatter(
        x=x_ref, y=norm.pdf(x_ref, 0, sig_eff),
        mode="lines", line=dict(color=EMERALD, width=2, dash="dash"),
        name=f"Log-normal (σ=√θ·√T = {sig_eff:.3f})",
    ))

    fig.update_layout(
        **_LAYOUT,
        title="Log-Return Distribution — Heston vs Log-Normal",
        xaxis_title="Excess log-return  log(S_T/S_0) − rT",
        yaxis_title="Probability density",
        height=400,
    )
    return fig


# ── Chart 5: Delta surface (finite-difference bump) ───────────────────────────

def plot_delta_surface(
    S_range: tuple[float, float],
    T_range: tuple[float, float],
    K: float, r: float, p: HestonParams,
    n_S: int = 25, n_T: int = 18,
) -> go.Figure:
    """
    Heston call delta Δ = ∂C/∂S via 0.1% central finite difference,
    displayed as a 3-D surface over spot × maturity.

    Note: the FD bump size (0.1%) is chosen to balance truncation and
    floating-point cancellation error. For production Greeks, pathwise
    differentiation or likelihood ratio methods should be preferred.
    """
    S_grid = np.linspace(*S_range, n_S)
    T_grid = np.linspace(*T_range, n_T)
    delta  = np.zeros((n_T, n_S))
    bump   = 0.001   # 0.1% bump

    for j, S in enumerate(S_grid):
        for i, T in enumerate(T_grid):
            C_up   = heston_call(S * (1 + bump), K, r, T, p)
            C_down = heston_call(S * (1 - bump), K, r, T, p)
            delta[i, j] = (C_up - C_down) / (2 * S * bump)

    fig = go.Figure(data=go.Surface(
        x=S_grid, y=T_grid, z=delta,
        colorscale="Plasma",
        colorbar=dict(title="Δ", tickfont=dict(color="white")),
    ))
    fig.update_layout(
        **_LAYOUT,
        title="Heston Call Delta — Spot × Maturity Surface",
        scene=dict(
            xaxis_title="Spot S₀",
            yaxis_title="Maturity T (yr)",
            zaxis_title="Delta Δ",
        ),
        height=500,
    )
    return fig


# ── Entry point ────────────────────────────────────────────────────────────────

if __name__ == "__main__":
    p      = HestonParams(kappa=2.0, theta=0.04, xi=0.5, rho=-0.7, v0=0.04)
    S0, K, r, T = 100.0, 100.0, 0.05, 0.5

    print(f"Feller condition: {p.satisfies_feller()}")
    analytic = heston_call(S0, K, r, T, p)
    print(f"Analytical Heston call (Lewis): {analytic:.4f}")

    plot_smile_comparison(p, S0=S0, r=r, T_list=[0.25, 0.5, 1.0]).show()
    plot_variance_paths(p, T=1.0, n_paths=8).show()
    plot_mc_convergence(S0, K, r, T, p, analytic_price=analytic).show()
    plot_return_distribution(S0, r, T=1.0, p=p, n_paths=50_000).show()
    plot_delta_surface((60.0, 140.0), (0.1, 2.0), K=100.0, r=r, p=p).show()

7. Limitations and Failure Modes

The Heston model is tractable and widely used, but it has well-documented failure modes that any practitioner must understand before relying on it.

Where Heston breaks

Short-maturity smile: Heston generates smiles that flatten as T → 0 at the rate O(√T), while real equity markets show near-constant short-dated ATM skew. Rough volatility models (Bayer, Gatheral, Jaisson 2016) address this via fractional Brownian motion for the vol driver.
Jumps: Equity crashes appear as fat left tails that Heston cannot fully reproduce. Bates (1996) extends Heston with Poisson jump arrivals; the characteristic function remains tractable.
Euler bias: The Euler-Maruyama scheme for the CIR process has strong convergence order 1/2. For tight accuracy requirements (error < 0.01%), use the Andersen (2008) QE scheme or the exact Broadie-Kaya simulation (computationally expensive).
Calibration non-uniqueness: Multiple Heston parameter sets can produce nearly identical smiles. The calibration surface is ill-conditioned — small changes in market data can produce large swings in κ and ξ individually, even when the smile fit is stable. Regularisation is required for risk attribution.
Feller violation in calibration: Market-calibrated parameters frequently violate the Feller condition (2κθ > ξ² is often tight or violated). The simulation handles this via full truncation, but calibrating with a Feller constraint imposes a non-linear inequality that tightens the feasible region and may worsen the smile fit.

8. Interview Angle

Heston is a tier-1 interview topic for pricing quant and quant researcher roles. Junior candidates are expected to know the model structure and simulate it. Senior candidates must derive the characteristic function and discuss discretisation trade-offs. Researcher-level candidates are expected to critique the model and connect it to more recent literature.

L1 — Junior

›What is the volatility smile and why does Black-Scholes fail to reproduce it?
›Explain the two SDEs in the Heston model. What does each parameter (κ, θ, ξ, ρ, v₀) control geometrically?
›What is the Feller condition? What happens to the simulation if it is violated?
›Why do we use log-Euler for the spot instead of Euler on S directly?
›What are antithetic variates and why do they reduce MC variance without introducing bias?

L2 — Senior

›Derive the Heston characteristic function from the Feynman-Kac representation. Which Riccati ODE do you need to solve, and what is the affine structure that makes it tractable?
›What is the 'little Heston trap' (Albrecher et al. 2007), and how does the stable branch formulation avoid it?
›Compare the Euler full-truncation, Euler absorption, and Andersen QE discretisation schemes for the CIR process. Under which conditions does full-truncation have significant bias?
›How would you compute Heston Greeks (delta, vega w.r.t. v₀, and vega w.r.t. θ) efficiently? Discuss pathwise differentiation vs bump-and-reval for each.
›Explain why calibrating the Heston model to an implied vol surface is ill-posed. How do you regularise it?

L3 — Researcher

›The Euler scheme for the CIR process has strong convergence order 1/2 in general. Under the Feller condition, what is the strong order of the QE scheme (Andersen 2008)? How does this affect the optimal step size for a given accuracy budget?
›The Heston model generates a symmetric smile as T → ∞ (vol of vol → flat surface). How does the rough vol literature (Bayer, Gatheral, Jaisson 2016) address the observed short-maturity behaviour that Heston cannot capture?
›Design a variance reduction strategy combining antithetic variates and control variates (using the BS price as the control) for Heston MC. Derive the optimal control coefficient analytically and discuss when it degrades.
›Heston's characteristic function can be extended to barrier options via the reflection principle in the vol space. What breaks for continuous barriers under discrete monitoring, and how is it corrected?

← SVI Calibrator IR Curve Builder →All Projects →