Derivadas en código

¿Con qué rapidez cambia?

La derivada de una función f en x es su tasa de cambio instantánea: la pendiente de la curva en ese punto. En finanzas casi toda sensibilidad es una derivada. El delta de una opción es cuánto se mueve su precio por unidad de movimiento del subyacente, su gamma es la rapidez con que cambia el propio delta, y la duración de un bono es una derivada del precio respecto al rendimiento. El cálculo da fórmulas exactas para funciones simples (la derivada de x**2 es 2*x), pero las funciones de valoración reales suelen ser demasiado complicadas o solo existen como código. Entonces estimamos la derivada numéricamente.

# 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

Diferencias finitas: adelante y centrada

La diferencia hacia adelante mira un pequeño paso h adelante: (f(x + h) - f(x)) / h. Su error decrece en proporción a h. La diferencia centrada mira un paso a cada lado: (f(x + h) - f(x - h)) / (2 * h). Los errores de ambos lados se cancelan, así que su error decrece como h al cuadrado y es exacta para cuadráticas. La segunda derivada (gamma) combina tres puntos: (f(x + h) - 2 * f(x) + f(x - h)) / h**2. Pasa la propia función como argumento (f), una habilidad de los temas de funciones.

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

Elegir el tamaño del paso

Un h más pequeño reduce el error de truncamiento de la fórmula, pero restar dos números casi iguales pierde dígitos: este error de redondeo crece al disminuir h, y en una segunda derivada crece más rápido porque además dividimos por h**2. Un h demasiado pequeño (como 1e-12) puede dar basura. Un buen compromiso es alrededor de 1e-5 para primeras derivadas con diferencias centradas y de 1e-4 para segundas derivadas. En el código de valoración real, «bump and reprice» es justo esto: volver a ejecutar el valorador con el spot subido y bajado una pequeña cantidad.

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: el método de Newton

Una derivada indica en qué dirección moverse para llegar a una raíz, un punto donde f(x) = 0. El método de Newton parte de una estimación x y la reemplaza repetidamente por x - f(x) / f'(x), siguiendo la tangente hasta cero. Converge muy rápido cerca de una raíz, por eso se usa para obtener una volatilidad implícita a partir de un precio de opción. Detente cuando abs(f(x)) baje de una tolerancia, rinde tras un máximo de iteraciones y lanza un error si la derivada es cero o el bucle no converge, en vez de devolver un número erróneo.

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

• Adelante (f(x+h) - f(x)) / h; centrada (f(x+h) - f(x-h)) / (2h); segunda (f(x+h) - 2f(x) + f(x-h)) / h**2.
• Ni muy grande ni muy pequeño: alrededor de 1e-5 para primeras y 1e-4 para segundas derivadas.
• Newton: x = x - f(x) / f'(x) hasta que abs(f(x)) sea mínimo; lanza un error si no converge.

A continuación: integrales en código, sumar muchas rebanadas finas para hallar áreas, totales y valores esperados.

Lab: qué se califica

• Estimar una pendiente mirando hacia adelante con forward_diff25 pts
• Estimar una pendiente desde ambos lados con central_diff25 pts
• Estimar la curvatura (gamma) con second_derivative25 pts
• Hallar dónde una función se anula con newton_root25 pts