Module 14 - ARCH/GARCH (v2)

Author

Kevin Rusu

library(fpp3)
library(fredr)
library(rugarch)

Part I: Model Equations and Connections

ARIMA Model

The ARIMA model stands for Autoregressive Integrated Moving Average. It has three components: an AR part that regresses the variable on its own lags, a differencing step that makes the series stationary, and an MA part that models the error as a linear combination of past forecast errors. The general ARIMA\((p,d,q)\) equation is:

\[\phi(B)(1 - B)^d y_t = c + \theta(B)\varepsilon_t\]

where \(\phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p\) is the AR polynomial, \(\theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q\) is the MA polynomial, \(B\) is the backshift operator, \(d\) is the number of differences applied, and \(\varepsilon_t \sim \text{WN}(0, \sigma^2)\) is a white noise error with constant variance.

ADL Model

The Autoregressive Distributed Lag model, or ADL\((p,q)\), extends the AR idea by bringing in current and lagged values of one or more external explanatory variables. The simplest ADL\((1,1)\) with one regressor \(x_t\) is:

\[y_t = \alpha + \phi_1 y_{t-1} + \beta_0 x_t + \beta_1 x_{t-1} + \varepsilon_t\]

The key difference from a pure ARIMA model is that ARIMA only uses the history of \(y_t\) itself, while ADL adds exogenous variables \(x_t\) and their lags. This is the same idea as the dynamic regression (ARIMAX) model from fpp3, except that ADL keeps the errors as white noise rather than allowing them to follow an ARMA process.

ARCH Model

The Autoregressive Conditional Heteroskedasticity model, or ARCH\((q)\), was introduced by Engle (1982). It models the variance of the error term as time-varying. The model has two equations:

Mean equation:

\[y_t = \mu + \varepsilon_t, \qquad \varepsilon_t = \sigma_t z_t, \qquad z_t \sim N(0,1)\]

Variance equation:

\[\sigma_t^2 = \omega + \alpha_1 \varepsilon_{t-1}^2 + \alpha_2 \varepsilon_{t-2}^2 + \cdots + \alpha_q \varepsilon_{t-q}^2\]

where \(\omega > 0\) and \(\alpha_i \geq 0\) to guarantee a positive variance. A large shock \(\varepsilon_{t-1}\) in either direction inflates \(\sigma_t^2\), producing volatility clustering.

GARCH Model

The Generalized ARCH model, or GARCH\((p,q)\), was introduced by Bollerslev (1986). It extends ARCH by adding lagged values of the conditional variance itself:

Mean equation:

\[y_t = \mu + \varepsilon_t, \qquad \varepsilon_t = \sigma_t z_t, \qquad z_t \sim N(0,1)\]

Variance equation:

\[\sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \varepsilon_{t-i}^2 + \sum_{j=1}^{p} \beta_j \sigma_{t-j}^2\]

where \(\omega > 0\), \(\alpha_i \geq 0\), \(\beta_j \geq 0\), and for covariance stationarity we require \(\sum \alpha_i + \sum \beta_j < 1\). The simplest and most common case is GARCH\((1,1)\):

\[\sigma_t^2 = \omega + \alpha_1 \varepsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2\]

Connections Between the Models

Model	Handles	Variance assumption
ARIMA	Conditional mean only	Constant (\(\sigma^2\))
ADL	Conditional mean + exogenous regressors	Constant (\(\sigma^2\))
ARCH	Simple mean; variance depends on past squared shocks	Time-varying
GARCH	Simple mean; variance depends on past shocks AND past variance	Time-varying, parsimonious

ARIMA and ADL model the conditional mean. ARCH and GARCH model the conditional variance. These ideas combine naturally: rugarch lets you set an ARMA structure for the mean equation and a GARCH structure for the variance equation in the same model.

Part II: GARCH Model on S&P 500 Daily Returns

Data

Daily S&P 500 closing levels from FRED (series SP500) from January 2015 to January 2025
converted levels to log returns (log returns are approximately stationary)

sp500_raw <- fredr(
  series_id         = "SP500",
  observation_start = as.Date("2015-01-01"),
  observation_end   = as.Date("2025-01-01"),
  frequency         = "d"
)

sp500 <- sp500_raw |>
  mutate(date = as.Date(date)) |>
  arrange(date) |>
  drop_na(value) |>
  mutate(log_return = difference(log(value))) |>
  drop_na(log_return) |>
  as_tsibble(index = date)

Time Series Plot

sp500 |>
  autoplot(log_return) +
  labs(
    title   = "S&P 500 Daily Log Returns (2015-2025)",
    y       = "Log Return",
    x       = NULL
  )

The plot shows clear volatility clustering: calm stretches alternate with bursts of large positive and negative returns (notably March 2020 and 2022). This is the primary motivation for using a GARCH model.

Model Specification: Equations with Greek Coefficients

I fit a GARCH(1,1) model with an AR(1) mean equation. Both \(p = 1\) and \(q = 1\) are greater than zero, satisfying the assignment requirement.

Mean equation

The mean of the return series follows an AR(1) process:

\[y_t = \mu + \phi_1 y_{t-1} + \varepsilon_t, \qquad \varepsilon_t = \sigma_t z_t, \qquad z_t \overset{\text{iid}}{\sim} N(0,1)\]

\(\mu\) is the unconditional mean return.
\(\phi_1\) is the AR(1) coefficient capturing first-order autocorrelation in returns.
\(\varepsilon_t\) is the shock, decomposed into a time-varying scale \(\sigma_t\) and a standardized noise \(z_t\).

Variance equation

\[\sigma_t^2 = \omega + \alpha_1 \varepsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2\]

