時系列ボラティリティ予測実践ガイド：GARCHモデルファミリーとPython実装

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt
from scipy import stats

# KOSPI 200 ETF（KODEX 200）またはS&P 500のデータ収集
spy = yf.download('SPY', start='2020-01-01', end='2026-03-01', auto_adjust=True)
returns = spy['Close'].pct_change().dropna() * 100  # パーセントリターン

# スタイライズドファクトの検証
print(f"Mean: {returns.mean():.4f}%")
print(f"Std Dev: {returns.std():.4f}%")
print(f"Skewness: {returns.skew():.4f}")
print(f"Kurtosis: {returns.kurtosis():.4f}")  # 超過尖度（正規分布=0）

# Jarque-Bera正規性検定
jb_stat, jb_pval = stats.jarque_bera(returns)
print(f"Jarque-Bera statistic: {jb_stat:.2f}, p-value: {jb_pval:.6f}")

# Ljung-Box検定：リターン vs 二乗リターン
from statsmodels.stats.diagnostic import acorr_ljungbox

lb_returns = acorr_ljungbox(returns, lags=10, return_df=True)
lb_squared = acorr_ljungbox(returns**2, lags=10, return_df=True)
print(f"\nReturns Ljung-Box p-value (lag=10): {lb_returns['lb_pvalue'].iloc[-1]:.4f}")
print(f"Returns^2 Ljung-Box p-value (lag=10): {lb_squared['lb_pvalue'].iloc[-1]:.6f}")
# 二乗リターンのp値が非常に小さければボラティリティクラスタリングの存在を示す

r_t = mu + epsilon_t
epsilon_t = sigma_t * z_t,  z_t ~ N(0,1)

sigma_t^2 = omega + alpha_1 * epsilon_{t-1}^2 + ... + alpha_q * epsilon_{t-q}^2

sigma_t^2 = omega + alpha * epsilon_{t-1}^2 + beta * sigma_{t-1}^2

sigma_long^2 = omega / (1 - alpha - beta)

from arch import arch_model
import yfinance as yf
import pandas as pd
import numpy as np

# データ準備
spy = yf.download('SPY', start='2018-01-01', end='2026-03-01', auto_adjust=True)
returns = spy['Close'].pct_change().dropna() * 100  # パーセントリターン

# Student-t分布によるGARCH(1,1)
garch11 = arch_model(
    returns,
    vol='Garch',
    p=1,          # GARCH次数（過去の分散）
    q=1,          # ARCH次数（過去のショック）
    mean='Constant',
    dist='t'      # Student-t分布（ファットテールを捕捉）
)
result = garch11.fit(disp='off')
print(result.summary())

# パラメータ抽出
omega = result.params['omega']
alpha = result.params['alpha[1]']
beta = result.params['beta[1]']
persistence = alpha + beta

print(f"\nomega: {omega:.6f}")
print(f"alpha: {alpha:.4f}")
print(f"beta: {beta:.4f}")
print(f"persistence (alpha+beta): {persistence:.4f}")

# 長期非条件付き分散と年率換算ボラティリティ
long_run_var = omega / (1 - alpha - beta)
annual_vol = np.sqrt(long_run_var * 252) / 100  # 小数比率に変換
print(f"Long-run annualized volatility: {annual_vol:.2%}")

# 条件付きボラティリティ時系列
cond_vol = result.conditional_volatility
print(f"\nConditional volatility - last 5 days:")
print(cond_vol.tail())

sigma_t^2 = omega + (alpha + gamma * I_{t-1}) * epsilon_{t-1}^2 + beta * sigma_{t-1}^2

# Student-t分布によるGJR-GARCH(1,1,1)
gjr = arch_model(
    returns,
    vol='Garch',
    p=1, o=1, q=1,   # o=1がGJRの非対称項
    mean='Constant',
    dist='t'
)
gjr_result = gjr.fit(disp='off')
print(gjr_result.summary())

gamma = gjr_result.params['gamma[1]']
alpha_gjr = gjr_result.params['alpha[1]']

print(f"\nalpha: {alpha_gjr:.4f}")
print(f"gamma: {gamma:.4f}")
print(f"Positive shock impact: {alpha_gjr:.4f}")
print(f"Negative shock impact: {alpha_gjr + gamma:.4f}")
print(f"Asymmetry ratio: {(alpha_gjr + gamma) / alpha_gjr:.2f}x")

