Finite-difference methods

Derivatives on a grid

A derivative is a limit, but a computer works with finite steps. A finite difference replaces f'(x) by the slope over a small step h: forward (f(x + h) - f(x)) / h, or central (f(x + h) - f(x - h)) / (2h), which is more accurate because the errors on the two sides cancel (its error shrinks like h**2, the forward one only like h). The same idea gives the second derivative, the curvature: (f(x + h) - 2 f(x) + f(x - h)) / h**2. For an option price as a function of the spot, that is the gamma.

import math

def forward(f, x, h):
    return (f(x + h) - f(x)) / h

def central(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

exact = math.cos(1.0)                               # the derivative of sin at 1
for h in (0.1, 0.01, 0.001):
    print(h, abs(forward(math.sin, 1.0, h) - exact), abs(central(math.sin, 1.0, h) - exact))

# Gamma of a call from the second difference of its Black-Scholes price
import math

def norm_cdf(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def bs_call(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))

def second_difference(f, x, h):
    return (f(x + h) - 2 * f(x) + f(x - h)) / h ** 2

price = lambda s: bs_call(s, 100, 0.05, 0.2, 1)
print(round(second_difference(price, 100, 0.5), 6))     # 0.018761
print(round(second_difference(price, 100, 0.01), 6))    # 0.018762

From a derivative to an equation

The point of the method is to solve partial differential equations. Black–Scholes is one: the option price changes in time according to its curvature in the spot. After a change of variables it becomes the heat equation, du/dt = D * d2u/dx2: heat spreads out, and the more curved the temperature profile, the faster it moves. Put the values on a grid, replace both derivatives by differences, and the equation becomes a rule to update every point from its neighbours.

# The explicit scheme: a spike of heat spreads to its neighbours, step by step
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

grid = np.zeros(9)
grid[4] = 1.0
for step in range(4):
    print(step, np.round(grid, 3))
    grid = explicit_step(grid, 0.25)

In code, alpha = D * dt / dx**2 ties the time step dt and the grid spacing dx together. One update, new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2]), moves every interior point at once with slices, and the two end points are the boundary conditions, held fixed here. The scheme is simple, and so is the catch.

Stability: the catch

The explicit scheme is stable only if alpha <= 0.5, that is D * dt / dx**2 <= 0.5. Beyond that, rounding errors are amplified at every step and the numbers explode in a sawtooth pattern. The condition is harsh: halving the grid spacing forces a time step four times smaller, so the cost grows fast. This is why professional code often uses implicit schemes, which solve a small linear system at each step and are stable for any dt; the explicit one is the simplest to write and the easiest to understand.

# The same spike with alpha = 0.4 (stable) and alpha = 0.6 (unstable) after 60 steps
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

for alpha in (0.4, 0.6):
    grid = np.zeros(21)
    grid[10] = 1.0
    for _ in range(60):
        grid = explicit_step(grid, alpha)
    print(alpha, "largest value", float(np.abs(grid).max()))       # 0.4: 0.0556, 0.6: about 3.1e7

• A finite difference replaces a derivative by values on a grid: central (f(x + h) - f(x - h)) / (2h) has error about h**2; the second difference (f(x + h) - 2 f(x) + f(x - h)) / h**2 gives curvature (gamma).
• PDEs become update rules on a grid: the explicit step of the heat equation is values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1]), written with slices.
• The explicit scheme needs D * dt / dx**2 <= 0.5 to stay stable; implicit schemes remove that limit at the price of solving a linear system per step.

Next up: covariance and correlation matrices, the first step of portfolio optimisation.

Lab · what's graded

• Estimate a curvature with second_difference25 pts
• Advance the grid one step with explicit_step25 pts
• Check the time step with is_stable25 pts
• March forward in time with solve_heat25 pts