Monte Carlo Pricing Framework in C++

Setup

The Pricing Problem

Before writing any code, fix the mathematical framework precisely.

We work on a filtered probability space $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \geq 0}, \mathbb{Q})$ where $\mathbb{Q}$ is the risk-neutral (pricing) measure. Under $\mathbb{Q}$, the underlying asset price $S_t$ follows the SDE:

$dS_t = \mu(t, S_t)\,dt + \sigma(t, S_t)\,dW_t^{\mathbb{Q}}$

The model is entirely determined by the pair of functions $\{(t,s) \mapsto \mu(t,s),\; (t,s) \mapsto \sigma(t,s)\}$. For Black-Scholes: $\mu(t,s) = rs$ and $\sigma(t,s) = \sigma_{\text{BS}} s$. For Dupire local volatility: $\mu(t,s) = rs$ and $\sigma(t,s) = \sigma_D(t,s) s$ where $\sigma_D$ is read from a calibrated surface.

The general pricing formula for a derivative with payoff functional $\psi$ over the path $(S_t)_{t \in [0,T]}$ is:

$\text{PV}(0) = \mathbb{E}^{\mathbb{Q}}\!\left[D(0,T)\cdot\psi\!\left((S_t)_{t\in[0,T]}\right)\right]$

where $D(0,T) = e^{-\int_0^T r(s)\,ds}$ is the stochastic discount factor. For constant $r$: $D(0,T) = e^{-rT}$.

The Monte Carlo approximation over $N$ independent paths:

$\text{PV}(0) \approx \frac{1}{N}\sum_{j=1}^{N} D(0,T)\cdot\psi\!\left(S^{(j)}_{t_0},\ldots,S^{(j)}_{t_m}\right)$

Three independent components are required:

1. N — the number of simulation paths. Determines statistical precision: standard error $\propto 1/\sqrt{N}$.
2. Payoff functor $\psi$ — a callable that maps a simulated path $(S_{t_0},\ldots,S_{t_m})$ to a real number. Separates the product specification from the simulation.
3. PathSimulator — generates the path $(S_{t_0}^{(j)},\ldots,S_{t_m}^{(j)})$ for each $j$. Encapsulates the discretisation scheme and the model.

Design Hierarchy

The object-oriented structure isolates each concern:

Model (abstract)
  ├── BlackScholesModel        μ = r·s,  σ = σ_BS·s
  └── DupireLocalVolModel      μ = r·s,  σ = σ_D(t,s)·s

PathSimulator (abstract)       owns a Model* via clone()
  ├── EulerPathSimulator       Euler-Maruyama discretisation
  └── MilsteinPathSimulator    Milstein correction

Payoff (abstract)              callable on std::vector<double>
  ├── CallPayoff               max(S_T - K, 0)
  ├── PutPayoff                max(K - S_T, 0)
  └── AsianCallPayoff          max(mean(S) - K, 0)

MonteCarlo                     engine; owns PathSimulator* via clone()

This hierarchy ensures that any combination of model, scheme, and payoff can be assembled without changing existing code — the Open/Closed Principle applied to a pricing context.

Conventions used throughout this module:

• Rate $r$: continuously compounded, annualised.
• Volatility $\sigma$: annualised, decimal (0.20 = 20%).
• Time $T$: in years.
• All prices in the same currency unit.

Theory

1. Model Abstract Class

The Model abstraction separates the SDE specification from the numerical scheme. Any class implementing drift_term and diffusion_term can be plugged into any PathSimulator without modification.

Two pure virtual methods are required:

• drift_term(double t, double s) — returns $\mu(t, s)$, the drift coefficient.
• diffusion_term(double t, double s) — returns $\sigma(t, s)$, the diffusion coefficient.

A clone() method returns a heap-allocated copy. This is the prototype pattern: the PathSimulator constructor takes a const Model& and immediately calls clone() to store its own owned copy. Without clone(), the simulator would hold a reference to an object that might go out of scope or be modified by the caller after construction.

Virtual destructor is mandatory whenever a class is intended to be used polymorphically through a base pointer. Omitting it causes undefined behaviour when a derived object is deleted through a Model*.

Black-Scholes model: $\mu(t, s) = r \cdot s, \qquad \sigma(t, s) = \sigma_{\text{BS}} \cdot s$

