Brownian Bridge™

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:

N — the number of simulation paths. Determines statistical precision: standard error $\propto 1/\sqrt{N}$ .
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.
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:

Mode	Order	What it measures
Strong	$\mathcal{O}(\sqrt{\Delta t})$	Pathwise accuracy: $\mathbb{E}[\\|S_T^{\text{exact}} - S_T^{\text{EM}}\\|]$
Weak	$\mathcal{O}(\Delta t)$	Functional accuracy: $

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:

Mode	Order	Improvement over Euler
Strong	$\mathcal{O}(\Delta t)$	One full order gained
Weak	$\mathcal{O}(\Delta t)$	Same as Euler

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

// ---------------------------------------------------------------------------
// Euler-Maruyama: S_{i+1} = S_i + μ(t_i, S_i)·Δt + σ(t_i, S_i)·ΔW_i
// Strong convergence: O(√Δt).  Weak convergence: O(Δt).
// ---------------------------------------------------------------------------
class EulerPathSimulator : public PathSimulator {
public:
    EulerPathSimulator(const Model& model,
                       const std::vector<double>& timePoints,
                       unsigned long long seed = 42ULL);

    std::vector<double> simulate_path(double s0) const override;
    EulerPathSimulator* clone() const override;
};

// ---------------------------------------------------------------------------
// Milstein: Euler + ½·σ·∂σ/∂s·[(ΔW)² − Δt]
// Strong convergence: O(Δt).  Requires σ differentiable in s.
// For GBM: ∂σ/∂s = σ_BS, so correction = ½·σ_BS²·S_i·[(ΔW)² − Δt].
// ---------------------------------------------------------------------------
class MilsteinPathSimulator : public PathSimulator {
public:
    MilsteinPathSimulator(const Model& model,
                          const std::vector<double>& timePoints,
                          unsigned long long seed = 42ULL);

    std::vector<double> simulate_path(double s0) const override;
    MilsteinPathSimulator* clone() const override;

private:
    // Finite-difference approximation of ∂σ/∂s at (t, s).
    // Used to make the Milstein scheme applicable to any Model,
    // not just Black-Scholes where ∂σ/∂s is known analytically.
    // Step h chosen to balance truncation and cancellation errors.
    static constexpr double FD_H = 1e-4;  // relative perturbation for FD
    double diff_derivative(double t, double s) const;
};

PathSimulator.cpp

// PathSimulator.cpp
#include "PathSimulator.h"
#include <cmath>
#include <stdexcept>

// ---------------------------------------------------------------------------
// PathSimulator base
// ---------------------------------------------------------------------------
PathSimulator::PathSimulator(const Model& model,
                             const std::vector<double>& timePoints,
                             unsigned long long seed)
    : _model(model.clone()),
      _timePoints(timePoints),
      _mt(seed),
      _normalDist(0.0, 1.0)
{
    if (timePoints.size() < 2)
        throw std::invalid_argument("PathSimulator: time grid must have at least 2 points");
    if (timePoints.front() != 0.0)
        throw std::invalid_argument("PathSimulator: first time point must be 0");
}

PathSimulator::~PathSimulator() {
    delete _model;
}

// ---------------------------------------------------------------------------
// EulerPathSimulator
// ---------------------------------------------------------------------------
EulerPathSimulator::EulerPathSimulator(const Model& model,
                                       const std::vector<double>& timePoints,
                                       unsigned long long seed)
    : PathSimulator(model, timePoints, seed) {}

std::vector<double> EulerPathSimulator::simulate_path(double s0) const {
    const std::size_t m = _timePoints.size();
    std::vector<double> path(m);
    path[0] = s0;

    for (std::size_t i = 0; i + 1 < m; ++i) {
        const double t  = _timePoints[i];
        const double dt = _timePoints[i + 1] - _timePoints[i];

        // Draw ΔW ~ N(0, Δt) = √Δt · Z, Z ~ N(0,1)
        const double dW = std::sqrt(dt) * _normalDist(_mt);

        const double s = path[i];
        // Euler-Maruyama: S_{i+1} = S_i + μ(t_i, S_i)·Δt + σ(t_i, S_i)·ΔW_i
        path[i + 1] = s
                    + _model->drift_term(t, s) * dt
                    + _model->diffusion_term(t, s) * dW;
    }
    return path;
}

