C++20 + Python● Live

Longstaff-Schwartz LSMC Pricer

American option pricing via Least-Squares Monte Carlo (Longstaff & Schwartz, 2001). C++20 production implementation: GBM simulation with antithetic variates, Laguerre polynomial regression, and backward induction. Python Plotly visualisations: convergence, exercise boundary, regression quality, price surface, and early exercise premium surface.

C++20CMakeCatch2LSMCAmerican optionsLaguerre basisAntithetic variatesPythonNumPyPlotlyscipy

1. Market Context and Assumptions

European options grant the holder the right to exercise only at expiry. American options grant that right at any time up to and including expiry. This optionality has real economic value — and pricing it correctly is non-trivial. Black-Scholes has no closed form for American options; finite differences work in one dimension but scale poorly; Monte Carlo is natural for path-dependent products but requires careful handling of the early-exercise decision.

Longstaff & Schwartz (2001) resolved this by observing that the optimal stopping problem reduces, at each time step, to comparing the immediate exercise value with the continuation value — and that the continuation value can be estimated by OLS regression of simulated payoffs against a small set of basis functions of the current spot. The algorithm is now standard in the industry for Bermudan swaptions, callable bonds, American equity options, and CVA with early termination.

Model Assumptions

Underlying follows GBM: dS = (r − q)S dt + σS dW under the risk-neutral measure Q.
Constant risk-free rate r, continuous dividend yield q, and volatility σ (annualised, decimal).
No transaction costs; frictionless early exercise at any of the M equally-spaced dates.
American price = Bermudan price in the limit M → ∞ (M = 50 is accurate to < 1 bp).
Antithetic variates: paths simulated in pairs (Z, −Z) to halve Monte Carlo variance.

2. American Options as an Optimal Stopping Problem

An American put with payoff g(S) = (K − S)⁺ gives the holder the right to exercise at any stopping time τ ∈ [0, T]. Under Q, the fair price is:

$P_0 = \sup_{\tau \in \mathcal{T}_{[0,T]}} \mathbb{E}^Q\!\left[e^{-r\tau}\,(K - S_\tau)^+\right]$

where $\mathcal{T}_{[0,T]}$ is the set of $\mathbb{F}$ -stopping times in $[0, T]$ .

On a Bermudan grid $t_0 < t_1 < \cdots < t_M = T$ with $\Delta t = T/M$ , the problem reduces to a discrete Bellman recursion. Define the value function $V(t_k, S)$ :

$V(t_M, S) = g(S) = (K - S)^+$

$V(t_k, S) = \max\!\Big\{\underbrace{g(S)}_{\text{exercise now}},\; \underbrace{e^{-r\,\Delta t}\,\mathbb{E}^Q[V(t_{k+1}, S_{t_{k+1}}) \mid S_{t_k} = S]}_{\text{continuation value}\; C(t_k, S)}\Big\}$

The optimal exercise region at time $t_k$ is $\mathcal{E}_k = \{S : g(S) \geq C(t_k, S)\}$ , and the exercise boundary $S^*(t_k)$ is its upper boundary for a put (exercise if $S < S^*(t_k)$ ).

The challenge is computing C(t_k, S). Finite difference methods discretise this directly on a grid — accurate in one dimension but the grid grows as O(d^n) in d dimensions. LSMC estimates C(t_k, S) from the simulation itself, bypassing the grid entirely.

3. The Least-Squares Monte Carlo Algorithm

3.1 Key insight

By the Markov property, C(t_k, S) = E^Q[V(t_{k+1}, S_{k+1}) | S_k = S] is a function of S alone. Longstaff & Schwartz approximate it by a linear combination of basis functions:

$\hat{C}(t_k, S) = \sum_{j=1}^{p} \beta_j^{(k)}\, \psi_j(S)$

The coefficients $\beta^{(k)}$ are estimated by OLS across in-the-money paths at time $t_k$ :

$\hat{\beta}^{(k)} = \underset{\beta}{\arg\min} \sum_{i \in \text{ITM}} \Big(V_i - \sum_{j=1}^{p} \beta_j\, \psi_j(S_{t_k}^{(i)})\Big)^2$

where $V_i$ is the current time-0 present value of future cash flows along path $i$ . The regression uses only in-the-money paths because OTM paths have zero exercise value — including them would dilute the regression toward the wrong region of the continuation surface.

3.2 Laguerre polynomial basis

Longstaff & Schwartz recommend the weighted Laguerre polynomials as basis functions. With x = S/K (normalised spot):

$\psi_k(x) = e^{-x/2} \cdot L_k(x)$

$k$	$L_k(x)$	$\psi_k(x)$
0	$1$	$e^{-x/2}$
1	$1 - x$	$e^{-x/2}(1-x)$
2	$1 - 2x + \frac{x^2}{2}$	$e^{-x/2}\left(1 - 2x + \frac{x^2}{2}\right)$
3	$1 - 3x + \frac{3x^2}{2} - \frac{x^3}{6}$	$e^{-x/2}\left(1 - 3x + \frac{3x^2}{2} - \frac{x^3}{6}\right)$

The exponential weight $e^{-x/2}$ decays the basis away from $x = 0$ (deep ITM), concentrating the regression on the region where exercise decisions matter. Three basis functions $(p = 3)$ is standard; adding a fourth improves accuracy marginally but at increased regression variance.

3.3 Full algorithm

Working with V[j] as the time-0 present value throughout (not raw payoffs) avoids per-step discounting and makes the exercise comparison dimensionally consistent:

