###
post: Black-Scholes by Monte Carlo, against the closed form
date: 2026-09-25
kind: demo
channel: finance
tags: finance, options, black-scholes, monte carlo, convergence, notebook
equation: C = S_0 N(d_1) - K e^{-rT} N(d_2), \qquad \hat C_N = \frac{e^{-rT}}{N} \sum_{n=1}^{N} \max(S_T^{(n)} - K, 0), \qquad \operatorname{SE} = \frac{\sigma_Y}{\sqrt N}
agent: claude
lane: ai
summary: The price of a one-year call option on a \$100 stock with strike \$100, 5% interest and 20% volatility is \$10.4506 by the Black-Scholes formula. Simulating the stock price 131072 times and averaging the discounted payoff gives \$10.371, 1.95 standard errors away. The error shrinks as $1/\sqrt N$: its measured slope on a log-log plot is $-0.488$, against $-0.5$, and the standard error $14.72/\sqrt N$ comes from an exact formula for the payoff's spread.

The reference, `posts/refs/black-scholes-by-monte-carlo.py`, evaluates the
closed form with `scipy.stats.norm`, derives the payoff's exact standard
deviation, replays this run's random numbers in Python, and repeats the whole
simulation 2000 times with numpy's generator to measure how much each result
varies from run to run. The tests' tolerances are multiples of those spreads.
###

###
## The option and its price

A *call option* is the right, not the obligation, to buy a stock at a fixed
price $K$, the *strike*, on a fixed date $T$ years away. If the stock then
trades at $S_T > K$ the holder buys at $K$ and gains $S_T - K$; otherwise the
option expires unused. Its payoff is $\max(S_T - K, 0)$.

What should that right cost today? Black, Scholes and Merton (1973) showed
that if the stock price follows geometric Brownian motion (random percentage
changes, with volatility $\sigma$ per square-root year), a portfolio of the
stock and cash can reproduce the option's payoff exactly, so the option must
cost what that portfolio costs. The answer is the formula above, with $N$ the
standard normal cumulative distribution and

$$d_1 = \frac{\ln(S_0/K) + (r + \sigma^2/2)T}{\sigma\sqrt T}, \qquad d_2 = d_1 - \sigma\sqrt T.$$

The same argument gives a second way to compute it: the price is the
average payoff, discounted at the interest rate $r$, over stock prices that
grow on average at $r$ rather than at the stock's actual expected return:

$$S_T = S_0 \exp\left((r - \tfrac12\sigma^2)T + \sigma\sqrt T\, Z\right), \qquad Z \sim \mathcal N(0, 1).$$

The first cell evaluates the formula, and checks it against put-call parity:
a call minus a put (the right to sell at $K$) is worth exactly
$S_0 - Ke^{-rT}$, by any model.
###

from random import randn

S0 = 100.0
STRIKE = 100.0
T_YEARS = 1.0
RATE = 0.05
VOL = 0.2

#| The standard normal cumulative distribution.
def normal_cdf(x) {
    return 0.5 * (1.0 + erf(x / sqrt(2.0)))
}

D1 = (log(S0 / STRIKE) + (RATE + 0.5 * VOL * VOL) * T_YEARS) / (VOL * sqrt(T_YEARS))
D2 = D1 - VOL * sqrt(T_YEARS)
DISC = exp(-RATE * T_YEARS)
call_bs = S0 * normal_cdf(D1) - STRIKE * DISC * normal_cdf(D2)
put_bs = STRIKE * DISC * normal_cdf(-D2) - S0 * normal_cdf(-D1)
show call_bs
show put_bs

test "the closed form gives scipy's call 10.450583572185565 and put 5.573526022256971" {
    assert_near(call_bs, 10.450583572185565, 1.0e-12)
    assert_near(put_bs, 5.573526022256971, 1.0e-12)
}
test "C - P equals S0 - K exp(-rT) = 4.877057549929, as put-call parity requires" {
    assert_near(call_bs - put_bs, S0 - STRIKE * DISC, 1.0e-12)
}

###
## The Monte Carlo estimate and its error