ln(sigma_t^2) = omega + alpha * (|z_{t-1}| - E|z_{t-1}|) + gamma * z_{t-1} + beta * ln(sigma_{t-1}^2)

# 歪みStudent-t分布によるEGARCH(1,1)
egarch = arch_model(
    returns,
    vol='EGARCH',
    p=1, q=1,
    mean='Constant',
    dist='skewt'   # 歪みStudent-t分布
)
egarch_result = egarch.fit(disp='off')
print(egarch_result.summary())

# EGARCHの非対称パラメータ確認
# gamma < 0 は負のショックがボラティリティをより増加させることを意味
gamma_egarch = egarch_result.params['gamma[1]']
print(f"\ngamma (asymmetry): {gamma_egarch:.4f}")
if gamma_egarch < 0:
    print("レバレッジ効果確認：負のショックがボラティリティをより増加させる")

# コアパッケージのインストール
pip install arch yfinance pandas numpy scipy matplotlib statsmodels

# オプションパッケージ
pip install seaborn quantstats

import numpy as np
import pandas as pd
import yfinance as yf
from arch import arch_model
from arch.unitroot import ADF, KPSS
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# --------------------------------------------------
# ステップ1：データ収集と前処理
# --------------------------------------------------
def prepare_returns(ticker: str, start: str, end: str) -> pd.Series:
    """リターンデータの収集と前処理"""
    data = yf.download(ticker, start=start, end=end, auto_adjust=True)
    prices = data['Close']
    log_returns = np.log(prices / prices.shift(1)).dropna() * 100
    return log_returns

returns = prepare_returns('SPY', '2018-01-01', '2026-03-01')

# --------------------------------------------------
# ステップ2：事前検定
# --------------------------------------------------
def pre_tests(returns: pd.Series) -> dict:
    """GARCHモデル適用前に必要な検定を実施"""
    results = {}

    # ADF検定：リターンが定常過程かどうかを確認
    adf = ADF(returns, lags=5)
    results['ADF_stat'] = adf.stat
    results['ADF_pvalue'] = adf.pvalue
    results['is_stationary'] = adf.pvalue < 0.05

    # ARCH効果検定（EngleのLM検定）
    from statsmodels.stats.diagnostic import het_arch
    lm_stat, lm_pvalue, _, _ = het_arch(returns, nlags=10)
    results['ARCH_LM_stat'] = lm_stat
    results['ARCH_LM_pvalue'] = lm_pvalue
    results['has_arch_effect'] = lm_pvalue < 0.05

    return results

tests = pre_tests(returns)
for key, val in tests.items():
    print(f"{key}: {val}")

# ADF p値が0.05未満（定常）+ ARCH LM p値が0.05未満（ARCH効果あり）
# 両条件を満たして初めてGARCHフィッティングに意味がある

# --------------------------------------------------
# ステップ3：複数モデルのフィッティング
# --------------------------------------------------
def fit_models(returns: pd.Series) -> dict:
    """複数のGARCHバリアントをフィッティングして結果を返す"""
    models = {}

    # 1) GARCH(1,1) - Student-t
    am = arch_model(returns, vol='Garch', p=1, q=1,
                    mean='Constant', dist='t')
    models['GARCH(1,1)-t'] = am.fit(disp='off')

    # 2) GJR-GARCH(1,1,1) - Student-t
    am = arch_model(returns, vol='Garch', p=1, o=1, q=1,
                    mean='Constant', dist='t')
    models['GJR-GARCH-t'] = am.fit(disp='off')

    # 3) EGARCH(1,1) - 歪みStudent-t
    am = arch_model(returns, vol='EGARCH', p=1, q=1,
                    mean='Constant', dist='skewt')
    models['EGARCH-skewt'] = am.fit(disp='off')

    # 4) GARCH(1,1) - 正規分布（比較用）
    am = arch_model(returns, vol='Garch', p=1, q=1,
                    mean='Constant', dist='normal')
    models['GARCH(1,1)-N'] = am.fit(disp='off')

    return models

models = fit_models(returns)

# --------------------------------------------------
# ステップ4：モデル比較（情報量基準）
# --------------------------------------------------
comparison = pd.DataFrame({
    name: {
        'AIC': res.aic,
        'BIC': res.bic,
        'Log-Likelihood': res.loglikelihood,
        'Params': res.num_params
    }
    for name, res in models.items()
}).T
comparison = comparison.sort_values('BIC')
print("\nModel comparison (sorted by BIC):")
print(comparison.to_string())