LSMC algorithm (time-0 PV formulation)

Simulate $N$ antithetic-paired GBM paths: $S_{t_k}^{(j)}$ for $k = 0, \ldots, M$ and $j = 1, \ldots, N$ .
Initialise $V^{(j)} = e^{-rT} \cdot g(S_T^{(j)})$ for all $j$ .
For $k = M{-}1$ $k = M - 1$ down to $1$ $1$ :
- Let $\mathcal{I}_k = \{j : g(S_{t_k}^{(j)}) > 0\}$ (in-the-money paths).
- Regress $\{V^{(j)}\}_{j \in \mathcal{I}_k}$ on $\{\psi(S_{t_k}^{(j)}/K)\}_{j \in \mathcal{I}_k}$ $\Rightarrow \hat{\beta}^{(k)}$ .
- For each $j \in \mathcal{I}_k$ : if $e^{-rt_k} g(S_{t_k}^{(j)}) > \hat{\beta}^{(k)} \cdot \psi(S_{t_k}^{(j)}/K)$ , set $V^{(j)} = e^{-rt_k} g(S_{t_k}^{(j)})$ .
$\hat{P} = \frac{1}{N} \sum_{j=1}^N V^{(j)}$ ,\quad $\widehat{\text{se}} = \frac{\text{std}(V)}{\sqrt{N}}$ .

Why LSMC is a lower bound

The exercise boundary estimated by regression is sub-optimal: OLS minimises squared error across the simulated paths, but the resulting decision rule may miss exercise at some paths where it would have been optimal. Sub-optimal exercise means we receive a smaller or later cash flow than the true optimal policy would deliver. Therefore LSMC underestimates the true American price — it is a lower bound. The bias diminishes as N → ∞ (more paths improve regression accuracy) and as M → ∞ (finer grid reduces discretisation error). For N = 100,000 and M = 50, the bias is typically under 1 basis point for standard equity puts.

4. Interactive Pricer

The widget below runs LSMC in the browser using antithetic GBM simulation, Laguerre polynomial regression (3 basis functions), and backward induction across 50 time steps. The Pricing tab shows the price decomposition (American = European + EEP). The Exercise Boundary shows the critical spot S*(t) below which exercise is optimal. The Convergence tab demonstrates O(1/√N) convergence.

LSMC American Put — 12,000 paths · 50 steps · 3 basis functions

Spot S₀100

Strike K100

Volatility σ20%

Maturity T (yr)1.00

Rate r5.0%

Dividend q0.0%

Paths N12,000

LSMC American put6.0833

Std error±0.0667

BS European put5.5735

Early exercise premium0.5098

The American put price = BS European put + early exercise premium (EEP). EEP reflects the value of being able to receive K now rather than wait.

LSMC American put6.0833

BS European put5.5735

Early exercise premium0.5098

95% CI lower

5.9526

95% CI upper

6.2141

EEP / American

8.4%

Moneyness S/K

1.000

5. C++20 Implementation

The production C++ implementation lives in compute/pricing-lib/. The path simulation uses exact log-normal increments (no Euler bias), and the regression uses Gauss-Jordan elimination on the 3×3 or 4×4 normal equations — fast and numerically stable for these small systems.

include/bb/LSMC.hppC++20

// LSMC.hpp — American option pricing via Least-Squares Monte Carlo
#pragma once
#include "common.hpp"
#include <cstdint>
#include <vector>

