Métodos de diferencias finitas

Derivadas en una malla

Una derivada es un límite, pero un ordenador trabaja con pasos finitos. Una diferencia finita sustituye f'(x) por la pendiente sobre un paso pequeño h: hacia delante (f(x + h) - f(x)) / h, o central (f(x + h) - f(x - h)) / (2h), más precisa porque los errores de los dos lados se compensan (su error baja como h**2, el de la progresiva solo como h). La misma idea da la segunda derivada, la curvatura: (f(x + h) - 2 f(x) + f(x - h)) / h**2. Para el precio de una opción en función del spot, es la 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

De una derivada a una ecuación

El objetivo del método es resolver ecuaciones en derivadas parciales. Black–Scholes es una: el precio de la opción cambia en el tiempo según su curvatura respecto al spot. Tras un cambio de variables se convierte en la ecuación del calor, du/dt = D * d2u/dx2: el calor se extiende, y cuanto más curvo es el perfil de temperatura, más rápido se mueve. Se ponen los valores en una malla, se sustituyen ambas derivadas por diferencias y la ecuación pasa a ser una regla para actualizar cada punto a partir de sus vecinos.

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

En código, alpha = D * dt / dx**2 une el paso de tiempo dt y el espaciado dx. Una actualización, new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2]), mueve todos los puntos interiores a la vez con rebanadas; los dos extremos son las condiciones de contorno, fijas aquí. El esquema es simple, y el inconveniente también.

La estabilidad: el inconveniente

El esquema explícito solo es estable si alpha <= 0.5, es decir D * dt / dx**2 <= 0.5. Más allá, los errores de redondeo se amplifican en cada paso y los números explotan en dientes de sierra. La condición es dura: reducir a la mitad el espaciado obliga a un paso de tiempo cuatro veces menor, así que el coste crece rápido. Por eso el código profesional suele usar esquemas implícitos, que resuelven un pequeño sistema lineal en cada paso y son estables para cualquier dt; el explícito es el más simple de escribir y entender.

# 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

• Una diferencia finita sustituye una derivada por valores en una malla: la central (f(x + h) - f(x - h)) / (2h) tiene error del orden de h**2; la segunda (f(x + h) - 2 f(x) + f(x - h)) / h**2 da la curvatura (gamma).
• Las EDP se vuelven reglas de actualización en una malla: el paso explícito de la ecuación del calor es values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1]), escrito con rebanadas.
• El esquema explícito exige D * dt / dx**2 <= 0.5 para ser estable; los implícitos eliminan ese límite a costa de resolver un sistema lineal por paso.

A continuación: las matrices de covarianza y correlación, el primer paso de la optimización de carteras.

Lab: qué se califica

• Estimar una curvatura con second_difference25 pts
• Avanzar la malla un paso con explicit_step25 pts
• Comprobar el paso de tiempo con is_stable25 pts
• Avanzar en el tiempo con solve_heat25 pts