GARCH Model Explained

Key Insights

GARCH models capture volatility clustering in financial time series where high-volatility periods cluster together, a phenomenon that simple constant-variance models miss entirely.
The GARCH(1,1) model—with just three parameters—effectively captures the conditional variance dynamics of most financial assets and serves as the industry standard for volatility forecasting.
Proper model validation through standardized residual diagnostics and ARCH-LM tests is critical; a well-fitted GARCH model should eliminate all remaining conditional heteroskedasticity from the residuals.

Introduction to Volatility Modeling

Volatility is the heartbeat of financial markets. It drives option pricing, risk management decisions, and portfolio allocation strategies. Yet most introductory time series courses assume constant variance—a dangerous simplification that ignores the reality of financial data.

Look at any stock return series and you’ll see volatility clustering: periods of high volatility followed by more high volatility, and calm periods that persist. The 2008 financial crisis, the 2020 COVID crash, and the 2022 tech selloff all demonstrate this phenomenon. Traditional models that assume constant variance are blind to these patterns.

This is where GARCH (Generalized Autoregressive Conditional Heteroskedasticity) models come in. They model the variance itself as a time-varying process, capturing the conditional nature of volatility.

import yfinance as yf
import numpy as np
import matplotlib.pyplot as plt

# Download S&P 500 data
data = yf.download('^GSPC', start='2015-01-01', end='2024-01-01')
returns = 100 * data['Adj Close'].pct_change().dropna()

# Plot returns showing volatility clustering
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(returns.index, returns.values)
ax.set_title('S&P 500 Daily Returns - Volatility Clustering', fontsize=14)
ax.set_ylabel('Returns (%)')
ax.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
plt.tight_layout()
plt.show()

The plot reveals distinct periods where volatility spikes and clusters—exactly what GARCH models are designed to capture.

From ARCH to GARCH

Robert Engle introduced the ARCH (Autoregressive Conditional Heteroskedasticity) model in 1982, earning him a Nobel Prize. The ARCH(q) model expresses today’s conditional variance as a function of past squared returns:

σ²ₜ = ω + α₁ε²ₜ₋₁ + α₂ε²ₜ₋₂ + … + αqε²ₜ₋q

The problem? ARCH models often require many lags to adequately capture volatility persistence, leading to parameter proliferation and estimation challenges.

Tim Bollerslev solved this in 1986 by introducing GARCH, which adds lagged conditional variance terms:

σ²ₜ = ω + α₁ε²ₜ₋₁ + … + αqε²ₜ₋q + β₁σ²ₜ₋₁ + … + βpσ²ₜ₋p

The GARCH(p,q) model has p lagged variance terms and q lagged squared residual terms. In practice, GARCH(1,1) works remarkably well for most financial series:

σ²ₜ = ω + αε²ₜ₋₁ + βσ²ₜ₋₁

Here’s a simple ARCH(1) implementation to understand the mechanics:

def arch1_simulate(n, omega, alpha, initial_variance=1.0):
    """Simulate ARCH(1) process"""
    variances = np.zeros(n)
    returns = np.zeros(n)
    variances[0] = initial_variance
    
    for t in range(n):
        if t > 0:
            variances[t] = omega + alpha * returns[t-1]**2
        returns[t] = np.sqrt(variances[t]) * np.random.normal()
    
    return returns, variances

# Simulate ARCH(1)
returns_arch, var_arch = arch1_simulate(1000, omega=0.1, alpha=0.7)

# GARCH(1,1) extends this by adding beta * variance[t-1]
# sigma_t^2 = omega + alpha * returns[t-1]^2 + beta * sigma[t-1]^2

The key difference: GARCH includes the lagged variance term (β), which creates long memory in volatility and reduces the number of parameters needed.

Implementing GARCH in Python

The arch library by Kevin Sheppard provides a robust implementation. Here’s a complete workflow:

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

# Download and prepare data
data = yf.download('SPY', start='2018-01-01', end='2024-01-01')
returns = 100 * data['Adj Close'].pct_change().dropna()

# Specify and fit GARCH(1,1) model
model = arch_model(returns, vol='Garch', p=1, q=1, dist='normal')
results = model.fit(disp='off')

# Display results
print(results.summary())

# Extract parameters
omega = results.params['omega']
alpha = results.params['alpha[1]']
beta = results.params['beta[1]']

print(f"\nEstimated Parameters:")
print(f"ω (omega): {omega:.6f}")
print(f"α (alpha): {alpha:.6f}")
print(f"β (beta):  {beta:.6f}")
print(f"Persistence (α + β): {alpha + beta:.6f}")

The persistence measure (α + β) is crucial. Values close to 1 indicate high volatility persistence—shocks to volatility decay slowly. Most equity indices show persistence around 0.95-0.99.

Model Diagnostics and Validation

A well-specified GARCH model should eliminate all conditional heteroskedasticity from the standardized residuals. Here’s how to verify:

from scipy import stats
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import acorr_ljungbox, het_arch

# Get standardized residuals
std_resid = results.std_resid

# Create diagnostic plots
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. Standardized residuals over time
axes[0, 0].plot(std_resid)
axes[0, 0].set_title('Standardized Residuals')
axes[0, 0].axhline(y=0, color='black', linestyle='--', linewidth=0.5)

# 2. Q-Q plot for normality
stats.probplot(std_resid.dropna(), dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot')

# 3. ACF of standardized residuals
plot_acf(std_resid.dropna(), lags=30, ax=axes[1, 0])
axes[1, 0].set_title('ACF of Standardized Residuals')

# 4. ACF of squared standardized residuals
plot_acf(std_resid.dropna()**2, lags=30, ax=axes[1, 1])
axes[1, 1].set_title('ACF of Squared Standardized Residuals')

plt.tight_layout()
plt.show()

# Statistical tests
# Ljung-Box test for autocorrelation in residuals
lb_test = acorr_ljungbox(std_resid.dropna(), lags=[10], return_df=True)
print("\nLjung-Box Test (residuals):")
print(lb_test)

# ARCH-LM test for remaining ARCH effects
arch_lm = het_arch(std_resid.dropna(), nlags=10)
print(f"\nARCH-LM Test:")
print(f"LM Statistic: {arch_lm[0]:.4f}")
print(f"P-value: {arch_lm[1]:.4f}")

The ARCH-LM test is critical. A p-value above 0.05 suggests no remaining ARCH effects—exactly what we want. If the test rejects, consider increasing the GARCH order or trying alternative specifications.

Forecasting Volatility

GARCH models excel at volatility forecasting, which is essential for VaR calculations, option pricing, and portfolio optimization:

# Generate forecasts
forecast_horizon = 30
forecasts = results.forecast(horizon=forecast_horizon, reindex=False)

# Extract variance forecasts
forecast_variance = forecasts.variance.values[-1, :]
forecast_volatility = np.sqrt(forecast_variance)

# Create forecast dates
last_date = returns.index[-1]
forecast_dates = pd.date_range(start=last_date, periods=forecast_horizon+1)[1:]

# Plot historical and forecasted volatility
fig, ax = plt.subplots(figsize=(14, 6))

# Historical conditional volatility
historical_vol = results.conditional_volatility
ax.plot(historical_vol.index[-252:], historical_vol.values[-252:], 
        label='Historical Volatility', color='blue')

# Forecasted volatility
ax.plot(forecast_dates, forecast_volatility, 
        label='Forecasted Volatility', color='red', linestyle='--')

# Add confidence bands (approximate)
upper_band = forecast_volatility * 1.2
lower_band = forecast_volatility * 0.8
ax.fill_between(forecast_dates, lower_band, upper_band, 
                 alpha=0.2, color='red', label='Forecast Band')

ax.set_title('GARCH(1,1) Volatility Forecast', fontsize=14)
ax.set_ylabel('Volatility (%)')
ax.legend()
plt.tight_layout()
plt.show()

Notice how GARCH forecasts mean-revert to the long-run variance. This is a fundamental property: short-term forecasts reflect recent volatility, while long-term forecasts converge to the unconditional variance.

Extensions and Practical Considerations

GARCH(1,1) assumes symmetric volatility response—positive and negative shocks of equal magnitude have the same impact. But equity markets exhibit leverage effects: negative returns increase volatility more than positive returns of the same size.

GJR-GARCH addresses this asymmetry:

# Compare GARCH(1,1) vs GJR-GARCH
model_garch = arch_model(returns, vol='Garch', p=1, q=1)
results_garch = model_garch.fit(disp='off')

model_gjr = arch_model(returns, vol='Garch', p=1, o=1, q=1)
results_gjr = model_gjr.fit(disp='off')

print("Model Comparison:")
print(f"GARCH(1,1) AIC: {results_garch.aic:.2f}")
print(f"GJR-GARCH AIC:  {results_gjr.aic:.2f}")

# Lower AIC indicates better fit
if results_gjr.aic < results_garch.aic:
    print("\nGJR-GARCH preferred - asymmetric effects present")
else:
    print("\nGARCH(1,1) sufficient - symmetric response")

Other useful extensions include:

EGARCH: Models log-variance, ensuring positivity without parameter constraints
TGARCH: Threshold GARCH with explicit asymmetry terms
FIGARCH: Captures fractional integration for extreme persistence

For model selection, use AIC or BIC. Start with GARCH(1,1)—it works for most applications. Only move to more complex specifications if diagnostics reveal problems or if information criteria clearly favor alternatives.

The key to successful GARCH modeling isn’t finding the most sophisticated variant. It’s understanding your data, validating your assumptions, and ensuring the model serves its intended purpose—whether that’s risk measurement, forecasting, or derivative pricing. Keep it simple until simplicity fails.