\(\omega > 0\) is the baseline (long-run) variance intercept.
\(\alpha_1 \geq 0\) is the ARCH effect: how much yesterday’s squared shock feeds into today’s variance.
\(\beta_1 \geq 0\) is the GARCH effect: how much yesterday’s conditional variance persists into today.

Model Estimation

returns_vec <- sp500$log_return

# Specify GARCH(1,1) with AR(1) mean, normal innovations
spec <- ugarchspec(
  mean.model     = list(armaOrder = c(1, 0), include.mean = TRUE),
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
  distribution.model = "norm"
)

fit <- ugarchfit(spec = spec, data = returns_vec)

# Full model output
fit


*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(1,0,0)
Distribution    : norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000856    0.000171   5.0062 0.000001
ar1    -0.052911    0.024001  -2.2046 0.027485
omega   0.000004    0.000003   1.2739 0.202694
alpha1  0.193011    0.023905   8.0741 0.000000
beta1   0.779551    0.029429  26.4893 0.000000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000856    0.001071  0.79952 0.423989
ar1    -0.052911    0.022728 -2.32799 0.019912
omega   0.000004    0.000034  0.11246 0.910461
alpha1  0.193011    0.160846  1.19998 0.230149
beta1   0.779551    0.273058  2.85490 0.004305

LogLikelihood : 7241.365 

Information Criteria
------------------------------------
                    
Akaike       -6.6511
Bayes        -6.6380
Shibata      -6.6511
Hannan-Quinn -6.6463

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.400  0.2367
Lag[2*(p+q)+(p+q)-1][2]     1.438  0.4725
Lag[4*(p+q)+(p+q)-1][5]     2.468  0.5755
d.o.f=1
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                    0.01728  0.8954
Lag[2*(p+q)+(p+q)-1][5]   0.73407  0.9163
Lag[4*(p+q)+(p+q)-1][9]   1.84909  0.9222
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]   0.05539 0.500 2.000  0.8139
ARCH Lag[5]   1.37764 1.440 1.667  0.6251
ARCH Lag[7]   2.02540 2.315 1.543  0.7119

Nyblom stability test
------------------------------------
Joint Statistic:  7.922
Individual Statistics:              
mu     0.06883
ar1    0.47752
omega  0.19459
alpha1 0.30488
beta1  0.75767

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.28 1.47 1.88
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                    t-value      prob sig
Sign Bias           3.39620 0.0006956 ***
Negative Sign Bias  0.44126 0.6590668    
Positive Sign Bias  0.06894 0.9450434    
Joint Effect       17.58373 0.0005359 ***


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     114.1    1.397e-15
2    30     120.3    4.338e-13
3    40     138.5    4.555e-13
4    50     145.2    1.833e-11


Elapsed time : 0.091398

coef(fit)

           mu           ar1         omega        alpha1         beta1 
 8.561239e-04 -5.291104e-02  3.877763e-06  1.930113e-01  7.795509e-01

Updated Equations with Numerical Coefficients

Mean equation:

\[y_t = 0.000856 - 0.0529 \, y_{t-1} + \varepsilon_t\]

Variance equation:

\[\sigma_t^2 = 0.00000388 + 0.1930 \, \varepsilon_{t-1}^2 + 0.7796 \, \sigma_{t-1}^2\]

Interpretation of Coefficients

alpha_hat <- coef(fit)["alpha1"]
beta_hat  <- coef(fit)["beta1"]
omega_hat <- coef(fit)["omega"]

persistence   <- alpha_hat + beta_hat
long_run_var  <- omega_hat / (1 - persistence)
long_run_vol  <- sqrt(long_run_var * 252)  # annualized (252 trading days)

cat("Persistence (alpha + beta):", round(persistence, 4), "\n")

Persistence (alpha + beta): 0.9726

cat("Long-run daily variance:   ", round(long_run_var, 8), "\n")

Long-run daily variance:    0.00014133

cat("Long-run annualized vol:   ", round(long_run_vol * 100, 2), "%\n")

Long-run annualized vol:    18.87 %

Volatility Clustering

The estimated \(\alpha_1 = 0.1930\) is positive and significant, confirming that a large return shock on day \(t-1\) raises the conditional variance on day \(t\), regardless of its sign. This is volatility clustering: big moves tend to be followed by more big moves, and quiet days tend to be followed by more quiet days. The return plot shows this clearly in the COVID-19 crash of March 2020 and the 2022 sell-off.

Persistence

Persistence is measured by \(\alpha_1 + \beta_1 = 0.1930 + 0.7796 = 0.9726\). A value this close to 1 means that shocks to variance dissipate very slowly. After a volatility spike, roughly 97% of the elevated variance carries into the next day, so it can take weeks or months to return to normal levels.

Long-Run Variance

Because persistence is less than 1, the process is covariance stationary and the variance has a finite unconditional level it gravitates toward:

\[\bar{\sigma}^2 = \frac{\omega}{1 - \alpha_1 - \beta_1} = \frac{0.00000388}{1 - 0.9726} = 0.0001416\]

Annualized, this implies a long-run volatility of about 18.9%, which is consistent with the historical average volatility of the S&P 500.

Mean Reversion

The mean reversion speed is governed by the gap \(1 - (\alpha_1 + \beta_1) = 0.0274\). Because this gap is small, reversion to the long-run variance is gradual rather than sharp. The model does not predict that volatility snaps back quickly after a crisis; it predicts a slow, steady decay back toward the 18.9% long-run level.

Summary Table

GARCH(1,1) Key Metrics
Concept	Value
ARCH effect (alpha1)	0.193
GARCH effect (beta1)	0.7796
Persistence (alpha + beta)	0.9726
Long-run daily variance	0.00014133
Long-run annualized volatility	18.87%
Mean reversion gap (1 - persistence)	0.0274