Derivadas em código

Com que rapidez muda?

A derivada de uma função f em x é sua taxa de variação instantânea: a inclinação da curva naquele ponto. Em finanças quase toda sensibilidade é uma derivada. O delta de uma opção é quanto seu preço se move por unidade de movimento do ativo subjacente, seu gama é a rapidez com que o próprio delta muda, e a duração de um título é uma derivada do preço em relação ao rendimento. O cálculo dá fórmulas exatas para funções simples (a derivada de x**2 é 2*x), mas funções de precificação reais costumam ser complicadas demais ou existir só como código. Então estimamos a derivada numericamente.

# A function is just a Python function
def f(x):
    return x * x

# Slope between two nearby points approaches the true slope (2x = 6 at x = 3)
x = 3.0
for h in (1.0, 0.1, 0.001):
    slope = (f(x + h) - f(x)) / h
    print(h, round(slope, 4))
# 1.0 7.0
# 0.1 6.1
# 0.001 6.001

Diferenças finitas: para a frente e central

A diferença para a frente olha um pequeno passo h adiante: (f(x + h) - f(x)) / h. Seu erro diminui em proporção a h. A diferença central olha um passo para cada lado: (f(x + h) - f(x - h)) / (2 * h). Os erros dos dois lados se cancelam, então seu erro diminui como h ao quadrado e ela é exata para quadráticas. A segunda derivada (gama) combina três pontos: (f(x + h) - 2 * f(x) + f(x - h)) / h**2. Passe a própria função como argumento (f), uma habilidade dos tópicos de funções.

import math

def forward_diff(f, x, h=1e-6):
    return (f(x + h) - f(x)) / h

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def second_derivative(f, x, h=1e-4):
    return (f(x + h) - 2 * f(x) + f(x - h)) / h**2

print(round(central_diff(math.sin, 0.0), 6))            # 1.0  (cos 0)
print(round(central_diff(lambda x: x**3, 2.0), 4))      # 12.0 (3 x^2)
print(round(second_derivative(lambda x: x**3, 2.0), 3)) # 12.0 (6 x)

# A quadratic: central is exact even for a big step, forward is not
print(central_diff(lambda x: x * x, 3.0, h=0.5))  # 6.0
print(forward_diff(lambda x: x * x, 3.0, h=0.5))  # 6.5

Escolher o tamanho do passo

Um h menor reduz o erro de truncamento da fórmula, mas subtrair dois números quase iguais perde dígitos: esse erro de arredondamento cresce quando h diminui, e numa segunda derivada cresce mais rápido porque também dividimos por h**2. Um h pequeno demais (como 1e-12) pode dar lixo. Um bom compromisso é cerca de 1e-5 para primeiras derivadas com diferenças centrais e cerca de 1e-4 para segundas derivadas. No código de precificação real, «bump and reprice» é exatamente isto: reexecutar o precificador com o spot aumentado e diminuído de uma pequena quantia.

import math

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

# True derivative of exp at 0 is 1; watch the error as h shrinks
for h in (1e-1, 1e-3, 1e-5, 1e-8, 1e-12):
    error = abs(central_diff(math.exp, 0.0, h) - 1.0)
    print(f"h={h:g} error={error:.1e}")
# h=0.1 error=1.7e-03
# h=0.001 error=1.7e-07
# h=1e-05 error=1.2e-11
# h=1e-08 error=5.3e-10     (rounding error takes over)
# h=1e-12 error=3.3e-05     (much worse: too small)

Usar derivadas: o método de Newton

Uma derivada indica em que direção ir para alcançar uma raiz, um ponto onde f(x) = 0. O método de Newton parte de um palpite x e o substitui repetidamente por x - f(x) / f'(x), seguindo a tangente até zero. Converge muito rápido perto de uma raiz, por isso é usado para obter uma volatilidade implícita a partir de um preço de opção. Pare quando abs(f(x)) ficar abaixo de uma tolerância, desista após um número máximo de iterações e gere um erro se a derivada for zero ou o laço não convergir, em vez de devolver um número errado.

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def newton_root(f, x0, tol=1e-10, max_iter=50):
    x = x0
    for _ in range(max_iter):
        if abs(f(x)) < tol:
            return x
        slope = central_diff(f, x)
        if slope == 0:
            raise ValueError("flat slope: Newton cannot continue")
        x = x - f(x) / slope
    raise ValueError("did not converge")

# The square root of 2 is the positive root of x^2 - 2
root = newton_root(lambda x: x * x - 2, 1.0)
print(round(root, 8))           # 1.41421356

# x^2 + 1 never reaches zero
try:
    newton_root(lambda x: x * x + 1, 1.0)
except ValueError:
    print("no root found")      # no root found

• Para a frente (f(x+h) - f(x)) / h; central (f(x+h) - f(x-h)) / (2h); segunda (f(x+h) - 2f(x) + f(x-h)) / h**2.
• Nem grande nem pequeno demais: cerca de 1e-5 para primeiras e 1e-4 para segundas derivadas.
• Newton: x = x - f(x) / f'(x) até abs(f(x)) ficar minúsculo; gere um erro se não convergir.

A seguir: integrais em código, somar muitas fatias finas para achar áreas, totais e valores esperados.

Lab: o que é avaliado

• Estimar uma inclinação olhando adiante com forward_diff25 pts
• Estimar uma inclinação dos dois lados com central_diff25 pts
• Estimar a curvatura (gama) com second_derivative25 pts
• Encontrar onde uma função se anula com newton_root25 pts