コードで学ぶ微分

どれくらいの速さで変わるか

関数 f の x における導関数は瞬間的な変化率で、その点での曲線の傾きです。金融の感応度はほとんどが微分です。オプションのデルタは原資産が1単位動いたときの価格変化、ガンマはデルタ自体の変化の速さ、債券のデュレーションは利回りに対する価格の微分です。微積分は単純な関数には正確な式を与えます（x**2 の導関数は 2*x）が、実際の価格関数は複雑すぎたり、コードとしてしか存在しなかったりします。その場合は導関数を数値的に推定します。

# A function is just a Python function
def f(x):
    return x * x

# Slope between two nearby points approaches the true slope (2x = 6 at x = 3)
x = 3.0
for h in (1.0, 0.1, 0.001):
    slope = (f(x + h) - f(x)) / h
    print(h, round(slope, 4))
# 1.0 7.0
# 0.1 6.1
# 0.001 6.001

有限差分：前進差分と中心差分

前進差分は少し先 h を見ます：(f(x + h) - f(x)) / h。誤差は h に比例して小さくなります。中心差分は両側に1歩ずつ見ます：(f(x + h) - f(x - h)) / (2 * h)。両側の誤差が打ち消し合うため、誤差は h の2乗で小さくなり、2次関数なら正確です。2階微分（ガンマ）は3点を組み合わせます：(f(x + h) - 2 * f(x) + f(x - h)) / h**2。関数そのものを引数（f）として渡します。これは関数のトピックで学んだ技術です。

import math

def forward_diff(f, x, h=1e-6):
    return (f(x + h) - f(x)) / h

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def second_derivative(f, x, h=1e-4):
    return (f(x + h) - 2 * f(x) + f(x - h)) / h**2

print(round(central_diff(math.sin, 0.0), 6))            # 1.0  (cos 0)
print(round(central_diff(lambda x: x**3, 2.0), 4))      # 12.0 (3 x^2)
print(round(second_derivative(lambda x: x**3, 2.0), 3)) # 12.0 (6 x)

# A quadratic: central is exact even for a big step, forward is not
print(central_diff(lambda x: x * x, 3.0, h=0.5))  # 6.0
print(forward_diff(lambda x: x * x, 3.0, h=0.5))  # 6.5

刻み幅の選び方

h を小さくすると式の打ち切り誤差は減りますが、ほぼ等しい2つの数を引くと桁落ちします。この丸め誤差は h が小さいほど増え、2階微分では h**2 でも割るためさらに速く増えます。h が小さすぎる（1e-12 など）と無意味な値になることもあります。中心差分の1階微分では約 1e-5、2階微分では約 1e-4 が妥当な折衷です。実際の価格計算コードの「バンプ・アンド・リプライス」はまさにこれで、スポットを少し上げ下げして価格計算をやり直します。

import math

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

# True derivative of exp at 0 is 1; watch the error as h shrinks
for h in (1e-1, 1e-3, 1e-5, 1e-8, 1e-12):
    error = abs(central_diff(math.exp, 0.0, h) - 1.0)
    print(f"h={h:g} error={error:.1e}")
# h=0.1 error=1.7e-03
# h=0.001 error=1.7e-07
# h=1e-05 error=1.2e-11
# h=1e-08 error=5.3e-10     (rounding error takes over)
# h=1e-12 error=3.3e-05     (much worse: too small)

微分の応用：ニュートン法

導関数は、根（f(x) = 0 となる点）に近づくにはどちらへ動けばよいかを教えてくれます。ニュートン法は推定値 x から始め、接線に沿ってゼロまで下りるように、x - f(x) / f'(x) で繰り返し置き換えます。根の近くでは非常に速く収束するため、オプション価格からインプライド・ボラティリティを逆算するのに使われます。abs(f(x)) が許容誤差を下回ったら止め、反復回数の上限で打ち切り、導関数が0だったり収束しなかったりしたら、誤った数値を返さずエラーを送出します。

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def newton_root(f, x0, tol=1e-10, max_iter=50):
    x = x0
    for _ in range(max_iter):
        if abs(f(x)) < tol:
            return x
        slope = central_diff(f, x)
        if slope == 0:
            raise ValueError("flat slope: Newton cannot continue")
        x = x - f(x) / slope
    raise ValueError("did not converge")

# The square root of 2 is the positive root of x^2 - 2
root = newton_root(lambda x: x * x - 2, 1.0)
print(round(root, 8))           # 1.41421356

# x^2 + 1 never reaches zero
try:
    newton_root(lambda x: x * x + 1, 1.0)
except ValueError:
    print("no root found")      # no root found

• 前進 (f(x+h) - f(x)) / h、中心 (f(x+h) - f(x-h)) / (2h)、2階 (f(x+h) - 2f(x) + f(x-h)) / h**2
• 大きすぎず小さすぎず：1階は約 1e-5、2階は約 1e-4
• ニュートン法：abs(f(x)) が十分小さくなるまで x = x - f(x) / f'(x)。収束しなければエラーにする

次はコードで積分です。多数の細い短冊を足し合わせて、面積・総量・期待値を求めます。

ラボ：採点される内容

• forward_diff で前向きに傾きを推定する25点
• central_diff で両側から傾きを推定する25点
• second_derivative で曲率（ガンマ）を推定する25点
• newton_root で関数がゼロになる点を見つける25点