GARCH(1,1)モデルの理論とPython実装：ボラティリティクラスタリングとARCH効果をLjung-Box検定で確認する

はじめに

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）を適用しました。

リターン自体は無相関（p=0.65）なのに、その2乗には極めて強い自己相関（p≈10⁻²⁰¹）があります。 対照実験として、GARCH効果を持たないi.i.d.正規乱数の2乗はp=0.70で無相関のままです。この違いこそが「ARCH効果」（分散の系列相関）であり、ARIMAのような線形・定分散モデルでは捉えられない構造です。GARCHモデルはこの分散自体の時間変化を明示的にモデル化します。

GARCH(1,1)モデル

リターン \(r_t\) を

とし、条件付き分散 \(\sigma_t^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に近いほど、ボラティリティショックの影響が長く持続します。以下、この節では「定常性条件」と「\(\alpha+\beta\) 」の意味を天下り的に受け入れるのではなく、式(2)から実際に導出します。

条件付き分散の可視化

導出に入る前に、式(1)(2)からシミュレートしたリターン系列が実際にどう振る舞うかを確認します。\(\omega=0.05\) , \(\alpha=0.10\) , \(\beta=0.85\) で2000サンプル生成し、真の条件付き分散 \(\sigma_t^2\) と、後述のMLEで推定したパラメータから再構成した条件付き分散 \(\hat\sigma_t^2\) を重ねて描画しました。

上段のリターン系列を見ると、\(t\approx1350\) 〜\(1450\) 付近で振幅が明らかに大きくなり、その前後の静かな期間と対照的です。下段はこれに対応する条件付き分散 \(\sigma_t^2\) の推移で、ショック直後に急上昇し、その後は指数的に減衰しながら無条件分散の水準へ回帰していく様子が見えます。真の \(\sigma_t^2\) （青線）とMLE推定量から再構成した \(\hat\sigma_t^2\) （赤破線）はほぼ完全に重なっており、後述のパラメータ推定が条件付き分散のダイナミクスを正しく捉えていることが視覚的にも確認できます。

定常性条件と無条件分散の導出

式(2)の再帰式が「定常」であるとはどういうことか、そして \(\alpha+\beta<1\) という条件がなぜ必要なのかを、天下り的にではなく導出します。

まず、\(\sigma_t^2\) が（弱）定常であれば、その期待値は時刻によらず一定の値 \(\sigma^2 \equiv E[\sigma_t^2]\) を取るはずです。ここで \(\varepsilon_{t-1} = \sigma_{t-1} z_{t-1}\) かつ \(z_{t-1}\) は \(\sigma_{t-1}^2\) （\(\mathcal{F}_{t-2}\) 可測）と独立に \(\mathcal{N}(0,1)\) に従うので、繰り返し期待値の法則（law of total expectation）より

（\(E[z_{t-1}^2]=1\) であることを使いました）。したがって式(2)の両辺の期待値を取ると

となります。定常性の仮定 \(E[\sigma_t^2]=E[\sigma_{t-1}^2]=\sigma^2\) を代入すると

が得られます。\(\omega>0\) を仮定しているので、この \(\sigma^2\) が有限な正の値として存在するためには \(1-\alpha-\beta>0\) 、すなわち \(\alpha+\beta<1\) が必要です。これが「定常性条件」の意味です。\(\alpha+\beta=1\) の場合は式(3)の右辺が発散し（有限な無条件分散が存在しない：IGARCH、後述のエッジケースで扱います）、\(\alpha+\beta>1\) の場合は \(\sigma^2<0\) となり無意味（そもそも定常な解が存在しない）です。

数値検証：無条件分散の理論値と実測値

式(3)の理論値 \(\omega/(1-\alpha-\beta)\) が、実際にシミュレートした系列とどれだけ一致するかを確認しました。

理論値1.0に対し、有限サンプル（\(n=2000\) ）の標本分散は0.954、条件付き分散パスの平均は0.970と、いずれも数パーセント以内で一致しています（残差は有限サンプルの標本誤差）。

ボラティリティ予測：平均回帰と半減期

式(3)の無条件分散は「\(t\to\infty\) での長期平均」ですが、実務で重要なのは「今日の情報 \(\mathcal{F}_t\) を条件として、\(h\) 期先の分散をどう予測するか」です。これも式(2)から直接導出できます。

\(\mathcal{F}_t\) を所与として1期先を予測すると、\(\varepsilon_t^2 = \sigma_t^2 z_t^2\) で \(z_t\) は \(\mathcal{F}_t\) と独立に \(E[z_t^2]=1\) なので \(E[\varepsilon_t^2 \mid \mathcal{F}_t] = \sigma_t^2\) です。したがって

ここで式(3)より \(\omega = \sigma^2(1-\alpha-\beta)\) を代入すると

という「無条件分散 \(\sigma^2\) への偏差が \((\alpha+\beta)\) 倍される」形になります。この関係を帰納法で \(h\) 期先に拡張すると

が成り立ちます（\(h=1\) で上式に一致し、\(h\to h+1\) の帰納ステップも同じ性質 \(E[\varepsilon_{t+h}^2\mid\mathcal{F}_t]=E[\sigma_{t+h}^2\mid\mathcal{F}_t]\) を使えば同様に示せます）。式(4)は、ボラティリティショックが幾何級数的に減衰しながら無条件分散へ平均回帰することを表しており、その減衰率が持続性 \(\alpha+\beta\) そのものであることを明示的に示しています。半減期（ショックの影響が半分になるまでの期間）は \((\alpha+\beta)^h = 1/2\) を解いて

で与えられます。\(\alpha+\beta=0.95\) のケースでは \(h_{1/2} = \ln(0.5)/\ln(0.95) \approx 13.5\) 期です。

数値検証：モンテカルロ法による平均回帰式の確認

式(4)をモンテカルロ法で検証しました。\(\sigma_t^2=5.0\) （無条件分散1.0から大きく乖離した状態）から2万本のパスをシミュレートし、\(h\) 期先の \(\sigma_{t+h}^2\) の標本平均と式(4)の理論値を比較しました。

モンテカルロ平均と理論値は、2万本のサンプル平均の標本誤差の範囲内（差はいずれも0.02未満）で一致しており、式(4)の平均回帰公式が正しいことが確認できました。\(h=14\) 付近で乖離 \(\sigma_t^2-\sigma^2=4\) がおよそ半分の2に近づいており、理論半減期 \(h_{1/2}\approx13.5\) と整合しています。

2乗残差のARMA(1,1)表現

GARCH(1,1)は、実は2乗残差 \(\varepsilon_t^2\) の系列に対するARMA(1,1)モデルと数学的に等価です。この事実は、https://yuhi-sa.github.io/posts/20260226_arima/1/ で扱ったARMAの枠組みとGARCHがどう結びつくかを明らかにするので、導出しておきます。

\(v_t \equiv \varepsilon_t^2 - \sigma_t^2\) と定義します。\(E[\varepsilon_t^2\mid\mathcal{F}_{t-1}]=\sigma_t^2\) なので、\(v_t\) は \(\mathcal{F}_{t-1}\) を条件として期待値0の誤差項（マルチンゲール差分）です。\(\sigma_t^2 = \varepsilon_t^2 - v_t\) をそのまま式(2)に代入すると

整理すると

となり、これはまさにARMA(1,1)モデルの形です（AR係数が \(\alpha+\beta\) 、MA係数が \(-\beta\) ）。この表現から、GARCH(1,1)の持続性 \(\alpha+\beta\) が2乗残差系列のAR係数と一致すること、そして \(\alpha+\beta \to 1\) が2乗残差系列における単位根（非定常性）に対応することが直接わかります。

数値検証：statsmodels ARIMAによる係数の復元

シミュレートしたリターン系列の2乗 \(r_t^2\) に対し、statsmodels.tsa.arima.model.ARIMA(order=(1,0,1)) でARMA(1,1)モデルを最尤推定し、式(5)の理論係数と比較しました。

MA係数は理論値 \(-\beta=-0.85\) とほぼ完全に一致（差0.0004）し、AR係数も理論値0.95に対し0.9649と近い値が復元できました（AR係数の方がやや誤差が大きいのは、2乗残差のARMA表現の誤差項 \(v_t\) が不均一分散かつ非正規であるため、通常のガウスARMA最尤推定は近似的なものになるからです）。それでも、GARCHの再帰式(2)が2乗残差のARMA(1,1)構造と対応するという理論的関係は数値的に明確に裏付けられました。

対数尤度と最尤推定

ガウス分布を仮定すると、対数尤度は

です。これは 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）を同一データに適用しました。

パラメータ・対数尤度とも高精度で一致しました（残差は最適化アルゴリズムの収束許容誤差レベル）。両者とも真値からはやや外れています（\(\alpha\) の推定値0.079 vs 真値0.10）が、これは2000サンプルという有限標本での推定誤差の範囲内です。\(\alpha+\beta \approx 0.96\) という高い持続性が正しく推定されており、シミュレーションで使った真の持続性 \(0.10+0.85=0.95\) に近い値が復元できています。

エッジケースと注意点

非定常性：IGARCHと発散

式(3)の導出で見た通り、\(\alpha+\beta \geq 1\) では有限な無条件分散が存在しません。特に \(\alpha+\beta=1\) ちょうどの場合は**IGARCH（Integrated GARCH）**と呼ばれ、2乗残差のARMA(1,1)表現（式(5)）でいうところの単位根過程に相当します。実務でよく使われる指数加重移動平均（EWMA、\(\lambda\) -減衰）による分散推定は、実はIGARCHの特殊ケース（\(\omega=0\) 、\(\beta=\lambda\) 、\(\alpha=1-\lambda\) ）に一致します。

\(\omega=0.05\) を固定し、\((\alpha,\beta)\) を変えて条件付き分散パスの挙動を数値的に確認しました。

定常ケースでは前半・後半で分散・平均ともほぼ変わりませんが、IGARCH境界（\(\alpha+\beta=1.00\) ）では後半の分散が前半の約16倍に増大し、平均も約3倍になっています。\(\alpha+\beta=1.05\) の発散域では、わずか2000サンプルで条件付き分散が天文学的な桁数まで爆発しました。これは式(4)の平均回帰公式で \((\alpha+\beta)^h\) が \(h\to\infty\) で収束しない（\(\alpha+\beta\geq1\) ）ことの直接的な帰結です。最尤推定を実装する際は、garch_negloglik 関数内で明示的に alpha + beta >= 1 を弾いている（境界条件のチェック）ことが、この発散を防ぐために不可欠です。

4次モーメント条件と裾の厚さ

GARCH(1,1)過程は、リターン自体がガウス分布に従うと仮定していても、周辺分布（無条件分布）はガウス分布より裾が厚くなる（超過尖度が正）という性質があります。これは、分散自体がランダムに変動することで、大きな \(\lvert\varepsilon_t\rvert\) が生じる確率が単純なガウス分布より高くなるためです。

ガウス分布の仮定の下で、GARCH(1,1)過程の4次モーメントが有限であるための条件（Bollerslev, 1986）は

であり、この条件下での超過尖度は

で与えられます。\(\alpha=0.10,\beta=0.85\) の場合、\(3\alpha^2+2\alpha\beta+\beta^2 = 0.9225 < 1\) で条件を満たし、理論超過尖度は \(\kappa = 0.7742\) です。これをモンテカルロ法（\(n=4000\) の系列を300本シミュレートし、各系列の標本超過尖度を平均）で検証したところ、**平均超過尖度は0.7529（標準誤差0.0356）**となり、理論値0.7742は標準誤差の6倍以内に収まりました。単一系列（\(n=2000\) ）の標本超過尖度は1.1745とやや理論値からずれていましたが、これは尖度のような高次モーメント推定量が単一標本ではばらつきが大きいためで、300回の反復平均を取ることで理論値への収束を確認できました。

レバレッジ効果：GARCH(1,1)の非対称性の限界

式(2)の再帰式は \(\varepsilon_{t-1}^2\) （2乗）にのみ依存するため、正のショックと負のショックが同じ大きさなら条件付き分散への影響も同じという対称性を仮定しています。しかし、株式市場では経験的に「悪いニュース（株価下落）の方が良いニュース（同じ大きさの株価上昇）よりもボラティリティを押し上げやすい」という非対称性（レバレッジ効果）が広く観測されています。標準的なGARCH(1,1)はこの非対称性を捉えられないという構造的な限界があり、これを補うために非対称項を加えたGJR-GARCH（\(\varepsilon_{t-1}<0\) のときに追加のペナルティ項を加える）やEGARCH（対数分散のモデル化により符号の効果を直接組み込む）といった拡張モデルが提案されています。

構造変化と機械学習との比較：最近の研究動向

GARCHモデルはパラメータ \((\omega,\alpha,\beta)\) が全期間で一定であることを仮定しますが、実際の金融市場では規制変更・金融危機・パンデミックなどによる**構造変化（structural break）**でボラティリティのレジーム自体が変わることがあります。Chung, Espinoza, & Quispe (2025)の Forecasting Financial Volatility Under Structural Breaks: A Comparative Study of GARCH Models and Deep Learning Techniques（Journal of Risk and Financial Management 18(9), 494, DOI: 10.3390/jrfm18090494）は、2000〜2024年のラテンアメリカ4か国の株価指数データを用い、修正ICSSアルゴリズムで構造変化点を検出した上でGARCHモデルとLSTM・CNNなどの深層学習モデルを比較しました。構造変化を無視するとボラティリティの持続性（\(\alpha+\beta\) ）が過大評価され予測精度が低下する一方、構造変化をレジーム分割としてGARCHに組み込んでも改善は限定的なケースが多く、中期〜長期の予測ホライズンでは深層学習モデルが非線形かつ時変のダイナミクスをより効果的に捉えて一貫してGARCH系モデルを上回ることが報告されています。この結果は、本記事で導出したGARCH(1,1)のようなパラメータ一定の古典的モデルが、構造変化が頻発する実データに対しては万能ではなく、レジームスイッチングGARCHやGARCH-深層学習ハイブリッドモデルとの併用が近年の実務的な潮流であることを示しています。

まとめ：ARIMAとGARCHの役割分担

https://yuhi-sa.github.io/posts/20260226_arima/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.
• Chung, V., Espinoza, J., & Quispe, R. (2025). Forecasting Financial Volatility Under Structural Breaks: A Comparative Study of GARCH Models and Deep Learning Techniques. Journal of Risk and Financial Management, 18(9), 494. DOI: 10.3390/jrfm18090494