価格はパーセントで変わる
算術ブラウン運動は、株価が5でも500でも同じドル額のノイズを足し、ゼロを下回ることもあります。株は違い、1%の動きはどの価格でも1%です。幾何ブラウン運動はランダムさを変化率に置きます：dS = mu * S * dt + sigma * S * dW。ブラック＝ショールズのモデルで、mu は年率の期待リターン、sigma は市場で提示されるボラティリティです。
import math
def gbm_step(s, mu, sigma, dt, z):
return s * math.exp((mu - 0.5 * sigma ** 2) * dt + sigma * math.sqrt(dt) * z)
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 0.0), 4)) # a calm day: only the (mu - sigma^2/2) drift acts
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 1.0), 4)) # a one-sigma up day
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, -50.0), 4)) # an absurd fall: still positive厳密解とマイナス2分の1シグマ二乗
方程を解くと S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t)) になります。なぜ奇妙な - sigma**2 / 2 があるのか。指数関数は凸なので、イェンセンの不等式によりノイズは exp(ノイズ) の平均を押し上げ、sigma**2 / 2 を引くとそれがちょうど打ち消されて、*期待*価格は率 mu で成長します。典型的な*経路*はそれより遅く mu - sigma**2 / 2 で伸びるため、株価の中央値は平均より低くなります。この項は次のトピックの伊藤の補題で導かれます。
import numpy as np
def gbm_path(s0, mu, sigma, dt, z):
log_moves = (mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z
return s0 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))
rng = np.random.default_rng(1)
z = rng.standard_normal(252)
path = gbm_path(100.0, 0.08, 0.2, 1 / 252, z)
print(path[:4].round(3), round(float(path[-1]), 3))多数の経路：平均と中央値
1年の経路を何千本もシミュレートして、終点がどこになるかを見てください。最終価格の*平均*は期待価格の公式 S0 * exp(mu * T) に近くなります。*中央値*はそれより低く S0 * exp((mu - sigma**2/2) * T) で、ボラティリティが高いと差は大きくなります。少数の驚くほど高い経路が平均を押し上げ、大半の経路は平均より下で終わります。この右に裾の長い対数正規分布が、長期の株価の形です。
import numpy as np
rng = np.random.default_rng(5)
s0, mu, sigma, T, steps, paths = 100.0, 0.08, 0.4, 1.0, 252, 20000
dt = T / steps
z = rng.standard_normal((paths, steps))
final = s0 * np.exp(np.sum((mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z, axis=1))
print(round(float(final.mean()), 1), round(s0 * np.exp(mu * T), 1)) # close to each other: the expected price
print(round(float(np.median(final)), 1), round(s0 * np.exp((mu - 0.5 * sigma ** 2) * T), 1)) # the median is lower
print(round(float(np.mean(final < s0)), 3)) # the share of paths that lose moneymuとsigmaを推定する
dt 年ごとにサンプリングした価格では、対数リターン r = diff(log(prices)) は、平均 (mu - sigma**2/2) * dt、標準偏差 sigma * sqrt(dt) の独立な正規分布になります。したがって sigma = std(r) / sqrt(dt)、mu = mean(r) / dt + sigma**2 / 2 です。ボラティリティの推定は数百点で十分ですが、ドリフトは悪名高くノイズが大きく、標準誤差は観測数ではなく総*期間*でしか縮みません。
import numpy as np
def estimate_gbm(prices, dt):
r = np.diff(np.log(prices))
sigma = r.std(ddof=1) / np.sqrt(dt)
return r.mean() / dt + 0.5 * sigma ** 2, sigma
rng = np.random.default_rng(8)
dt = 1 / 252
for years in (1, 10, 100):
n = int(years * 252)
log_moves = (0.10 - 0.5 * 0.25 ** 2) * dt + 0.25 * np.sqrt(dt) * rng.standard_normal(n)
prices = 100 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))
mu, sigma = estimate_gbm(prices, dt)
print(years, "years:", round(float(mu), 3), round(float(sigma), 3)) # true mu 0.10 and sigma 0.25- GBM：
S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t))。変化率がランダムで価格は正のまま - sigma**2 / 2はexpの凸性を打ち消す。期待価格はS0 * exp(mu * t)だが、中央値の経路はmu - sigma**2 / 2で成長する- 対数リターンから推定する：
sigma = std(r) / sqrt(dt)、mu = mean(r) / dt + sigma**2 / 2。ボラティリティは精度よく、ドリフトは粗くしか推定できない
次は伊藤の補題です。- sigma**2 / 2 がどこから来るのかを説明する連鎖律です。
ラボ：採点される内容
gbm_stepで厳密解の1ステップを進める25点gbm_pathで価格経路をシミュレートする25点gbm_expectedで全ての将来を平均する25点estimate_gbmでパラメータを読み戻す25点