GARCH(1,1) in Python: Modeling Volatility Clustering and Confirming the ARCH Effect with the Ljung-Box Test

Introduction

The ARIMA models covered in https://yuhi-sa.github.io/en/posts/20260226_arima/1/ and https://yuhi-sa.github.io/en/posts/20260715_arima_state_space_kalman/1/ capture a time series’ mean structure (autoregression, moving average), but not the phenomenon where variance itself changes over time — volatility clustering, where turbulent periods are followed by more turbulence and calm periods stay calm. This article models that variance clustering as a conditional-variance recursion, the GARCH (Generalized AutoRegressive Conditional Heteroskedasticity) model, first confirming its necessity with a Ljung-Box test, then implementing it from scratch.

Volatility Clustering and the ARCH Effect

A phenomenon commonly observed in financial time series: returns themselves show no autocorrelation, but squared returns (a proxy for variance) show strong autocorrelation. To confirm this directly, we applied the Ljung-Box test (lag 10) to a return series simulated from a GARCH(1,1) process.

The returns themselves are uncorrelated (p=0.65), yet their squares show extremely strong autocorrelation (p≈10⁻²⁰¹). As a control, squared i.i.d. Gaussian noise (no GARCH effect) stays uncorrelated at p=0.70. This distinction is exactly the “ARCH effect” (serial correlation in variance) — a structure that linear, constant-variance models like ARIMA cannot capture. GARCH models this time-varying nature of the variance explicitly.

The GARCH(1,1) Model

Let the return \(r_t\) be

with the conditional variance \(\sigma_t^2\) following the recursion

where \(\omega > 0\) , \(\alpha, \beta \geq 0\) , and \(\alpha + \beta < 1\) (the stationarity condition). Equation (2) says, intuitively: “today’s variance \(\sigma_t^2\) is a weighted combination of yesterday’s shock magnitude \(\varepsilon_{t-1}^2\) (the shock’s impact, weighted by \(\alpha\) ) and yesterday’s variance \(\sigma_{t-1}^2\) (persistence of the variance level itself, weighted by \(\beta\) ).” The closer \(\alpha + \beta\) is to 1, the longer volatility shocks persist.

Log-Likelihood and Maximum Likelihood Estimation

Under a Gaussian assumption, the log-likelihood is

This has the same form as the log-likelihood derived from the Kalman filter’s prediction-error decomposition in https://yuhi-sa.github.io/en/posts/20260715_arima_state_space_kalman/1/ (the difference: instead of the state-space model’s covariance \(P_{t|t-1}\) , GARCH’s own recursion updates the conditional variance \(\sigma_t^2\) directly).

Python Implementation

import numpy as np
from scipy.optimize import minimize

# Simulate a GARCH(1,1) process (to generate verification data)
def simulate_garch(omega, alpha, beta, n, rng):
    eps, sigma2 = np.zeros(n), np.zeros(n)
    sigma2[0] = omega / (1 - alpha - beta)  # initialize at the unconditional variance
    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 estimates the initial variance \(\sigma_0^2\) as an exponentially weighted average (decay 0.94) of the first 75 samples, matching the convention used by the arch library (initializing with the plain sample variance slightly distorts the likelihood computation early on).

Numerical Experiment: Comparison Against the arch Library

We simulated 2000 samples of GARCH(1,1) returns with true parameters \(\omega=0.05\) , \(\alpha=0.10\) , \(\beta=0.85\) , then fit both the from-scratch MLE above and arch.arch_model (vol="Garch", p=1, q=1) to the same data.

Parameters and log-likelihood agree closely (the residual differences are at the level of optimizer convergence tolerance). Both estimates deviate somewhat from the true values (estimated \(\alpha \approx 0.079\) vs. true 0.10), which is ordinary finite-sample estimation error at \(N=2000\) . The high persistence \(\alpha+\beta \approx 0.96\) is correctly recovered, close to the true persistence used in the simulation, \(0.10+0.85=0.95\) .

Summary: How ARIMA and GARCH Divide the Work

https://yuhi-sa.github.io/en/posts/20260226_arima/1/ and this article are complementary.

In practice, combining both into an ARIMA-GARCH model (ARIMA for the mean, GARCH for the residual variance) is the standard starting point for financial time-series analysis.

Related Articles

• ARIMA and SARIMA: Theory and Python Implementation - The foundation for modeling a time series’ mean structure; GARCH complements it on the variance side.
• ARMA/ARIMA in State-Space Form: Kalman-Filter Maximum Likelihood Estimation - Explains the same “prediction-error decomposition” framework behind this article’s GARCH likelihood, in the ARIMA context.
• Autocorrelation Function: Theory and Python Implementation - Foundational concepts behind this article’s Ljung-Box test and the ACF.
• PACF and AR Order Identification - The AR order-selection approach; similar information criteria apply when choosing GARCH’s own order \((p,q)\) .
• Time Series Anomaly Detection: From Statistical Methods to Kalman Filters - GARCH’s conditional variance can also drive dynamic thresholds for anomaly detection.

References

• 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.