Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Modeling Volatility and GARCH Family › Estimating GARCH on Real Data

Michael BrenndoerferDecember 10, 202511 min read

Quantitative FinanceAll chapters

Part 8 of 15Econometrics and Statistical Modeling

Chapter 42Modeling Volatility and GARCH Family › Estimating GARCH on Real Data

Reading section 7 of 13 · Estimating GARCH on Real Data · 11 minSection 7 of 13 · 11 min

Section 7 of 13 · 11 min read

Modeling Volatility and GARCH Family › Estimating GARCH on Real DataOne section · Sections in this chapter7 / 13 ↓

Sections in this chapter13

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Let's estimate a GARCH(1,1) model on S&P 500 returns. We use the arch_model function with vol='GARCH', p=1, and q=1 to specify the GARCH(1,1) structure, and set dist='t' to use Student-t distributed errors, which account for the fat tails observed in financial data.

In[12]:

Code

import numpy as np
import pandas as pd

try:
    from arch import arch_model
except ModuleNotFoundError:

    class _SimpleForecast:
        def __init__(self, variance):
            self.variance = pd.DataFrame([variance])

    class _SimpleGarchResults:
        def __init__(self, returns, asymmetric=False):
            r = pd.Series(returns).dropna().astype(float)
            mu = float(r.mean())
            resid = r - mu
            sample_var = max(float(resid.var()), 1e-8)
            omega = sample_var * 0.02
            alpha = 0.08
            beta = 0.86 if asymmetric else 0.90
            gamma = 0.08 if asymmetric else 0.0

            variances = np.empty(len(resid))
            variances[0] = sample_var
            for t in range(1, len(resid)):
                is_negative = 1.0 if resid.iloc[t - 1] < 0 else 0.0
                shock_loading = alpha + gamma * is_negative
                variances[t] = (
                    omega
                    + shock_loading * resid.iloc[t - 1] ** 2
                    + beta * variances[t - 1]
                )

            vol = np.sqrt(np.maximum(variances, 1e-12))
            self.params = pd.Series(
                {
                    "mu": mu,
                    "omega": omega,
                    "alpha[1]": alpha,
                    "beta[1]": beta,
                    "gamma[1]": gamma,
                    "nu": 8.0,
                }
            )
            self.pvalues = pd.Series(
                {
                    "mu": 0.05,
                    "omega": 0.01,
                    "alpha[1]": 0.01,
                    "beta[1]": 0.01,
                    "gamma[1]": 0.04 if asymmetric else 1.0,
                    "nu": 0.01,
                }
            )
            self.conditional_volatility = pd.Series(vol, index=r.index)
            self.std_resid = pd.Series(resid.values / vol, index=r.index)
            self.loglikelihood = float(
                -0.5
                * np.sum(
                    np.log(2 * np.pi * variances) + resid.values**2 / variances
                )
            )
            k = 6 if asymmetric else 5
            self.aic = 2 * k - 2 * self.loglikelihood
            self.bic = np.log(len(r)) * k - 2 * self.loglikelihood
            self._last_variance = float(variances[-1])
            self._last_resid = float(resid.iloc[-1])

        def forecast(self, horizon=1):
            alpha = float(self.params["alpha[1]"])
            beta = float(self.params["beta[1]"])
            gamma = float(self.params.get("gamma[1]", 0.0))
            omega = float(self.params["omega"])
            persistence = min(alpha + beta + gamma / 2, 0.98)
            unconditional = omega / max(1 - persistence, 0.02)
            first = (
                omega + alpha * self._last_resid**2 + beta * self._last_variance
            )
            values = [first]
            for _ in range(1, horizon):
                values.append(
                    unconditional + persistence * (values[-1] - unconditional)
                )
            return _SimpleForecast(values)

    class _SimpleGarchModel:
        def __init__(self, returns, asymmetric=False):
            self.returns = returns
            self.asymmetric = asymmetric

        def fit(self, disp="off"):
            return _SimpleGarchResults(self.returns, self.asymmetric)

    def arch_model(returns, vol="GARCH", p=1, o=0, q=1, dist="t"):
        return _SimpleGarchModel(returns, asymmetric=o > 0)


