Notebooks · posts/black-scholes-by-monte-carlo

Black-Scholes by Monte Carlo, against the closed form

Adam Zettel · September 25, 2026 · written with an agent

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/ N$ : its measured slope on a log-log plot is $- 0.488$ , against $- 0.5$ , and the standard error $14.72/ N$ comes from an exact formula for the payoff's spread.

C = S_{0} N (d_{1}) - K e^{- r T} N (d_{2}), C^_{N} = \frac{e ^{- r T}}{N} n = 1 \sum N max (S_{T}^{(n)} - K, 0), SE = \frac{σ _{Y}}{N}

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.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

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 $σ$ 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 + σ ^{2} /2 ) T}{σ T}, d_{2} = d_{1} - σ 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 ((r - \frac{1}{2} σ^{2}) T + σ T Z), Z \sim 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} - K e^{- r T}$ , by any model.

[1]lines 47–74

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^{- r T} 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 $σ_{Y} / N$ , the standard error, where $σ_{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.

$σ_{Y}$ has a closed form too. The payoff's second moment needs $E [S_{T}^{2} 1_{S_{T} > K}] = S_{0}^{2} e^{(2 r + σ^{2}) T} N (d_{1} + σ T)$ , $E [S_{T} 1_{S_{T} > K}] = S_{0} e^{r T} N (d_{1})$ and $P (S_{T} > K) = N (d_{2})$ , and $σ_{Y}^{2} = E [Y^{2}] - C^{2}$ . Here $σ_{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.

[2]lines 100–189

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)

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 $lo g (RMS)$ against $lo g N$ for $N = 16, \dots, 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$ .

[3]lines 207–250

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")

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 $σ_{Y}$ rather than the exponent.

[4]lines 259–260

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