Finite-Differenzen-Verfahren

Ableitungen auf einem Gitter

Eine Ableitung ist ein Grenzwert, ein Computer arbeitet aber mit endlichen Schritten. Eine finite Differenz ersetzt f'(x) durch die Steigung über einen kleinen Schritt h: vorwärts (f(x + h) - f(x)) / h oder zentral (f(x + h) - f(x - h)) / (2h), genauer, weil sich die Fehler beider Seiten aufheben (der Fehler schrumpft wie h**2, beim Vorwärtsquotienten nur wie h). Dieselbe Idee liefert die zweite Ableitung, die Krümmung: (f(x + h) - 2 f(x) + f(x - h)) / h**2. Für einen Optionspreis als Funktion des Kassakurses ist das das 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

Von der Ableitung zur Gleichung

Der Zweck des Verfahrens ist das Lösen partieller Differentialgleichungen. Black-Scholes ist eine: Der Optionspreis ändert sich in der Zeit gemäß seiner Krümmung im Kassakurs. Nach einer Variablentransformation wird daraus die Wärmeleitungsgleichung du/dt = D * d2u/dx2: Wärme breitet sich aus, und je stärker gekrümmt das Temperaturprofil, desto schneller bewegt es sich. Man legt die Werte auf ein Gitter, ersetzt beide Ableitungen durch Differenzen, und die Gleichung wird zur Regel, jeden Punkt aus seinen Nachbarn zu aktualisieren.

# 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)

Im Code verbindet alpha = D * dt / dx**2 den Zeitschritt dt und den Gitterabstand dx. Ein Update, new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2]), bewegt mit Slices alle inneren Punkte auf einmal; die beiden Enden sind die Randbedingungen, hier fest. Das Schema ist einfach, und der Haken ebenso.

Stabilität: der Haken

Das explizite Schema ist nur stabil, wenn alpha <= 0.5, also D * dt / dx**2 <= 0.5. Darüber werden Rundungsfehler in jedem Schritt verstärkt, und die Zahlen explodieren im Sägezahnmuster. Die Bedingung ist hart: Halbiert man den Gitterabstand, muss der Zeitschritt viermal kleiner werden, die Kosten wachsen schnell. Deshalb nutzt professioneller Code oft implizite Schemata, die in jedem Schritt ein kleines lineares System lösen und für jedes dt stabil sind; das explizite ist am einfachsten zu schreiben und zu verstehen.

# 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

• Eine finite Differenz ersetzt eine Ableitung durch Werte auf einem Gitter: die zentrale (f(x + h) - f(x - h)) / (2h) hat Fehler etwa h**2; die zweite Differenz (f(x + h) - 2 f(x) + f(x - h)) / h**2 liefert die Krümmung (Gamma).
• PDEs werden zu Update-Regeln auf einem Gitter: Der explizite Schritt der Wärmeleitungsgleichung ist values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1]), mit Slices geschrieben.
• Das explizite Schema braucht D * dt / dx**2 <= 0.5, um stabil zu bleiben; implizite Schemata heben die Grenze auf, kosten aber ein lineares System pro Schritt.

Als Nächstes: Kovarianz- und Korrelationsmatrizen, der erste Schritt der Portfoliooptimierung.

Lab: Was bewertet wird

• Mit second_difference eine Krümmung schätzen25 Pkt.
• Das Gitter mit explicit_step um einen Schritt weiterführen25 Pkt.
• Mit is_stable den Zeitschritt prüfen25 Pkt.
• Mit solve_heat in der Zeit vorwärtsrechnen25 Pkt.