Intégrales en code

Additionner de fines tranches

Une intégrale mesure l'aire sous une courbe : on additionne la valeur de f sur de nombreux petits pas. En finance, les intégrales sont partout, déguisées : le gain cumulé d'une position dans le temps, la valeur actuelle d'un flux continu, un gain espéré comme moyenne du gain pondérée par une densité de probabilité, ou l'exposition totale d'un portefeuille sur une période. Peu de fonctions réelles ont une formule propre, donc on approche l'aire par des tranches. Découpez l'intervalle [a, b] en n tranches égales de largeur h = (b - a) / n et additionnez leurs aires.

# 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

De meilleures tranches : trapèzes et Simpson

Les rectangles à gauche sont grossiers : chaque tranche fait comme si la courbe était plate. La règle des trapèzes relie les deux extrémités de chaque tranche par une droite, donc l'aire d'une tranche est h * (f(gauche) + f(droite)) / 2 ; en les additionnant, les points intérieurs comptent deux fois et les deux extrémités une fois : h * (f(a) / 2 + f(a + h) + ... + f(b - h) + f(b) / 2). Elle est exacte pour les droites. La règle de Simpson ajuste une parabole sur chaque triplet de points : h / 3 * (f(a) + 4 * f(points impairs) + 2 * f(points intérieurs pairs) + f(b)), ce qui exige un n pair, et elle est exacte pour les cubiques. Pour des fonctions lisses, Simpson est en général bien plus précise pour le même travail.

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)

Intégrer des données que l'on possède déjà

Il n'y a souvent pas de fonction, seulement des mesures à intervalles réguliers : une liste d'expositions échantillonnées chaque jour. La règle des trapèzes s'applique directement à de tels échantillons. Avec des valeurs ys espacées de dt, chaque paire voisine forme un trapèze d'aire dt * (ys[i] + ys[i + 1]) / 2, et la somme approche l'aire sous la courbe échantillonnée, par exemple le total d'expositions-jours. Avec moins de deux échantillons il n'y a pas d'intervalle, donc l'aire vaut 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

• Tranches de largeur h = (b - a) / n ; rectangles à gauche h * sum(f(a + i*h)) pour i < n.
• Trapèzes : les extrémités comptent pour moitié ; Simpson : poids 1, 4, 2, 4, ..., 4, 1 fois h / 3, n pair.
• Pour des données échantillonnées, dt * (premier / 2 + valeurs du milieu + dernier / 2) ; testez sur des intégrales connues.

À suivre : vecteurs et produits scalaires, la brique de base des rendements de portefeuille, des corrélations et, plus tard, de NumPy.

Lab : ce qui est noté

• Approcher une aire par des rectangles avec left_riemann25 pts
• Faire mieux avec des trapèzes grâce à trapezoid25 pts
• Ajuster des paraboles avec simpson25 pts
• Intégrer des données mesurées avec trapezoid_samples25 pts