Finite differences, Monte Carlo, and quasi-Monte Carlo — the three families that price every derivative for which a closed form doesn't exist (i.e. nearly all of them).
Closed-form prices exist for European options under Black–Scholes — and almost nothing else. Everything else is a numerical method. Wilmott groups them into three:
Discretise the pricing PDE on a grid in $(S, t)$ space. Step backwards from the terminal payoff to today. Great when there's a PDE and the dimension is low (1–3).
Discretise the SDE for $S$. Simulate paths, average the payoff, discount. Error $\propto 1/\sqrt{N}$ regardless of dimension — the high-dimensional default choice.
When the price reduces to an integral — e.g. via a Fourier transform of the characteristic function for Heston or VG — just integrate it. The narrowest scope, but cheapest where it applies.
Almost every production pricer is some combination of these three. The fun is figuring out which one suits the contract: path-dependence steers you toward MC; early-exercise steers you toward FD; closed-form characteristic functions invite Fourier integration.
This deck covers chapters 27, 28, 29 and 30 of Wilmott. The interactive widget on slide 10 lets you price a European call by either Monte Carlo or explicit finite differences and watch the convergence (the $1/\sqrt N$ MC envelope and the FD grid error) against the closed-form Black–Scholes value.
The Black–Scholes PDE for option price $V(S, t)$:
$$\frac{\partial V}{\partial t} + \tfrac{1}{2} \sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + r S \frac{\partial V}{\partial S} - r V = 0.$$
Lay down a grid: $S_i = i \cdot \delta S$ for $i = 0, \ldots, I$ and $t_k = k \cdot \delta t$ for $k = 0, \ldots, N$. Three time-stepping schemes:
Use values at $t_{k+1}$ to compute $t_k$. Each new value is an explicit weighted average of three neighbours.
Pro: trivial to code.
Con: conditionally stable — $\delta t$ must shrink as $\delta S$ shrinks.
Solve a tridiagonal linear system at each time step.
Pro: unconditionally stable.
Con: needs a linear solver; first-order in time.
The average of explicit and implicit.
Pro: stable, second-order in both $\delta t$ and $\delta S$.
Con: can produce oscillations near non-smooth payoffs.
For the textbook Black–Scholes PDE, the workhorse in research is CN; the workhorse in production is often implicit with a small CN smoother near the strike (the "Rannacher start" trick) to suppress oscillations.
Three flavours of finite-difference approximation. Let $V_i^k$ denote $V(S_i, t_k)$.
$$\frac{\partial V}{\partial S} \approx \frac{V_{i+1}^k - V_i^k}{\delta S}.$$
First-order accurate. Asymmetric — biased toward the right neighbour.
$$\frac{\partial V}{\partial S} \approx \frac{V_i^k - V_{i-1}^k}{\delta S}.$$
First-order accurate. Asymmetric the other way.
$$\frac{\partial V}{\partial S} \approx \frac{V_{i+1}^k - V_{i-1}^k}{2 \delta S}.$$
Second-order accurate. The default for Δ — clean and symmetric.
For the second derivative (Gamma):
$$\frac{\partial^2 V}{\partial S^2} \approx \frac{V_{i+1}^k - 2 V_i^k + V_{i-1}^k}{(\delta S)^2}.$$
And for the time derivative, either a forward or backward step depending on the scheme:
$$\frac{\partial V}{\partial t} \approx \frac{V_i^{k+1} - V_i^k}{\delta t}.$$
The grid has boundaries; the PDE doesn't. You impose Dirichlet ($V = 0$ at $S = 0$ for a call) or asymptotic ($V \to S - K e^{-rT}$ as $S \to \infty$) conditions at the edges. The corners matter — a wrong boundary can corrupt the whole price.
Insert the finite differences into the PDE and rearrange:
$$V_i^k = a_i V_{i-1}^{k+1} + b_i V_i^{k+1} + c_i V_{i+1}^{k+1},$$
with coefficients
$$a_i = \tfrac{1}{2}\delta t (\sigma^2 i^2 - r i),\;\; b_i = 1 - \delta t (\sigma^2 i^2 + r),\;\; c_i = \tfrac{1}{2}\delta t (\sigma^2 i^2 + r i).$$
Step backwards from the terminal payoff $V_i^N = \max(S_i - K, 0)$. At each step, every interior $V_i^k$ is computed from three neighbours at the next time slice.
For all $a_i, b_i, c_i$ to be non-negative and the scheme to be stable:
$$\delta t \le \frac{1}{\sigma^2 I^2 + r}.$$
Halving $\delta S$ forces $\delta t$ to drop by a factor of four. The total cost is $O(I \cdot N) = O(I^3)$.
Oscillations that double in amplitude every step. The price blows up in a few iterations. Useful pedagogically — you can see it happen on a small grid.
Practical version: use a log-transformed grid in $x = \log S$, which kills the $i^2$ in the coefficients and gives a stable scheme with $\delta t = O(\delta x^2)$. The CN scheme on the log-grid is the standard implementation.
An American option can be exercised at any time before $T$. At every grid node, the holder compares two values:
The American option is worth the larger of the two:
$$V_i^k = \max\!\big(V_\text{cont}, \, \text{intrinsic}(S_i)\big).$$
At each $t_k$ there's a critical $S^*(t_k)$ above (call) or below (put) which exercise is optimal. The "max" step traces this boundary out for free — it's the kink in the value function.
The same idea works for implicit and CN schemes, but the "max" step has to be combined with the linear-solve step — this becomes a linear complementarity problem, and is the cleanest formulation of the American pricing PDE.
Backward induction trivially compares "hold" vs. "exercise" at each grid node. Monte Carlo runs forward in time and doesn't see future paths from a given state — the Longstaff–Schwartz trick (slide 08) is the way to retrofit it.
The risk-neutral price of a European-style option is
$$V_0 = e^{-rT}\, \mathbb{E}^\mathbb{Q}[\text{payoff}(S_T)].$$
MC estimates the expectation by simulating paths under $dS = r S\, dt + \sigma S\, dW$ and averaging the discounted payoff:
$$\hat V_0 = e^{-rT} \cdot \frac{1}{N} \sum_{n=1}^N \text{payoff}(S_T^{(n)}).$$
For lognormal $S$, you don't need to discretise: $S_T = S_0 \exp\!\big((r - \tfrac{1}{2}\sigma^2)T + \sigma \sqrt T\, Z\big)$ with $Z \sim \mathcal{N}(0,1)$.
The 95% confidence half-width is $\approx 1.96\, \hat\sigma / \sqrt{N}$ where $\hat\sigma$ is the sample std of the discounted payoffs. Always report this alongside the price.
From two i.i.d. uniform $(0,1)$ variables $U_1, U_2$, produce two i.i.d. standard normals:
$$Z_1 = \sqrt{-2 \ln U_1}\, \cos(2 \pi U_2),$$
$$Z_2 = \sqrt{-2 \ln U_1}\, \sin(2 \pi U_2).$$
Or the polar version (Marsaglia), which avoids the trig functions.
For every path driven by $Z$, also run a "mirror" path driven by $-Z$. The two estimators are negatively correlated, so their average has lower variance. Free (almost) — just one extra payoff per path.
Use a related instrument with a known price $g$ as a control:
$$\hat V_\text{CV} = \hat V - \beta\, (\hat g - g).$$
Choose $\beta$ to minimise variance. Spectacular when $V$ and $g$ are highly correlated (e.g. pricing an Asian via a geometric Asian).
All three are variance reduction techniques — they don't change the asymptotic $1/\sqrt N$ rate, but they shrink the constant in front. A 4× variance reduction is equivalent to 4× the paths at the same cost.
For a vanilla European, plain MC gets within 1% of BS at $N \approx 10^4$–$10^5$. With antithetic and a control variate, you can get there at $N \approx 10^3$.
An American option's value depends on the optimal exercise policy. To use MC, we need to know — on a forward-simulated path — whether to exercise at each step. But that requires knowing the continuation value, which is the very expectation we're trying to compute.
Longstaff and Schwartz (2001) solved this with a regression on basis functions:
The conditional expectation $\mathbb{E}[C^{k+1} \mid S^k]$ is a function of $S^k$. We approximate that function by least-squares regression on a few basis functions. As $N$ and the number of basis functions grow, the approximation converges to the true conditional expectation.
You're effectively learning the exercise boundary from the simulated paths themselves.
LS is the dominant method for American/Bermudan options when the state space is multi-dimensional and FD is impractical. It's the everyday workhorse for Bermudan-swaption pricers in fixed-income desks.
Random samples are great for variance estimates but they cluster — you get holes and clumps that look uniform only in the limit. Low-discrepancy (quasi-random) sequences fill space more evenly, by construction.
Uses the radical-inverse function in different prime bases for each dimension. Quick to implement; degrades in high dimensions because higher primes give noticeable correlation.
Built from a set of "direction numbers" via bitwise XOR. Excellent space-filling up to thousands of dimensions, especially with the modern Joe–Kuo direction numbers. The industry standard for QMC.
Error rate up to $O((\log N)^d / N)$ — essentially $1/N$ for moderate dimensions, vs. $1/\sqrt N$ for MC. For smooth payoffs in dimension < 50 this is a huge speedup.
The catch: error bounds rely on integrand smoothness, and confidence intervals are harder to get than for plain MC. The standard fix is randomised QMC — scramble or shift the Sobol sequence to recover a probabilistic error estimate while keeping most of the variance reduction.
Pricing a European option in 1D: probably overkill. Pricing a 30-dimensional CMS spread option with mostly smooth payoff: a 10×–100× speedup over plain MC is realistic. The crossover depends heavily on payoff smoothness.
A European call. Pick $S_0$, $K$, $T$, $\sigma$, $r$, $N$ (paths for MC / grid resolution for FD), and a method. For MC you'll see ~30 sample paths and the payoff histogram; for FD you'll see a heatmap of the option-value surface $V(S, t)$. The metric grid reports the estimated price, the standard error or grid error, and the closed-form Black–Scholes value.
Toggle between MC and FD to see the trade-off. The MC histogram shows how the bulk of paths finish out-of-the-money (payoff = 0) and a long right tail carries the value — this is why MC for far-out-of-the-money options is so noisy. The FD heatmap is the entire solution $V(S, t)$ — one solve gives you the price for every spot.