Finite-Difference Methods

Every method in this module so far has approached pricing sideways — a tree that mimics the stock's randomness, a simulation that samples it. There is a more direct route: the option's value V(S,t) already satisfies a known equation, the Black–Scholes PDE. Instead of solving it with pen and paper — only possible for the tidiest payoffs — put it on a grid of stock prices and times, replace every derivative with a difference of neighbouring grid values, and march the whole price surface backward from expiry to today. This is finite differences: the numerical-methods workhorse of computational physics, aimed squarely at the equation this course has been building around all along. We won't grind through the full derivation here — the general recipe for turning derivatives into differences is exactly the same one used to solve the heat equation — but we will get concrete enough to run it, and see clearly when it beats the alternatives and when it doesn't.

The grid, and the explicit update rule

Chop the stock-price axis into M steps of width \Delta S (from 0 up to some S_{\max} comfortably above anything the option cares about) and time into N steps of width \Delta t, working backward from expiry (where the payoff is known exactly) toward today. Write f_{i,j} for the option's value at time step i and stock-price node j (so S = j\,\Delta S). Substituting central/forward differences for the PDE's derivatives and solving for the values one step earlier in time gives the classic explicit update:

The update. f_{i,j} = a_j\,f_{i+1,j-1} + b_j\,f_{i+1,j} + c_j\,f_{i+1,j+1}, with coefficients depending only on the node j and the model parameters, a_j = \tfrac12\Delta t\bigl(\sigma^2 j^2 - rj\bigr), \quad b_j = 1 - \Delta t\bigl(\sigma^2 j^2 + r\bigr), \quad c_j = \tfrac12\Delta t\bigl(\sigma^2 j^2 + rj\bigr).
Boundary rows. At S = 0 and S = S_{\max} the PDE degenerates, so those two edges of the grid are set directly from the payoff's known behaviour there (e.g. a put is worth Ke^{-r\tau} at S=0 and essentially 0 at S_{\max}), rather than from the update rule.
It is conditionally stable. Exactly as with the heat equation, too large a \Delta t relative to \Delta S^2 makes the scheme blow up into oscillations instead of converging — the same family of stability limits, now dressed in the Black–Scholes PDE's coefficients.

Read the update rule as a weighted blend, just like the heat-equation stencil: each new value is built from the three neighbouring values one time-step later (closer to expiry), weighted by a_j, b_j, c_j. That is the entire engine — everything else is bookkeeping.

Early exercise is one more line, again

This should look familiar. Exactly as on a binomial tree, pricing an American option on the grid needs only one extra step: after computing each f_{i,j} from the update rule, clamp it against the intrinsic value,

f_{i,j} \leftarrow \max\bigl(f_{i,j},\ \text{intrinsic}(j\,\Delta S)\bigr),

applied at every node, every time step — the grid's version of the tree's \max(\text{exercise},\ \text{continue}) rule. Run the same coefficients, same boundaries, same backward march, with this one clamp added, and out comes an American price on the very same grid used for the European one. Below, both are computed with the classic textbook example — a stock at S_0=50, struck at K=50, with r=10\%, \sigma=40\%, and T = 5 months to expiry:

// Explicit finite-difference pricer for a put, on the Black-Scholes PDE grid, // European and American, following the classic textbook coefficients. function erf(x: number): number { const sign = x < 0 ? -1 : 1; x = Math.abs(x); const a1 = 0.254829592, a2 = -0.284496736, a3 = 1.421413741; const a4 = -1.453152027, a5 = 1.061405429, p = 0.3275911; const t = 1 / (1 + p * x); const y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-x * x); return sign * y; } function normCdf(x: number): number { return 0.5 * (1 + erf(x / Math.sqrt(2))); } function bsPut(S0: number, K: number, r: number, sig: number, T: number): number { const d1 = (Math.log(S0 / K) + (r + 0.5 * sig * sig) * T) / (sig * Math.sqrt(T)); const d2 = d1 - sig * Math.sqrt(T); return K * Math.exp(-r * T) * normCdf(-d2) - S0 * normCdf(-d1); } function fdPut( S0: number, K: number, r: number, sig: number, T: number, Smax: number, M: number, N: number, american: boolean, ): number { const dS = Smax / M; const dt = T / N; let f: number[] = new Array(M + 1); for (let j = 0; j <= M; j++) f[j] = Math.max(K - j * dS, 0); // payoff at expiry const a: number[] = [], b: number[] = [], c: number[] = []; for (let j = 0; j <= M; j++) { a[j] = 0.5 * dt * (sig * sig * j * j - r * j); b[j] = 1 - dt * (sig * sig * j * j + r); c[j] = 0.5 * dt * (sig * sig * j * j + r * j); } for (let i = N - 1; i >= 0; i--) { const fNew: number[] = new Array(M + 1); for (let j = 1; j < M; j++) { fNew[j] = a[j] * f[j - 1] + b[j] * f[j] + c[j] * f[j + 1]; if (american) fNew[j] = Math.max(fNew[j], K - j * dS); // clamp to intrinsic } fNew[0] = K; // boundary: S = 0, a put is worth K (ignoring discounting at this resolution) fNew[M] = 0; // boundary: S = Smax, the put is worthless f = fNew; } const j0 = S0 / dS; // read off the grid at S0 (falls exactly on a node here) return f[Math.round(j0)]; } const S0 = 50, K = 50, r = 0.10, sig = 0.40, T = 5 / 12; console.log(`Black-Scholes European put (closed form) = ${bsPut(S0, K, r, sig, T).toFixed(4)}`); const Smax = 100, M = 20, N = 100; console.log(`Finite-difference European put = ${fdPut(S0, K, r, sig, T, Smax, M, N, false).toFixed(4)}`); console.log(`Finite-difference American put = ${fdPut(S0, K, r, sig, T, Smax, M, N, true).toFixed(4)}`);