def diagnose_model(result, model_name: str = ""):
    """GARCHモデルのフィット品質を診断"""
    std_resid = result.std_resid
    print(f"=== {model_name} Diagnostics ===")

    # 1. 標準化残差の記述統計
    print(f"Mean: {std_resid.mean():.4f} (expected: 0)")
    print(f"Std Dev: {std_resid.std():.4f} (expected: 1)")
    print(f"Skewness: {std_resid.skew():.4f}")
    print(f"Kurtosis: {std_resid.kurtosis():.4f}")

    # 2. Ljung-Box検定：残差の自己相関
    from statsmodels.stats.diagnostic import acorr_ljungbox
    lb_resid = acorr_ljungbox(std_resid, lags=10, return_df=True)
    lb_sq = acorr_ljungbox(std_resid**2, lags=10, return_df=True)

    print(f"\nLjung-Box (residuals, lag=10) p-value: "
          f"{lb_resid['lb_pvalue'].iloc[-1]:.4f}")
    print(f"Ljung-Box (residuals^2, lag=10) p-value: "
          f"{lb_sq['lb_pvalue'].iloc[-1]:.4f}")

    # p値 > 0.05 は自己相関なし = モデルが構造を十分に捕捉

    # 3. ARCH-LM検定：残差のARCH効果
    from statsmodels.stats.diagnostic import het_arch
    lm_stat, lm_pval, _, _ = het_arch(std_resid, nlags=10)
    print(f"ARCH-LM (residual effect) p-value: {lm_pval:.4f}")

    # p値 > 0.05 は残差にARCH効果なし = 分散構造を十分に捕捉

    # 4. 正規性 vs t分布のフィット
    jb_stat, jb_pval = stats.jarque_bera(std_resid)
    print(f"Jarque-Bera p-value: {jb_pval:.6f}")

    print()
    return {
        'lb_resid_pval': lb_resid['lb_pvalue'].iloc[-1],
        'lb_sq_pval': lb_sq['lb_pvalue'].iloc[-1],
        'arch_lm_pval': lm_pval,
        'jb_pval': jb_pval
    }

# 最良モデルの診断
best_model_name = comparison.index[0]  # BIC最小のモデル
best_result = models[best_model_name]
diag = diagnose_model(best_result, best_model_name)

def calculate_garch_var(
    returns: pd.Series,
    model_result,
    confidence_levels: list = [0.01, 0.05],
    investment: float = 1_000_000
) -> pd.DataFrame:
    """
    GARCHモデルに基づくVaRとCVaR（ES）の計算
    """
    cond_vol = model_result.conditional_volatility
    mu = model_result.params.get('mu', 0)

    dist = model_result.model.distribution
    params = model_result.params

    records = []
    for alpha in confidence_levels:
        if 'nu' in params.index:
            nu = params['nu']
            q = stats.t.ppf(alpha, df=nu)
        else:
            q = stats.norm.ppf(alpha)

        var_pct = -(mu + cond_vol * q) / 100
        var_amount = var_pct * investment

        if 'nu' in params.index:
            nu = params['nu']
            pdf_q = stats.t.pdf(q, df=nu)
            cvar_factor = -(
                (nu + q**2) / (nu - 1)
            ) * pdf_q / alpha
        else:
            pdf_q = stats.norm.pdf(q)
            cvar_factor = -pdf_q / alpha

        cvar_pct = -(mu + cond_vol * cvar_factor) / 100
        cvar_amount = cvar_pct * investment

        for date in cond_vol.index:
            records.append({
                'date': date,
                'alpha': alpha,
                'confidence': f"{(1-alpha)*100:.0f}%",
                'VaR_pct': var_pct.loc[date],
                'VaR_amount': var_amount.loc[date],
                'CVaR_pct': cvar_pct.loc[date],
                'CVaR_amount': cvar_amount.loc[date],
            })

    return pd.DataFrame(records)

var_df = calculate_garch_var(returns, best_result)

recent = var_df[var_df['alpha'] == 0.01].tail(5)
print("直近5日間の99% VaR（投資額100万ウォン基準）:")
print(recent[['date', 'confidence', 'VaR_pct', 'VaR_amount',
              'CVaR_pct', 'CVaR_amount']].to_string(index=False))

