Integrale im Code

Dünne Streifen aufaddieren

Ein Integral misst die Fläche unter einer Kurve: Man addiert den Wert von f über viele winzige Schritte. In der Finanzwelt stecken Integrale überall verkleidet: der kumulierte Gewinn einer Position über die Zeit, der Barwert eines kontinuierlichen Zahlungsstroms, ein erwarteter Payoff als mit einer Wahrscheinlichkeitsdichte gewichteter Mittelwert oder das gesamte Exposure eines Portfolios über einen Zeitraum. Wenige echte Integranden haben eine saubere Formel, also nähern wir die Fläche mit Streifen an. Teilen Sie das Intervall [a, b] in n gleiche Streifen der Breite h = (b - a) / n und addieren Sie deren Flächen.

# 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

Bessere Streifen: Trapeze und Simpson

Linke Rechtecke sind grob: Jeder Streifen tut so, als wäre die Kurve flach. Die Trapezregel verbindet die beiden Enden jedes Streifens mit einer Geraden, die Fläche eines Streifens ist also h * (f(links) + f(rechts)) / 2; beim Addieren zählen die inneren Punkte doppelt und die beiden Endpunkte einfach: h * (f(a) / 2 + f(a + h) + ... + f(b - h) + f(b) / 2). Sie ist für Geraden exakt. Die Simpson-Regel legt durch je drei Punkte eine Parabel: h / 3 * (f(a) + 4 * f(ungerade Punkte) + 2 * f(gerade innere Punkte) + f(b)), was ein gerades n braucht, und ist für kubische Funktionen exakt. Bei glatten Funktionen ist Simpson meist bei gleichem Aufwand weit genauer.

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)

Bereits vorhandene Daten integrieren

Oft gibt es keine Funktion, nur Messungen zu regelmäßigen Zeiten: eine Liste täglich erfasster Exposures. Die Trapezregel funktioniert direkt auf solchen Stichproben. Bei Werten ys im Abstand dt bildet jedes benachbarte Paar ein Trapez der Fläche dt * (ys[i] + ys[i + 1]) / 2, und die Summe nähert die Fläche unter der abgetasteten Kurve an, etwa die gesamten Exposure-Tage. Bei weniger als zwei Werten gibt es kein Intervall, die Fläche ist also 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

• Streifen der Breite h = (b - a) / n; linke Rechtecke h * sum(f(a + i*h)) für i < n.
• Trapez: die Enden zählen halb; Simpson: Gewichte 1, 4, 2, 4, ..., 4, 1 mal h / 3, n gerade.
• Für Messreihen dt * (erster / 2 + mittlere Werte + letzter / 2); an bekannten Integralen testen.

Als Nächstes: Vektoren und Skalarprodukte, der Baustein für Portfoliorenditen, Korrelationen und später NumPy.

Lab: Was bewertet wird

• Mit left_riemann eine Fläche durch Rechtecke annähern25 Pkt.
• Mit trapezoid und Trapezen genauer werden25 Pkt.
• Mit simpson Parabeln anpassen25 Pkt.
• Mit trapezoid_samples Messdaten integrieren25 Pkt.