The finite-difference European price lands close to the exact Black–Scholes value — the residual gap is pure discretization error from the coarse 20\times100 grid, shrinking as M and N grow. And the American price comes out higher than both, by exactly the kind of early-exercise premium the tree found in the previous page — the same economic fact, arrived at from the PDE side instead of the tree side.

Explicit vs implicit — and finite differences vs Monte Carlo

The explicit scheme above is cheap per step but conditionally stable, exactly like the heat equation's FTCS scheme. Evaluate the right-hand side at the new (unknown) time level instead, and you get the implicit scheme — unconditionally stable at the cost of solving a (tridiagonal) linear system at every time step — or its accuracy-boosting cousin, Crank–Nicolson, which averages the explicit and implicit updates and is the industry-standard choice in practice. None of that changes the picture painted above: whichever scheme fills in the grid, the method's real strength and real limitation both come from the grid itself.

Method	Many correlated assets	Early exercise	Greeks
Binomial tree	Poor — nodes explode per extra asset	Native (max at each node)	Rough, from nearby nodes
Finite differences	Poor — grid cells explode per dimension	Native (clamp at each node)	Direct — read off neighbouring grid values
Monte Carlo	Excellent — cost barely changes with dimension	Awkward — needs extra machinery (e.g. regression-based methods)	Needs extra work (bump-and-reprice, or pathwise estimators)

Finite differences win precisely where Monte Carlo struggles: American-style early exercise falls out of the same backward march used for the European price, and because the grid holds V(S,t) at every node simultaneously, the Greeks — Delta, Gamma — are just finite differences of neighbouring grid values, essentially free. Monte Carlo wins precisely where finite differences struggle: a basket option on ten correlated stocks needs a ten-dimensional grid, and a grid with even a modest 50 points per axis in ten dimensions has 50^{10} cells — the same curse of dimensionality that haunts every grid method, exactly as it does in Monte Carlo Methods. Simulation's error stays at a flat 1/\sqrt N no matter how many underlyings the basket option has — which is exactly why the desk keeps both tools on hand, and picks the one that fits the payoff in front of them.

Once the backward march reaches today's row, you don't just have V(S_0, 0) — you have V at every nearby stock price on that same row, because the grid computed them all simultaneously as a side effect of filling in the surface. Delta and Gamma are derivatives of V with respect to S — so just take finite differences of the neighbouring grid values you already have, no extra pricing runs required:

\Delta \approx \frac{f_{0,j+1} - f_{0,j-1}}{2\,\Delta S}, \qquad \Gamma \approx \frac{f_{0,j+1} - 2f_{0,j} + f_{0,j-1}}{\Delta S^2}.

Compare that with Monte Carlo, where getting a Greek generally means re-running the whole simulation with a slightly bumped input and differencing two noisy price estimates against each other — expensive, and the noise in each price estimate makes the noise in the difference even worse. This is the concrete version of the table above: a grid method doesn't just price the option, it happens to compute its whole neighbourhood of sensitivities along the way.

The explicit scheme's stability condition doesn't vanish just because the equation changed. Refine \Delta S without also shrinking \Delta t by enough, and the explicit Black–Scholes scheme explodes exactly the way the explicit heat-equation scheme does — the numbers are different, but the failure mode (growing, alternating-sign oscillations from amplifying the grid's roughest mode) is the same disease with the same cure: a smaller \Delta t, or switch to an implicit scheme.
The early-exercise clamp must be applied at every node, every time step — the same mistake as on the tree (checking it only near expiry) silently throws away most of the early-exercise premium.
A grid is not a free lunch in high dimensions. The temptation to "just add one more stock-price axis" for a two-asset payoff is reasonable — a 2-D grid is still very tractable. Try it for five or ten correlated assets and the cell count becomes astronomical long before you've written the second axis; that is precisely the point at which the desk reaches for Monte Carlo instead, not because finite differences got harder to code, but because they became computationally impossible.