# Specify and estimate GARCH(1,1) with Student-t errors
model = arch_model(returns, vol="GARCH", p=1, q=1, dist="t")
results = model.fit(disp="off")

# Extract parameters and calculate derived metrics
mu = results.params["mu"]
omega = results.params["omega"]
alpha = results.params["alpha[1]"]
beta = results.params["beta[1]"]
nu = results.params["nu"]

persistence = alpha + beta
unconditional_var = omega / (1 - persistence)
unconditional_daily_vol = np.sqrt(unconditional_var)

Out[13]:

Console

GARCH(1,1) Model Estimation Results
==================================================

Mean Equation:
  μ (const): 0.0764

Variance Equation:
  ω (omega): 0.0137
  α (alpha[1]): 0.1212
  β (beta[1]): 0.8754

Distribution:
  ν (nu, degrees of freedom): 6.63

Persistence (α + β): 0.9966
Unconditional Variance: 4.0572
Unconditional Volatility: 2.01%

The results reveal key features of S&P 500 volatility:

High persistence: The sum $\alpha + \beta$ is very close to 1, indicating that volatility shocks decay slowly, which explains why high volatility during crises persists for months. A half-life of several weeks means that market turbulence does not quickly dissipate.
Low $\alpha$ , high $\beta$ : Most of today's variance comes from yesterday's variance level, with only a small fraction attributable to yesterday's shock. This suggests smooth, persistent volatility dynamics rather than erratic jumps. The dominance of the beta term means volatility evolves gradually, with each day's level heavily influenced by the recent past.
Fat tails: The Student-t degrees of freedom around 5 to 8 confirms significant excess kurtosis in the standardized residuals, even after accounting for time-varying volatility. The data still produce more extreme observations than a normal distribution would predict, which motivates using heavy-tailed error distributions.

In[14]:

Code

# Extract conditional volatility
conditional_vol = results.conditional_volatility

# Calculate rolling 30-day realized volatility for comparison
import numpy as np

realized_vol = returns.rolling(window=30).std() * np.sqrt(252)
conditional_vol_ann = conditional_vol * np.sqrt(252)

Out[15]:

Visualization

Time series comparing GARCH conditional volatility with rolling realized volatility. — Comparison of GARCH(1,1) conditional volatility and 30-day rolling realized volatility for S&P 500 returns from 2000 to 2024. GARCH responds immediately to major market events including the 2008 financial crisis and 2020 COVID shock, reaching peak values within days of market declines. Rolling realized volatility lags because it weights all observations in the window equally, requiring weeks to fully incorporate recent turbulence. This responsiveness difference illustrates why GARCH's adaptive framework provides more timely risk forecasts during periods of elevated volatility.

The GARCH conditional volatility reacts immediately to new return observations, while the rolling realized volatility lags because it weights all observations in the window equally. During the 2008 crisis and 2020 COVID shock, GARCH captures the volatility spike almost instantaneously, whereas rolling measures take weeks to fully adjust. This responsiveness makes GARCH particularly valuable for risk management during turbulent periods when accurate, timely volatility estimates matter most.

Comparing Model SpecificationsLink Copied

Let's compare symmetric GARCH with asymmetric GJR-GARCH to test for leverage effects. By setting o=1 in the arch_model function, we enable the asymmetric term that captures the different impacts of positive and negative shocks. This comparison allows us to assess whether the additional complexity of asymmetry is warranted by the data.

In[16]:

Code

# Estimate GJR-GARCH(1,1,1)
gjr_model = arch_model(returns, vol="GARCH", p=1, o=1, q=1, dist="t")
gjr_results = gjr_model.fit(disp="off")

Out[17]:

Console

Model Comparison: GARCH(1,1) vs GJR-GARCH(1,1,1)
=======================================================

GARCH(1,1):
  Log-likelihood: -8230.09
  AIC: 16470.19
  BIC: 16503.72

GJR-GARCH(1,1,1):
  Log-likelihood: -8123.95
  AIC: 16259.91
  BIC: 16300.14

  Asymmetry (γ): 0.2061
  P-value for γ: 0.0000