EulerPathSimulator* EulerPathSimulator::clone() const {
    // Construct from the stored model clone; reuse same time grid.
    // A new RNG seed is used so the clone produces independent paths.
    return new EulerPathSimulator(*_model, _timePoints);
}

// ---------------------------------------------------------------------------
// MilsteinPathSimulator
// ---------------------------------------------------------------------------
MilsteinPathSimulator::MilsteinPathSimulator(const Model& model,
                                             const std::vector<double>& timePoints,
                                             unsigned long long seed)
    : PathSimulator(model, timePoints, seed) {}

double MilsteinPathSimulator::diff_derivative(double t, double s) const {
    // Central finite difference: ∂σ/∂s ≈ [σ(t, s+h) − σ(t, s−h)] / (2h)
    // Use relative step to preserve accuracy across different price scales.
    const double h = FD_H * s;
    return (_model->diffusion_term(t, s + h) - _model->diffusion_term(t, s - h))
           / (2.0 * h);
}

std::vector<double> MilsteinPathSimulator::simulate_path(double s0) const {
    const std::size_t m = _timePoints.size();
    std::vector<double> path(m);
    path[0] = s0;

    for (std::size_t i = 0; i + 1 < m; ++i) {
        const double t  = _timePoints[i];
        const double dt = _timePoints[i + 1] - _timePoints[i];
        const double Z  = _normalDist(_mt);
        const double dW = std::sqrt(dt) * Z;

        const double s     = path[i];
        const double mu    = _model->drift_term(t, s);
        const double sigma = _model->diffusion_term(t, s);
        const double dsigma_ds = diff_derivative(t, s);

        // Milstein update:
        // S_{i+1} = S_i + μ·Δt + σ·ΔW + ½·σ·(∂σ/∂s)·[(ΔW)² − Δt]
        const double milstein_correction =
            0.5 * sigma * dsigma_ds * (dW * dW - dt);

        path[i + 1] = s + mu * dt + sigma * dW + milstein_correction;
    }
    return path;
}

MilsteinPathSimulator* MilsteinPathSimulator::clone() const {
    return new MilsteinPathSimulator(*_model, _timePoints);
}

Payoff.h

// Payoff.h
// Abstract payoff functor and concrete vanilla implementations.
// A payoff callable takes the full simulated path and returns the
// undiscounted payoff. Discounting is applied by the MonteCarlo engine.
#pragma once

#include <vector>
#include <algorithm>  // std::max

class Payoff {
public:
    // path: simulated asset prices at each time grid point.
    // Returns the undiscounted payoff at expiry (or along the path for exotics).
    virtual double operator()(const std::vector<double>& path) const = 0;

    virtual Payoff* clone() const = 0;
    virtual ~Payoff() = default;
};

// European call: max(S_T − K, 0)
class CallPayoff : public Payoff {
public:
    explicit CallPayoff(double strike) : _strike(strike) {}
    double operator()(const std::vector<double>& path) const override {
        // path.back() = S_T
        return std::max(path.back() - _strike, 0.0);
    }
    CallPayoff* clone() const override { return new CallPayoff(_strike); }
private:
    double _strike;
};

// European put: max(K − S_T, 0)
class PutPayoff : public Payoff {
public:
    explicit PutPayoff(double strike) : _strike(strike) {}
    double operator()(const std::vector<double>& path) const override {
        return std::max(_strike - path.back(), 0.0);
    }
    PutPayoff* clone() const override { return new PutPayoff(_strike); }
private:
    double _strike;
};

// Asian arithmetic call: max(mean(S_{t_1},...,S_{t_m}) − K, 0)
// Arithmetic mean of all path values including S_0; for a path starting at t_0=0.
class AsianCallPayoff : public Payoff {
public:
    explicit AsianCallPayoff(double strike) : _strike(strike) {}
    double operator()(const std::vector<double>& path) const override;
    AsianCallPayoff* clone() const override { return new AsianCallPayoff(_strike); }
private:
    double _strike;
};

Payoff.cpp

// Payoff.cpp
#include "Payoff.h"
#include <numeric>  // std::accumulate

double AsianCallPayoff::operator()(const std::vector<double>& path) const {
    if (path.empty()) return 0.0;
    // Arithmetic mean over the entire path including S_0
    const double mean = std::accumulate(path.begin(), path.end(), 0.0)
                        / static_cast<double>(path.size());
    return std::max(mean - _strike, 0.0);
}