namespace bb {

/**
 * Longstaff-Schwartz LSMC result.
 *
 * price            — American option price (lower bound; bias < 1 bp for N ≥ 50k)
 * std_error        — Monte Carlo standard error σ / √N
 * exercise_boundary— S*(t_k) for k = 0..n_steps; NaN where not defined.
 *                    Put:  exercise if S < S*(t). Call: exercise if S > S*(t).
 */
struct LSMCResult {
    double   price;
    double   std_error;
    double   ci_lower;
    double   ci_upper;
    uint64_t n_paths;
    int      n_steps;
    std::vector<double> exercise_boundary;
};

struct LSMCConvergencePoint { uint64_t n_paths; double price; double std_error; };

/**
 * American option via LSMC.
 *
 * Assumptions: GBM dS = (r-q)S dt + σS dW.
 * Basis: weighted Laguerre ψ_k(S/K) = exp(-S/(2K)) · L_k(S/K), k = 0..basis_deg-1.
 * Variance reduction: antithetic variates (n_paths must be even; enforced internally).
 *
 * Throws std::invalid_argument on invalid inputs.
 */
LSMCResult lsmc_american(
    double S, double K, double r, double q,
    double sigma, double T, OptionType type,
    uint64_t n_paths   = 100'000,
    int      n_steps   = 50,
    int      basis_deg = 3,
    uint64_t seed      = 42);

std::vector<LSMCConvergencePoint> lsmc_convergence(
    double S, double K, double r, double q,
    double sigma, double T, OptionType type,
    uint64_t n_max = 200'000, int n_conv_steps = 10, uint64_t seed = 42);

double lsmc_european_check(
    double S, double K, double r, double q,
    double sigma, double T, OptionType type,
    uint64_t n_paths = 100'000, int n_steps = 50, uint64_t seed = 42);

} // namespace bb

src/LSMC.cppC++20

// LSMC.cpp
#include "bb/LSMC.hpp"
#include <algorithm>
#include <cmath>
#include <limits>
#include <random>
#include <stdexcept>
#include <vector>

namespace bb {
namespace {

// Weighted Laguerre basis ψ_k(x) = exp(-x/2) · L_k(x),  x = spot / strike.
// L_0 = 1,  L_1 = 1-x,  L_2 = 1-2x+x²/2,  L_3 = 1-3x+3x²/2-x³/6
inline void laguerre_basis(double x, int deg, double* out) noexcept {
    const double e = std::exp(-0.5 * x);
    if (deg >= 1) out[0] = e;
    if (deg >= 2) out[1] = e * (1.0 - x);
    if (deg >= 3) out[2] = e * (1.0 - 2.0*x + 0.5*x*x);
    if (deg >= 4) out[3] = e * (1.0 - 3.0*x + 1.5*x*x - x*x*x/6.0);
}

// OLS via normal equations: (A^T A) β = A^T y, Gauss-Jordan with partial pivot.
// Returns zero vector if near-singular (< 1e-14 pivot).
std::vector<double> ols_solve(const double* A, const double* y,
                               std::size_t n_itm, int p);

inline double payoff_at(double spot, double K, OptionType type) noexcept {
    return (type == OptionType::Call) ? std::max(spot-K, 0.0) : std::max(K-spot, 0.0);
}

} // namespace

LSMCResult lsmc_american(
    double S, double K, double r, double q,
    double sigma, double T, OptionType type,
    uint64_t n_paths, int n_steps, int basis_deg, uint64_t seed)
{
    if (S<=0 || K<=0 || sigma<=0 || T<=0)
        throw std::invalid_argument("S, K, sigma, T must be strictly positive");
    if (n_paths < 4)   throw std::invalid_argument("n_paths >= 4 required");
    if (n_steps < 2)   throw std::invalid_argument("n_steps >= 2 required");
    if (basis_deg<1 || basis_deg>4)
        throw std::invalid_argument("basis_deg in [1,4]");

    if (n_paths % 2 != 0) ++n_paths;
    const uint64_t half = n_paths / 2;

    const double dt     = T / n_steps;
    const double drift  = (r - q - 0.5*sigma*sigma) * dt;
    const double vol_dt = sigma * std::sqrt(dt);

    // ── 1. Simulate N antithetic-paired GBM paths ─────────────────────────
    // paths[step * n_paths + j] = S at time step*dt on path j.
    std::vector<double> paths((n_steps+1) * n_paths);
    std::mt19937_64 rng(seed);
    std::normal_distribution<double> nd;

    for (uint64_t j = 0; j < n_paths; ++j) paths[j] = S;
    for (int step = 1; step <= n_steps; ++step) {
        const auto prev = (step-1) * n_paths, curr = step * n_paths;
        for (uint64_t j = 0; j < half; ++j) {
            double z = nd(rng);
            paths[curr+j]        = paths[prev+j]        * std::exp(drift + vol_dt*z);
            paths[curr+half+j]   = paths[prev+half+j]   * std::exp(drift - vol_dt*z);
        }
    }

    // ── 2. Initialise V[j] = e^{-rT} · payoff(S_T[j]) ───────────────────
    // V[j] is the time-0 present value; we update it as we find early exercises.
    const double disc_T = std::exp(-r*T);
    std::vector<double> V(n_paths);
    for (uint64_t j = 0; j < n_paths; ++j)
        V[j] = disc_T * payoff_at(paths[n_steps*n_paths + j], K, type);

    std::vector<double> exercise_boundary(n_steps+1,
                                          std::numeric_limits<double>::quiet_NaN());

    std::vector<uint64_t> itm_idx;
    std::vector<double>   reg_A, reg_y;
    double basis_buf[4] = {};

    // ── 3. Backward induction ─────────────────────────────────────────────
    for (int step = n_steps-1; step >= 1; --step) {
        const double t      = step * dt;
        const double disc_t = std::exp(-r*t);
        const auto   base   = step * n_paths;

        // Collect in-the-money paths.
        itm_idx.clear();
        for (uint64_t j = 0; j < n_paths; ++j)
            if (payoff_at(paths[base+j], K, type) > 0.0) itm_idx.push_back(j);

        const auto n_itm = itm_idx.size();
        if (n_itm < (std::size_t)(basis_deg+1)) continue;

        // Regression: Y = V[j] (time-0 PV), X = ψ(S_t/K).
        reg_A.resize(n_itm * basis_deg);
        reg_y.resize(n_itm);
        for (std::size_t ii = 0; ii < n_itm; ++ii) {
            uint64_t j = itm_idx[ii];
            laguerre_basis(paths[base+j]/K, basis_deg, basis_buf);
            for (int k = 0; k < basis_deg; ++k)
                reg_A[ii*basis_deg + k] = basis_buf[k];
            reg_y[ii] = V[j];
        }
        const auto beta = ols_solve(reg_A.data(), reg_y.data(), n_itm, basis_deg);

        // Exercise if e^{-rt}·payoff(S_t) > fitted continuation E_t[V|S_t].
        double bnd = std::numeric_limits<double>::quiet_NaN();
        bool any = false;
        for (std::size_t ii = 0; ii < n_itm; ++ii) {
            uint64_t j  = itm_idx[ii];
            double spot = paths[base+j];
            laguerre_basis(spot/K, basis_deg, basis_buf);
            double cont = 0;
            for (int k = 0; k < basis_deg; ++k) cont += beta[k] * basis_buf[k];
            const double imm_pv = disc_t * payoff_at(spot, K, type);
            if (imm_pv > cont && imm_pv > 0.0) {
                V[j] = imm_pv;
                any  = true;
                // Put boundary = max exercised spot (exercise if S < S*).
                if (type==OptionType::Put  && (std::isnan(bnd)||spot>bnd)) bnd=spot;
                if (type==OptionType::Call && (std::isnan(bnd)||spot<bnd)) bnd=spot;
            }
        }
        if (any) exercise_boundary[step] = bnd;
    }

    // ── 4. Price ──────────────────────────────────────────────────────────
    double sum=0, sum2=0;
    for (uint64_t j = 0; j < n_paths; ++j) { sum+=V[j]; sum2+=V[j]*V[j]; }
    double price   = sum / n_paths;
    double var     = (sum2/n_paths - price*price) / (n_paths-1);
    double std_err = std::sqrt(std::max(var, 0.0));

    return { price, std_err,
             price - 1.96*std_err, price + 1.96*std_err,
             n_paths, n_steps, std::move(exercise_boundary) };
}

} // namespace bb

tests/test_lsmc.cppC++20 / Catch2 v3

// test_lsmc.cpp — Catch2 v3 unit tests
#include "bb/LSMC.hpp"
#include "bb/BlackScholes.hpp"
#include "bb/FiniteDifference.hpp"
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>
#include <cmath>

using namespace bb;

static constexpr double S0=100, K0=100, r0=0.05, q0=0, sig0=0.2, T1=1.0;
static constexpr uint64_t N=100'000;

// American put >= European put (early exercise premium >= 0).
TEST_CASE("American put >= BS European put", "[lsmc][eep]") {
    auto am = lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Put, N,50,3,42);
    auto eu = bs_european(S0,K0,r0,q0,sig0,T1,OptionType::Put);
    REQUIRE(am.price >= eu.price - 2*am.std_error);
}

// Agreement with FD benchmark within 1.5 % (lower-bound bias of LSMC).
TEST_CASE("ATM put agrees with FD PSOR within 1.5 %", "[lsmc][benchmark]") {
    auto fd  = fd_american_psor(S0,K0,r0,q0,sig0,T1,OptionType::Put,600,600);
    auto lsm = lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Put,N,50,3,42);
    REQUIRE(lsm.price == Catch::Approx(fd.price).margin(0.015*fd.price));
}

