コードで学ぶ積分

細い短冊を足し合わせる

積分は曲線の下の面積を測るもので、多数の小さな刻みで f の値を足し合わせます。金融では、積分はさまざまな形で現れます。時間をかけたポジションの累積損益、連続キャッシュフローの現在価値、確率密度で重み付けした期待ペイオフ、一定期間のポートフォリオの総エクスポージャーなどです。きれいな式を持つ被積分関数は多くないので、面積を短冊で近似します。区間 [a, b] を幅 h = (b - a) / n の n 個の等しい短冊に分け、面積を足します。

# 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

より良い短冊：台形と Simpson

左端の長方形は粗い近似で、各短冊で曲線が平らだと見なします。台形公式は各短冊の両端を直線で結ぶので、短冊の面積は h * (f(左) + f(右)) / 2 になり、足し合わせると内側の点は2回、両端は1回数えられます：h * (f(a) / 2 + f(a + h) + ... + f(b - h) + f(b) / 2)。直線に対して正確です。Simpson の公式は3点ごとに放物線を当てはめ、h / 3 * (f(a) + 4 * f(奇数番目の点) + 2 * f(内側の偶数番目の点) + f(b)) となり、n は偶数が必要で、3次関数に対して正確です。滑らかな関数では、同じ計算量でもたいてい Simpson のほうがはるかに高精度です。

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)

手元のデータを積分する

関数がなく、一定間隔の測定値だけがあることもよくあります。たとえば毎日サンプリングしたエクスポージャーのリストです。台形公式はそのようなサンプルにそのまま使えます。間隔 dt の値 ys では、隣り合う各ペアが面積 dt * (ys[i] + ys[i + 1]) / 2 の台形になり、その合計がサンプルされた曲線の下の面積（たとえばエクスポージャー日数の合計）を近似します。サンプルが2つ未満なら区間がないので面積は 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

• 幅 h = (b - a) / n の短冊。左端長方形は i < n について h * sum(f(a + i*h))
• 台形：両端は半分に数える。Simpson：重み 1, 4, 2, 4, ..., 4, 1 に h / 3 を掛ける。n は偶数
• サンプルデータには dt * (最初 / 2 + 中間の値 + 最後 / 2) を使い、積分が分かっている例でテストする

次はベクトルと内積です。ポートフォリオのリターン、相関、そして後に学ぶ NumPy の基本要素です。

ラボ：採点される内容

• left_riemann で長方形による面積の近似を行う25点
• trapezoid で台形を使ってより正確に近似する25点
• simpson で放物線を当てはめる25点
• trapezoid_samples で測定データを積分する25点