시계열 변동성 예측 실전 가이드: 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  # 백분율 수익률

# Stylized Facts 확인
print(f"평균: {returns.mean():.4f}%")
print(f"표준편차: {returns.std():.4f}%")
print(f"왜도(Skewness): {returns.skew():.4f}")
print(f"첨도(Kurtosis): {returns.kurtosis():.4f}")  # 정규분포 = 0 (excess kurtosis)

# Jarque-Bera 정규성 검정
jb_stat, jb_pval = stats.jarque_bera(returns)
print(f"Jarque-Bera 통계량: {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"\n수익률 Ljung-Box p-value (lag=10): {lb_returns['lb_pvalue'].iloc[-1]:.4f}")
print(f"수익률^2 Ljung-Box p-value (lag=10): {lb_squared['lb_pvalue'].iloc[-1]:.6f}")
# 수익률 제곱의 p-value가 매우 작으면 변동성 군집 존재

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  # 백분율 수익률

# GARCH(1,1) with Student-t 분포
garch11 = arch_model(
    returns,
    vol='Garch',
    p=1,          # GARCH 차수 (과거 분산)
    q=1,          # ARCH 차수 (과거 충격)
    mean='Constant',
    dist='t'      # Student-t 분포 (fat tails 포착)
)
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"장기 연율화 변동성: {annual_vol:.2%}")

# 조건부 변동성 시계열
cond_vol = result.conditional_volatility
print(f"\n조건부 변동성 - 최근 5일:")
print(cond_vol.tail())

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

# GJR-GARCH(1,1,1) with Student-t
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"양의 충격 영향: {alpha_gjr:.4f}")
print(f"음의 충격 영향: {alpha_gjr + gamma:.4f}")
print(f"비대칭 비율: {(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)

# EGARCH(1,1) with Skewed Student-t
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 (비대칭): {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-value < 0.05 (정상) + ARCH LM p-value < 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) - Skewed 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) - Normal (비교용)
    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("\n모델 비교 (BIC 기준 정렬):")
print(comparison.to_string())

def diagnose_model(result, model_name: str = ""):
    """GARCH 모델의 적합도를 진단한다."""
    std_resid = result.std_resid
    print(f"=== {model_name} 진단 ===")

    # 1. 표준화 잔차의 기술 통계
    print(f"평균: {std_resid.mean():.4f} (기대값: 0)")
    print(f"표준편차: {std_resid.std():.4f} (기대값: 1)")
    print(f"왜도: {std_resid.skew():.4f}")
    print(f"첨도: {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 (잔차, lag=10) p-value: "
          f"{lb_resid['lb_pvalue'].iloc[-1]:.4f}")
    print(f"Ljung-Box (잔차^2, lag=10) p-value: "
          f"{lb_sq['lb_pvalue'].iloc[-1]:.4f}")

    # p-value > 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 (잔여 효과) p-value: {lm_pval:.4f}")

    # p-value > 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)를 계산한다.

    Parameters
    ----------
    returns : pd.Series
        수익률 시계열 (백분율)
    model_result : arch result
        적합된 GARCH 모델 결과
    confidence_levels : list
        VaR 신뢰 수준 (0.01 = 99%, 0.05 = 95%)
    investment : float
        투자 금액 (원 단위)

    Returns
    -------
    pd.DataFrame
        날짜별 VaR, CVaR 값
    """
    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:
        # 분포의 분위수 (quantile)
        # Student-t의 경우 자유도를 반영
        if 'nu' in params.index:
            nu = params['nu']
            q = stats.t.ppf(alpha, df=nu)
        else:
            q = stats.norm.ppf(alpha)

        # VaR = -(mu + sigma * quantile)
        var_pct = -(mu + cond_vol * q) / 100
        var_amount = var_pct * investment

        # CVaR (Conditional VaR / Expected Shortfall)
        # E[X | X < VaR] 계산
        if 'nu' in params.index:
            nu = params['nu']
            # Student-t의 CVaR 공식
            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 계산 실행
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:
    """
    Rolling window 방식의 VaR 백테스팅을 수행한다.

    Parameters
    ----------
    returns : pd.Series
        수익률 시계열 (백분율)
    window : int
        추정 윈도우 크기 (거래일)
    alpha : float
        VaR 신뢰 수준

    Returns
    -------
    dict
        백테스팅 결과 (초과 횟수, 비율, Kupiec 검정 등)
    """
    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)

            # 1일 앞 예측
            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)

            # VaR 계산
            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

    # Kupiec POF 검정 (우도비 검정)
    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-value: {kupiec_pval:.4f}")
    print(f"  모델 수용: {'Yes' if result['model_accepted'] else 'No'}")

    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. 수익률 스케일 확인 (0.01 단위가 아닌 백분율로)
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-19 이후 데이터만 사용
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')