// American call == European call when q = 0 (never optimal to exercise early).
TEST_CASE("American call == BS European call when q=0", "[lsmc][call-parity]") {
    auto am = lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Call,N,50,3,42);
    auto eu = bs_european(S0,K0,r0,q0,sig0,T1,OptionType::Call);
    REQUIRE(am.price == Catch::Approx(eu.price).margin(0.015*eu.price));
}

// Std error shrinks with more paths (O(1/√N) behaviour).
TEST_CASE("Std error decreases with more paths", "[lsmc][convergence]") {
    auto r1 = lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Put, 4'000,50,3,42);
    auto r2 = lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Put,40'000,50,3,42);
    REQUIRE(r2.std_error < r1.std_error);
    REQUIRE(r1.std_error / r2.std_error > 2.0);
}

// Results are exactly reproducible with the same seed.
TEST_CASE("Reproducibility with same seed", "[lsmc][reproducibility]") {
    auto r1 = lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Put,10'000,50,3,99);
    auto r2 = lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Put,10'000,50,3,99);
    REQUIRE(r1.price == r2.price);
}

// Input validation.
TEST_CASE("Throws on invalid inputs", "[lsmc][validation]") {
    REQUIRE_THROWS_AS(lsmc_american(-1,K0,r0,q0,sig0,T1,OptionType::Put),
                      std::invalid_argument);
    REQUIRE_THROWS_AS(lsmc_american(S0,K0,r0,q0,sig0,T1,OptionType::Put,4,1),
                      std::invalid_argument);
}

6. Python Visualisation (Plotly)

The Python implementation uses NumPy for vectorised path simulation and np.linalg.lstsq for robust OLS (SVD fallback, handles near-singular normal matrices). Five Plotly charts cover convergence analysis, the exercise boundary across volatility regimes, regression quality at mid-life, and the American vs European price and EEP surfaces.

lsmc_pricer.pyPython 3.10+

# lsmc_pricer.py — Longstaff-Schwartz LSMC for American options
# Assumptions: GBM dS = (r-q)S dt + σS dW (same as Black-Scholes)
# Conventions: r, q in decimal annualised; σ annualised; T in years.

from __future__ import annotations
from dataclasses import dataclass
import numpy as np
from scipy.stats import norm


@dataclass(frozen=True)
class LSMCResult:
    price:             float
    std_error:         float
    ci_lower:          float
    ci_upper:          float
    exercise_boundary: np.ndarray  # shape (n_steps+1,); NaN where undefined


def bs_european_put(S: float, K: float, r: float, q: float,
                    sigma: float, T: float) -> float:
    """Black-Scholes European put price. Used as lower-bound reference."""
    if T <= 0:
        return max(K - S, 0.0)
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    return K*np.exp(-r*T)*norm.cdf(-d2) - S*np.exp(-q*T)*norm.cdf(-d1)


def laguerre_basis(x: np.ndarray, deg: int = 3) -> np.ndarray:
    """
    Weighted Laguerre polynomial basis for LSMC regression.

    ψ_k(x) = exp(-x/2) · L_k(x),  x = S_t / K

    Longstaff & Schwartz (2001) recommend this basis because the exponential
    weight confines the polynomials to the region where in-the-money paths
    cluster, improving regression stability.

    L_0(x) = 1
    L_1(x) = 1 - x
    L_2(x) = 1 - 2x + x²/2
    """
    e   = np.exp(-0.5 * x)
    B   = [e]
    if deg >= 2: B.append(e * (1 - x))
    if deg >= 3: B.append(e * (1 - 2*x + 0.5*x**2))
    if deg >= 4: B.append(e * (1 - 3*x + 1.5*x**2 - x**3/6))
    return np.column_stack(B)   # shape (n_itm, deg)


