Debugging Monte Carlo Option Prices Against Black–Scholes

Benchmark a European call using exact risk-neutral terminal prices, then compare the simulation gap with its sampling standard error. With spot and strike both

100

, one year remaining, volatility

20\%

, and a continuously compounded rate of

4\%

, Black–Scholes gives

9.9251

. If the simulation instead grows the stock with a physical drift of

10\%

, the discounted payoff can average around

14

. That gap reflects a different distribution.

Test the simulated stock before testing the option

For $Z\sim N(0,1)$ , simulate

S_T=S_0\exp\left[(r-\sigma^2/2)T+\sigma\sqrt{T}Z\right].

The log drift includes $-\sigma^2/2$ because the exponential of a mean-zero normal has mean above one. The discounted payoff is $Y=e^{-rT}(S_T-K)^+$ . With independent draws, estimate the price by $\bar Y$ and its standard error by $s_Y/\sqrt n$ . Compare the analytic price with $\bar Y\pm1.96s_Y/\sqrt n$ as an approximate diagnostic interval. Even a correct implementation can miss that interval in about one run out of twenty under the usual normal approximation.

Interpret the price difference in standard errors

Use 100,000 independent normals with seed 96 and identical draws for the two drift settings. The risk-neutral simulation yields about $9.8375$ with standard error $0.0454$ . The difference from $9.9251$ is about $1.93$ standard errors. The physical-drift version yields $13.9797$ with standard error $0.0539$ , more than $75$ standard errors away from Black–Scholes. Increasing the number of simulations reduces random error around each simulation's own population mean, so it helps distinguish the two cases.

With correctly scaled volatility, the standard deviation of simulated log returns should also be close to

20\%

over one year. If you swap the exact formula for Euler time steps, run a time-step convergence check at roughly the same sampling error.

Use the same normal draws when comparing the two drift settings. Then changes in the estimated value can be traced to the model input rather than to a different random sample. Check the simulated discounted mean stock against

S_0

first. If it is materially off by more than its sampling uncertainty, inspect the risk-neutral drift, time units and discounting. If the stock check passes but the option check fails, next test the payoff on hand-picked stock outcomes below, at and above strike. These tests localise the source of the discrepancy before you increase the simulation size.

A reproducible paired simulation

import numpy as np
from scipy.stats import norm

spot, strike, maturity, rate, vol = 100, 100, 1, 0.04, 0.20
n = 100_000
z = np.random.default_rng(96).standard_normal(n)
discount = np.exp(-rate * maturity)
d1 = (np.log(spot / strike) + (rate + vol**2 / 2) * maturity) / (vol * np.sqrt(maturity))
d2 = d1 - vol * np.sqrt(maturity)
reference = spot * norm.cdf(d1) - strike * discount * norm.cdf(d2)

for drift in (rate, 0.10):
    terminal = spot * np.exp((drift - vol**2 / 2) * maturity + vol * np.sqrt(maturity) * z)
    pv = discount * np.maximum(terminal - strike, 0)
    stderr = pv.std(ddof=1) / np.sqrt(n)
    print(drift, pv.mean(), stderr, (pv.mean() - reference) / stderr)
    print("discounted terminal stock", discount * terminal.mean())

Check the analytic value independently with the options pricer before measuring the Monte Carlo error.

Debugging Monte Carlo Option Prices Against Black–Scholes

Test the simulated stock before testing the option

Interpret the price difference in standard errors

A reproducible paired simulation

Change the drift and volatility of geometric Brownian motion

References and next steps