MonteCarlo.h

// MonteCarlo.h
// Monte Carlo pricing engine.
// Owns a PathSimulator (via clone); does not own the Payoff (caller manages lifetime).
#pragma once

#include "PathSimulator.h"
#include "Payoff.h"
#include <cmath>  // std::sqrt

class MonteCarlo {
public:
    // s0:        initial spot price
    // discount:  discount factor D(0,T) = exp(-r·T)
    // simulator: the path generator — cloned internally
    MonteCarlo(double s0, double discount, const PathSimulator& simulator);

    MonteCarlo(const MonteCarlo& other);
    MonteCarlo& operator=(const MonteCarlo& other);
    ~MonteCarlo();

    // price: estimate E[payoff] and multiply by discount.
    // Returns the Monte Carlo price estimate.
    double price(int N, const Payoff& payoff) const;

    // price_with_stderr: additionally computes the standard error.
    // se_out: output parameter for the standard error estimate.
    double price_with_stderr(int N, const Payoff& payoff, double& se_out) const;

private:
    double           _s0;
    double           _discount;
    PathSimulator*   _simulator;  // heap-allocated clone; owned
};

MonteCarlo.cpp

// MonteCarlo.cpp
#include "MonteCarlo.h"
#include <cmath>
#include <stdexcept>

MonteCarlo::MonteCarlo(double s0, double discount, const PathSimulator& simulator)
    : _s0(s0), _discount(discount), _simulator(simulator.clone()) {}

MonteCarlo::MonteCarlo(const MonteCarlo& other)
    : _s0(other._s0),
      _discount(other._discount),
      _simulator(other._simulator->clone()) {}

MonteCarlo& MonteCarlo::operator=(const MonteCarlo& other) {
    if (this != &other) {
        delete _simulator;
        _s0       = other._s0;
        _discount = other._discount;
        _simulator = other._simulator->clone();
    }
    return *this;
}

MonteCarlo::~MonteCarlo() {
    delete _simulator;
}

double MonteCarlo::price(int N, const Payoff& payoff) const {
    if (N <= 0)
        throw std::invalid_argument("MonteCarlo::price: N must be positive");
    double sum = 0.0;
    for (int j = 0; j < N; ++j) {
        const auto path = _simulator->simulate_path(_s0);
        sum += payoff(path);
    }
    return _discount * sum / static_cast<double>(N);
}

double MonteCarlo::price_with_stderr(int N, const Payoff& payoff, double& se_out) const {
    if (N <= 1)
        throw std::invalid_argument("MonteCarlo::price_with_stderr: N must be > 1");

    double sum  = 0.0;
    double sum2 = 0.0;
    for (int j = 0; j < N; ++j) {
        const auto path = _simulator->simulate_path(_s0);
        const double v  = payoff(path);
        sum  += v;
        sum2 += v * v;
    }
    const double dN   = static_cast<double>(N);
    const double mean = sum / dN;
    // Sample variance of the payoff
    const double var  = (sum2 / dN - mean * mean) * (dN / (dN - 1.0));
    // Standard error of the MC estimate (before discount)
    se_out = _discount * std::sqrt(var / dN);
    return _discount * mean;
}

main.cpp

// main.cpp
// European call: MC price (Euler + Milstein) vs Black-Scholes analytic.
// Compile (single-file convenience build):
//   g++ -std=c++17 -Wall -Wextra -Wpedantic -Werror -O2 \
//       Model.cpp PathSimulator.cpp Payoff.cpp MonteCarlo.cpp main.cpp \
//       -o mc_pricer
// Run: ./mc_pricer
#include "Model.h"
#include "PathSimulator.h"
#include "Payoff.h"
#include "MonteCarlo.h"

#include <iostream>
#include <iomanip>
#include <cmath>
#include <numeric>
#include <vector>