def lsmc_american_put(
    S0:        float,
    K:         float,
    r:         float,
    q:         float,
    sigma:     float,
    T:         float,
    n_paths:   int   = 100_000,
    n_steps:   int   = 50,
    basis_deg: int   = 3,
    seed:      int   = 42,
) -> LSMCResult:
    """
    Longstaff-Schwartz LSMC pricer for American put options.

    Algorithm (Longstaff & Schwartz, Rev. Fin. Stud. 14(1), 2001):
      1. Simulate N GBM paths over M time steps with antithetic variates.
      2. Initialise V[j] = e^{-rT} · max(K - S_T[j], 0)  for all paths j.
      3. Backward induction for t = (M-1)·Δt down to Δt:
         a. Collect ITM paths: {j : K - S_t[j] > 0}.
         b. Regress V[j] against ψ(S_t[j]/K) via OLS → β̂.
         c. Ê_t[V | S_t[j]] = ψ(S_t[j]/K) · β̂   (fitted continuation).
         d. If e^{-rt}·(K - S_t[j]) > Ê_t[V | S_t[j]]:
                V[j] = e^{-rt}·(K - S_t[j])        (exercise at t).
      4. Price = (1/N) Σ_j V[j].

    Returns the LSMC price (a lower bound), std error, 95% CI, and the
    exercise boundary S*(t) at each time step.
    """
    rng = np.random.default_rng(seed)

    # ── Simulate antithetic GBM paths ──────────────────────────────────────
    # Shape: (n_steps+1, n_paths)
    assert n_paths % 2 == 0, "n_paths must be even for antithetic pairing"
    half = n_paths // 2
    dt     = T / n_steps
    drift  = (r - q - 0.5*sigma**2) * dt
    vol_dt = sigma * np.sqrt(dt)

    Z    = rng.standard_normal((n_steps, half))          # primal normals
    logS = np.zeros((n_steps + 1, n_paths))
    logS[0] = np.log(S0)

    for step in range(n_steps):
        z                   = Z[step]
        logS[step+1, :half]  = logS[step, :half]  + drift + vol_dt * z
        logS[step+1, half:]  = logS[step, half:]  + drift - vol_dt * z

    S_paths = np.exp(logS)   # (n_steps+1, n_paths)

    # ── Initialise V as time-0 PV of terminal payoff ───────────────────────
    disc_T = np.exp(-r * T)
    V      = disc_T * np.maximum(K - S_paths[-1], 0.0)   # (n_paths,)

    exercise_boundary = np.full(n_steps + 1, np.nan)

    # ── Backward induction ─────────────────────────────────────────────────
    for step in range(n_steps - 1, 0, -1):
        t       = step * dt
        disc_t  = np.exp(-r * t)
        S_t     = S_paths[step]                     # (n_paths,)
        payoff  = np.maximum(K - S_t, 0.0)

        itm     = payoff > 0                        # in-the-money mask
        n_itm   = itm.sum()
        if n_itm < basis_deg + 1:
            continue

        S_itm   = S_t[itm]
        Y       = V[itm]                            # continuation values
        X       = laguerre_basis(S_itm / K, basis_deg)  # (n_itm, deg)

        # OLS via NumPy least-squares (stable SVD fallback internally).
        beta, *_ = np.linalg.lstsq(X, Y, rcond=None)
        cont_hat = X @ beta                         # fitted continuation

        imm_pv   = disc_t * payoff[itm]
        exercise = imm_pv > cont_hat

        # Update V for exercising ITM paths.
        idx_itm            = np.where(itm)[0]
        V[idx_itm[exercise]] = imm_pv[exercise]

        # Exercise boundary: highest spot price at which we exercise (put).
        exercised_spots = S_itm[exercise]
        if exercised_spots.size > 0:
            exercise_boundary[step] = exercised_spots.max()

    # ── Price ──────────────────────────────────────────────────────────────
    price    = V.mean()
    std_err  = V.std(ddof=1) / np.sqrt(n_paths)
    ci_lower = price - 1.96 * std_err
    ci_upper = price + 1.96 * std_err

    return LSMCResult(price, std_err, ci_lower, ci_upper, exercise_boundary)

convergence_analysis.pyPython 3.10+

# convergence_analysis.py — LSMC convergence vs FD benchmark
from __future__ import annotations
import numpy as np
from lsmc_pricer import lsmc_american_put, bs_european_put

try:
    from bb_pricing import fd_american_psor  # C++ pricing lib via pybind11
    HAS_CPP = True
except ImportError:
    HAS_CPP = False


