Monte Carlo Methods in Finance

Risk-neutral valuation (Quant 03) says an option's price today is the discounted expectation of its payoff under the risk-neutral measure \(\mathbb{Q}\): \(V_0 = e^{-rT}\,\mathbb{E}^{\mathbb{Q}}[\,\text{payoff}\,]\). For a European call on a single lognormal stock that expectation has a closed form — the Black–Scholes formula. But the moment the payoff depends on the whole path (an Asian average, a barrier that may knock out), or on many assets at once (a basket, a worst-of), the integral becomes high-dimensional and no formula survives. Monte Carlo replaces the integral with an average over simulated futures.

The justification is the law of large numbers: draw \(N\) independent payoff samples \(Y_1, \dots, Y_N\) with mean \(\mu = \mathbb{E}[Y]\), and their sample mean converges to \(\mu\). That sample mean, discounted, is our price estimate:

How wrong can \(\hat{V}_N\) be? The central limit theorem turns the spread of the samples into an error bar. With \(\sigma_Y\) the standard deviation of the (discounted) payoff, the estimator's own standard deviation — its standard error — is

Why bother, when grids and trees exist? Because their cost explodes with dimension. A finite-difference PDE solver or a binomial tree (Quant 02) is excellent for one or two underlyings but suffers the curse of dimensionality: work grows exponentially in the number of state variables. Monte Carlo's \(1/\sqrt{N}\) error is independent of dimension — the same hundred-thousand paths price a one-asset call and a fifty-asset basket. Above three or four dimensions, simulation is the only game in town.

To average payoffs you must first generate scenarios — and they must be drawn under \(\mathbb{Q}\), where the underlying drifts at the riskless rate \(r\) rather than its real-world drift \(\mu\). For a stock following geometric Brownian motion (Quant 03), one terminal price needs exactly one Gaussian draw, because the exact solution of the SDE is known in closed form:

When you do need the whole path — for an Asian average or a barrier check — you discretize the SDE over a time grid \(0 = t_0 < t_1 < \cdots < t_m = T\). The log-Euler (exact-in-log) scheme steps the log-price with one independent Gaussian per step:

The whole recipe for a European price is then four lines: draw \(Z\), map to \(S_T\), apply the payoff, discount the mean. The first Python cell does exactly this and lays the result alongside the closed form so you can see the agreement — and the standard error that quantifies it.

# Monte-Carlo price a European call vs Black-Scholes, with a standard error
import numpy as np, math
rng = np.random.default_rng(0)

def N(x):  # standard normal CDF via the error function
    return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))
def bs_call(S, K, r, sig, T):
    d1 = (math.log(S/K) + (r + 0.5*sig*sig)*T) / (sig*math.sqrt(T))
    d2 = d1 - sig*math.sqrt(T)
    return S*N(d1) - K*math.exp(-r*T)*N(d2)

S0, K, r, sig, T = 100.0, 105.0, 0.05, 0.20, 1.0
M = 200_000
Z  = rng.standard_normal(M)
ST = S0 * np.exp((r - 0.5*sig*sig)*T + sig*math.sqrt(T)*Z)   # EQ Q5.3, risk-neutral
disc_payoff = math.exp(-r*T) * np.maximum(ST - K, 0.0)        # discounted payoffs

price = disc_payoff.mean()                                    # EQ Q5.1
se    = disc_payoff.std(ddof=1) / math.sqrt(M)                # EQ Q5.2
exact = bs_call(S0, K, r, sig, T)

print(f"Monte-Carlo call : {price:.4f}   +/- {1.96*se:.4f}  (95% CI)")
print(f"Black-Scholes    : {exact:.4f}")
print(f"gap / SE         : {(price-exact)/se:+.2f}  (should sit within ~2)")
print(f"exact lies in CI : {abs(price-exact) <= 1.96*se}")

Because the error is \(\sigma_Y/\sqrt{N}\), there are only two ways to make a Monte Carlo estimate more accurate: run more paths (linear cost, square-root payoff) or shrink \(\sigma_Y\) itself. The second is almost free and is what separates a textbook simulation from a production one. Two techniques dominate.

Antithetic variates

The Gaussian is symmetric: if \(Z\) is a valid draw, so is \(-Z\), with the same probability. So pair every path driven by \(Z\) with a mirror path driven by \(-Z\), and average the two payoffs into one sample. The pairs cost one extra (cheap) evaluation but reuse the same random number, and — for a payoff that moves monotonically with \(Z\) — the two halves are negatively correlated, which cancels noise:

Control variates

If you can't reduce a payoff's spread, subtract off a correlated quantity whose expectation you already know. Let \(X\) be a control with known mean \(\mathbb{E}[X]\) (for an Asian option, the stock itself, or a geometric-average Asian which prices in closed form). Form a corrected estimator:

These compose. Production pricers stack antithetics, a strong control variate, and quasi-random (low-discrepancy Sobol) sequences, which replace pseudo-random points with a deliberately even sweep of the unit cube and can lift the convergence rate toward \(1/N\) in low effective dimension. The discipline is always the same: never grow \(N\) when you can shrink \(\sigma_Y\).

# Antithetic variates: measure the variance reduction vs plain Monte Carlo
import numpy as np, math
rng = np.random.default_rng(1)

S0, K, r, sig, T = 100.0, 105.0, 0.05, 0.20, 1.0
M = 100_000                                   # plain uses M draws; antithetic M/2 pairs
drift, vol, disc = (r - 0.5*sig*sig)*T, sig*math.sqrt(T), math.exp(-r*T)
def call_payoff(Z):
    return disc * np.maximum(S0*np.exp(drift + vol*Z) - K, 0.0)

# plain MC: M independent draws
Yp = call_payoff(rng.standard_normal(M))
price_plain, se_plain = Yp.mean(), Yp.std(ddof=1)/math.sqrt(M)

# antithetic: M/2 shocks, each paired with its mirror -Z, averaged into one sample
Zh = rng.standard_normal(M//2)
Ya = 0.5 * (call_payoff(Zh) + call_payoff(-Zh))   # M/2 antithetic samples
price_anti, se_anti = Ya.mean(), Ya.std(ddof=1)/math.sqrt(M//2)

rho = np.corrcoef(call_payoff(Zh), call_payoff(-Zh))[0, 1]
print(f"plain MC     : {price_plain:.4f}  SE {se_plain:.4f}")
print(f"antithetic   : {price_anti:.4f}  SE {se_anti:.4f}")
print(f"mirror corr  : {rho:+.3f}   (negative -> variance falls)")
print(f"variance cut : {(se_plain/se_anti)**2:.2f}x   for the same {M:,} evaluations")

Monte Carlo earns its keep where closed forms stop. Two workhorse exotics show why. Both depend not on \(S_T\) alone but on the trajectory the price took to get there — so you must simulate the whole path (EQ Q5.4), not just the endpoint.

An Asian option pays on the average price over a set of monitoring dates, not the final price. Averaging damps both the payoff's volatility and any last-minute manipulation, which is why Asians are popular in commodity and FX markets:

A barrier option activates or extinguishes if the price touches a level \(B\) at any point. A down-and-out call, say, pays the usual call payoff unless the path ever falls to \(B\), in which case it is worthless:

The pattern generalizes: lookbacks (payoff on the running max/min), cliquets (sums of capped period returns), autocallables (early-redemption ladders), and worst-of baskets across many correlated underlyings are all priced by the same loop — simulate paths, apply the payoff functional, discount the mean. The payoff code changes; the engine does not. That uniformity is exactly why simulation became the lingua franca of the exotics desk.

A price is only half the job — desks hedge with the Greeks (Quant 03), and those derivatives must come out of the same simulation. The naïve approach, finite differences (bump-and-revalue), reprices at \(S_0\) and \(S_0+h\) and takes the slope. It is universal but treacherous:

Two cleaner families avoid the bump entirely. The pathwise (infinitesimal perturbation) estimator differentiates the payoff itself — for a call, \(\partial_{S_0}\max(S_T-K,0) = \mathbf{1}\{S_T>K\}\,S_T/S_0\) — giving an unbiased delta in a single pass, valid whenever the payoff is Lipschitz (it fails at the kink of a digital). The likelihood-ratio method instead differentiates the probability density, moving the derivative onto a smooth weight so it works even for discontinuous payoffs, at the cost of higher variance. Modern pricers increasingly use adjoint algorithmic differentiation (AAD), which computes all Greeks in one reverse pass at a cost independent of the number of parameters — the same trick that powers backpropagation in deep learning.

1. Boyle, P. (1977). Options: A Monte Carlo Approach. Journal of Financial Economics 4(3) — the paper that brought simulation to option pricing.
2. Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer — the definitive reference for variance reduction, path simulation, and Greeks.
3. Boyle, P., Broadie, M. & Glasserman, P. (1997). Monte Carlo Methods for Security Pricing. Journal of Economic Dynamics and Control 21 — the canonical survey of MC techniques in finance.
4. Longstaff, F. & Schwartz, E. (2001). Valuing American Options by Simulation: A Simple Least-Squares Approach. Review of Financial Studies 14(1) — least-squares Monte Carlo for early-exercise payoffs.
5. Giles, M. & Glasserman, P. (2006). Smoking Adjoints: Fast Monte Carlo Greeks. RISK Magazine — adjoint algorithmic differentiation for computing all Greeks in one reverse pass.
6. Kloeden, P. & Platen, E. (1992). Numerical Solution of Stochastic Differential Equations. Springer — Euler–Maruyama and higher-order path-discretization schemes.