// ---------------------------------------------------------------------------
// Black-Scholes closed-form call price (Abramowitz & Stegun N(x) approximation)
// Assumptions: European, no dividends, constant r and σ, continuous compounding.
// ---------------------------------------------------------------------------
static double norm_cdf(double x) {
    // Horner form of A&S 26.2.17, max |ε| < 7.5e-8
    const double a1 =  0.319381530;
    const double a2 = -0.356563782;
    const double a3 =  1.781477937;
    const double a4 = -1.821255978;
    const double a5 =  1.330274429;
    const double k  = 1.0 / (1.0 + 0.2316419 * std::abs(x));
    const double poly = k*(a1 + k*(a2 + k*(a3 + k*(a4 + k*a5))));
    const double n    = 1.0 - (1.0 / std::sqrt(2.0 * M_PI)) * std::exp(-0.5*x*x) * poly;
    return x >= 0.0 ? n : 1.0 - n;
}

static double bs_call(double S, double K, double r, double sigma, double T) {
    const double sqrtT = sigma * std::sqrt(T);
    const double d1 = (std::log(S / K) + (r + 0.5*sigma*sigma)*T) / sqrtT;
    const double d2 = d1 - sqrtT;
    return S * norm_cdf(d1) - K * std::exp(-r*T) * norm_cdf(d2);
}

int main() {
    // -----------------------------------------------------------------------
    // Parameters: ATM call, T=1, standard inputs for textbook comparison
    // -----------------------------------------------------------------------
    const double S0    = 100.0;  // spot price
    const double K     = 100.0;  // strike
    const double r     = 0.05;   // continuously compounded risk-free rate
    const double sigma = 0.20;   // annualised vol (decimal)
    const double T     = 1.0;    // expiry in years
    const int    N     = 100000; // number of MC paths
    const int    m     = 252;    // time steps (daily monitoring)

    // Time grid: t_0=0, t_1=T/m, ..., t_m=T
    std::vector<double> timeGrid(m + 1);
    for (int i = 0; i <= m; ++i)
        timeGrid[i] = T * i / m;

    // -----------------------------------------------------------------------
    // Analytic reference
    // -----------------------------------------------------------------------
    const double analytic = bs_call(S0, K, r, sigma, T);
    const double discount = std::exp(-r * T);

    // -----------------------------------------------------------------------
    // Monte Carlo: Euler-Maruyama
    // -----------------------------------------------------------------------
    BlackScholesModel  bsModel(r, sigma);
    EulerPathSimulator eulerSim(bsModel, timeGrid, /*seed=*/1234ULL);
    CallPayoff         callPO(K);
    MonteCarlo         mcEuler(S0, discount, eulerSim);

    double se_euler = 0.0;
    const double priceEuler = mcEuler.price_with_stderr(N, callPO, se_euler);

    // -----------------------------------------------------------------------
    // Monte Carlo: Milstein
    // -----------------------------------------------------------------------
    MilsteinPathSimulator milsteinSim(bsModel, timeGrid, /*seed=*/5678ULL);
    MonteCarlo            mcMilstein(S0, discount, milsteinSim);

    double se_milstein = 0.0;
    const double priceMilstein = mcMilstein.price_with_stderr(N, callPO, se_milstein);

    // -----------------------------------------------------------------------
    // Report
    // -----------------------------------------------------------------------
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "=== European Call: S=100, K=100, r=5%, σ=20%, T=1y ===\n\n";
    std::cout << "Black-Scholes analytic:     " << analytic << "\n\n";
    std::cout << "Euler-Maruyama MC:\n";
    std::cout << "  Price:      " << priceEuler  << "\n";
    std::cout << "  Std error:  " << se_euler    << "\n";
    std::cout << "  95% CI:     ["
              << priceEuler - 1.96*se_euler << ", "
              << priceEuler + 1.96*se_euler << "]\n";
    std::cout << "  In CI?      "
              << ((analytic >= priceEuler - 1.96*se_euler &&
                   analytic <= priceEuler + 1.96*se_euler) ? "YES" : "NO")
              << "\n\n";
    std::cout << "Milstein MC:\n";
    std::cout << "  Price:      " << priceMilstein << "\n";
    std::cout << "  Std error:  " << se_milstein   << "\n";
    std::cout << "  95% CI:     ["
              << priceMilstein - 1.96*se_milstein << ", "
              << priceMilstein + 1.96*se_milstein << "]\n";
    std::cout << "  In CI?      "
              << ((analytic >= priceMilstein - 1.96*se_milstein &&
                   analytic <= priceMilstein + 1.96*se_milstein) ? "YES" : "NO")
              << "\n\n";
    std::cout << "Paths: " << N << "   Time steps: " << m << "\n";
    return 0;
}

Validation

Test 1: Euler vs. Black-Scholes Closed Form

Inputs: $S_0 = 100$ , $K = 100$ , $r = 0.05$ , $\sigma = 0.20$ , $T = 1$ year.

Black-Scholes analytic call price:

$d_1 = \frac{\ln(1) + (0.05 + 0.02) \times 1}{0.20 \times 1} = \frac{0.07}{0.20} = 0.3500$ $d_2 = 0.3500 - 0.20 = 0.1500$ $C_{\text{BS}} = 100 \cdot N(0.3500) - 100 e^{-0.05} \cdot N(0.1500) \approx 10.4506$

With $N = 100{,}000$ paths and $m = 252$ steps, the Euler MC price should sit inside the 95% confidence interval $\bar{C} \pm 1.96\,\hat{\text{SE}}$ approximately 95% of the time. The standard error for this problem is approximately $\hat{\text{SE}} \approx 0.03$ – $0.05$ , giving a 95% CI half-width of about $\pm 0.06$ – $\pm 0.10$ .

Test 2: Milstein Pathwise Accuracy

For GBM with $m = 10$ steps (deliberately coarse to expose scheme error), compute the root mean square deviation between simulated final values and the exact solution $S_T = S_0 \exp\!\left[(r - \tfrac{1}{2}\sigma^2)T + \sigma W_T\right]$ over $N = 10{,}000$ paths with the same $Z_i$ draws. Milstein RMSE should be smaller than Euler RMSE by approximately a factor of $\sqrt{m}$ for moderate step sizes — consistent with the one-order improvement in strong convergence rate.

Test 3: Asian Payoff Check

For an Asian call with flat drift ( $r = 0$ , $\sigma = 0$ , deterministic path), the payoff $\max(\bar{S} - K, 0)$ equals $\max(S_0 - K, 0)$ exactly. For $S_0 = 110$ , $K = 100$ , the price should be $10.0$ exactly with zero variance across all paths. This validates the AsianCallPayoff implementation independently of the scheme.

Test 4: Put-Call Parity

Run the MC pricer for both a call and a put with identical parameters. Verify:

$\hat{C}_{\text{MC}} - \hat{P}_{\text{MC}} \approx S_0 - K e^{-rT}$

This is a model-independent check that cannot pass if either the call or put payoff is incorrectly implemented.

Limitations and Pitfalls

Seed Management

A fixed seed (e.g., 42) produces reproducible paths — mandatory for unit tests. However, two instances of PathSimulator constructed with the same seed will produce identical Brownian paths, which is correct for single-threaded use but breaks variance estimation across parallel workers. In production: initialise each thread's generator with a seed derived from std::random_device or a seed sequence that ensures independence. The C++ standard does not guarantee std::random_device is truly random on all platforms — check your target OS.

Euler Instability for Non-Linear Diffusion

For models where $\sigma(t,s)$ grows faster than linearly in $s$ (e.g., CEV with $\beta > 1$ , or rough volatility models), Euler can produce negative values of $S_{t_{i+1}}$ . The full-truncation scheme replaces $s$ with $\max(s, 0)$ in both drift and diffusion evaluations before each step. This prevents non-physical negative prices at the cost of introducing a small bias near zero.

Path Memory

Returning std::vector<double> from simulate_path allocates on the heap per path. For vanilla European options, the payoff depends only on $S_T$ — there is no need to store the full path. Storing and discarding 252 doubles per path for $N = 10^6$ paths involves $252 \times 10^6 \times 8 \approx 2$ GB of allocation and deallocation. For European options, pass a running maximum, running sum, or terminal value instead. For genuinely path-dependent products (Asian, barrier), the full path is required — but stream the payoff computation immediately after simulation so at most one path resides in memory at a time.

Convergence Rate for Discontinuous Payoffs

The weak convergence rate $\mathcal{O}(\Delta t)$ assumes a smooth test function $f$ . For digital options ( $\psi = \mathbf{1}_{S_T > K}$ ) and barrier options, the payoff has a discontinuity. Near the barrier or strike, the weak error has an additional contribution of order $\mathcal{O}(\sqrt{\Delta t})$ from the discrete monitoring. This is not a scheme deficiency — it is a fundamental feature of discrete approximation to continuous monitoring. Brownian bridge interpolation can partially recover the convergence rate for barrier options.