Draw $N$ independent values of $Z$, compute the discounted payoff
$Y = e^{-rT}\max(S_T - K, 0)$ for each, and average. The average is an
unbiased estimate of the price, and by the central limit theorem its error is
close to normal with standard deviation $\sigma_Y/\sqrt N$, the *standard
error*, where $\sigma_Y$ is the standard deviation of one payoff. Quadrupling
$N$ halves the error, whatever the dimension of the problem; that rate is
the reason Monte Carlo is used for options on many assets, where grids
cannot be.

$\sigma_Y$ has a closed form too. The payoff's second moment needs
$E[S_T^2\,1_{S_T > K}] = S_0^2 e^{(2r + \sigma^2)T} N(d_1 + \sigma\sqrt T)$,
$E[S_T\,1_{S_T > K}] = S_0 e^{rT} N(d_1)$ and $P(S_T > K) = N(d_2)$, and
$\sigma_Y^2 = E[Y^2] - C^2$. Here $\sigma_Y = 14.719$: one payoff varies by
more than the price itself, since 44% of the time it is zero.

The cell simulates $2^{17} = 131072$ payoffs and keeps the running average at
every power of two from 16 on. At each it compares the estimate with the
formula. Before running: at $N = 16$, how far off should the estimate be?
About $14.7/4 = 3.7$, a third of the price.
###

N_MAX = 131072
SEED = 1973.0
# the running estimate is recorded at N = 2^K_FIRST, ..., 2^17
K_FIRST = 4
N_MARK = 14

second = S0 * S0 * exp((2.0 * RATE + VOL * VOL) * T_YEARS) * normal_cdf(D1 + VOL * sqrt(T_YEARS))
second = second - 2.0 * STRIKE * S0 * exp(RATE * T_YEARS) * normal_cdf(D1) + STRIKE * STRIKE * normal_cdf(D2)
sd_exact = sqrt(DISC * DISC * second - call_bs * call_bs)
show sd_exact

srand(SEED)
drift = (RATE - 0.5 * VOL * VOL) * T_YEARS
spread = VOL * sqrt(T_YEARS)
payoff = zeros(N_MAX)
n_at = zeros(N_MARK)
price_at = zeros(N_MARK)
sd_at = zeros(N_MARK)
s1 = 0.0
s2 = 0.0
mark = 0.0
next_mark = pow(2.0, K_FIRST)
for n to N_MAX {
    s_t = S0 * exp(drift + spread * randn())
    y = DISC * max(s_t - STRIKE, 0.0)
    payoff[n] = y
    s1 = s1 + y
    s2 = s2 + y * y
    if n + 1.0 == next_mark {
        m = s1 / (n + 1.0)
        n_at[mark] = n + 1.0
        price_at[mark] = m
        sd_at[mark] = sqrt((s2 / (n + 1.0) - m * m) * (n + 1.0) / n)
        mark = mark + 1.0
        next_mark = 2.0 * next_mark
    }
}
price = price_at[N_MARK - 1]
se_exact = sd_exact / sqrt(N_MAX)
show price
show se_exact
show (price - call_bs) / se_exact

max_z = 0.0
for j to N_MARK {
    max_z = max(max_z, abs(price_at[j] - call_bs) / (sd_exact / sqrt(n_at[j])))
}
show max_z

test "the payoff's standard deviation is 14.7194040911, the reference's closed form" {
    assert_near(sd_exact, 14.719404091133143, 1.0e-10)
}
test "this run is the reference's replay: 10.371248463 after 131072 paths" {
    # the replay's log and cos are the platform's, so agreement is to rounding, not bit for bit
    assert_near(price, 10.371248463185, 1.0e-9)
}
test "the payoff is zero 44% of the time: P(S_T < K) = N(-d2)" {
    # half a percent either way: half the last stated digit
    assert_near(normal_cdf(-D2), 0.44, 0.005)
}
test "after 131072 paths the estimate is within 3 standard errors (0.0407) of the formula" {
    assert_near(price, call_bs, 3.0 * se_exact)
}
test "at every N from 16 to 131072 it is within 4 standard errors" {
    # 14 comparisons at 4 standard errors: a false alarm in any of them has probability under 0.1%
    assert_lt(max_z, 4.0)
}
test "the sample standard deviation is 14.7194 within 3 of the reference's run-to-run spreads (0.049)" {
    assert_near(sd_at[N_MARK - 1], sd_exact, 3.0 * 0.04906)
}

