Finite-difference methods for the BS PDE

Finite-difference methods solve PDEs numerically by approximating derivatives with differences on a grid. In finance: solving the Black-Scholes PDE for vanilla and exotic options, American-option pricing via PSOR, multi-factor PDEs for hybrid products. FD methods complement Monte Carlo: better for low-dimensional path-dependent options, worse for high-dimensional basket options.

The Black-Scholes PDE — discretise

Discretise (S, t) on a rectangular grid: S_i = i · ΔS for i = 0, 1, ..., I; t_j = j · Δt for j = 0, 1, ..., J. Let V_{i,j} ≈ V(S_i, t_j). Replace derivatives by finite differences.

Finite-difference approximations

∂V/∂t ≈ (V_{i,j+1} - V_{i,j}) / Δt    (forward)
∂V/∂S ≈ (V_{i+1,j} - V_{i-1,j}) / (2 ΔS)    (central)
∂²V/∂S² ≈ (V_{i+1,j} - 2 V_{i,j} + V_{i-1,j}) / (ΔS)²

Central differences are O(Δx²) accurate; forward/backward are O(Δx). Always use central where possible.

Explicit scheme

Use forward time-difference and central space-differences:

V_{i,j+1} = V_{i,j} + Δt [r S_i (V_{i+1,j} - V_{i-1,j})/(2ΔS) + (σ²S_i²/2)(V_{i+1,j} - 2V_{i,j} + V_{i-1,j})/(ΔS)² - r V_{i,j}]

Trivially easy: explicit forward-propagation. Stability requires Δt small relative to (ΔS)² — the CFL-style condition. With reasonable grids, this often forces tiny time steps and slow simulation.

Implicit scheme

Use backward time-difference: V_{i,j+1} - V_{i,j} = Δt · L V_{i,j+1}, where L is the spatial differential operator. Now V_{i,j+1} appears on both sides — must solve a linear system at each time step.

(I - Δt L) V_{j+1} = V_j

Unconditionally stable — Δt can be large. Costs more per step (tridiagonal solve) but takes far fewer steps.

Crank-Nicolson — the production default

Average the explicit and implicit schemes. Second-order accurate in both time and space. Unconditionally stable. The default for BS-PDE solving in production.

(I - Δt/2 · L) V_{j+1} = (I + Δt/2 · L) V_j

Boundary conditions

• Terminal: V_{i,J} = payoff(S_i) — the option's payoff at expiry.
• Lower boundary (S = 0): V(0, t) = K · exp(-r(T-t)) for puts, 0 for calls.
• Upper boundary (S = S_max): V(S_max, t) ≈ S_max - K · exp(-r(T-t)) for calls, 0 for puts. Choose S_max large enough that this asymptote is accurate.

Stretched grids

Concentrate grid points near the strike where gamma is large. Logarithmic or sinh-stretched grids can be 2-5× more accurate than uniform grids for the same number of nodes.

American options via PSOR

American options can be exercised early. The PDE becomes a linear complementarity problem (LCP):

max(L V, payoff - V) = 0

Either the PDE is satisfied (no exercise) or the option equals its intrinsic value (exercise optimal). Solved by Projected Successive Over-Relaxation (PSOR): take a tentative implicit step, project values up to the payoff if they fall below. Iterate to convergence. Modern alternatives: penalty methods, the LCP formulation solved by Brennan-Schwartz tridiagonal.

Multi-factor PDEs

For two state variables (e.g., stock and stochastic vol in Heston), the PDE is 2D and requires alternating-direction implicit (ADI) methods. ADI splits each step into two: implicit in one variable, explicit in the other, alternating. Standard in production for two-factor models. Three factors and more: PDE becomes very expensive; switch to Monte Carlo.