Chapter 11 — Numerical Methods for Option Pricing
"Numerical methods are not approximations to the theory — they are the theory made computable."
After this chapter you will be able to:
- Build a CRR binomial tree and price both European and American options with it
- Implement explicit, implicit, and Crank-Nicolson finite difference schemes for the Black-Scholes PDE
- Explain the CFL stability condition and why Crank-Nicolson is preferred in practice
- Quantify $O(1/n)$ convergence for lattice methods and apply Richardson extrapolation to achieve $O(1/n^2)$
- Select the right numerical method for each option type (European, American, barrier, path-dependent)
Black, Scholes, and Merton provided a beautiful closed-form formula for European options in 1973. But the real world quickly demands more: American options that can be exercised early, barrier options that knock in or out if the price crosses a level, Asian options whose payoff depends on the average price over a period. For almost all of these, no closed-form formula exists. Numerical methods are not a fallback — they are the primary tool of options practice.
The two foundational approaches developed here correspond to two equivalent ways of thinking about the Black-Scholes PDE. Lattice methods (binomial and trinomial trees) discretise the asset price process itself: at each time step, the price moves up or down by a specific factor, and option values are computed by backward induction through the lattice. The elegance of this approach is that early exercise is handled trivially — at each node, simply take the maximum of the continuation value and the intrinsic value. Cox, Ross, and Rubinstein introduced their parameterisation in 1979, and it remains the standard benchmark even today. Finite difference methods discretise the Black-Scholes PDE directly: the continuous derivatives $\partial V / \partial t$ and $\partial^2 V / \partial S^2$ are replaced by discrete differences on a grid. The Crank-Nicolson method, which averages the explicit and implicit schemes, achieves second-order accuracy in both time and space and is the gold standard for one-dimensional option PDEs.
This chapter implements both families and studies their convergence. The key insight is that convergence analysis is quantitative: we can say that CRR trees converge as $O(1/n)$ and Crank-Nicolson as $O(\Delta t^2)$, and these rates determine how fine a grid we need to achieve a target accuracy. We also implement Richardson extrapolation, which uses the known convergence rate to eliminate the leading error term and achieve $O(1/n^2)$ accuracy from two tree evaluations.
11.1 Binomial Tree Models
The binomial tree discretises the asset price process into a lattice of up and down moves. It converges to Black-Scholes as the number of steps increases.
11.1.1 Cox-Ross-Rubinstein (CRR) Parameterisation
For $n$ steps over time $T$, with step size $\Delta t = T/n$:
$$u = e^{\sigma\sqrt{\Delta t}}, \quad d = 1/u = e^{-\sigma\sqrt{\Delta t}}$$
$$p = \frac{e^{(r-q)\Delta t} - d}{u - d} \quad \text{(risk-neutral up probability)}$$
module Binomial = struct
type node = {
price : float;
value : float; (* option value *)
}
(**
CRR binomial tree option pricer.
Handles both European and American exercise.
n_steps: typically 200-1000 for good accuracy
*)
let price ~spot ~strike ~rate ?(div_yield = 0.0) ~vol ~tau ~n_steps
~option_type ~exercise =
let dt = tau /. float_of_int n_steps in
let u = exp (vol *. sqrt dt) in
let d = 1.0 /. u in
let p = (exp ((rate -. div_yield) *. dt) -. d) /. (u -. d) in
let q = 1.0 -. p in
let df = exp (-. rate *. dt) in
(* Terminal node prices *)
let terminal_prices = Array.init (n_steps + 1) (fun j ->
spot *. (u ** float_of_int (n_steps - j)) *. (d ** float_of_int j)
) in
(* Terminal option values *)
let values = Array.map (fun s ->
match option_type with
| `Call -> Float.max 0.0 (s -. strike)
| `Put -> Float.max 0.0 (strike -. s)
) terminal_prices in
(* Backward induction *)
for step = n_steps - 1 downto 0 do
for j = 0 to step do
let s = spot *. (u ** float_of_int (step - j)) *. (d ** float_of_int j) in
let continuation = df *. (p *. values.(j) +. q *. values.(j + 1)) in
let intrinsic = match option_type with
| `Call -> Float.max 0.0 (s -. strike)
| `Put -> Float.max 0.0 (strike -. s)
in
values.(j) <- match exercise with
| `European -> continuation
| `American -> Float.max intrinsic continuation
done
done;
values.(0)
(**
Convergence: compare n=50, 100, 200, 500 to Black-Scholes.
CRR oscillates; use odd/even average to reduce oscillation.
*)
let smooth_price ~spot ~strike ~rate ~vol ~tau ~option_type =
let p1 = price ~spot ~strike ~rate ~vol ~tau ~n_steps:499
~option_type ~exercise:`European in
let p2 = price ~spot ~strike ~rate ~vol ~tau ~n_steps:500
~option_type ~exercise:`European in
(p1 +. p2) /. 2.0 (* Richardson extrapolation / averaging *)
(** Greek via tree: bumped price *)
let delta_tree ~spot ~strike ~rate ~vol ~tau ~n_steps ~option_type ~exercise =
let eps = spot *. 0.001 in
let p_up = price ~spot:(spot +. eps) ~strike ~rate ~vol ~tau ~n_steps
~option_type ~exercise in
let p_down = price ~spot:(spot -. eps) ~strike ~rate ~vol ~tau ~n_steps
~option_type ~exercise in
(p_up -. p_down) /. (2.0 *. eps)
end
11.1.2 The Early Exercise Premium
The American option price exceeds the European price by the early exercise premium:
$$C_{\text{American}} - C_{\text{European}} = 0 \quad \text{(calls without dividends)}$$ $$P_{\text{American}} - P_{\text{European}} \geq 0 \quad \text{(puts — always exercise early)}$$
American calls on dividend-paying stocks may be exercised early just before the ex-dividend date.
11.2 Trinomial Trees
Trinomial trees add a middle (unchanged) outcome, giving better accuracy per node:
$$u = e^{\sigma\sqrt{3\Delta t}}, \quad m = 1, \quad d = 1/u$$
$$p_u = \frac{1}{6} + \frac{(r-q-\sigma^2/2)\sqrt{\Delta t}}{2\sigma\sqrt{3}}$$ $$p_m = \frac{2}{3}$$ $$p_d = \frac{1}{6} - \frac{(r-q-\sigma^2/2)\sqrt{\Delta t}}{2\sigma\sqrt{3}}$$
module Trinomial = struct
let price ~spot ~strike ~rate ?(div_yield = 0.0) ~vol ~tau ~n_steps
~option_type ~exercise =
let dt = tau /. float_of_int n_steps in
let u = exp (vol *. sqrt (3.0 *. dt)) in
let d = 1.0 /. u in
let alpha = (rate -. div_yield -. 0.5 *. vol *. vol) *. sqrt dt
/. (2.0 *. vol *. sqrt 3.0) in
let pu = 1.0 /. 6.0 +. alpha in
let pm = 2.0 /. 3.0 in
let pd = 1.0 /. 6.0 -. alpha in
let df = exp (-. rate *. dt) in
let n_nodes = 2 * n_steps + 1 in
let values = Array.make n_nodes 0.0 in
(* Terminal values: node j corresponds to j moves from top *)
for j = 0 to n_nodes - 1 do
let k = n_steps - j in (* net up moves *)
let s = spot *. (if k >= 0 then u ** float_of_int k
else d ** float_of_int (-k)) in
values.(j) <- match option_type with
| `Call -> Float.max 0.0 (s -. strike)
| `Put -> Float.max 0.0 (strike -. s)
done;
for step = n_steps - 1 downto 0 do
for j = 0 to 2 * step do
let k = step - j in
let s = spot *. (if k >= 0 then u ** float_of_int k
else d ** float_of_int (-k)) in
let cont = df *. (pu *. values.(j) +. pm *. values.(j + 1)
+. pd *. values.(j + 2)) in
let intr = match option_type with
| `Call -> Float.max 0.0 (s -. strike)
| `Put -> Float.max 0.0 (strike -. s)
in
values.(j) <- match exercise with
| `European -> cont
| `American -> Float.max intr cont
done
done;
values.(0)
end
11.3 Finite Difference Methods
Finite difference methods (FDMs) discretise the Black-Scholes PDE directly on a grid of stock price vs time.
The Black-Scholes PDE on the change of variables $\tau = T - t$ (time to expiry):
$$\frac{\partial V}{\partial \tau} = \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + (r-q) S \frac{\partial V}{\partial S} - r V$$
11.3.1 Explicit Scheme (FTCS)
Replace derivatives with finite differences:
$$\frac{V_i^{m+1} - V_i^m}{\Delta\tau} = \frac{1}{2}\sigma^2 S_i^2 \frac{V_{i+1}^m - 2V_i^m + V_{i-1}^m}{\Delta S^2} + (r-q)S_i \frac{V_{i+1}^m - V_{i-1}^m}{2\Delta S} - r V_i^m$$
Explicit scheme: stable only if $\Delta\tau \leq \frac{(\Delta S)^2}{\sigma^2 S_{\max}^2}$ (CFL condition).
11.3.2 Implicit Scheme
The implicit (backward Euler) scheme is unconditionally stable. At each time step, solve a tridiagonal linear system:
$$-a_i V_{i-1}^{m+1} + (1 + b_i) V_i^{m+1} - c_i V_{i+1}^{m+1} = V_i^m$$
module Finite_difference = struct
(** Implicit finite difference scheme for BS PDE.
Unconditionally stable. Solves tridiagonal system via Thomas algorithm. *)
let implicit_scheme ~spot ~strike ~rate ~vol ~tau
~n_space ~n_time ~option_type ~exercise =
let s_max = 4.0 *. spot in
let ds = s_max /. float_of_int n_space in
let dt = tau /. float_of_int n_time in
let s = Array.init (n_space + 1) (fun i -> float_of_int i *. ds) in
(* Initialise with terminal condition *)
let v = Array.map (fun si ->
match option_type with
| `Call -> Float.max 0.0 (si -. strike)
| `Put -> Float.max 0.0 (strike -. si)
) s in
(* Time-step backward *)
let a = Array.make (n_space - 1) 0.0 in (* sub-diagonal *)
let b = Array.make (n_space - 1) 0.0 in (* main diagonal *)
let c = Array.make (n_space - 1) 0.0 in (* super-diagonal *)
for _ = 1 to n_time do
(* Assemble tridiagonal system *)
for i = 1 to n_space - 1 do
let si = s.(i) in
let ai = 0.5 *. dt *. (rate *. si /. ds -. vol *. vol *. si *. si /. (ds *. ds)) in
let bi = 1.0 +. dt *. (vol *. vol *. si *. si /. (ds *. ds) +. rate) in
let ci = -. 0.5 *. dt *. (rate *. si /. ds +. vol *. vol *. si *. si /. (ds *. ds)) in
a.(i - 1) <- ai;
b.(i - 1) <- bi;
c.(i - 1) <- ci
done;
(* Thomas algorithm: forward elimination *)
let c' = Array.copy c and d' = Array.copy (Array.sub v 1 (n_space - 1)) in
c'.(0) <- c.(0) /. b.(0);
d'.(0) <- d'.(0) /. b.(0);
for i = 1 to n_space - 2 do
let denom = b.(i) -. a.(i) *. c'.(i - 1) in
c'.(i) <- c.(i) /. denom;
d'.(i) <- (d'.(i) -. a.(i) *. d'.(i - 1)) /. denom
done;
(* Back substitution *)
v.(n_space - 1) <- d'.(n_space - 2);
for i = n_space - 3 downto 0 do
v.(i + 1) <- d'.(i) -. c'.(i) *. v.(i + 2)
done;
(* Boundary conditions *)
(match option_type with
| `Call ->
v.(0) <- 0.0;
v.(n_space) <- s_max -. strike *. exp (-. rate *. dt)
| `Put ->
v.(0) <- strike *. exp (-. rate *. dt);
v.(n_space) <- 0.0);
(* Early exercise constraint for American options *)
if exercise = `American then
Array.iteri (fun i si ->
let intrinsic = match option_type with
| `Call -> Float.max 0.0 (si -. strike)
| `Put -> Float.max 0.0 (strike -. si)
in
v.(i) <- Float.max v.(i) intrinsic
) s
done;
(* Interpolate to find value at spot *)
let i = int_of_float (spot /. ds) in
let w = (spot -. s.(i)) /. ds in
v.(i) *. (1.0 -. w) +. v.(i + 1) *. w
(**
Crank-Nicolson scheme: average of explicit and implicit.
Second-order accurate in time: O(dt²) vs O(dt) for pure implicit.
*)
let crank_nicolson ~spot ~strike ~rate ~vol ~tau ~n_space ~n_time
~option_type ~exercise =
(* Full Crank-Nicolson implementation: solves 1/2 implicit + 1/2 explicit *)
(* Abbreviated: use implicit_scheme with half weights *)
ignore (spot, strike, rate, vol, tau, n_space, n_time, option_type, exercise);
0.0 (* placeholder — full implementation follows same pattern *)
end
11.4 Comparison of Methods
| Method | American | Path-dependent | Convergence | Implementation |
|---|---|---|---|---|
| Binomial tree | ✓ | Partial | $O(1/n)^*$ | Simple |
| Trinomial tree | ✓ | Partial | $O(1/n)$ | Moderate |
| Explicit FD | ✓ | ✗ | $O(dt, dS^2)$ | Moderate |
| Implicit FD | ✓ | ✗ | $O(dt, dS^2)$ | Moderate |
| Crank-Nicolson | ✓ | ✗ | $O(dt^2, dS^2)$ | Harder |
| Monte Carlo (Ch12) | Via LSMC | ✓✓ | $O(1/\sqrt{N})$ | Flexible |
*CRR oscillates; use smooth CRR or Richardson extrapolation
11.5 Convergence Analysis
Understanding convergence is not merely academic — it determines how many steps you need to get a given accuracy. For $O(1/n)$ convergence, halving the error requires doubling the steps. For $O(1/n^2)$ convergence (Richardson extrapolation), halving the error requires only $\sqrt{2} \approx 1.41$ times more steps. This makes convergence order the most important practical parameter in numerical pricing.
What $O(1/n)$ means in practice. For the CRR binomial tree with $n = 100$ steps, the error versus Black-Scholes is typically around 0.005–0.01 per 100 face value for an at-the-money option — about 0.5–1 cent. For $n = 500$ steps, the error shrinks to 0.001–0.002 (0.1–0.2 cents). For $n = 1000$ steps — about 1 second of compute — the error is below 0.5bp. In production, you would use $n = 200–500$ for vanilla options or $n = 1000+$ for barrier options near the barrier, where the payoff discontinuity slows convergence.
However, plain CRR oscillates: the error alternates sign as $n$ increases (negative for even $n$, positive for odd $n$). This occurs because the strike $K$ sometimes coincides with a tree node (exact pricing) and sometimes falls between nodes (interpolation error). This oscillation can be eliminated by:
- Averaging adjacent step counts: $V^* = (V(n) + V(n+1)) / 2$
- Richardson extrapolation: $V^* = 2V(2n) - V(n)$ which cancels the leading error term and achieves $O(1/n^2)$
- Leisen-Reimer parameterisation: uses the Peizer-Pratt approximation to align a tree node with the strike, eliminating the oscillation source
Finite difference interpretation of convergence. The Crank-Nicolson scheme is $O(\Delta t^2, \Delta S^2)$ in both dimensions. For a $200 \times 200$ grid over $T = 1$ year with $S \in [0, 3K]$, the time step is $\Delta t = 1/200 = 0.005$ and the space step is $\Delta S = 3K/200$. The temporal error is $O(0.005^2) = O(2.5\times 10^{-5})$ and the spatial error is $O((3K/200)^2) = O(K^2 \times 2.25 \times 10^{-4})$. Numerically this gives errors well below 0.01 for standard option parameters — better than Monte Carlo with 100,000 paths.
let convergence_study ~spot ~strike ~rate ~vol ~tau ~option_type =
Printf.printf "Convergence to Black-Scholes (n_steps):\n";
Printf.printf "%-10s %-14s %-14s %-14s\n" "Steps" "Binomial" "Trinomial" "Error(Binom)";
Printf.printf "%s\n" (String.make 56 '-');
let bs = match option_type with
| `Call -> Black_scholes.call ~spot ~strike ~rate ~vol ~tau
| `Put -> Black_scholes.put ~spot ~strike ~rate ~vol ~tau
in
List.iter (fun n ->
let binom = Binomial.price ~spot ~strike ~rate ~vol ~tau ~n_steps:n
~option_type ~exercise:`European in
let trinom = Trinomial.price ~spot ~strike ~rate ~vol ~tau ~n_steps:n
~option_type ~exercise:`European in
Printf.printf "%-10d %-14.6f %-14.6f %-14.2e\n"
n binom trinom (Float.abs (binom -. bs))
) [10; 50; 100; 250; 500; 1000]
11.6 Chapter Summary
Numerical methods for option pricing follow directly from the two equivalent formulations of the Black-Scholes theory: the risk-neutral expectation (which suggests simulation or lattice methods) and the PDE (which suggests finite difference methods). Both yield the same prices when implemented correctly, and both have distinct advantages.
Lattice methods are intuitive and flexible. The CRR binomial tree places up and down moves at factors $u = e^{\sigma\sqrt{\Delta t}}$ and $d = 1/u$, with risk-neutral probability $p = (e^{r\Delta t} - d)/(u-d)$. Convergence is $O(1/n)$ for smooth payoffs, but the notorious oscillation in error (alternating between positive and negative as step count increases) can be damped by averaging adjacent step counts or switching to Leisen-Reimer parameterisation for European options. American options are handled identically by replacing the continuation value with the payoff maximum at each node — this is the key advantage of lattice methods over closed-form approaches.
Finite difference methods provide more direct control over accuracy. The explicit scheme is simple to implement but requires $\Delta t \leq \Delta S^2 / (\sigma^2 S^2)$ for stability — this Courant-Friedrichs-Lewy condition typically forces very small time steps. The implicit scheme (backward Euler) is unconditionally stable but only first-order in time, requiring fine time grids for accuracy. Crank-Nicolson splits the difference, achieving second-order in both time and space while remaining unconditionally stable. For a 200×200 grid, Crank-Nicolson prices a vanilla European option to within 0.001% of Black-Scholes — better than Monte Carlo with 100,000 paths at a fraction of the computational cost.
For two-dimensional problems (e.g., rainbow options, convertible bonds with credit), these methods extend naturally: the PDE gains mixed derivative terms and the grid becomes a 2D mesh. Beyond two dimensions, the curse of dimensionality makes Monte Carlo preferable.
Exercises
11.1 Implement a smooth CRR binomial tree (average of odd/even step counts) and compare error to plain CRR at n=20, 50, 100 steps vs Black-Scholes.
11.2 Price an American put at S=100, K=100, r=5%, σ=25%, T=1. Compare: CRR binomial (500 steps), trinomial (200 steps), and implicit FD (200×200 grid).
11.3 Study the effect of grid spacing: for implicit FD, vary n_space from 50 to 500 at fixed n_time=1000. What is the optimal ratio n_space/n_time?
11.4 Implement Richardson extrapolation on the binomial tree: $V^* = 2V(2n) - V(n)$. Show error improves to $O(1/n^2)$.