Simulation pricing with variance reduction

Overview

Monte Carlo prices by the law of large numbers. Simulate the terminal spot S_T under the risk-neutral measure many times, take the payoff at each path, discount, and average. The price is unbiased for any finite number of paths; what improves with more paths is the precision, measured by the standard error, which falls as 1/√N.

The implementation is built around three questions a reader should be able to answer from the code: how do you get to the terminal spot, how do you reduce the variance, and how do you know the number you printed is anywhere near the truth?

Getting to the terminal spot

The dynamics under the risk-neutral measure are geometric Brownian motion:

There are three ways to integrate this to maturity, and the implementation gives you all three so you can see what each one costs:

Exact sampling

Geometric Brownian motion is integrable, so the terminal spot can be drawn directly:

One normal draw per path, no time-stepping, no discretisation error. Any bias in an exact-sampling price is pure sampling noise — the kind that the standard error measures. The implementation uses this as the default scheme because it is the one that makes the fewest assumptions.

inline double terminal_spot(McScheme scheme, double s0, double drift, double sig,
                            double T, int steps, const std::vector<double>& z) {
    if (scheme == McScheme::Exact)
        return s0 * std::exp((drift - 0.5*sig*sig)*T + sig*std::sqrt(T)*z[0]);
    // ... Euler / Milstein below
}

Euler–Maruyama

The simplest discretisation of an SDE. Replace the differential with a finite step:

Weak order 1 — the expected value converges at O(Δt). For a European option the bias is small at reasonable step counts, but it is there, and it is the kind of bias that disappears if you forget to check it. The implementation clamps negative spot values to zero, because Euler can step a positive process negative; the comment in the source notes that this clamping is a real source of bias, not a formality.

Milstein

Adds the Stratonovich correction term to the Euler step:

Strong order 1 — the pathwise error converges at O(Δt), better than Euler’s O(√Δt). For a European option the weak-order bias is similar to Euler’s at the same step count, but Milstein is the correct generalisation when you need pathwise accuracy (e.g. for path-dependent contracts). The implementation includes it so the three schemes form a complete picture: exact (no error), Euler (small weak bias), Milstein (strong-order correct).

Variance reduction

The standard error of a plain Monte Carlo estimate is σ_payoff / √N, and σ_payoff for an option payoff can be large relative to the price. Two techniques are used here, both optional and both independent:

Both are applied in the same draw function, so a path with both techniques produces one adjusted payoff that carries the benefit of each. The demo main.cpp reports the factor by which the standard error shrinks when both are turned on — around 7.5× at the same path count.

Control variate

The control variate is the undiscounted terminal spot S_T. Its expectation under the risk-neutral measure is the forward price S₀ · exp((r-q)T), which is known exactly. The terminal spot is highly correlated with the option payoff (especially for a call), so subtracting a β-weighted deviation of S_T - forward from the payoff removes much of the variance without biasing the result:

The β is estimated on a small pilot run and then held fixed for the main simulation. Estimating β on the same paths it corrects would bias the result; the implementation uses a separate pilot of up to min(paths/10, 20000) paths for this, so the correction stays independent of the sample it is applied to.

if (cfg.control_variate) {
    const int64_t pilot = std::min<int64_t>(cfg.paths/10, 20000);
    // ... estimate beta from pilot ...
}
// main loop
sum += payoff - beta * (control - forward);

The forward used in the control is the risk-neutral expectation of the terminal spot:

Standard error and 95% interval

Every result carries a standard error, computed from the sample variance of the adjusted payoffs with Bessel’s correction:

The discount factor is applied after the mean and standard error are computed, so they are stated in present-value terms. The 95% interval is price ± 1.959964 · standard error, and it is reported on every result so a claim of agreement with the analytic price can be checked rather than eyeballed.

struct McResult {
    double price;
    double standard_error;
    int64_t paths;
    [[nodiscard]] double ci_low()  const noexcept { return price - 1.959964 * standard_error; }
    [[nodiscard]] double ci_high() const noexcept { return price + 1.959964 * standard_error; }
};

The standard error falls as 1/√N, which the test suite checks by measuring it at {2000, 8000, 32000, 128000} paths and fitting the observed order — expected 0.5, observed in the range 0.4 to 0.6.

American exercise

Monte Carlo cannot price an American option by forward simulation alone. An American contract requires knowing the optimal stopping rule — the value of exercising now versus continuing — and forward simulation gives you the latter only by a regression against future paths (Longstaff–Schwartz). The implementation rejects American exercise outright rather than returning a number that looks like a price:

if (opt.exercise == Exercise::American)
    throw std::invalid_argument(
        "monte_carlo: American exercise needs a regression method "
        "(Longstaff-Schwartz); forward simulation alone cannot price it");

Source