Integrais em código

Somar fatias finas

Uma integral mede a área sob uma curva: some o valor de f em muitos passos minúsculos. Em finanças as integrais estão por toda parte, disfarçadas: o lucro acumulado de uma posição ao longo do tempo, o valor presente de um fluxo contínuo, um retorno esperado como a média do retorno ponderada por uma densidade de probabilidade, ou a exposição total de uma carteira num período. Poucas funções reais têm uma fórmula limpa, então aproximamos a área com fatias. Divida o intervalo [a, b] em n fatias iguais de largura h = (b - a) / n e some as á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

Fatias melhores: trapézios e Simpson

Retângulos à esquerda são grosseiros: cada fatia finge que a curva é plana. A regra do trapézio liga as duas pontas de cada fatia por uma reta, então a área de uma fatia é h * (f(esquerda) + f(direita)) / 2; ao somá-las, os pontos internos contam duas vezes e as pontas uma: h * (f(a) / 2 + f(a + h) + ... + f(b - h) + f(b) / 2). É exata para retas. A regra de Simpson ajusta uma parábola a cada três pontos: h / 3 * (f(a) + 4 * f(pontos ímpares) + 2 * f(pontos internos pares) + f(b)), o que exige um n par, e é exata para cúbicas. Para funções suaves, Simpson costuma ser muito mais precisa com o mesmo trabalho.

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 dados que você já tem

Muitas vezes não há função, só medições em tempos regulares: uma lista de exposições amostradas todo dia. A regra do trapézio funciona diretamente nessas amostras. Com valores ys espaçados de dt, cada par vizinho forma um trapézio de área dt * (ys[i] + ys[i + 1]) / 2, e a soma aproxima a área sob a curva amostrada, por exemplo o total de exposição-dias. Com menos de duas amostras não há intervalo, então a área é 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

• Fatias de largura h = (b - a) / n; retângulos à esquerda h * sum(f(a + i*h)) para i < n.
• Trapézio: as pontas contam metade; Simpson: pesos 1, 4, 2, 4, ..., 4, 1 vezes h / 3, n par.
• Para dados amostrados use dt * (primeiro / 2 + valores do meio + último / 2); teste em integrais conhecidas.

A seguir: vetores e produtos escalares, o bloco básico dos retornos de carteira, das correlações e, mais tarde, do NumPy.

Lab: o que é avaliado

• Aproximar uma área com retângulos usando left_riemann25 pts
• Melhorar com trapézios usando trapezoid25 pts
• Ajustar parábolas com simpson25 pts
• Integrar dados medidos com trapezoid_samples25 pts