差分法

格子上の微分

微分は極限ですが、コンピュータは有限の刻みで計算します。差分は f'(x) を小さな刻み h の傾きで置き換えます。前進差分 (f(x + h) - f(x)) / h、あるいは中心差分 (f(x + h) - f(x - h)) / (2h) で、後者は両側の誤差が打ち消し合うのでより正確です（誤差は h**2 で縮み、前進差分は h でしか縮みません）。同じ考え方で2階微分、つまり曲率が得られます：(f(x + h) - 2 f(x) + f(x - h)) / h**2。オプション価格を現物価格の関数と見れば、これがガンマです。

import math

def forward(f, x, h):
    return (f(x + h) - f(x)) / h

def central(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

exact = math.cos(1.0)                               # the derivative of sin at 1
for h in (0.1, 0.01, 0.001):
    print(h, abs(forward(math.sin, 1.0, h) - exact), abs(central(math.sin, 1.0, h) - exact))

# Gamma of a call from the second difference of its Black-Scholes price
import math

def norm_cdf(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def bs_call(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))

def second_difference(f, x, h):
    return (f(x + h) - 2 * f(x) + f(x - h)) / h ** 2

price = lambda s: bs_call(s, 100, 0.05, 0.2, 1)
print(round(second_difference(price, 100, 0.5), 6))     # 0.018761
print(round(second_difference(price, 100, 0.01), 6))    # 0.018762

微分から方程式へ

この方法の目的は偏微分方程式を解くことです。ブラック＝ショールズもその1つで、オプション価格は現物価格に関する曲率に応じて時間とともに変わります。変数変換すると熱方程式 du/dt = D * d2u/dx2 になります。熱は広がり、温度分布の曲がりが強いほど速く動きます。値を格子に置き、2つの微分を差分に置き換えると、方程式は各点を隣の点から更新する規則になります。

# The explicit scheme: a spike of heat spreads to its neighbours, step by step
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

grid = np.zeros(9)
grid[4] = 1.0
for step in range(4):
    print(step, np.round(grid, 3))
    grid = explicit_step(grid, 0.25)

コードでは alpha = D * dt / dx**2 が時間刻み dt と格子間隔 dx を結びつけます。1回の更新 new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2]) は、スライスですべての内部点を一度に動かし、両端は境界条件でここでは固定です。手法は単純で、落とし穴も単純です。

安定性：落とし穴

陽解法が安定なのは alpha <= 0.5、つまり D * dt / dx**2 <= 0.5 のときだけです。それを超えると丸め誤差が毎ステップ増幅され、数値はのこぎり状に発散します。条件は厳しく、格子間隔を半分にすると時間刻みは4分の1にしなければならず、コストは急増します。そのため実務のコードでは、各ステップで小さな連立一次方程式を解き、どんな dt でも安定な陰解法がよく使われます。陽解法は書くのも理解するのもいちばん簡単です。

# The same spike with alpha = 0.4 (stable) and alpha = 0.6 (unstable) after 60 steps
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

for alpha in (0.4, 0.6):
    grid = np.zeros(21)
    grid[10] = 1.0
    for _ in range(60):
        grid = explicit_step(grid, alpha)
    print(alpha, "largest value", float(np.abs(grid).max()))       # 0.4: 0.0556, 0.6: about 3.1e7

• 差分は微分を格子上の値で置き換える。中心差分 (f(x + h) - f(x - h)) / (2h) の誤差は h**2 程度、2階差分 (f(x + h) - 2 f(x) + f(x - h)) / h**2 は曲率（ガンマ）を与える
• 偏微分方程式は格子上の更新規則になる。熱方程式の陽的ステップは values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1]) で、スライスで書く
• 陽解法が安定であるには D * dt / dx**2 <= 0.5 が必要。陰解法はこの制限をなくすが、各ステップで連立一次方程式を解く代償がある

次は共分散行列と相関行列です。ポートフォリオ最適化の第一歩です。

ラボ：採点される内容

• second_difference で曲率を推定する25点
• explicit_step で格子を1ステップ進める25点
• is_stable で時間刻みを確認する25点
• solve_heat で時間を進めながら解く25点