Finite Difference Methods: Solving PDEs on a Grid
Chapter 77 is where we stop talking about pricing theory and start actually building a pricer. This is the chapter where the Black-Scholes equation stops being a formula on paper and becomes code on a screen. If you have ever used the binomial model, congratulations, you already know the basic idea. Finite differences are just a more powerful, more flexible version of the same concept.
Wilmott is clear about his preference. He uses finite-difference methods about 75% of the time in practice. Monte Carlo gets 20%. Explicit formulas get the rest, and those are almost always just plain Black-Scholes for calls and puts. The binomial method? He has used it seriously exactly once, and that was more for modeling insight than numerical accuracy.
The Grid
In the binomial model, you have a tree that branches out from the current stock price. In finite differences, you replace that tree with a rectangular grid. One axis is the stock price S, the other is time t. The grid has evenly spaced nodes in both directions.
The time step is delta_t. The asset step is delta_S. Your grid covers stock prices from 0 up to some maximum (usually three or four times the strike price) and time from now to expiry. Each node on the grid will eventually hold the option value at that stock price and time.
Why a regular grid instead of a tree? Because going from a differential equation to a difference equation is much easier on a regular grid. And because there are many ways to improve finite-difference methods that the binomial model simply cannot match.
Approximating the Derivatives
The Black-Scholes equation has three types of derivatives: theta (time derivative), delta (first S derivative), and gamma (second S derivative). We need to approximate each one using the grid values.
Theta
The time derivative is the simplest. If you know the option value at time t and at time t - delta_t, then:
$$\theta \approx \frac{V_i^{k+1} - V_i^k}{\delta t}$$
The superscript k is the time index (k increases as we go backward from expiry). This is accurate to O(delta_t). Not great, but sufficient for now.
Delta
For the first S derivative, you have three choices:
Forward difference: Use the current point and the one above it. Error: O(delta_S).
Backward difference: Use the current point and the one below it. Error: O(delta_S).
Central difference: Use the points above and below:
$$\Delta \approx \frac{V_{i+1}^k - V_{i-1}^k}{2\delta S}$$
Error: O(delta_S^2). Much better. The symmetry causes the first-order error terms to cancel. Central difference is almost always the right choice.
Sometimes at the boundaries of your grid (i=0 or i=I), you cannot use a central difference because you are missing a neighbor. In that case, you can build a three-point one-sided difference that still achieves O(delta_S^2) accuracy:
$$\Delta \approx \frac{-3V_i^k + 4V_{i+1}^k - V_{i+2}^k}{2\delta S}$$
This uses points on one side only but gets the same order of accuracy as the central difference.
Gamma
The second derivative is estimated by taking the difference of two deltas. Use the forward difference as an estimate at S + delta_S/2, the backward difference at S - delta_S/2, and take their difference:
$$\Gamma \approx \frac{V_{i+1}^k - 2V_i^k + V_{i-1}^k}{\delta S^2}$$
Also accurate to O(delta_S^2). This is the standard gamma approximation and you will use it everywhere.
Final Conditions (Where You Start)
In finance, PDEs are solved backward in time. You start at expiry and work toward the present. The “final condition” is the payoff function:
$$V_i^0 = \text{Payoff}(S_i)$$
For a call: max(S_i - E, 0). For a put: max(E - S_i, 0). This gets the scheme started.
Boundary Conditions (The Edges)
At the edges of your grid (S = 0 and S = S_max), you need boundary conditions. These depend on the contract.
At S = 0: For a call, the value is always zero. For a put, V(0,t) = E * exp(-r(T-t)). More generally, since both drift and diffusion switch off at S = 0, the payoff is guaranteed: V(0, t) = Payoff(0) * exp(-r(T-t)).
At S = S_max: For a put, the value is zero for large S. For a call, it approaches S - E*exp(-r(T-t)). The most useful general boundary condition is:
$$\frac{\partial^2 V}{\partial S^2} = 0 \quad \text{(gamma equals zero)}$$
In finite-difference form: V_I = 2*V_{I-1} - V_{I-2}. This works for almost any contract whose payoff is at most linear in S for large S. It is independent of the specific contract, which means you do not need to change your code when you switch from pricing calls to puts.
Barriers: If the option has a barrier at S_u, set V = 0 (or the rebate value) at that level. Ideally, make your grid spacing so the barrier falls exactly on a grid point. If it does not, use a fictitious point and linear interpolation to maintain O(delta_S^2) accuracy. Never just set V to zero at the nearest grid point. That kills your accuracy.
The Explicit Method
Now we put it all together. Substitute the finite-difference approximations into the Black-Scholes equation:
$$V_i^{k+1} = A_i \cdot V_{i-1}^k + B_i \cdot V_i^k + C_i \cdot V_{i+1}^k$$
where A, B, and C are coefficients that depend on the volatility, interest rate, dividend yield, and the grid spacing. This is the explicit method. You know everything on the right side (from the current time step), so you can directly calculate the left side (the next time step backward).
In code, the core loop looks like this (pseudocode):
For k = 1 to number_of_time_steps:
For i = 1 to I-1:
V_new[i] = A[i]*V_old[i-1] + B[i]*V_old[i] + C[i]*V_old[i+1]
Apply boundary conditions at i=0 and i=I
V_old = V_new
That is it. The entire explicit finite-difference pricer is a nested loop with a few lines per iteration.
For the Black-Scholes equation specifically, the coefficients work out to:
$$A_i = \delta t \left( \frac{1}{2}\sigma^2 i^2 - \frac{1}{2}(r-D)i \right)$$ $$B_i = 1 - \delta t \left( \sigma^2 i^2 + r \right)$$ $$C_i = \delta t \left( \frac{1}{2}\sigma^2 i^2 + \frac{1}{2}(r-D)i \right)$$
where i is the grid index (so S = i * delta_S).
Stability: The Catch
The explicit method is simple, but it has a serious limitation. It only works if the time step is small enough relative to the asset step. The stability condition is approximately:
$$\delta t < \frac{\delta S^2}{2a}$$
where a is the diffusion coefficient. For Black-Scholes, this becomes:
$$\delta t < \frac{1}{\sigma^2 I^2}$$
where I is the number of asset steps. This is a severe constraint. If you want better accuracy by halving delta_S (doubling I), you must quarter delta_t. The computation time goes up by a factor of eight, but accuracy only improves by a factor of four.
There is also a constraint on the first-derivative term that requires the “convection” to be small relative to the “diffusion.” For Black-Scholes this is usually not a problem unless volatility is very small. For interest rate models where volatility can be tiny, this constraint can bite.
The good news: if you violate the stability condition, you will know immediately. The instability produces wild oscillations that are impossible to miss. You will never get a quietly wrong answer from an unstable explicit scheme.
Convergence
The global error in the explicit method is O(delta_t, delta_S^2). This comes from accumulating the local truncation error O(delta_t^2, delta_t * delta_S^2) over O(1/delta_t) time steps.
This is the same accuracy you get from the binomial model. But the finite-difference framework gives you a clear path to doing better: use an implicit method (next chapter) and the error improves to O(delta_t^2, delta_S^2) with no stability constraint. The binomial method does not offer this upgrade path.
The Takeaway
The explicit finite-difference method is the simplest way to solve the Black-Scholes equation (and many other PDEs) numerically. The concept is: approximate derivatives using grid values, substitute into the equation, and march backward in time.
You need three things:
- A final condition (the payoff)
- Boundary conditions (what happens at the edges)
- The time-stepping formula (the A, B, C coefficients)
The method is easy to code, easy to debug, and produces correct results as long as you respect the stability constraint. It handles different contracts by just changing the final and boundary conditions. It handles different models by changing the coefficients.
The main limitation is speed. The stability constraint forces small time steps, which makes computation expensive for high accuracy. This is where implicit methods come in, which we will cover next.
If you are new to numerical methods in finance, start here. Code the explicit method for a European call. Check it against the Black-Scholes formula. Then try a put. Then try an American option (add three lines). Then try a barrier. Each step teaches you something new about how the math becomes code.
Previous post: Numerical Methods: How Computers Actually Price Derivatives
Next post: Advanced Finite Differences: Implicit, Crank-Nicolson, and More