This is GBM (Geometric Brownian Motion). The drift is linear in $s$; the diffusion is also linear in $s$. For constant coefficients, the SDE has the exact solution $S_T = S_0 \exp\!\left[(r - \tfrac{1}{2}\sigma^2)T + \sigma W_T\right]$, which can be simulated in a single step. We retain the general scheme for generality.

Dupire local volatility model: $\mu(t, s) = r \cdot s, \qquad \sigma(t, s) = \sigma_D(t, s) \cdot s$

where $\sigma_D(t,s)$ is the local volatility function extracted from the implied vol surface (see Module 6). The model matches all market call prices by construction; the function $\sigma_D$ is not constant, so the exact solution above is not available and discretisation is mandatory.

2. Euler-Maruyama Scheme

Discretise $[0,T]$ uniformly: $t_0 = 0 < t_1 < \cdots < t_m = T$, with step $\Delta t = T/m$.

The Euler-Maruyama update at step $i$:

$S_{t_{i+1}} = S_{t_i} + \mu(t_i, S_{t_i})\,\Delta t + \sigma(t_i, S_{t_i})\,\Delta W_i$

where $\Delta W_i = W_{t_{i+1}} - W_{t_i} \sim \mathcal{N}(0, \Delta t)$. In practice: $\Delta W_i = \sqrt{\Delta t}\, Z_i$ with $Z_i \overset{\text{iid}}{\sim} \mathcal{N}(0,1)$.

Convergence properties:

For option pricing, weak convergence governs the bias in the price estimate. For smooth payoffs (vanilla calls, puts), the pricing error shrinks as $\mathcal{O}(\Delta t)$. For discontinuous payoffs (digital options, barrier options), the effective weak rate degrades.

Assumptions for convergence: $\mu$ and $\sigma$ satisfy a global Lipschitz condition in $s$ and at most linear growth. GBM satisfies both ($\mu = rs$, $\sigma = \sigma s$ are globally Lipschitz). CEV with $\beta > 1$ violates linear growth and can cause Euler to produce negative prices — see Limitations.

3. Milstein Scheme

The Milstein scheme adds one term to the Euler update, derived by applying Itô's lemma to the diffusion coefficient $\sigma(t,S_t)$ itself:

$S_{t_{i+1}} = S_{t_i} + \mu(t_i, S_{t_i})\,\Delta t + \sigma(t_i, S_{t_i})\,\Delta W_i + \frac{1}{2}\sigma(t_i, S_{t_i})\,\frac{\partial \sigma}{\partial s}(t_i, S_{t_i})\left[(\Delta W_i)^2 - \Delta t\right]$

The additional term $\tfrac{1}{2}\sigma \cdot \partial_s\sigma \cdot [(\Delta W)^2 - \Delta t]$ captures the curvature of the diffusion coefficient with respect to the state — the contribution of quadratic variation to the path that Euler ignores.

For GBM: $\sigma(t,s) = \sigma s$, so $\partial\sigma/\partial s = \sigma$. The correction term becomes:

$\frac{1}{2}\sigma^2 S_{t_i}\left[(\Delta W_i)^2 - \Delta t\right]$

Convergence:

Milstein's strong order $\mathcal{O}(\Delta t)$ means the pathwise error is squared relative to Euler at the same step size. This matters for path-dependent options where the quality of each individual path, not just the average, affects the result.

Requirement: $\sigma$ must be differentiable in $s$. For models where $\partial_s\sigma$ does not exist (e.g., Heston where the variance process has a square root diffusion), Milstein requires a transformation or an alternative formulation.

4. Random Number Generator

Wrap std::mt19937_64 (Mersenne Twister, 64-bit state) seeded at construction. Produce standard normal samples via std::normal_distribution<double>.

The generator is a value member of PathSimulator — not a pointer, not a global. This ensures:

• No heap allocation per path for random number generation.
• Each PathSimulator instance has independent state (important if multiple simulators are constructed for parallel execution).
• Reproducibility: seeding with a fixed value yields identical paths across runs, which is mandatory for unit tests.

Production note: A fixed seed makes tests deterministic but all instances seeded identically will produce identical sequences. For parallel Monte Carlo, use a seed sequence (std::seed_seq) derived from an independent source per thread.

5. PathSimulator Architecture

PathSimulator holds:

• A time grid std::vector<double> _timePoints storing $(t_0, t_1, \ldots, t_m)$.
• A const Model* _model — a heap-allocated clone of the model passed at construction.
• A std::mt19937_64 _mt — the random number generator.
• A std::normal_distribution<double> _normalDist.

