Búsqueda de raíces: bisección y Newton

Resolver f(x) = 0

Muchas preguntas cuantitativas se reducen a hallar la x donde una función se anula: el rendimiento que valora un bono, la volatilidad que iguala una opción, el tipo que hace cero el valor de un proyecto. Una raíz de f es una x con f(x) = 0. Normalmente no hay fórmula, así que un programa se va acercando. Dos métodos hacen casi todo el trabajo. La bisección es lenta pero segura; el método de Newton es rápido pero requiere cuidado. Ambos reciben la función misma como argumento, algo que Python permite porque las funciones son 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

Bisección: lenta y segura

Si f es continua y f(a) y f(b) tienen signos opuestos, hay una raíz entre ambos. Se evalúa el punto medio: una de las dos mitades sigue con signos opuestos en sus extremos, así que la raíz queda atrapada allí. Cada paso divide el intervalo por dos, sea cual sea f, de modo que el error tras n pasos es (b - a) / 2**n y se puede saber de antemano cuántos pasos hacen falta: ceil(log2((b - a) / tol)). Esa garantía, sin necesidad de derivada, hace de la bisección la opción segura por defecto.

# 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 la tangente

El método de Newton sustituye f por su tangente en la estimación actual y salta adonde esa recta cruza cero: x_new = x - f(x) / df(x). Cerca de una raíz simple el error se *eleva aproximadamente al cuadrado* en cada paso, de modo que el número de cifras correctas se duplica: de 1 a 2, 4, 8, 16 en cuatro pasos. Necesita la derivada df (para la volatilidad implícita es la vega, que ya tienes) y un buen punto 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))

Cuándo falla Newton

La tangente es una imagen local. Si la derivada es cero o diminuta, el salto es enorme o indefinido; si el inicio está lejos de la raíz, las iteraciones pueden alejarse u oscilar (para atan(x) desde x = 2 crecen sin límite). Un solucionador robusto protege la división, limita las iteraciones y lanza un error en lugar de devolver un número falso. Muchas bibliotecas parten de un intervalo, usan Newton dentro y vuelven a la bisección cuando un paso se sale.

• Una raíz de f es una x con f(x) = 0; los solucionadores toman f (y df) como argumentos función.
• Bisección: quedarse con la mitad donde cambia el signo; ceil(log2((b - a) / tol)) pasos, garantizados, sin derivada.
• Newton: x - f(x) / df(x); las cifras se duplican en cada paso cerca de una raíz simple, pero necesita df, un buen inicio y salvaguardas (derivada cero, no convergencia).

A continuación: los árboles binomiales, valorar una opción haciendo subir y bajar la acción y trabajando hacia atrás.

Lab: qué se califica

• Atrapar una raíz dividiendo a la mitad con bisection25 pts
• Contar las divisiones a la mitad con bisection_steps25 pts
• Seguir la tangente una vez con newton_step25 pts
• Converger a una raíz con newton25 pts