The Clone Pattern vs. `std::unique_ptr`

The raw pointer _model managed by delete _model in the destructor is correct but fragile: if a constructor throws after _model is assigned, the destructor may not run (depending on where the exception occurs), causing a leak. Production code should replace const Model* _model with std::unique_ptr<Model> and initialise it in the constructor initialiser list: _model(model.clone()). The clone() return type would then be std::unique_ptr<Model>. This module uses raw pointers to match the course project's pedagogical intent; understand the RAII alternative for production code.

Interview Angle

Level 1 (Junior Quant Developer)

"What is the Euler-Maruyama scheme? What is its strong convergence order?"

Expected: the SDE $dS_t = \mu\,dt + \sigma\,dW_t$ is discretised as $S_{i+1} = S_i + \mu(t_i, S_i)\Delta t + \sigma(t_i, S_i)\sqrt{\Delta t}\,Z_i$ where $Z_i \sim \mathcal{N}(0,1)$ i.i.d. Strong convergence order is $\mathcal{O}(\sqrt{\Delta t})$ : the expected absolute pathwise error between the discrete approximation and the true solution goes to zero as $\sqrt{\Delta t}$ . Weak order is $\mathcal{O}(\Delta t)$ , governing the accuracy of price estimates for smooth payoffs.

"How do you estimate the standard error of a Monte Carlo price?"

Expected: run $N$ paths, collect discounted payoffs $\psi_1, \ldots, \psi_N$ . Compute sample mean $\bar{\psi} = \frac{1}{N}\sum \psi_j$ and sample variance $\hat{\sigma}^2 = \frac{1}{N-1}\sum(\psi_j - \bar{\psi})^2$ . Standard error is $\hat{\sigma}/\sqrt{N}$ . A 95% confidence interval is $\bar{\psi} \pm 1.96\,\hat{\sigma}/\sqrt{N}$ . For the price to be meaningful, report both the estimate and this interval.

Level 2 (Mid-level / Pricing Quant)

"What does the Milstein correction add over Euler, and why does it improve strong convergence?"

Expected: Euler applies Itô's formula only to first order. Milstein expands $\sigma(t, S_t)$ using Itô's lemma, capturing the contribution of the quadratic variation $d\langle W\rangle_t = dt$ . The added term $\tfrac{1}{2}\sigma \cdot \partial_s\sigma \cdot [(\Delta W)^2 - \Delta t]$ accounts for the curvature of the diffusion coefficient in $s$ . This raises the strong convergence order from $\mathcal{O}(\sqrt{\Delta t})$ to $\mathcal{O}(\Delta t)$ — meaning pathwise errors decrease twice as fast as $\Delta t \to 0$ .

"Why does PathSimulator use model.clone() rather than storing const Model& _model?"

Expected: the Model argument in the constructor may be a temporary, a stack variable, or an object that the caller subsequently modifies or destroys. Storing a reference to it would leave _simulator with a dangling or aliased reference. clone() creates a heap-allocated independent copy that the PathSimulator owns for its entire lifetime, regardless of what happens to the original. This is the prototype design pattern: copying by virtual dispatch without knowing the concrete type.

Level 3 (Senior / Quant Researcher)

"For a barrier option with daily monitoring, why does Euler's convergence deteriorate compared to a vanilla call?"

Expected: weak convergence theory requires the test function $f$ to be smooth (at least twice continuously differentiable). For a vanilla call, $\psi(S_T) = \max(S_T - K, 0)$ is convex but only Lipschitz, and convergence still holds at the expected rate in practice. For a daily-monitored barrier option, the payoff depends on whether any $S_{t_i}$ exceeds the barrier — a discontinuous functional. Near the barrier, the probability of crossing in one discrete step is $\mathcal{O}(\sqrt{\Delta t})$ , contributing a bias term of order $\mathcal{O}(\sqrt{\Delta t})$ on top of the scheme's own $\mathcal{O}(\Delta t)$ weak error. The dominant error is thus $\mathcal{O}(\sqrt{\Delta t})$ , as for strong convergence. Brownian bridge correction introduces a continuous-time crossing probability within each step, partially recovering the $\mathcal{O}(\Delta t)$ rate.