f(x) = 0 を解く
多くの定量的な問いは、関数がゼロを横切る x を見つけることに行き着きます。債券を値付けする利回り、オプションに合うボラティリティ、プロジェクトの価値をゼロにする金利などです。f の根とは f(x) = 0 となる x のことです。たいてい公式はないので、プログラムが少しずつ近づきます。ほぼすべてを担う方法は2つで、二分法は遅いが確実、ニュートン法は速いが注意が要ります。どちらも関数そのものを引数に取り、Python では関数が値なのでそれができます。
def bisection(f, a, b, tol=1e-10):
if f(a) * f(b) > 0:
raise ValueError("f(a) and f(b) must have opposite signs")
while b - a > tol:
mid = (a + b) / 2
if f(a) * f(mid) <= 0: # the root is in the left half
b = mid
else:
a = mid
return (a + b) / 2
root = bisection(lambda x: x * x - 2, 0, 2)
print(root) # 1.414213562355144
print(root * root) # 1.9999999999492266二分法：遅いが確実
f が連続で f(a) と f(b) の符号が逆なら、その間に根があります。中点を評価すると、2つの半分のどちらかは端点の符号がまだ逆なので、根はそこに閉じ込められます。f がどんな形でも1回ごとに区間は半分になるので、n 回後の誤差は (b - a) / 2**n で、必要な回数は事前に ceil(log2((b - a) / tol)) と分かります。導関数が要らないこの保証が、二分法を安全な既定の選択にしています。
# How many halvings does a given tolerance cost? About 3.3 per decimal digit
import math
def bisection_steps(a, b, tol):
if b - a <= tol:
return 0
return math.ceil(math.log2((b - a) / tol))
for tol in (1e-2, 1e-4, 1e-8, 1e-12):
print(tol, bisection_steps(0, 2, tol)) # 1e-2: 8, 1e-4: 15, 1e-8: 28, 1e-12: 41ニュートン法：接線をたどる
ニュートン法は、現在の推定値で f を接線に置き換え、その直線がゼロを横切る点へ飛びます：x_new = x - f(x) / df(x)。単根の近くでは誤差が1回ごとにほぼ*2乗*されるので、正しい桁数は1から2、4、8、16と4回で倍々に増えます。導関数 df（インプライド・ボラティリティならすでにあるベガ）と良い初期値が必要です。
# Newton for sqrt(2): print every iterate and the error that remains
import math
def newton_step(f, df, x):
return x - f(x) / df(x)
f = lambda x: x * x - 2
df = lambda x: 2 * x
x = 1.0
for step in range(6):
print(step, x, abs(x - math.sqrt(2)))
x = newton_step(f, df, x)# Newton on the implied volatility of a call, using vega as the derivative
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def norm_pdf(x):
return math.exp(-0.5 * x * x) / math.sqrt(2 * math.pi)
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 vega(spot, strike, rate, sigma, t):
d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
return spot * norm_pdf(d1) * math.sqrt(t)
market = 8.0
sigma = 0.5 # a rough start
for step in range(1, 6):
error = bs_call(100, 100, 0.05, sigma, 0.5) - market
sigma -= error / vega(100, 100, 0.05, sigma, 0.5)
print(step, round(sigma, 10))ニュートン法が失敗するとき
接線は局所的な像です。導関数がゼロか非常に小さければ、ジャンプは巨大になるか定義できません。出発点が根から遠いと、反復は遠くへ飛んだり巡回したりします（atan(x) を x = 2 から始めると限りなく増えます）。そのため堅牢なソルバーは割り算を守り、反復回数に上限を設け、誤った数を返さずにエラーを出します。多くのライブラリは区間を挟んだ状態から始め、その内側でニュートン法を使い、はみ出したら二分法に切り替えます。
fの根はf(x) = 0となるx。ソルバーはf（とdf）を関数の引数として受け取る- 二分法：符号が変わる側の半分を残す。
ceil(log2((b - a) / tol))回で保証され、導関数は不要 - ニュートン法：
x - f(x) / df(x)。単根の近くでは1回ごとに桁が倍になるが、df・良い初期値・ガード（導関数ゼロ、非収束）が必要
次は二項ツリーです。株価を上下に動かし、後ろ向きにたどってオプションを値付けします。
ラボ：採点される内容
bisectionで区間を半分にして根を挟み込む25点bisection_stepsで半分にする回数を数える25点newton_stepで接線を1回たどる25点newtonで根に収束させる25点