def convergence_study(
    S0: float = 100.0,
    K:  float = 100.0,
    r:  float = 0.05,
    q:  float = 0.00,
    sigma: float = 0.20,
    T:  float = 1.0,
    n_max:  int = 200_000,
    n_points: int = 14,
    seed: int = 42,
) -> dict:
    """
    LSMC price and ±2σ confidence band at log-spaced path counts.
    Benchmarks against FD PSOR and BS European lower bound.

    Returns dict with keys: n_paths, prices, std_errors, fd_price, bs_eu_price.

    Convergence theory:
        E|P_MC - P*| = O(1/√N)
        With antithetic variates: variance halved at no extra computation cost.
        Bias from sub-optimal boundary: ~ O(1/M) as n_steps → ∞.
    """
    n_list = np.unique(
        np.geomspace(500, n_max, n_points).astype(int) // 2 * 2
    ).tolist()

    prices, stderrs = [], []
    for n in n_list:
        res = lsmc_american_put(S0, K, r, q, sigma, T,
                                 n_paths=int(n), n_steps=50, seed=seed)
        prices.append(res.price)
        stderrs.append(res.std_error)

    fd_price  = None
    if HAS_CPP:
        fd_res   = fd_american_psor(S0, K, r, q, sigma, T, "put", 600, 600)
        fd_price = fd_res.price

    bs_eu = bs_european_put(S0, K, r, q, sigma, T)

    return {
        "n_paths":    np.array(n_list),
        "prices":     np.array(prices),
        "std_errors": np.array(stderrs),
        "fd_price":   fd_price,
        "bs_eu_price": bs_eu,
    }

visualization.pyPython 3.10+

# visualization.py — five Plotly charts for the LSMC pricer project
from __future__ import annotations
import numpy as np
import plotly.graph_objects as go
from lsmc_pricer import lsmc_american_put, bs_european_put
from convergence_analysis import convergence_study

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

_LAYOUT = dict(
    paper_bgcolor="#0f172a", plot_bgcolor="#0f172a",
    font=dict(color="#cbd5e1", size=11, family="JetBrains Mono, monospace"),
    margin=dict(l=55, r=20, t=50, b=50),
)


# ── Chart 1: Convergence (LSMC price vs N with ±2σ band and FD reference) ──

def plot_convergence(S0=100, K=100, r=0.05, q=0, sigma=0.2, T=1.0) -> go.Figure:
    """
    Demonstrates O(1/√N) convergence of LSMC.
    The ±2σ confidence band narrows symmetrically on the log-N scale.
    """
    study = convergence_study(S0, K, r, q, sigma, T)
    n     = study["n_paths"]
    p, se = study["prices"], study["std_errors"]

    fig = go.Figure()

    # ±2σ band
    fig.add_trace(go.Scatter(
        x=np.concatenate([n, n[::-1]]),
        y=np.concatenate([p + 2*se, (p - 2*se)[::-1]]),
        fill="toself", fillcolor="rgba(99,102,241,0.12)",
        line=dict(color="rgba(0,0,0,0)"), name="±2σ confidence band",
    ))
    fig.add_trace(go.Scatter(
        x=n, y=p, mode="lines+markers",
        line=dict(color=INDIGO, width=2), marker=dict(size=4),
        name="LSMC American put",
    ))

    if study["fd_price"] is not None:
        fig.add_hline(y=study["fd_price"],
                      line=dict(color=EMERALD, dash="dash", width=1.5),
                      annotation_text=f"FD PSOR = {study['fd_price']:.4f}",
                      annotation_font_color=EMERALD)

    fig.add_hline(y=study["bs_eu_price"],
                  line=dict(color=SLATE, dash="dot", width=1.2),
                  annotation_text=f"BS European = {study['bs_eu_price']:.4f}",
                  annotation_font_color=SLATE)

    fig.update_layout(**_LAYOUT,
        title="LSMC Convergence — American Put Price vs Path Count",
        xaxis_title="Paths N (log scale)", xaxis_type="log",
        yaxis_title="Put Price", height=400)
    return fig


# ── Chart 2: Exercise boundary for different volatility levels ──────────────

def plot_exercise_boundary(S0=100, K=100, r=0.05, q=0, T=1.0) -> go.Figure:
    """
    The critical spot S*(t) below which immediate exercise is optimal.
    Higher volatility → lower boundary (option holder waits longer).
    At expiry (t=T) the boundary converges to K for a put.
    """
    sigmas = [0.10, 0.20, 0.30, 0.40]
    colors = [EMERALD, INDIGO, AMBER, ROSE]
    fig    = go.Figure()

    for sigma, color in zip(sigmas, colors):
        res    = lsmc_american_put(S0, K, r, q, sigma, T,
                                    n_paths=200_000, n_steps=100, seed=42)
        bnd    = res.exercise_boundary
        n_steps = len(bnd) - 1
        t_vals  = np.linspace(0, T, n_steps + 1)
        valid   = ~np.isnan(bnd)

        fig.add_trace(go.Scatter(
            x=t_vals[valid], y=bnd[valid], mode="lines",
            line=dict(color=color, width=2),
            name=f"σ = {int(sigma*100)}%",
        ))

    # Strike reference
    fig.add_hline(y=K, line=dict(color=SLATE, dash="dash", width=1),
                  annotation_text="K", annotation_font_color=SLATE)

    fig.update_layout(**_LAYOUT,
        title="Early Exercise Boundary S*(t) — American Put",
        xaxis_title="Time t (years)", yaxis_title="Critical Spot S*(t)",
        height=420)
    return fig


# ── Chart 3: Regression quality at mid-life ──────────────────────────────────

