Integrales en código

Sumar rebanadas finas

Una integral mide el área bajo una curva: suma el valor de f sobre muchos pasos diminutos. En finanzas las integrales están en todas partes, disfrazadas: la ganancia acumulada de una posición en el tiempo, el valor actual de un flujo continuo, un pago esperado como el promedio del pago ponderado por una densidad de probabilidad o la exposición total de una cartera en un periodo. Pocas funciones reales tienen una fórmula limpia, así que aproximamos el área con rebanadas. Divide el intervalo [a, b] en n rebanadas iguales de ancho h = (b - a) / n y suma las áreas.

# Area under x^2 on [0, 1] is exactly 1/3. Left rectangles, n slices:
def left_riemann(f, a, b, n):
    h = (b - a) / n
    return sum(f(a + i * h) for i in range(n)) * h

square = lambda x: x * x
for n in (10, 100, 1000):
    print(n, round(left_riemann(square, 0, 1, n), 6))
# 10 0.285
# 100 0.32835
# 1000 0.332834

Mejores rebanadas: trapecios y Simpson

Los rectángulos por la izquierda son toscos: cada rebanada finge que la curva es plana. La regla del trapecio une los dos extremos de cada rebanada con una recta, así que el área de una rebanada es h * (f(izquierda) + f(derecha)) / 2; al sumarlas, los puntos interiores cuentan dos veces y los extremos una: h * (f(a) / 2 + f(a + h) + ... + f(b - h) + f(b) / 2). Es exacta para rectas. La regla de Simpson ajusta una parábola a cada tres puntos: h / 3 * (f(a) + 4 * f(puntos impares) + 2 * f(puntos interiores pares) + f(b)), lo que exige un n par, y es exacta para cúbicas. Para funciones suaves, Simpson suele ser mucho más precisa con el mismo trabajo.

def trapezoid(f, a, b, n):
    h = (b - a) / n
    return h * (f(a) / 2 + sum(f(a + i * h) for i in range(1, n)) + f(b) / 2)

def simpson(f, a, b, n):
    if n % 2:
        raise ValueError("n must be even")
    h = (b - a) / n
    odd = sum(f(a + i * h) for i in range(1, n, 2))
    even = sum(f(a + i * h) for i in range(2, n, 2))
    return h / 3 * (f(a) + 4 * odd + 2 * even + f(b))

square = lambda x: x * x
for n in (10, 100, 1000):
    print(n, round(trapezoid(square, 0, 1, n), 6), round(simpson(square, 0, 1, n), 6))
# 10 0.335 0.333333
# 100 0.33335 0.333333
# 1000 0.333334 0.333333

print(simpson(lambda x: x**3, 0, 2, 4))   # 4.0  (exact for cubics)

Integrar datos que ya tienes

A menudo no hay función, solo mediciones a intervalos regulares: una lista de exposiciones muestreadas cada día. La regla del trapecio funciona directamente con esas muestras. Con valores ys separados dt, cada par vecino forma un trapecio de área dt * (ys[i] + ys[i + 1]) / 2, y la suma aproxima el área bajo la curva muestreada, por ejemplo el total de exposición-días. Con menos de dos muestras no hay intervalo, así que el área es 0.0.

def trapezoid_samples(ys, dt):
    if len(ys) < 2:
        return 0.0
    return dt * (ys[0] / 2 + sum(ys[1:-1]) + ys[-1] / 2)

# Exposure sampled daily: 0, 2, 4, 2, 0 million
exposure = [0, 2, 4, 2, 0]
print(trapezoid_samples(exposure, 1.0))   # 8.0  exposure-days
print(trapezoid_samples([5.0], 1.0))      # 0.0

import math

# An expected value is an integral: E[X^2] for a standard normal is 1
def density(x):
    return math.exp(-x * x / 2) / math.sqrt(2 * math.pi)

def simpson(f, a, b, n):
    h = (b - a) / n
    odd = sum(f(a + i * h) for i in range(1, n, 2))
    even = sum(f(a + i * h) for i in range(2, n, 2))
    return h / 3 * (f(a) + 4 * odd + 2 * even + f(b))

second_moment = simpson(lambda x: x * x * density(x), -8, 8, 400)
print(round(second_moment, 6))                   # 1.0

# Probability that X is within one standard deviation: about 68%
print(round(simpson(density, -1, 1, 100), 4))    # 0.6827

• Rebanadas de ancho h = (b - a) / n; rectángulos por la izquierda h * sum(f(a + i*h)) para i < n.
• Trapecio: los extremos cuentan la mitad; Simpson: pesos 1, 4, 2, 4, ..., 4, 1 por h / 3, n par.
• Para datos muestreados usa dt * (primero / 2 + valores intermedios + último / 2); prueba con integrales conocidas.

A continuación: vectores y productos escalares, la pieza básica de los rendimientos de cartera, las correlaciones y, más adelante, NumPy.

Lab: qué se califica

• Aproximar un área con rectángulos usando left_riemann25 pts
• Mejorar con trapecios usando trapezoid25 pts
• Ajustar parábolas con simpson25 pts
• Integrar datos medidos con trapezoid_samples25 pts