The constructor signature is PathSimulator(const Model& model, const std::vector<double>& timePoints). It calls model.clone() to store _model. The destructor deletes _model.

The pure virtual method simulate_path(double s0) returns std::vector<double> — the simulated path values at each time point. Derived classes implement the specific scheme.

6. Monte Carlo Engine

The MonteCarlo class stores an initial spot _s0, a discount factor _discount, and a PathSimulator* _simulator (owned via clone). Its pricing method:

price(N, payoff):
    sum = 0
    for j = 1..N:
        path = _simulator->simulate_path(_s0)
        sum += payoff(path)
    return _discount * sum / N

Standard error of the estimate: $\hat{\text{SE}} = \hat{\sigma}_{\psi} / \sqrt{N}$, where $\hat{\sigma}_{\psi}^2 = \frac{1}{N-1}\sum_{j=1}^N [\psi_j - \bar{\psi}]^2$ is the sample variance of the discounted payoffs. A 95% confidence interval: $\bar{\psi} \pm 1.96\,\hat{\text{SE}}$.

Implementation

All files target C++17 with -Wall -Wextra -Wpedantic -Werror.

Model.h

// Model.h
// Abstract base class for SDE drift and diffusion coefficients.
// Compile: part of a multi-file project; see CMakeLists.txt or main.cpp build comment.
#pragma once

class Model {
public:
    // Pure virtual: drift coefficient μ(t, s) in dS = μ dt + σ dW
    virtual double drift_term(double t, double s) const = 0;

    // Pure virtual: diffusion coefficient σ(t, s)
    virtual double diffusion_term(double t, double s) const = 0;

    // Prototype pattern: returns a heap-allocated copy of the concrete model.
    // PathSimulator calls this so it owns an independent copy of the model.
    virtual Model* clone() const = 0;

    // Virtual destructor: required for correct cleanup through a Model* pointer.
    virtual ~Model() = default;
};

// ---------------------------------------------------------------------------
// Black-Scholes model: dS = r·S dt + σ·S dW
// ---------------------------------------------------------------------------
class BlackScholesModel : public Model {
public:
    // rate: continuously compounded, annualised risk-free rate
    // sigma: annualised volatility (decimal: 0.20 = 20%)
    BlackScholesModel(double rate, double sigma);

    double drift_term(double t, double s) const override;
    double diffusion_term(double t, double s) const override;
    BlackScholesModel* clone() const override;

private:
    double _rate;
    double _sigma;
};

Model.cpp

// Model.cpp
#include "Model.h"

BlackScholesModel::BlackScholesModel(double rate, double sigma)
    : _rate(rate), _sigma(sigma) {}

double BlackScholesModel::drift_term(double /*t*/, double s) const {
    // μ(t, s) = r·s  — linear drift; risk-neutral measure removes risk premium
    return _rate * s;
}

double BlackScholesModel::diffusion_term(double /*t*/, double s) const {
    // σ(t, s) = σ_BS·s  — proportional diffusion gives lognormal paths
    return _sigma * s;
}

BlackScholesModel* BlackScholesModel::clone() const {
    return new BlackScholesModel(_rate, _sigma);
}

PathSimulator.h

// PathSimulator.h
// Abstract base for path discretisation schemes.
// Concrete classes implement Euler-Maruyama and Milstein.
#pragma once

#include "Model.h"
#include <vector>
#include <random>

class PathSimulator {
public:
    // model: the SDE model — a clone is stored internally.
    // timePoints: the time grid {t_0=0, t_1, ..., t_m=T}.
    // seed: RNG seed for reproducibility; default 42 for unit tests.
    PathSimulator(const Model& model,
                  const std::vector<double>& timePoints,
                  unsigned long long seed = 42ULL);

    // Simulate one path from s0 and return the path at each time point.
    // Returns std::vector<double> of size timePoints.size().
    virtual std::vector<double> simulate_path(double s0) const = 0;

    virtual PathSimulator* clone() const = 0;
    virtual ~PathSimulator();

protected:
    const Model*         _model;       // heap-allocated clone; owned
    std::vector<double>  _timePoints;  // t_0, t_1, ..., t_m

    // Mutable: simulate_path is logically const (same distribution) but
    // advances the RNG state. Declared mutable to allow const simulate_path.
    mutable std::mt19937_64                    _mt;
    mutable std::normal_distribution<double>   _normalDist;
};