ボラティリティは集中する
日次リターンはほとんど予測できませんが、その大きさは予測できます。今日大きく動けば明日も大きく動く可能性が高く、穏やかな日は続けて現れます。このボラティリティ・クラスタリングは金融で最も確立された事実の一つで、単純なVaRの背後にある一定のボラティリティという仮定を崩します。GARCH（一般化自己回帰条件付き不均一分散）は条件付き分散を直接モデル化するので、リスクやオプション価格が市場のムードに追従できます。
import numpy as np
rng = np.random.default_rng(4)
n = 3000
omega, alpha, beta = 2e-6, 0.10, 0.85
var = np.empty(n)
returns = np.empty(n)
var[0] = omega / (1 - alpha - beta)
returns[0] = np.sqrt(var[0]) * rng.normal()
for t in range(1, n):
var[t] = omega + alpha * returns[t - 1] ** 2 + beta * var[t - 1]
returns[t] = np.sqrt(var[t]) * rng.normal()
# Clustering: the size of a return today predicts the size tomorrow, though the sign is a coin flip
print(round(float(np.corrcoef(returns[:-1], returns[1:])[0, 1]), 3)) # about 0
print(round(float(np.corrcoef(np.abs(returns[:-1]), np.abs(returns[1:]))[0, 1]), 3)) # clearly positiveGARCH(1,1)の漸化式
モデルのパラメータは3つで、var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1] です。omega は下限、alpha は前日のショックへの反応、beta は前日の分散の持続性です。各日が前日に依存するため、初期分散から出発してループで系列を作ります。日次の株式では alpha が0.05〜0.1、beta が0.85〜0.95程度が典型で、ショックはゆっくり薄れます。
import numpy as np
def garch_variance(returns, omega, alpha, beta, initial_var):
var = np.empty(len(returns))
var[0] = initial_var
for t in range(1, len(returns)):
var[t] = omega + alpha * returns[t - 1] ** 2 + beta * var[t - 1]
return var
returns = np.array([0.01, -0.02, 0.015, 0.0, -0.03, 0.01, 0.0, 0.0])
var = garch_variance(returns, 1e-6, 0.08, 0.9, 1e-4)
print(np.round(np.sqrt(var) * 100, 3)) # daily volatility in %: it jumps after the -3% day, then fades落ち着く水準と予測
alpha + beta < 1 なら分散は平均回帰し、長期水準 omega / (1 - alpha - beta) に落ち着きます。1以上なら落ち着きません。予測では、明日の分散は今日のリターンと分散から分かります。その先は期待される二乗リターンが分散に等しいので、各ステップは omega + (alpha + beta) * 直前 となり、予測は今日の水準から長期水準へ、持続性 alpha + beta で決まる速さで近づきます。
def long_run_variance(omega, alpha, beta):
if alpha + beta >= 1:
raise ValueError("alpha + beta must be below 1")
return omega / (1 - alpha - beta)
def garch_forecast(omega, alpha, beta, last_var, last_return, steps):
out = [omega + alpha * last_return ** 2 + beta * last_var]
for _ in range(steps - 1):
out.append(omega + (alpha + beta) * out[-1])
return out
omega, alpha, beta = 1e-6, 0.08, 0.9
print(long_run_variance(omega, alpha, beta) ** 0.5) # about 0.00707, 0.71% a day
forecast = garch_forecast(omega, alpha, beta, 4e-4, -0.03, 60) # after a -3% day in a volatile spell
print([round(v ** 0.5 * 100, 3) for v in (forecast[0], forecast[9], forecast[29], forecast[59])])パラメータの選び方：尤度
パラメータは未知で、データから最尤法で推定します。リターンが分散 var[t] の正規分布に従うなら、1日の対数尤度は -0.5 * (log(2 * pi * var[t]) + r[t]**2 / var[t]) です。これを合計すると、(omega, alpha, beta) の選び方を採点する1つの数が得られ、最良のパラメータは観測されたリターンを最も起こりやすくします。オプティマイザ（例：scipy.optimize.minimize）が多くの組を試し、最高のものを残します。
# Score three parameter choices on data simulated with alpha = 0.10, beta = 0.85
import numpy as np
rng = np.random.default_rng(4)
n = 3000
true = (2e-6, 0.10, 0.85)
var = np.empty(n); returns = np.empty(n)
var[0] = true[0] / (1 - true[1] - true[2]); returns[0] = np.sqrt(var[0]) * rng.normal()
for t in range(1, n):
var[t] = true[0] + true[1] * returns[t - 1] ** 2 + true[2] * var[t - 1]
returns[t] = np.sqrt(var[t]) * rng.normal()
def loglik(params):
omega, alpha, beta = params
v = np.empty(n)
v[0] = returns.var()
for t in range(1, n):
v[t] = omega + alpha * returns[t - 1] ** 2 + beta * v[t - 1]
return float(np.sum(-0.5 * (np.log(2 * np.pi * v) + returns ** 2 / v)))
for params in ((2e-6, 0.10, 0.85), (2e-6, 0.30, 0.50), (1e-5, 0.0, 0.0)):
print(params, round(loglik(params), 1)) # the true parameters score highest- ボラティリティは集中する。GARCH(1,1)は分散をループで作る：
var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1] alpha + beta < 1ならomega / (1 - alpha - beta)に戻る。予測はomega + (alpha + beta) * 直前を使う- パラメータは最尤法
sum(-0.5 * (log(2 * pi * var) + r**2 / var))で選ぶ。裾にはファットテールのショックを使う
次はブラウン運動です。オプション価格モデルの基礎にある、連続時間のランダムウォークです。
ラボ：採点される内容
garch_varianceで漸化式を回す25点long_run_varianceで落ち着く水準を求める25点garch_forecastで次の数日を予測する25点garch_loglikで当てはまりを採点する25点