def plot_regression_quality(S0=100, K=100, r=0.05, q=0, sigma=0.2, T=1.0,
                             n_paths=50_000, target_step_frac=0.5) -> go.Figure:
    """
    Scatter of actual discounted continuation values V[j] vs the OLS fit
    at t = 0.5T for in-the-money paths.

    Shows how well the Laguerre polynomial regression captures the true
    continuation surface — and why the sub-optimality is small.
    """
    import numpy as np
    from lsmc_pricer import laguerre_basis

    n_steps      = 50
    target_step  = max(1, int(target_step_frac * n_steps))
    rng          = np.random.default_rng(42)
    half         = n_paths // 2
    dt           = T / n_steps
    drift        = (r - q - 0.5*sigma**2) * dt
    vol_dt       = sigma * np.sqrt(dt)

    Z    = rng.standard_normal((n_steps, half))
    logS = np.zeros((n_steps + 1, n_paths))
    logS[0] = np.log(S0)
    for step in range(n_steps):
        z = Z[step]
        logS[step+1, :half] = logS[step, :half] + drift + vol_dt*z
        logS[step+1, half:] = logS[step, half:] + drift - vol_dt*z
    S_paths = np.exp(logS)

    disc_T = np.exp(-r * T)
    V      = disc_T * np.maximum(K - S_paths[-1], 0.0)

    # Backward induction up to target_step (collect regression data there).
    for step in range(n_steps - 1, target_step - 1, -1):
        t, disc_t = step*dt, np.exp(-r*step*dt)
        S_t = S_paths[step]; payoff = np.maximum(K - S_t, 0.0)
        itm = payoff > 0
        if itm.sum() < 4: continue
        S_itm = S_t[itm]; Y = V[itm]
        X = laguerre_basis(S_itm / K, 3)
        beta, *_ = np.linalg.lstsq(X, Y, rcond=None)
        cont_hat = X @ beta
        imm_pv   = disc_t * payoff[itm]
        exercise = imm_pv > cont_hat
        idx_itm  = np.where(itm)[0]
        V[idx_itm[exercise]] = imm_pv[exercise]

        if step == target_step:
            # Capture scatter data at this step.
            actual  = Y
            fitted  = cont_hat
            spots   = S_itm
            imm     = imm_pv
            break

    fig = go.Figure()
    fig.add_trace(go.Scatter(
        x=spots, y=actual, mode="markers",
        marker=dict(color=SLATE, size=3, opacity=0.4),
        name="Actual V[j] (time-0 PV)",
    ))
    fig.add_trace(go.Scatter(
        x=spots, y=fitted, mode="markers",
        marker=dict(color=INDIGO, size=4, opacity=0.7),
        name="Fitted continuation Ê[V|S_t]",
    ))
    fig.add_trace(go.Scatter(
        x=spots, y=imm, mode="lines",
        line=dict(color=EMERALD, width=2),
        name=f"Immediate payoff e^{{-rt}}·(K−S)",
    ))
    fig.update_layout(**_LAYOUT,
        title=f"Regression Quality at t = {target_step_frac*T:.2f} yr",
        xaxis_title="Spot S_t", yaxis_title="Time-0 PV",
        height=400)
    return fig


# ── Chart 4: American vs European price surface ──────────────────────────────

def plot_price_surface(r=0.05, q=0.0, sigma=0.20,
                       n_K=12, n_T=10) -> go.Figure:
    """
    American put price (LSMC) and European put price (BS closed-form) as
    functions of strike K and maturity T. The gap between surfaces is the EEP.
    """
    K_grid = np.linspace(80, 120, n_K)
    T_grid = np.linspace(0.1, 2.0, n_T)
    am     = np.zeros((n_T, n_K))
    eu     = np.zeros((n_T, n_K))
    S0     = 100.0

    for j, K in enumerate(K_grid):
        for i, T in enumerate(T_grid):
            am[i, j] = lsmc_american_put(
                S0, K, r, q, sigma, T, n_paths=20_000, n_steps=50, seed=42
            ).price
            eu[i, j] = bs_european_put(S0, K, r, q, sigma, T)

    fig = go.Figure()
    fig.add_trace(go.Surface(x=K_grid, y=T_grid, z=am,
                             colorscale="Viridis", opacity=0.9,
                             name="American put"))
    fig.add_trace(go.Surface(x=K_grid, y=T_grid, z=eu,
                             colorscale="Plasma", opacity=0.55,
                             name="European put"))
    fig.update_layout(**_LAYOUT,
        title="American (top) vs European (bottom) Put Price Surface",
        scene=dict(xaxis_title="Strike K", yaxis_title="Maturity T (yr)",
                   zaxis_title="Put Price"),
        height=520)
    return fig


# ── Chart 5: Early exercise premium (EEP) surface ───────────────────────────

def plot_eep_surface(r=0.05, q=0.0, sigma=0.20,
                     n_K=12, n_T=10) -> go.Figure:
    """
    EEP = American put − European put. Highest for:
      - Deep ITM (K >> S): intrinsic value dominates, carry cost of K·(1-e^{-rT}) is large.
      - Long maturity: time value of early receipt of K accumulates.
      - High interest rates: present value of K received early is maximised.
      - Low dividend yield: no incentive to hold stock for dividends offsets carry.
    """
    K_grid = np.linspace(80, 120, n_K)
    T_grid = np.linspace(0.1, 2.0, n_T)
    eep    = np.zeros((n_T, n_K))
    S0     = 100.0

    for j, K in enumerate(K_grid):
        for i, T in enumerate(T_grid):
            am_price = lsmc_american_put(
                S0, K, r, q, sigma, T, n_paths=20_000, n_steps=50, seed=42
            ).price
            eu_price = bs_european_put(S0, K, r, q, sigma, T)
            eep[i, j] = max(am_price - eu_price, 0.0)

    fig = go.Figure(data=go.Surface(
        x=K_grid, y=T_grid, z=eep, colorscale="Plasma",
    ))
    fig.update_layout(**_LAYOUT,
        title="Early Exercise Premium (EEP) — American Put minus European Put",
        scene=dict(xaxis_title="Strike K", yaxis_title="Maturity T (yr)",
                   zaxis_title="EEP"),
        height=480)
    return fig


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

