Monte Carlo Simulation for Option Pricing
Trees and grids both fight the same enemy as they scale up: add a second underlying asset, and a binomial
tree's node count and a finite-difference grid's cell count both explode — two dimensions of stock price
instead of one, then three, then ten for a basket of correlated names. Neither method survives that
growth. Monte Carlo simulation sidesteps the problem entirely, because it never builds a grid at all — it
just plays the risk-neutral stock
process forward, many times, on a computer, and averages what happens. It is how a desk
prices the payoffs that a tree or a grid simply cannot reach: baskets on ten names, payoffs that depend on
an entire price path rather than just the final value. This page builds the recipe from the discounted
expectation you already know from
American
options on a tree, runs it live, and is honest about where it is overkill and where nothing
else will do.
The recipe: simulate, discount, average
Risk-neutral pricing says V_0 = e^{-rT}\,\mathbb{E}_{\mathbb{Q}}[\text{payoff}].
An expectation is just an average over enough random draws — so simulate the expectation instead
of solving for it. Under \mathbb{Q} the stock follows risk-neutral geometric
Brownian motion, dS_t = rS_t\,dt + \sigma S_t\,d\tilde W_t, whose solution gives
the terminal price directly in closed form:
S_T = S_0\,\exp\!\Bigl(\bigl(r - \tfrac12\sigma^2\bigr)T + \sigma\sqrt{T}\,Z\Bigr), \qquad Z \sim N(0,1).
- Draw N independent standard normals Z_1,\dots,Z_N
and turn each into a simulated terminal price S_T^{(i)} using the formula
above.
- Evaluate the payoff on each simulated path, e.g. \max(S_T^{(i)}-K, 0) for
a call, and average: \hat V_0 = e^{-rT}\cdot\frac{1}{N}\sum_i \text{payoff}^{(i)}.
- The estimate's standard error shrinks like \sigma_{\text{payoff}}/\sqrt{N}
— the same slow-but-dimension-blind law from
Monte Carlo Methods.
Crucially, that error rate does not depend on how many underlying assets the payoff involves.
For a single European call this is spectacular overkill — Black–Scholes gives the exact answer instantly,
with no sampling error at all. The method earns its keep the moment the payoff needs the whole
path, not just the endpoint. But it is worth seeing it nail the easy case first, so you know exactly
what "correct" looks like before trusting it on a hard one.
// Monte Carlo pricer for a European call, sampling S_T directly from the
// closed-form risk-neutral GBM solution (no need to step through time for
// a payoff that only depends on the terminal price).
function mulberry32(seed: number) {
let a = seed >>> 0;
return function (): number {
a = (a + 0x6d2b79f5) | 0;
let t = Math.imul(a ^ (a >>> 15), 1 | a);
t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
};
}
const rand = mulberry32(2024);
function gauss(): number {
// Box-Muller: two uniforms in -> one standard normal out.
const u1 = Math.max(rand(), 1e-12);
const u2 = rand();
return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
}
const S0 = 100, K = 100, r = 0.05, sig = 0.2, T = 1;
// Reference: the exact Black-Scholes price for these inputs is 10.4506.
function mcCall(N: number): { price: number; se: number } {
let sum = 0;
let sumSq = 0;
for (let i = 0; i < N; i++) {
const Z = gauss();
const ST = S0 * Math.exp((r - 0.5 * sig * sig) * T + sig * Math.sqrt(T) * Z);
const payoff = Math.max(ST - K, 0);
sum += payoff;
sumSq += payoff * payoff;
}
const mean = sum / N;
const variance = sumSq / N - mean * mean;
const price = Math.exp(-r * T) * mean;
const se = Math.exp(-r * T) * Math.sqrt(variance / N);
return { price, se };
}
for (const N of [1000, 10000, 50000, 200000]) {
const { price, se } = mcCall(N);
console.log(
`N=${N.toString().padStart(6)} price=${price.toFixed(4)} ±${se.toFixed(4)} (Black-Scholes = 10.4506)`,
);
}
Watch the price wander toward 10.4506 as N grows,
and the \pm standard-error column shrink alongside it — but only by a factor of
about \sqrt{10}\approx3.16 each time N is multiplied
by ten. Change the seed and re-run: the price moves, but always within roughly one standard error of
10.4506.
Where terminal sampling breaks down: path-dependent payoffs
Sampling S_T alone worked because a European call's payoff only cares about the
stock's value at the very end. Many real contracts don't have that luxury. An
Asian
option, for instance, pays off on the average stock price over the option's life —
there is no way to know that average without walking through the whole path. Trees can do this only by
exploding the state space (the tree must remember the running average at every node, destroying
recombination); Monte Carlo barely notices, because it was already generating a full path's worth of
randomness one time-step at a time — we just need to step through it instead of jumping straight to
T:
// Monte Carlo pricer for an Asian call (payoff on the ARITHMETIC AVERAGE
// stock price along the path) -- this needs the full path, not just S_T.
function mulberry32(seed: number) {
let a = seed >>> 0;
return function (): number {
a = (a + 0x6d2b79f5) | 0;
let t = Math.imul(a ^ (a >>> 15), 1 | a);
t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
};
}
const rand = mulberry32(777);
function gauss(): number {
const u1 = Math.max(rand(), 1e-12);
const u2 = rand();
return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
}
const S0 = 100, K = 100, r = 0.05, sig = 0.2, T = 1;
const steps = 50;
const dt = T / steps;
function asianPathPayoff(): number {
let S = S0;
let sum = 0;
for (let i = 0; i < steps; i++) {
const Z = gauss();
// Step the risk-neutral GBM forward by dt, then record the price.
S = S * Math.exp((r - 0.5 * sig * sig) * dt + sig * Math.sqrt(dt) * Z);
sum += S;
}
const average = sum / steps;
return Math.max(average - K, 0);
}
function mcAsian(N: number): number {
let sum = 0;
for (let i = 0; i < N; i++) sum += asianPathPayoff();
return Math.exp(-r * T) * (sum / N);
}
for (const N of [2000, 10000, 50000]) {
console.log(`N=${N.toString().padStart(5)} Asian call price ≈ ${mcAsian(N).toFixed(4)}`);
}
console.log("(No closed form exists for the arithmetic-average Asian call --");
console.log(" simulation is not a shortcut here, it's the only general method.)");
Notice the estimate settles well below the plain European call's 10.45
— averaging the path damps the payoff's volatility, since an extreme terminal price gets diluted by all
the ordinary prices that came before it. That is a genuine, economically meaningful difference between
two contracts, and Monte Carlo captured it just by keeping the full path instead of jumping to the end.
They don't just buy more compute — they make each path work harder. Antithetic variates
pairs every draw Z with its mirror image -Z: since a
path and its antithetic partner are negatively correlated, averaging the two payoffs cancels out some of
the noise for free, at essentially no extra cost (one extra path per pair, reusing the same random
numbers). Control variates is cleverer still: price a similar instrument that
does have a closed form (e.g. use the geometric-average Asian option, which has one, to correct
the arithmetic-average Asian option, which doesn't) on the very same simulated paths, and subtract off the
known error. Both tricks shrink the effective \sigma_{\text{payoff}} in the
1/\sqrt{N} law rather than fighting the law itself — routinely buying a
five- to fifty-fold reduction in the paths needed for a given accuracy, which is the difference between a
pricing run that finishes in a coffee break and one that doesn't finish before the market closes.
-
More paths fixes variance, not bias. Simulating under the real-world drift
\mu instead of r, or forgetting the discount
factor e^{-rT} altogether, are systematic errors — running a
billion paths of the wrong model just gives you an extremely precise wrong answer. Only genuine sampling
noise shrinks with N; a wrong assumption does not.
-
Too few time steps is its own bias, separate from too few paths. The Asian-option code
above steps through 50 discrete points to approximate a continuous average.
Cut that to 5 steps and every path becomes a cruder approximation of the true
continuous path — a discretization error that more paths at the same coarse step count cannot
fix, only more (and finer) time steps can.
-
The reported standard error is an estimate of precision, not of correctness. A tight
\pm band around a systematically wrong number (see the first point) looks
reassuring and is completely misleading — always sanity-check a Monte Carlo price against a closed form
or a tree on a case where one is available, exactly as this page did for the European call.