はじめに
https://yuhi-sa.github.io/posts/20260226_arima/1/ と https://yuhi-sa.github.io/posts/20260715_arima_state_space_kalman/1/ で扱ったARIMAモデルは、時系列の平均の構造（自己回帰・移動平均）をモデル化しますが、分散が時間とともに変動するという現象（ボラティリティクラスタリング：荒い値動きの後には荒い値動きが続き、静かな期間には静かな値動きが続く）は捉えられません。本記事では、この分散のクラスタリングを条件付き分散の再帰式としてモデル化するGARCH（Generalized AutoRegressive Conditional Heteroskedasticity）モデルを、Ljung-Box検定でその必要性を確認した上で、フルスクラッチ実装します。
ボラティリティクラスタリングとARCH効果
金融時系列などでよく観測される現象として、リターン自体には自己相関がないのに、リターンの2乗（≒分散の代理指標）には強い自己相関があるというものがあります。これを直接確認するため、GARCH(1,1)過程からシミュレートしたリターン系列に対し、Ljung-Box検定（ラグ10）を適用しました。
| 対象 | Ljung-Box統計量 | p値 | 自己相関 |
|---|---|---|---|
| リターン自体 | 7.77 | 0.651 | なし（想定通り） |
| リターンの2乗 | 966.80 | 2.6×10⁻²⁰¹ | 強く存在 |
| （対照）i.i.d.正規乱数の2乗 | 7.23 | 0.704 | なし |
リターン自体は無相関（p=0.65）なのに、その2乗には極めて強い自己相関（p≈10⁻²⁰¹）があります。 対照実験として、GARCH効果を持たないi.i.d.正規乱数の2乗はp=0.70で無相関のままです。この違いこそが「ARCH効果」（分散の系列相関）であり、ARIMAのような線形・定分散モデルでは捉えられない構造です。GARCHモデルはこの分散自体の時間変化を明示的にモデル化します。
GARCH(1,1)モデル
リターン \(r_t\) を
\[ r_t = \mu + \varepsilon_t, \qquad \varepsilon_t = \sigma_t z_t, \qquad z_t \sim \mathcal{N}(0,1) \tag{1} \]とし、条件付き分散 \(\sigma_t^2\) が次の再帰式に従うとします。
\[ \sigma_t^2 = \omega + \alpha \varepsilon_{t-1}^2 + \beta \sigma_{t-1}^2 \tag{2} \]\(\omega > 0\) 、\(\alpha, \beta \geq 0\) 、\(\alpha + \beta < 1\) （定常性条件）です。式(2)の直感的な意味は、「今日の分散 \(\sigma_t^2\) は、昨日の予測誤差の大きさ \(\varepsilon_{t-1}^2\) （ショックの影響、\(\alpha\) ）と昨日の分散 \(\sigma_{t-1}^2\) （分散自体の持続性、\(\beta\) ）の加重和で決まる」というものです。\(\alpha + \beta\) が1に近いほど、ボラティリティショックの影響が長く持続します。
対数尤度と最尤推定
ガウス分布を仮定すると、対数尤度は
\[ \log L(\omega,\alpha,\beta) = -\frac{1}{2}\sum_{t=1}^{N}\left[\log(2\pi\sigma_t^2) + \frac{\varepsilon_t^2}{\sigma_t^2}\right] \tag{3} \]です。これは https://yuhi-sa.github.io/posts/20260715_arima_state_space_kalman/1/ で導出したカルマンフィルタの予測誤差分解による対数尤度と同じ形をしています（違いは、状態空間モデルの共分散 \(P_{t|t-1}\) の代わりに、GARCH自身の再帰式で条件付き分散 \(\sigma_t^2\) を直接更新する点です）。
Python実装
import numpy as np
from scipy.optimize import minimize
# GARCH(1,1)過程のシミュレーション（検証用データ生成）
def simulate_garch(omega, alpha, beta, n, rng):
eps, sigma2 = np.zeros(n), np.zeros(n)
sigma2[0] = omega / (1 - alpha - beta) # 無条件分散で初期化
z = rng.standard_normal(n)
eps[0] = np.sqrt(sigma2[0]) * z[0]
for t in range(1, n):
sigma2[t] = omega + alpha * eps[t - 1] ** 2 + beta * sigma2[t - 1]
eps[t] = np.sqrt(sigma2[t]) * z[t]
return eps
rng = np.random.default_rng(7)
returns = simulate_garch(omega=0.05, alpha=0.10, beta=0.85, n=2000, rng=rng)
def arch_style_backcast(r, decay=0.94, tau=75):
tau = min(tau, len(r))
w = decay ** np.arange(tau)
w /= w.sum()
return np.sum(w * r[:tau] ** 2)
def garch_negloglik(params, r, sigma2_0):
omega, alpha, beta = params
if omega <= 0 or alpha < 0 or beta < 0 or alpha + beta >= 1:
return 1e10
N = len(r)
sigma2 = np.zeros(N)
sigma2[0] = sigma2_0
for t in range(1, N):
sigma2[t] = omega + alpha * r[t - 1] ** 2 + beta * sigma2[t - 1]
ll = -0.5 * np.sum(np.log(2 * np.pi * sigma2) + r**2 / sigma2)
return -ll
sigma2_0 = arch_style_backcast(returns)
res = minimize(
garch_negloglik, x0=[0.1, 0.05, 0.8], args=(returns, sigma2_0),
method="L-BFGS-B", bounds=[(1e-6, None), (0, 1), (0, 1)],
options={"maxiter": 10000, "ftol": 1e-12, "gtol": 1e-10},
)
omega_hat, alpha_hat, beta_hat = res.x
arch_style_backcast は、初期分散 \(\sigma_0^2\)
を先頭75サンプルの指数加重平均（減衰率0.94）で推定するもので、arch ライブラリの慣例に合わせています（単純なサンプル分散で初期化すると、序盤の尤度計算がわずかに歪みます）。
数値実験：archライブラリとの比較
真のパラメータ \(\omega=0.05\)
, \(\alpha=0.10\)
, \(\beta=0.85\)
でGARCH(1,1)リターンを2000サンプル生成し、上記のフルスクラッチMLEと arch.arch_model（vol="Garch", p=1, q=1）を同一データに適用しました。
| パラメータ | フルスクラッチ | archライブラリ | 差分 |
|---|---|---|---|
| \(\omega\) | 0.036171 | 0.035934 | 2.4×10⁻⁴ |
| \(\alpha\) | 0.078660 | 0.078473 | 1.9×10⁻⁴ |
| \(\beta\) | 0.881839 | 0.882264 | 4.3×10⁻⁴ |
| 対数尤度 | -2678.305 | -2678.354 | 0.05 |
パラメータ・対数尤度とも高精度で一致しました（残差は最適化アルゴリズムの収束許容誤差レベル）。両者とも真値からはやや外れています（\(\alpha\) の推定値0.079 vs 真値0.10）が、これは2000サンプルという有限標本での推定誤差の範囲内です。\(\alpha+\beta \approx 0.96\) という高い持続性が正しく推定されており、シミュレーションで使った真の持続性 \(0.10+0.85=0.95\) に近い値が復元できています。
まとめ：ARIMAとGARCHの役割分担
https://yuhi-sa.github.io/posts/20260226_arima/1/ と本記事は補完関係にあります。
| モデル | 対象 | 捉える構造 |
|---|---|---|
| ARIMA | 平均 \(E[r_t]\) | 自己回帰・移動平均（水準の予測可能性） |
| GARCH | 条件付き分散 \(\text{Var}[r_t \mid \mathcal{F}_{t-1}]\) | ボラティリティクラスタリング（リスクの予測可能性） |
実務では両者を組み合わせた ARIMA-GARCH モデル（平均をARIMAで、残差の分散をGARCHでモデル化）が金融時系列分析の標準的な出発点になります。
おすすめ書籍
ARIMAだけでなくGARCHファミリー（ARCH/GARCH/EGARCH等）も体系的に解説されており、本記事の理論的背景を補完する定番書です。
※ 上記は Amazon アソシエイトのリンクです。
関連記事
- ARIMA・SARIMAの理論とPython実装 - 時系列の平均構造をモデル化する基礎です。GARCHは分散側を補完します。
- ARMA/ARIMAの状態空間表現とカルマンフィルタによる最尤推定 - 本記事のGARCH尤度と同じ「予測誤差分解」の枠組みを、ARIMAの文脈で解説しています。
- 自己相関関数の理論とPython実装 - 本記事のLjung-Box検定・ACFの基礎となる概念です。
- PACF・ARモデル次数同定 - AR次数同定の考え方。GARCHの次数 \((p,q)\) 選定にも類似の情報量規準が使えます。
- 時系列データの異常検知：統計的手法からカルマンフィルタまで - GARCHの条件付き分散は、異常検知における動的な閾値設定にも応用できます。
参考文献
- Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31(3), 307-327.
- Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987-1007.
- Sheppard, K., et al. ARCH: Autoregressive Conditional Heteroskedasticity Models in Python (arch library documentation).
- Ljung, G. M., & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297-303.