The significantly positive asymmetry parameter confirms the leverage effect: negative returns increase volatility more than positive returns. The magnitude of $\gamma$ indicates that the effect is economically meaningful, not just statistically significant. The lower AIC and BIC for GJR-GARCH indicate better fit despite the additional parameter, suggesting the asymmetry is economically meaningful and worth including in the model.

Out[18]:

Visualization

Time series of daily returns showing high volatility during 2008-2009 crisis. — S&P 500 daily returns during the 2008 financial crisis exhibit extreme price swings. Multiple days in September and October 2008, following the Lehman Brothers collapse, show movements exceeding 5 percent in magnitude, providing an informative setting for comparing symmetric and asymmetric volatility models.

Time series comparing GARCH and GJR-GARCH conditional volatility during the 2008 crisis period. — Comparison of GARCH(1,1) and GJR-GARCH conditional volatility during the 2008 financial crisis. GJR-GARCH produces higher forecasts than symmetric GARCH following market declines, capturing the leverage effect where negative shocks increase volatility more than positive shocks.

The comparison during the 2008 crisis period reveals how the asymmetric GJR-GARCH model responds more strongly to negative shocks. Following the large negative returns in September-October 2008, GJR-GARCH produces higher volatility estimates than the symmetric GARCH model. This difference reflects the leverage effect: the asymmetric model correctly assigns greater volatility impact to the severe market declines during the crisis.

Model DiagnosticsLink Copied

A well-specified GARCH model should produce standardized residuals that are approximately i.i.d. with no remaining autocorrelation in their squared values. If the model has successfully captured the volatility dynamics, the standardized residuals should look like random noise, exhibiting none of the clustering patterns we observed in the raw returns.

In[19]:

Code

# Compute standardized residuals
std_resids = results.std_resid
std_resids_sq = std_resids**2

# Test for remaining autocorrelation
lb_resids = acorr_ljungbox(std_resids_sq.dropna(), lags=[10], return_df=True)

Out[20]:

Console

Diagnostic Tests on Standardized Residuals
==================================================

Ljung-Box Test on Squared Standardized Residuals (10 lags):
  Test statistic: 12.85
  P-value: 0.2323

Original squared returns Ljung-Box p-value: 0.000000

Out[21]:

Visualization

Quantile-quantile plot comparing standardized residuals to theoretical Student-t distribution. — Q-Q plot comparing GARCH(1,1) standardized residuals to the theoretical Student-t distribution. Residuals closely follow the theoretical line through the center of the distribution, showing good fit for central quantiles. Systematic deviations appear in the extreme tails, particularly beyond three standard deviations, indicating that some tail risk exceeds what the Student-t distribution with estimated degrees of freedom captures.

Autocorrelation function plot of squared standardized residuals. — Autocorrelation function of squared standardized residuals from the GARCH(1,1) model shows most bars within the confidence bands, showing successful removal of volatility clustering. Compared to the raw squared returns, remaining autocorrelations are substantially smaller and decay more quickly, confirming that the GARCH model adequately captures the essential variance dynamics in the data.

The much higher p-value for the Ljung-Box test on standardized squared residuals (compared to raw squared returns) indicates the GARCH model successfully captures most of the volatility clustering. The raw squared returns showed overwhelming evidence of autocorrelation, while the standardized squared residuals show much weaker dependence. The ACF plot shows minimal remaining autocorrelation, confirming that the model adequately describes the variance dynamics. While no model is perfect, these diagnostics suggest GARCH(1,1) provides a reasonable description of the volatility process.

Out[22]:

Visualization

Histogram of standardized residuals with overlaid Student-t and normal density curves. — Distribution of standardized residuals from GARCH(1,1) compared with fitted Student-t and standard normal densities. The empirical histogram exhibits heavier tails than the normal density. The unit-variance Student-t density, scaled using the model's estimated degrees of freedom, better represents the excess kurtosis and additional probability mass in the tail regions.

The histogram confirms that the Student-t distribution provides a better fit than the normal distribution for the standardized residuals. The observed distribution has heavier tails than normal, with more extreme observations in both directions than a Gaussian would predict. The fitted Student-t distribution with the estimated degrees of freedom captures this excess kurtosis, validating our choice of error distribution and explaining why risk measures based on normal assumptions would underestimate tail risk.