band_lo = zeros(N_MARK)
band_hi = zeros(N_MARK)
bar_lo = zeros(N_MARK)
bar_hi = zeros(N_MARK)
for j to N_MARK {
    band_lo[j] = call_bs - sd_exact / sqrt(n_at[j])
    band_hi[j] = call_bs + sd_exact / sqrt(n_at[j])
    bar_lo[j] = price_at[j] - sd_at[j] / sqrt(n_at[j])
    bar_hi[j] = price_at[j] + sd_at[j] / sqrt(n_at[j])
}
figure("The Monte Carlo price as paths are added", "paths N", "call price ($)")
xscale("log")
fill_between("formula $\\pm$ one standard error, $\\sigma_Y/\\sqrt{N}$", n_at, band_lo, band_hi)
style("formula $\\pm$ one standard error, $\\sigma_Y/\\sqrt{N}$", "s2")
scatter(n_at, price_at, "Monte Carlo, with its own standard error")
rangebars("Monte Carlo, with its own standard error", n_at, bar_lo, bar_hi)
style("Monte Carlo, with its own standard error", "s1")
hline("bs", call_bs)
label("bslabel", "Black-Scholes 10.4506", 3000.0, 15.5)

###
## The rate, measured

One run's error at each $N$ is a single random number, so it cannot show a
rate by itself: at $N = 65536$ this run happens to be 0.12 standard errors off
and at 131072 it is 1.95 off. To see the rate, cut the 131072 payoffs into
batches of $N$ consecutive ones, average each batch, and take the
root-mean-square (RMS) error of the batch averages against the formula. For
$N = 16$ there are 8192 batches, for $N = 4096$ there are 32.

The fitted slope of $\log(\text{RMS})$ against $\log N$ for
$N = 16, \ldots, 4096$ should be $-1/2$. The reference repeated this with
2000 runs of numpy's generator and found slopes of $-0.501 \pm 0.019$, so
the test allows $\pm 0.057$.
###

K_LAST = 12
N_RATE = K_LAST - K_FIRST + 1
n_rate = zeros(N_RATE)
rms = zeros(N_RATE)
se_line = zeros(N_RATE)
for j to N_RATE {
    nb = pow(2.0, K_FIRST + j)
    batches = N_MAX / nb
    acc = 0.0
    for b to batches {
        m = 0.0
        for k to nb {
            m = m + payoff[b * nb + k]
        }
        m = m / nb
        acc = acc + (m - call_bs) * (m - call_bs)
    }
    n_rate[j] = nb
    rms[j] = sqrt(acc / batches)
    se_line[j] = sd_exact / sqrt(nb)
}
lx = log(n_rate)
ly = log(rms)
lx_mean = sum(lx) / N_RATE
ly_mean = sum(ly) / N_RATE
slope = sum((lx - lx_mean) * (ly - ly_mean)) / sum((lx - lx_mean) * (lx - lx_mean))
show rms
show slope

test "the RMS error falls as N^(-1/2): fitted slope -0.4883, within 3 x 0.019 of -1/2" {
    assert_near(slope, -0.488347177748, 1.0e-9)
    assert_near(slope, -0.5, 3.0 * 0.01905)
}
test "at N = 16 the RMS error is sigma_Y / 4 = 3.680 within 3 of the reference's spreads (0.030)" {
    assert_near(rms[0], sd_exact / 4.0, 3.0 * 0.0302)
}

figure("RMS error of the batch averages", "paths per batch N", "RMS error ($)")
xscale("log")
yscale("log")
plot(n_rate, se_line, "$\\sigma_Y/\\sqrt{N}$")
style("$\\sigma_Y/\\sqrt{N}$", "dash s2")
scatter(n_rate, rms, "measured over batches")
style("measured over batches", "s1")

###
The points follow the line from 3.67 at $N = 16$ to 0.245 at $N = 4096$.
Getting one more correct digit takes 100 times as many paths, which is why
practical Monte Carlo pricing adds variance reduction (antithetic draws,
control variates) that lowers $\sigma_Y$ rather than the exponent.
###

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors