Busca de raízes: bissecção e Newton

Resolver f(x) = 0

Muitas perguntas quantitativas se resumem a achar o x em que uma função se anula: o rendimento que precifica um título, a volatilidade que casa com uma opção, a taxa que zera o valor de um projeto. Uma raiz de f é um x com f(x) = 0. Em geral não há fórmula, então um programa vai se aproximando. Dois métodos fazem quase todo o trabalho. A bisseção é lenta, mas certa; o método de Newton é rápido, mas exige cuidado. Ambos recebem a própria função como argumento, o que o Python permite porque funções são valores.

def bisection(f, a, b, tol=1e-10):
    if f(a) * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    while b - a > tol:
        mid = (a + b) / 2
        if f(a) * f(mid) <= 0:       # the root is in the left half
            b = mid
        else:
            a = mid
    return (a + b) / 2

root = bisection(lambda x: x * x - 2, 0, 2)
print(root)                          # 1.414213562355144
print(root * root)                   # 1.9999999999492266

Bisseção: lenta e segura

Se f é contínua e f(a) e f(b) têm sinais opostos, há uma raiz entre eles. Avalia-se o ponto médio: uma das metades ainda tem sinais opostos nas pontas, então a raiz fica presa ali. Cada passo divide o intervalo por dois, seja qual for f, de modo que o erro após n passos é (b - a) / 2**n e dá para saber de antemão quantos passos são necessários: ceil(log2((b - a) / tol)). Essa garantia, sem precisar de derivada, é o que faz da bisseção a escolha segura por padrão.

# How many halvings does a given tolerance cost? About 3.3 per decimal digit
import math

def bisection_steps(a, b, tol):
    if b - a <= tol:
        return 0
    return math.ceil(math.log2((b - a) / tol))

for tol in (1e-2, 1e-4, 1e-8, 1e-12):
    print(tol, bisection_steps(0, 2, tol))       # 1e-2: 8, 1e-4: 15, 1e-8: 28, 1e-12: 41

Newton: seguir a tangente

O método de Newton substitui f pela sua tangente no palpite atual e salta para onde essa reta cruza zero: x_new = x - f(x) / df(x). Perto de uma raiz simples o erro é aproximadamente *elevado ao quadrado* a cada passo, de modo que o número de dígitos corretos dobra: de 1 para 2, 4, 8, 16 em quatro passos. Exige a derivada df (na volatilidade implícita é o vega, que você já tem) e um bom ponto de partida.

# Newton for sqrt(2): print every iterate and the error that remains
import math

def newton_step(f, df, x):
    return x - f(x) / df(x)

f = lambda x: x * x - 2
df = lambda x: 2 * x
x = 1.0
for step in range(6):
    print(step, x, abs(x - math.sqrt(2)))
    x = newton_step(f, df, x)

# Newton on the implied volatility of a call, using vega as the derivative
import math

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

def norm_pdf(x):
    return math.exp(-0.5 * x * x) / math.sqrt(2 * math.pi)

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 vega(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_pdf(d1) * math.sqrt(t)

market = 8.0
sigma = 0.5                                              # a rough start
for step in range(1, 6):
    error = bs_call(100, 100, 0.05, sigma, 0.5) - market
    sigma -= error / vega(100, 100, 0.05, sigma, 0.5)
    print(step, round(sigma, 10))

Quando o Newton falha

A tangente é uma imagem local. Se a derivada é zero ou minúscula, o salto é enorme ou indefinido; se o início está longe da raiz, as iteradas podem se afastar ou oscilar (para atan(x) a partir de x = 2 crescem sem limite). Um solucionador robusto protege a divisão, limita as iterações e levanta um erro em vez de devolver um número errado. Muitas bibliotecas partem de um intervalo, usam Newton dentro dele e voltam à bisseção quando um passo sai.

• Uma raiz de f é um x com f(x) = 0; os solucionadores recebem f (e df) como argumentos função.
• Bisseção: ficar com a metade em que o sinal muda; ceil(log2((b - a) / tol)) passos, garantidos, sem derivada.
• Newton: x - f(x) / df(x); os dígitos dobram a cada passo perto de uma raiz simples, mas exige df, um bom início e salvaguardas (derivada zero, não convergência).

A seguir: as árvores binomiais, precificar uma opção fazendo a ação subir e descer e trabalhando de trás para frente.

Lab: o que é avaliado

• Encurralar uma raiz dividindo ao meio com bisection25 pts
• Contar as divisões ao meio com bisection_steps25 pts
• Seguir a tangente uma vez com newton_step25 pts
• Convergir para uma raiz com newton25 pts