Métodos de diferenças finitas

Derivadas em uma grade

Uma derivada é um limite, mas o computador trabalha com passos finitos. Uma diferença finita substitui f'(x) pela inclinação sobre um passo pequeno h: progressiva (f(x + h) - f(x)) / h ou central (f(x + h) - f(x - h)) / (2h), mais precisa porque os erros dos dois lados se compensam (o erro cai como h**2, o da progressiva só como h). A mesma ideia dá a segunda derivada, a curvatura: (f(x + h) - 2 f(x) + f(x - h)) / h**2. Para o preço de uma opção em função do spot, isso é o gama.

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 uma derivada a uma equação

O objetivo do método é resolver equações diferenciais parciais. Black–Scholes é uma: o preço da opção muda no tempo conforme sua curvatura em relação ao spot. Após uma mudança de variáveis ela vira a equação do calor, du/dt = D * d2u/dx2: o calor se espalha, e quanto mais curvo o perfil de temperatura, mais rápido ele se move. Colocam-se os valores numa grade, substituem-se as duas derivadas por diferenças e a equação vira uma regra para atualizar cada ponto a partir de seus vizinhos.

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

No código, alpha = D * dt / dx**2 liga o passo de tempo dt e o espaçamento dx. Uma atualização, new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2]), move todos os pontos internos de uma vez com fatias; as duas pontas são as condições de contorno, fixas aqui. O esquema é simples, e a armadilha também.

A estabilidade: a armadilha

O esquema explícito só é estável se alpha <= 0.5, isto é, D * dt / dx**2 <= 0.5. Além disso, os erros de arredondamento são amplificados a cada passo e os números explodem em dente de serra. A condição é dura: reduzir o espaçamento à metade exige um passo de tempo quatro vezes menor, e o custo cresce rápido. Por isso o código profissional costuma usar esquemas implícitos, que resolvem um pequeno sistema linear a cada passo e são estáveis para qualquer dt; o explícito é o mais simples de escrever e 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

• Uma diferença finita substitui uma derivada por valores numa grade: a central (f(x + h) - f(x - h)) / (2h) tem erro da ordem de h**2; a segunda (f(x + h) - 2 f(x) + f(x - h)) / h**2 dá a curvatura (gama).
• As EDPs viram regras de atualização numa grade: o passo explícito da equação do calor é values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1]), escrito com fatias.
• O esquema explícito exige D * dt / dx**2 <= 0.5 para ser estável; os implícitos removem esse limite ao custo de resolver um sistema linear por passo.

A seguir: as matrizes de covariância e correlação, o primeiro passo da otimização de carteiras.

Lab: o que é avaliado

• Estimar uma curvatura com second_difference25 pts
• Avançar a grade um passo com explicit_step25 pts
• Verificar o passo de tempo com is_stable25 pts
• Avançar no tempo com solve_heat25 pts