Lesson 8.2: Finite differences and PDE methods
Learn to price options by solving the Black-Scholes PDE on a grid: explicit, implicit and Crank-Nicolson schemes, stability, American options and Greeks.
By the end of this lesson you will be able to set up the Black–Scholes equation on a grid with its boundary and initial conditions, write the explicit, implicit and Crank–Nicolson finite-difference schemes, say why the explicit scheme can blow up and what step size it needs, check the order of convergence and use Richardson extrapolation, price an American option with a projection step, and read the Greeks off the grid. You will show that a Crank–Nicolson grid of 400 price nodes and 400 time steps prices the six-month call at 7.1573 against the exact 7.1559 and that the error falls by four when the grid is doubled, that an explicit scheme with 100 price nodes needs at least 200 time steps and 400 nodes need at least 3,200, that Richardson extrapolation of two grids gives 7.15588, and that the grid's delta, gamma and theta agree with the formulas to three or four digits. You should have read Lessons 6.2, 7.2 and 7.4. The mathematics is the Black–Scholes equation, Taylor series and a tridiagonal linear system. This course is education about how markets and research work, not investment advice.
Lesson 7.4 noted that a binomial tree is a crude finite-difference scheme for the Black–Scholes equation. This lesson takes that idea seriously. Instead of simulating random paths, as Monte Carlo does, a finite-difference method solves the pricing equation directly, on a grid of prices and times. For problems with one or two sources of risk, it is usually the most accurate and fastest method, it delivers the Greeks almost for free, and it handles American options naturally. Its weakness is the curse of dimensionality: the grid grows exponentially with the number of risk factors, which is where Monte Carlo takes over.
Setting up the problem
The Black–Scholes equation of Lesson 7.2 is
Vt + ½ σ² S² VSS + r S VS − r V = 0
with a terminal condition at expiry, V(S, T) = payoff(S). It is convenient to turn time around. Let τ = T − t, the time remaining. Then Vτ = −Vt and the equation becomes
Vτ = ½ σ² S² VSS + r S VS − r V
with an initial condition at τ = 0: V(S, 0) = max(S − K, 0) for a call. The solution then evolves forward in τ from the payoff to today, like heat spreading in a rod. We need three more things.
A grid: prices Sj = j × ΔS for j = 0 to M, up to a maximum price Smax, and times τn = n × Δt for n = 0 to N, with Δt = T ÷ N.
Boundary conditions at the edges of the price range. For a call, V(0, τ) = 0, since a share worth nothing will not end above the strike, and V(Smax, τ) = Smax − K e^(−rτ), since a call far in the money is worth the share less the discounted strike. For a put, V(0, τ) = K e^(−rτ) and V(Smax, τ) = 0.
Difference formulas for the derivatives, from Taylor series: VS ≈ (V{j+1} − V{j−1}) ÷ (2ΔS) and VSS ≈ (V{j+1} − 2Vj + V{j−1}) ÷ ΔS². Substituting Sj = jΔS, the ΔS cancels from the coefficients and the equation at node j becomes
dVj/dτ = aj V{j−1} + bj Vj + cj V{j+1}
with aj = ½(σ² j² − r j), bj = −(σ² j² + r) and cj = ½(σ² j² + r j). Everything is a linear combination of three neighbouring nodes, so the whole grid moves forward by a tridiagonal operator. The coefficients involve only the node number j, not the size of the step.
Three ways to step forward in time
The remaining choice is how to advance from time level n to n + 1. Write V^n for the vector of values at all price nodes at time τn, and L for the operator above. The general θ-scheme is
V^{n+1} − V^n = Δt × [ θ × L V^{n+1} + (1 − θ) × L V^n ]
Scheme · θ · What it does · Accuracy · StabilityExplicit (forward Euler) · 0 · Computes V^{n+1} directly from known values · First order in time · Only if Δt is small enough
Implicit (backward Euler) · 1 · Solves a tridiagonal system for V^{n+1} · First order in time · Always stable
Crank–Nicolson · ½ · Averages the two · Second order in time · Always stable
The explicit scheme is the simplest: each new value is a weighted average of three old ones, with no equation to solve. The implicit schemes need a tridiagonal solve at every step, which costs only about M operations with the Thomas algorithm, so they are cheap. Crank–Nicolson is the standard choice, since it is second-order accurate in both space and time.
Stability: why the explicit scheme can blow up
An unstable scheme does not just give an inaccurate answer. Errors grow at every step, and the numbers become astronomically large. For the explicit scheme on the Black–Scholes equation, the coefficient of the middle node is positive only if Δt is small. The condition, at the largest node j = M where the diffusion term is biggest, is approximately
Δt ≤ 1 ÷ (σ² M²)
where M = Smax ÷ ΔS is the number of price nodes. It couples the time step to the square of the number of price nodes. Halving the price spacing means quartering the time step, so the cost grows with the cube of the grid size in price.
We tested it for the six-month call with σ = 20 per cent, so Δt ≤ 25 ÷ M² per year in these units, giving at least N = T σ² M² time steps.
Price nodes M · Time steps N · Stability limit on N · Result100 · 100 · ≥ 200 · Unstable: the largest value on the grid reached 7.3 × 10³³
100 · 200 · ≥ 200 · Stable: price 7.1817
200 · 100 · ≥ 800 · Unstable: the price at 100 came out at 3.2 × 10¹⁹
400 · 1,600 · ≥ 3,200 · Unstable: the result was not a number
400 · 3,200 · ≥ 3,200 · Stable: price 7.1575
The first row shows the danger: the price at S = 100 looked normal, 7.185, because the blow-up was still confined to the nodes near the top boundary, but it was there, and a longer run or a finer grid would have carried it down to the price you read. A check on the largest value on the grid, and on the stability condition, is a basic precaution. For this reason the explicit scheme is rarely used in practice, except for teaching and for problems where the stability limit is mild.
Accuracy: second-order convergence and Richardson extrapolation
The implicit and Crank–Nicolson schemes are stable for any step, so we can choose the grid by accuracy alone. We solved the six-month call with Smax = 300 and the same number of time steps as price nodes, and compared with the exact price of 7.1559.
Grid (price nodes and time steps) · Crank–Nicolson price · Error · Implicit price · Error20 · 7.8190 · +0.6631 · ·
50 · 7.2493 · +0.0934 · ·
100 · 7.1782 · +0.0223 · 7.1711 · +0.0152
200 · 7.1616 · +0.0057 · 7.1580 · +0.0021
400 · 7.1573 · +0.0014 · 7.1555 · −0.0004
800 · 7.1563 · +0.0004 · ·
Each time the grid is doubled, the Crank–Nicolson error falls by about a factor of four: 0.0223, 0.0057, 0.0014, 0.0004. That is second-order convergence: the error is proportional to the square of the grid spacing. It is much better than the first-order convergence of a binomial tree (Lesson 7.4), where doubling the steps only halved the error. Here the implicit scheme happens to be more accurate at the same grid, since Crank–Nicolson is sensitive to the kink in the payoff at the strike, which excites small oscillations. A common remedy, Rannacher smoothing, replaces the first few Crank–Nicolson steps by implicit ones.
A known order of convergence is a gift: it lets you remove most of the error. If the error is proportional to h², then the fine-grid price Vf and the coarse-grid price Vc (with twice the spacing) satisfy Vc − Vexact ≈ 4 × (Vf − Vexact), which gives the Richardson extrapolation
Vexact ≈ Vf + (Vf − Vc) ÷ 3
Using the 200 and 400 grids, the extrapolated price is 7.1573 + (7.1573 − 7.1616) ÷ 3 = 7.15588, against the exact 7.15590: the error is 0.00002, better than the 800-node grid, at a fraction of the cost.
The choice of domain and grid
The maximum price Smax must be far enough from the strike that the boundary condition is accurate. A common choice is three to five times the strike, or the strike times exp(4σ√T), several standard deviations above it. Too small a domain pollutes the price near the strike, and too large wastes nodes.
A uniform grid wastes nodes. The solution changes fastest near the strike, where the payoff has its kink, and slowly far away. A non-uniform grid, with the nodes concentrated near the strike, gives the same accuracy with far fewer.
Change the variable. Working with x = ln S turns the coefficients of the equation into constants: Vτ = ½σ² Vxx + (r − σ²/2) Vx − rV. The equation on this scale is simpler, and its boundaries sit at plus and minus several standard deviations.
Different boundary conditions may be better than a fixed value: for a call, assuming VSS = 0 at Smax, that is, a payoff that is linear far away, is accurate and needs no formula.
American options on a grid
An American option can be exercised at any time, so its value can never fall below its intrinsic value. On the grid this is a simple addition: after solving for each time step, project the solution onto the constraint: Vj = max(Vj, payoffj). The projection can be done in the simplest way, as above, or more accurately with the PSOR method (projected successive over-relaxation) or a penalty term. It is the finite-difference counterpart of the early-exercise check in the tree of Lesson 7.4.
We priced the one-year at-the-money put of Lesson 7.4 (share 100, strike 100, 6 per cent interest, 20 per cent volatility) on a Crank–Nicolson grid of 400 nodes and 400 steps.
· Finite differences · CheckEuropean put · 5.1669 · Black–Scholes 5.1660
American put · 5.7978 · 500-step binomial tree 5.7974
American put, share at 60, strike 100 · 40.00 · Intrinsic value 40: exercise at once
The two methods agree to better than 0.001 for the American put. The grid also gives the early-exercise boundary, the price below which the put should be exercised at each time, as the set of nodes where the projection is active.
The Greeks come with the grid
Because the grid holds the option value at every price, the Greeks come directly from the neighbouring nodes, with no extra runs:
Delta = (V{j+1} − V{j−1}) ÷ (2ΔS).
Gamma = (V{j+1} − 2Vj + V{j−1}) ÷ ΔS².
Theta = −(change in V as τ increases), read from successive time levels.
To get the grid to include the spot exactly, choose Smax and M so that the spot is a node. We used Smax = 400 and 400 nodes, so ΔS = 1 and the spot of 100 is node 100, with 400 time steps:
· Grid · ExactPrice · 7.1524 · 7.1559
Delta · 0.61125 · 0.61135
Gamma · 0.027123 · 0.027103
Theta per year · −8.663 · −8.659
The Greeks agree to three or four significant digits. Gamma, a second derivative, is the least accurate, as always. A tree or a Monte Carlo simulation has to work much harder for the same accuracy, and this is one of the main reasons that finite-difference methods are used in production for single-factor problems, in particular for local-volatility models (Lesson 8.4), where the volatility depends on price and time but the same grid and the same tridiagonal solve apply unchanged.
Beyond one dimension
A problem with two state variables, such as a share and its stochastic volatility, needs a two-dimensional grid with M × M nodes, and the matrix at each step is no longer tridiagonal. Alternating direction implicit (ADI) methods split the step into a sequence of tridiagonal solves, one per dimension, and handle two or three factors efficiently. Beyond three factors, the number of nodes grows exponentially, the curse of dimensionality, and Monte Carlo, whose cost is not exponential in the dimension, becomes the better tool.
· Binomial tree · Finite differences · Monte CarloBest for · Teaching, American options in one factor · One to three factors, American options, fast Greeks · Many factors, path-dependent payoffs
Error · First order · Second order (Crank–Nicolson) · 1 ÷ √N, statistical
American options · Natural · Natural (projection) · Needs regression (Longstaff–Schwartz)
Path dependence · Difficult · Needs extra state variables · Natural
Greeks · Noisy · Free from the grid · Need extra work
Dimension · One or two · Up to about three · Any
The code below builds the grid and solves it with the three schemes, tests the stability limit, prices the American put and reads the Greeks.
python
import numpy as np
from scipy.stats import norm
from scipy.linalg import solvebanded
from math import exp, log, sqrt
S0, K, r, sigma, T = 100.0, 100.0, 0.06, 0.20, 0.5
def bs(S, K, r, sig, T, kind="call"):
d1 = (log(S / K) + (r + 0.5 sig 2) T) / (sig sqrt(T))
d2 = d1 - sig sqrt(T)
if kind == "call":
return S norm.cdf(d1) - K exp(-r T) norm.cdf(d2)
return K exp(-r T) norm.cdf(-d2) - S norm.cdf(-d1)
def fd(M, N, theta=0.5, Smax=300.0, kind="call", american=False, T=T, K=K, r=r, sig=sigma, explicit=False):
"""Solve Vtau = 0.5 sig^2 S^2 VSS + r S VS - r V on M+1 price nodes with N time steps."""
dS, dt = Smax / M, T / N
S = np.linspace(0, Smax, M + 1)
j = np.arange(1, M)
payoff = np.maximum(S - K, 0) if kind == "call" else np.maximum(K - S, 0)
V = payoff.copy()
lo = 0.5 (sig 2 j 2 - r j) # coefficients of V[j-1], V[j], V[j+1]
di = -(sig 2 j 2 + r)
up = 0.5 (sig 2 j 2 + r j)
for n in range(N):
tau = (n + 1) dt
if kind == "call":
bl, bu = 0.0, Smax - K exp(-r tau)
else:
bl, bu = K exp(-r tau), 0.0
if explicit:
new = V.copy()
new[1:-1] = V[1:-1] + dt (lo V[:-2] + di V[1:-1] + up V[2:])
new[0], new[-1] = bl, bu
V = new
continue
rhs = V[1:-1] + (1 - theta) dt (lo V[:-2] + di V[1:-1] + up V[2:])
rhs[0] += theta dt lo[0] bl
rhs[-1] += theta dt up[-1] bu
ab = np.zeros((3, M - 1)) # banded storage of the tridiagonal matrix
ab[0, 1:] = -theta dt up[:-1]
ab[1, :] = 1 - theta dt di
ab[2, :-1] = -theta dt lo[1:]
new = solvebanded((1, 1), ab, rhs)
if american:
new = np.maximum(new, payoff[1:-1]) # projection onto the exercise constraint
V = np.concatenate([[bl], new, [bu]])
return S, V
def price(M, N, kw):
S, V = fd(M, N, kw)
return float(np.interp(S0, S, V))
exact = bs(S0, K, r, sigma, T)
print("exact", round(exact, 5))
print("Crank-Nicolson:")
for M in (20, 50, 100, 200, 400, 800):
p = price(M, M)
print(f" grid {M:>3}: price {p:.5f} | error {p - exact:+.5f}")
print("implicit (theta = 1):")
for M in (100, 200, 400):
p = price(M, M, theta=1.0)
print(f" grid {M:>3}: price {p:.5f} | error {p - exact:+.5f}")
p200, p400 = price(200, 200), price(400, 400)
print("Richardson extrapolation of 200 and 400:", round(p400 + (p400 - p200) / 3, 5))
print("explicit scheme, stability limit N >= T sigma^2 M^2:")
with np.errstate(all="ignore"):
for M, N in ((100, 100), (100, 200), (200, 100), (400, 1600), (400, 3200)):
S, V = fd(M, N, explicit=True)
need = T sigma 2 M 2
print(f" M {M}, N {N} (need at least {need:.0f}): price at 100 = {np.interp(S0, S, V):.5g} | largest value on the grid = {np.nanmax(np.abs(V)):.3g}")
print("American put (one year):")
ps, pk, pr, pT = 100, 100, 0.06, 1.0
print(" European put on the grid", round(price(400, 400, kind="put", T=pT), 4), "| Black-Scholes", round(bs(100, 100, 0.06, 0.2, 1.0, "put"), 4),
"| American put", round(price(400, 400, kind="put", american=True, T=pT), 4))
S, V = fd(400, 400, kind="put", american=True, T=pT)
print(" American put, share at 60, strike 100:", round(float(np.interp(60, S, V)), 3))
Greeks from a grid in which the spot is a node: Smax = 400, 400 nodes, 400 steps
S, V = fd(400, 400, Smax=400.0)
dS, i = S[1] - S[0], 100
delta = (V[i + 1] - V[i - 1]) / (2 dS)
gamma = (V[i + 1] - 2 V[i] + V[i - 1]) / dS 2
, V1 = fd(400, 400, Smax=400.0, T=T + 0.0005)
, V2 = fd(400, 400, Smax=400.0, T=T - 0.0005)
thetayear = -(V1[i] - V2[i]) / 0.001
d1 = (r + 0.5 sigma 2) T / (sigma sqrt(T))
print("grid price", round(V[i], 4), "delta", round(delta, 5), "gamma", round(gamma, 6), "theta", round(thetayear, 3))
print("exact ", round(exact, 4), " ", round(norm.cdf(d1), 5), " ", round(norm.pdf(d1) / (S0 sigma sqrt(T)), 6), " -8.659")
Common mistakes
Using the explicit scheme without checking stability. The condition Δt ≤ 1 ÷ (σ²M²) links the steps, and a blow-up can hide near the boundary.
Choosing a domain that is too small. The boundary condition then contaminates the price at the spot.
Forgetting the boundary conditions' form for puts and calls. They are different, and an error there shifts the whole solution.
Not checking convergence. Run two or three grids and look at the pattern of errors.
Reading Greeks at a node that is not the spot. Interpolate, or choose the grid so the spot is a node.
Applying Crank–Nicolson to a payoff with a kink without smoothing. Small oscillations appear in the Greeks, and a few implicit steps cure them.
Using a uniform grid where a non-uniform grid would be much better. Concentrate nodes near the strike.
Applying the American projection in the wrong place. It must be applied after each time step, to the solved values.
Using finite differences for many factors. The grid grows exponentially, so use Monte Carlo.
Key takeaways
Turning time round (τ = T − t) makes the Black–Scholes equation Vτ = ½σ²S²VSS + rSVS − rV, with the payoff as the initial condition and boundary conditions at S = 0 and Smax. On the grid, each node depends on its two neighbours, giving a tridiagonal system.
The θ-scheme covers the explicit (θ = 0), implicit (θ = 1) and Crank–Nicolson (θ = ½) methods. Implicit and Crank–Nicolson are stable for any step and need a cheap tridiagonal solve.
The explicit scheme is stable only if Δt ≤ 1 ÷ (σ²M²): 100 price nodes need at least 200 steps for six months at 20 per cent volatility and 400 nodes need 3,200. Violating it gave values of 7.3 × 10³³ and 3.2 × 10¹⁹.
Crank–Nicolson converges at second order: errors of 0.0223, 0.0057, 0.0014 and 0.0004 as the grid doubles from 100 to 800. Richardson extrapolation of the 200 and 400 grids gave 7.15588 against the exact 7.15590.
American options need only a projection, V = max(V, payoff), after each step: the one-year put was 5.7978 against 5.7974 from a tree, and a deep in-the-money put was worth its intrinsic value of 40.
The Greeks come from neighbouring nodes: delta 0.61125, gamma 0.027123 and theta −8.663 against the exact 0.61135, 0.027103 and −8.659.
Finite differences are best for one to three factors and for local-volatility models. Beyond that, the grid grows exponentially and Monte Carlo is the tool.
Check yourself
1. A grid has 250 price nodes. What is the fewest time steps an explicit scheme needs for a one-year option with 25 per cent volatility?
N ≥ T σ² M² = 1 × 0.0625 × 62,500 = 3,906 steps.
2. Crank–Nicolson prices an option at 5.482 on a grid of 100 and 5.446 on a grid of 200. What is the extrapolated price?
Vf + (Vf − Vc) ÷ 3 = 5.446 + (5.446 − 5.482) ÷ 3 = 5.446 − 0.012 = 5.434.
3. Why does the American projection step produce the early-exercise boundary?
The nodes where the solved value falls below the intrinsic value are replaced by the intrinsic value, so they are the prices at which exercising is optimal. The edge of that set at each time is the boundary.
4. A call on a grid with a maximum price of 150 and a strike of 100 is priced for a share at 100 with σ = 40 per cent and a year to go. Is that a sensible domain?
No. The one-standard-deviation range of the log price is 0.4, so a price of 150 is only about one standard deviation above the spot. The boundary condition is inaccurate there. A maximum of about 400 or more is needed.
5. A colleague says doubling the number of nodes in a Crank–Nicolson run should halve the error, as in a tree. Is that right?
No. Crank–Nicolson is second order in the grid spacing, so doubling the grid reduces the error by about a factor of four, and the extrapolation formula uses that.
Try this
Set aside forty-five minutes. You need Python and the code in this lesson.
1. Run the code. Check the convergence tables, the explicit stability results, the American puts and the Greeks.
2. Plot the error of Crank–Nicolson against the grid size on a log-log scale. Is the slope about −2?
3. Change Smax from 300 to 200 and to 600 with the grid at 400. How does the price at 100 change?
4. Replace the first four Crank–Nicolson steps by implicit steps (Rannacher smoothing) and see how the gamma near the strike changes.
5. Price an American put with a dividend yield of 3 per cent by lowering the drift in the equation to r − q. How does the early-exercise premium change?
6. Use the stability simulator to find, for your own volatility and grid, the time step an explicit scheme requires, and compare the cost with a Crank–Nicolson run at the same accuracy.
7. Write a short note: why you would use finite differences rather than Monte Carlo for an American put on one share, and rather than a tree.
A good result is a convergence table with errors falling by four per doubling, an American price that matches a tree and Greeks that match the formulas. The next lesson turns from equities to interest rates: bond mathematics, yield curves and models of the short rate.
Read more
Daniel Duffy's Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach is a thorough guide for practitioners. Paul Wilmott, Sam Howison and Jeff Dewynne's The Mathematics of Financial Derivatives explains the equation and the numerical schemes at a gentle pace, and Yves Achdou and Olivier Pironneau's Computational Methods for Option Pricing is a more mathematical treatment. John Crank and Phyllis Nicolson's 1947 paper "A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat-Conduction Type" is the origin of the scheme. R. Rannacher's 1984 paper on finite element solution of diffusion problems with irregular data explains the smoothing, and Tavella and Randall's Pricing Financial Instruments: The Finite Difference Method gives non-uniform grids and many examples.
Tags: Quant Trader Finite differences PDE methods Crank-Nicolson American options