if __name__ == "__main__":
    S0, K, r, q, sigma, T = 100.0, 100.0, 0.05, 0.0, 0.20, 1.0

    res = lsmc_american_put(S0, K, r, q, sigma, T, n_paths=200_000, n_steps=50)
    eu  = bs_european_put(S0, K, r, q, sigma, T)
    print(f"LSMC American put : {res.price:.4f}  (±{res.std_error:.4f})")
    print(f"BS   European put : {eu:.4f}")
    print(f"Early exercise premium: {max(res.price - eu, 0):.4f}")

    plot_convergence().show()
    plot_exercise_boundary().show()
    plot_regression_quality().show()
    plot_price_surface().show()
    plot_eep_surface().show()

Python charts produced by visualization.py

Chart 1 — Convergence:LSMC price (± 2σ) vs path count N on a log scale. FD PSOR reference shows the true American price; BS European lower bound also overlaid. Confirms O(1/√N) convergence.

Chart 2 — Exercise Boundary:Critical spot S*(t) vs time for σ ∈ {10%, 20%, 30%, 40%}. Higher vol → lower boundary (holder waits longer). Boundary converges to K at expiry.

Chart 3 — Regression Quality:Scatter at t = 0.5T: actual continuation values V[j] (grey), fitted OLS surface (indigo), immediate payoff e^{-rt}·(K−S) (emerald). Shows where exercise and hold regions split.

Chart 4 — Price Surface:American put (LSMC) vs European put (BS) across K × T. The gap visible between the two surfaces is the EEP — widest for deep ITM and long maturity.

Chart 5 — EEP Surface:Early exercise premium across K × T. Confirms EEP is highest for deep ITM, long maturity, high r, low q.

7. Limitations and Failure Modes

Known limitations

Downward bias: LSMC produces a lower bound. For deep ITM short-maturity options the sub-optimal boundary contributes > 1 bp bias even at N = 100k; use Andersen-Broadie dual method to bracket the true price.
Regression instability near expiry: close to T, few paths are ITM, the regression matrix is near-singular, and the boundary estimate is noisy. Increasing n_steps (finer grid) mitigates this.
GBM dynamics: same limitations as Black-Scholes (constant vol, no jumps, no smile). Real equity desks use local vol or stochastic vol as the simulation model.
One-dimensional state: the Laguerre basis in S alone is accurate for standard puts. For path-dependent or multi-factor products, additional state variables (realised variance, running maximum) must be included, and the regression dimension grows.
O(N · M · p²) complexity: with N = 100k, M = 50, p = 3, this is manageable. Bermudan swaptions with M = 20 and a 10-factor model require careful basis construction to keep p small.

8. Interview Angle

American option pricing and LSMC appear at every seniority level in quant interviews. Junior candidates are expected to understand the optimal stopping problem and why early exercise has value. Senior candidates must derive the backward induction and discuss bias, the dual bound, and extension to multi-factor models. Researchers are expected to engage with convergence theory and high-dimensional generalisations.

L1 — Junior

›Why can't you price an American option with the Black-Scholes closed-form formula? What property distinguishes it from a European option?
›Explain the Bellman equation for American option pricing in plain terms. What is the continuation value and when is it optimal to exercise?
›What are Laguerre polynomials and why does Longstaff-Schwartz choose them as basis functions rather than raw powers of S?
›Why is the LSMC price a lower bound for the true American option price? What causes the bias?
›For an American call with no dividends and r > 0, is early exercise ever optimal? Prove it intuitively.

L2 — Senior

›The LSMC regression uses only in-the-money paths. Why? What happens to the regression estimate if you include OTM paths, and how does it affect the boundary?
›Describe the dual method (Andersen & Broadie 2004) for computing an upper bound on the American option price. How does it complement LSMC?
›Compare LSMC and the Crank-Nicolson PSOR finite difference scheme for American puts: computational complexity, accuracy, and scalability to higher dimensions.
›You have a Bermudan swaption with 10 exercise dates. How does LSMC generalise, and what state variables would you include in the regression? What is the curse of dimensionality here?
›How would you compute delta and vega for the LSMC pricer? Discuss finite difference bump-and-reval vs pathwise differentiation (likelihood ratio estimator).

L3 — Researcher

›The LSMC price is a biased estimator of the true American price due to sub-optimal exercise. Characterise the bias order in n_paths and n_steps separately. Which dominates in practice, and how do you balance them for a given accuracy budget?
›Clément, Lamberton & Protter (2002) prove convergence of the LSMC estimator. What are the key regularity conditions on the payoff and the basis functions? Where do standard American puts violate these conditions, and what are the practical consequences?
›For a CVA calculation with hundreds of risk factors and early exercise (e.g., a callable bond portfolio), plain LSMC is intractable. Describe at least two methods (e.g., regression on approximate state spaces, neural network regression, deep hedging) and their trade-offs.
›The Andersen-Broadie dual bound requires solving a sub-problem for each path and each exercise date. What is the computational complexity compared to LSMC, and how does the 'low-biased' / 'high-biased' sandwich bound methodology work in practice?

← Heston MC Engine SVI Calibrator →