def backtest_var(
    returns: pd.Series,
    model_type: str = 'Garch',
    p: int = 1,
    o: int = 0,
    q: int = 1,
    dist: str = 't',
    window: int = 1000,
    alpha: float = 0.01
) -> dict:
    """
    ローリングウィンドウVaRバックテストの実行
    """
    violations = []
    var_series = []
    actual_series = []

    n_forecasts = len(returns) - window
    print(f"バックテスト開始：{n_forecasts}日分の予測")

    for i in range(window, len(returns)):
        train = returns.iloc[i-window:i]
        actual = returns.iloc[i]

        try:
            am = arch_model(train, vol=model_type,
                           p=p, o=o, q=q,
                           mean='Constant', dist=dist)
            res = am.fit(disp='off', show_warning=False)

            forecast = res.forecast(horizon=1)
            mean_f = forecast.mean.iloc[-1, 0]
            var_f = forecast.variance.iloc[-1, 0]
            vol_f = np.sqrt(var_f)

            if 'nu' in res.params.index:
                nu = res.params['nu']
                quantile = stats.t.ppf(alpha, df=nu)
            else:
                quantile = stats.norm.ppf(alpha)

            var_value = -(mean_f + vol_f * quantile)

            is_violation = actual < -var_value

            violations.append(is_violation)
            var_series.append(var_value)
            actual_series.append(actual)
        except Exception:
            violations.append(False)
            var_series.append(np.nan)
            actual_series.append(actual)

    violations = np.array(violations)
    n_violations = violations.sum()
    violation_rate = n_violations / len(violations)
    expected_rate = alpha

    n = len(violations)
    x = n_violations

    if x == 0 or x == n:
        kupiec_stat = np.nan
        kupiec_pval = np.nan
    else:
        lr = -2 * (
            np.log((1 - alpha)**(n - x) * alpha**x)
            - np.log((1 - x/n)**(n - x) * (x/n)**x)
        )
        kupiec_pval = 1 - stats.chi2.cdf(lr, df=1)
        kupiec_stat = lr

    result = {
        'n_forecasts': n,
        'n_violations': int(n_violations),
        'violation_rate': violation_rate,
        'expected_rate': expected_rate,
        'kupiec_stat': kupiec_stat,
        'kupiec_pval': kupiec_pval,
        'model_accepted': kupiec_pval > 0.05 if not np.isnan(kupiec_pval) else None,
    }

    print(f"\nバックテスト結果:")
    print(f"  期待超過率: {expected_rate:.2%}")
    print(f"  実際超過率: {violation_rate:.2%}")
    print(f"  超過回数: {n_violations}/{n}")
    print(f"  Kupiec LR統計量: {kupiec_stat:.4f}")
    print(f"  Kupiec p値: {kupiec_pval:.4f}")
    print(f"  モデル採択: {'はい' if result['model_accepted'] else 'いいえ'}")

    return result

bt_result = backtest_var(
    returns,
    model_type='Garch',
    p=1, o=1, q=1,   # GJR-GARCH
    dist='t',
    window=1000,
    alpha=0.01
)

ConvergenceWarning: Maximum number of iterations has been exceeded.

# 1. リターンのスケールを確認（パーセンテージであるべき、小数ではない）
returns_pct = returns * 100  # 小数 -> パーセンテージ

# 2. 異なる最適化方法を試す
result = model.fit(
    disp='off',
    options={'maxiter': 500},
    starting_values=np.array([0.01, 0.05, 0.90])  # omega, alpha, beta
)

# 推定期間を短縮（構造変化後のデータのみ使用）
# 例：COVID後のデータのみ使用
returns_post_covid = returns['2020-07-01':]
am = arch_model(returns_post_covid, vol='Garch', p=1, q=1,
                mean='Constant', dist='t')
result = am.fit(disp='off')
persistence = result.params['alpha[1]'] + result.params['beta[1]']
print(f"Persistence: {persistence:.4f}")

# 自由度を固定
am = arch_model(returns, vol='Garch', p=1, q=1,
                mean='Constant', dist='t')
# 自由度を固定（例：nu=5）
constraints = {'nu': 5.0}
result = am.fit(disp='off')
# またはGED分布を代替として使用
am_ged = arch_model(returns, vol='Garch', p=1, q=1,
                    mean='Constant', dist='ged')
result_ged = am_ged.fit(disp='off')

# 悪い例：正規分布仮定
bad_model = arch_model(returns, vol='Garch', p=1, q=1,
                       mean='Constant', dist='normal')

# 良い例：Student-t分布をデフォルトに
good_model = arch_model(returns, vol='Garch', p=1, q=1,
                        mean='Constant', dist='t')