library(lubridate)
library(zoo)
library(reshape2)
library(car)
library(tseries)
source("common/functions.r", chdir = TRUE)
library(tsgarch)
library(xts)
library(ggplot2)18 Univariate volatility
Volatility is not observed. It has to be inferred from the returns it generated, and every risk number in the rest of the book rests on that inference. This chapter specifies and estimates the main models and inspects the conditional volatility they imply inside the estimation sample. The out-of-sample forecasts used for risk calculations come later, in Chapter 26.
The mathematics of univariate volatility models are discussed in Chapter 2 of Financial Risk Forecasting.
We use maximum likelihood methods for estimating volatility models and implement them with tsgarch in R, arch in Python and ARCHModels in Julia. The tsgarch package is part of the tsmodels ecosystem developed by Alexios Galanos, the same author as the earlier rugarch package. We use tsgarch because it estimates GARCH models cleanly and keeps the code readable.
NotePerformance considerations
The tsgarch package covers many GARCH variants but is slower than simpler alternatives. For applications requiring many repeated estimations (such as rolling window backtests), see the speed comparison in Chapter 26, which benchmarks various GARCH implementations across R, Python and Julia.
This chapter uses a long data sample, and all the estimation works without problems. In Chapter 20, we will see what can go wrong with volatility models and how to fix them.
18.1 Data and libraries
The GARCH package in each language, plus the diagnostic and plotting libraries the chapter uses later.
import numpy as np
import pandas as pd
import sys
sys.path.insert(0, 'common')
from functions import ProcessRawData
from arch import arch_model
from plotnine import ggplot, aes, geom_line, labs, theme_minimalusing CSV, DataFrames, Dates, Statistics
using ARCHModels
using TidierPlots, CairoMakie
using Distributions
using StatsBase
using HypothesisTests
using Random
include("common/functions.jl");data = ProcessRawData()
dates = as.Date(as.character(data$sp500$date), format="%Y%m%d")
y = xts(data$sp500$y, order.by = dates)
y = na.omit(y)
y = y - mean(y)
y = y * 100data = ProcessRawData()
y = data['sp500']['y'].values
y = y[~np.isnan(y)]
y = y - np.mean(y)
y = y * 100data = ProcessRawData();
y = collect(skipmissing(data["sp500"].y));
y = y .- mean(y);
y = y .* 100;18.1.1 What to do with the mean?
We can either de-mean the returns — subtract the sample mean — or estimate the mean along with the other parameters.
However, since the mean is quite small, we can usually ignore it if the sample size is large. Chapter 20 shows how that can fail.
To estimate the mean alongside the volatility parameters, include a constant in the mean equation. Here is a brief example:
spec_mean = garch_modelspec(
y = y,
model = "garch",
order = c(1, 1),
constant = TRUE,
distribution = "norm"
)
fit_mean = estimate(spec_mean)
coef(fit_mean)model_mean = arch_model(y, vol='GARCH', p=1, q=1, dist='normal', mean='Constant')
fit_mean = model_mean.fit(disp='off')
print(fit_mean.params)fit_mean = fit(GARCH{1,1}, y; meanspec=Intercept)
coef(fit_mean)The estimated mean (\(\Mean\)) is typically small and statistically insignificant for daily returns. Where it does matter, it enters as a constant in the mean equation and is estimated jointly with the variance parameters.
18.1.2 Scaling returns
Daily returns are typically expressed in decimal form, where a 1% return is represented as 0.01. These small values cause numerical problems for the optimisation algorithms used in GARCH estimation. When the scale of the data is very small the variance and its derivatives become tiny, the gradients are computed near machine precision and become unreliable, and the optimiser either stalls or stops at the values it started from.
Parameters such as \(\ARCHcoeff\) and \(\GARCHcoeff\) are then reported at the values the optimiser started from, 0.1 and 0.88 for example, rather than at the values that maximise the likelihood. Those numbers are starting values and not bounds — tsgarch and ARCHModels both restrict \(\ARCHcoeff\) and \(\GARCHcoeff\) to the unit interval.
The standard solution is to express returns in percentage form by multiplying by 100, so that a 1% return becomes 1.0 rather than 0.01.
The rescaling affects the interpretation of \(\GARCHconst\) (the constant in the variance equation), which scales with the square of the data. If returns are multiplied by 100, \(\GARCHconst\) is multiplied by 10,000. The \(\ARCHcoeff\) and \(\GARCHcoeff\) parameters, being unitless persistence measures, remain unchanged.
We multiply returns by 100 throughout this chapter to ensure reliable estimation across R, Python and Julia.
18.2 EWMA
The exponentially weighted moving average (EWMA) model sets the conditional variance to a weighted average of past squared returns, where the weight on an observation declines exponentially as it recedes into the past. Unlike more complex GARCH models, EWMA requires only one parameter (\(\EWMAdecay\)) and needs no optimisation, which makes it a useful benchmark for comparison.
The EWMA volatility forecast is given by:
\[\Vol_t^2 = \EWMAdecay \Vol_{t-1}^2 + (1-\EWMAdecay) \CompoundReturns_{t-1}^2\]
where \(\Vol_t^2\) is the conditional variance at time \(t\), \(\CompoundReturns_{t-1}\) is the return at time \(t-1\) and \(\EWMAdecay\) is the decay parameter (typically around 0.94 for daily data).
Nothing is estimated here. We fix \(\EWMAdecay\) at 0.94, start the recursion from the sample variance and run it forward, then plot the returns with ±2 standard deviation bands.
EWMA = vector(length=length(y))
lambda = 0.94
EWMA[1] = var(y)
for (i in 2:length(y)){
EWMA[i] =
lambda * EWMA[i-1]+
(1-lambda) * y[i-1] ^2
}
EWMA=sqrt(EWMA)
par(mar=c(2,3.5,1,0))
matplot(
cbind(y,2*EWMA,-2*EWMA),
type='l',
lty=1,
col=c("black","red","red"),
bty='l',
ylab="",
main="S&P 500 returns with ±2 EWMA sd",
las=1
)lam = 0.94
EWMA = np.zeros(len(y))
EWMA[0] = np.var(y, ddof=1)
for i in range(1, len(y)):
EWMA[i] = lam * EWMA[i-1] + (1-lam) * y[i-1]**2
EWMA = np.sqrt(EWMA)
df_ewma = pd.DataFrame({
't': range(len(y)),
'returns': y,
'upper': 2*EWMA,
'lower': -2*EWMA
})
p = (ggplot(df_ewma, aes(x='t'))
+ geom_line(aes(y='returns'), color='black', size=0.5)
+ geom_line(aes(y='upper'), color='red', size=0.8)
+ geom_line(aes(y='lower'), color='red', size=0.8)
+ labs(title='S&P 500 returns with ±2 EWMA sd', x='', y='')
+ theme_minimal())
p.save("_figs/vol_ewma_py.png", width=10, height=4, dpi=100)
lambda = 0.94;
EWMA_jl = zeros(length(y));
EWMA_jl[1] = var(y);
for i in 2:length(y)
EWMA_jl[i] = lambda * EWMA_jl[i-1] + (1-lambda) * y[i-1]^2
end
EWMA_jl = sqrt.(EWMA_jl);
df_ewma = DataFrame(t = 1:length(y), y = y, upper = 2*EWMA_jl, lower = -2*EWMA_jl);
p = ggplot(df_ewma) +
geom_line(@aes(x = t, y = y), color = "black", linewidth = 0.3) +
geom_line(@aes(x = t, y = upper), color = "red", linewidth = 0.5) +
geom_line(@aes(x = t, y = lower), color = "red", linewidth = 0.5) +
labs(title = "S&P 500 returns with ±2 EWMA sd", x = "", y = "") +
theme_minimal();
ggsave("_figs/vol_ewma_jl.png", p);
EWMA is a useful baseline, but volatility never reverts under it, so a shock persists indefinitely, and one decay parameter is all it has. GARCH fixes both. Volatility reverts towards a long-run level, and the lag structure and innovation distribution are ours to choose.
18.3 GARCH
GARCH models let volatility depend on past volatility (the GARCH component) and past squared returns (the ARCH component), whilst accommodating different distributional assumptions for the innovations.
The GARCH(1,1) conditional variance is
\[\Vol_t^2 = \GARCHconst + \ARCHcoeff \CompoundReturns_{t-1}^2 + \GARCHcoeff \Vol_{t-1}^2\]
where \(\ARCHcoeff\) is the weight on the most recent squared return, which is how news enters, and \(\GARCHcoeff\) is the weight on the previous conditional variance, which is how the model remembers. The parameters must satisfy \(\GARCHconst > 0\), \(\ARCHcoeff \geq 0\) and \(\GARCHcoeff \geq 0\) for the variance to stay positive, and \(\ARCHcoeff + \GARCHcoeff < 1\) for it to be covariance stationary. The unconditional variance is then
\[\Var(\CompoundReturns_t) = \frac{\GARCHconst}{1 - \ARCHcoeff - \GARCHcoeff}\]
and volatility reverts towards it, which is what EWMA lacks. The sum \(\ARCHcoeff + \GARCHcoeff\) measures persistence and is usually close to one for daily equity returns.
We implement GARCH models using the tsgarch package in R, the arch package in Python and the ARCHModels package in Julia. All three support specification, estimation and forecasting for the common GARCH variants, including the asymmetric models that capture the leverage effect in equity returns.
18.3.1 Setup
In all three packages, we first specify the model and then estimate it.
To estimate a univariate GARCH model with tsgarch:
- Create a model specification using
garch_modelspec()- this includes the model type, order, distribution and data - Fit the specified model using
estimate()
The main arguments to garch_modelspec() are:
model: The GARCH model type (“garch”, “egarch”, “gjrgarch”, “aparch”, “fgarch”, “cgarch”, “igarch”, “ewma”)order: A vector c(p, q) specifying GARCH(p,q) orderdistribution: The innovation distribution (“norm”, “std”, “snorm”, “sstd”, “ged”, “sged”, “nig”, “gh”, “jsu”, “ghst”)constant: Whether to include a mean (default FALSE)
To estimate a univariate GARCH model with arch:
- Create a model using
arch_model()- this specifies the model type, order, distribution and data - Fit the model using
.fit()
The main arguments to arch_model() are:
vol: The volatility model (“GARCH”, “ARCH”, “EGARCH”, “FIGARCH”, “APARCH”, “HARCH”)p,q: The GARCH order parameterso: The asymmetric order (for GJR-GARCH, use o=1)dist: The distribution (“normal”, “t”, “skewt”, “ged”)mean: The mean model (“Zero”, “Constant”, “AR”, etc.)
To estimate a univariate GARCH model with ARCHModels:
- Fit the model directly using
fit(GARCH{p,q}, data)with optional distribution specification
The main arguments are:
- Model type:
GARCH{p,q},EGARCH{o,p,q},TGARCH{o,p,q}(GJR-GARCH) dist: The distribution (StdNormal,StdT,StdSkewT,StdGED)meanspec: The mean specification (NoIntercept,Intercept)
18.3.2 Sample size considerations
Models with more parameters require more data. Insufficient sample sizes cause parameter instability, convergence failures and unreliable standard errors.
As a rule of thumb, a GARCH(1,1) needs at least 500 observations, while a GARCH-t needs around 3,000. Asymmetric models with skewed distributions require even larger samples due to their additional parameters. We see an example of estimation failure from insufficient data in Chapter 20.
18.3.3 The models
Now we estimate several GARCH models and see what each assumption changes.
We will put the results into a list (R), dictionary (Python) or Dict (Julia) so we can automatically process them afterwards.
Results = list()Results = {}Results_jl = Dict{String, Any}();18.3.3.1 Gaussian GARCH(1,1)
GARCH(1,1) builds persistence into volatility. High volatility tends to be followed by high volatility.
spec.1 = garch_modelspec(
y = y,
model = "garch",
order = c(1, 1),
distribution = "norm"
)
Results$GARCH11 = estimate(spec.1)model = arch_model(y, vol='GARCH', p=1, q=1, dist='normal', mean='Zero')
Results['GARCH11'] = model.fit(disp='off')Results_jl["GARCH11"] = fit(GARCH{1,1}, y; meanspec=NoIntercept);18.3.3.1.1 Printing and plotting
The fitted model stores the main results, and we can pull them out directly:
# View coefficients
p = coef(Results$GARCH11)
cat(sprintf("omega: %.6f\nalpha: %.6f\nbeta: %.6f\n", p["omega"], p["alpha1"], p["beta1"]))omega: 0.018975
alpha: 0.105270
beta: 0.879347# View log-likelihood
cat("Log-likelihood:", round(as.numeric(logLik(Results$GARCH11)), 2), "\n")Log-likelihood: -11923.08# View fitted conditional volatility
cat(round(head(as.numeric(Results$GARCH11$sigma)), 4), "\n")1.1432 1.0849 1.0671 1.0621 1.0146 1.0394# View coefficients
p = Results['GARCH11'].params
print(f"omega: {p['omega']:.6f}")
print(f"alpha: {p['alpha[1]']:.6f}")
print(f"beta: {p['beta[1]']:.6f}")omega: 0.018981
alpha: 0.105314
beta: 0.879324# View log-likelihood
print(f"Log-likelihood: {Results['GARCH11'].loglikelihood:.2f}")Log-likelihood: -11923.25# View fitted conditional volatility
print(np.round(Results['GARCH11'].conditional_volatility[:6], 4))[1.0038 0.9559 0.9525 0.9616 0.9222 0.9607]# View coefficients
using Printf
c = coef(Results_jl["GARCH11"])
# ARCHModels coefficient order is (omega, beta, alpha)
@printf("omega: %.6f\nalpha: %.6f\nbeta: %.6f\n", c[1], c[3], c[2])omega: 0.018976
alpha: 0.105267
beta: 0.879347# View log-likelihood
using Printf
@printf("Log-likelihood: %.2f\n", loglikelihood(Results_jl["GARCH11"]))Log-likelihood: -11923.08# View fitted conditional volatility
println(round.(volatilities(Results_jl["GARCH11"])[1:6], digits=4))[1.1437, 1.0854, 1.0675, 1.0625, 1.015, 1.0397]We can also plot the estimated conditional volatility:
df_sigma = data.frame(t = 1:length(sigma(Results$GARCH11)), sigma = as.numeric(sigma(Results$GARCH11)))
ggplot(df_sigma, aes(x = t, y = sigma)) +
geom_line() +
labs(title = "GARCH(1,1) Conditional Volatility", x = "", y = "") +
theme_minimal()
df_vol = pd.DataFrame({
't': range(len(Results['GARCH11'].conditional_volatility)),
'sigma': Results['GARCH11'].conditional_volatility
})
p = (ggplot(df_vol, aes(x='t', y='sigma'))
+ geom_line()
+ labs(title='GARCH(1,1) Conditional Volatility', x='', y='')
+ theme_minimal())
p.save("_figs/vol_garch_sigma_py.png", width=10, height=4, dpi=100)
df_vol = DataFrame(t = 1:length(volatilities(Results_jl["GARCH11"])), sigma = volatilities(Results_jl["GARCH11"]));
p = ggplot(df_vol, @aes(x = t, y = sigma)) +
geom_line(linewidth = 0.5) +
labs(title = "GARCH(1,1) Conditional Volatility", x = "", y = "") +
theme_minimal();
ggsave("_figs/vol_garch_sigma_jl.png", p);
18.3.3.1.2 Comparing estimates
The three packages fit the same specification to the same data, so the estimates should agree to several decimals. The table below puts them side by side.
|Parameter | R| Python| Julia|
|:---------|--------:|--------:|--------:|
|omega | 0.018975| 0.020245| 0.020272|
|alpha | 0.105270| 0.006221| 0.006314|
|beta | 0.879347| 0.895360| 0.895225|18.3.3.2 GARCH-t(1,1)
We now allow for heavier tails by using Student-t distributed innovations instead of normal. This better captures the extreme returns commonly observed in financial markets whilst maintaining the same volatility dynamics.
We can see that the coefficients now include shape (R), nu (Python) or the fourth coefficient (Julia), which is the estimated degrees of freedom:
spec.2 = garch_modelspec(
y = y,
model = "garch",
order = c(1, 1),
distribution = "std"
)
Results$tGARCH11 = estimate(spec.2)
p = coef(Results$tGARCH11)
cat(sprintf("omega: %.6f\nalpha: %.6f\nbeta: %.6f\nshape: %.6f\n", p["omega"], p["alpha1"], p["beta1"], p["shape"]))omega: 0.011564
alpha: 0.099805
beta: 0.894969
shape: 6.300055model = arch_model(y, vol='GARCH', p=1, q=1, dist='t', mean='Zero')
Results['tGARCH11'] = model.fit(disp='off')
p = Results['tGARCH11'].params
print(f"omega: {p['omega']:.6f}")
print(f"alpha: {p['alpha[1]']:.6f}")
print(f"beta: {p['beta[1]']:.6f}")
print(f"nu: {p['nu']:.6f}")omega: 0.011563
alpha: 0.099827
beta: 0.894965
nu: 6.301738using Printf
Results_jl["tGARCH11"] = fit(GARCH{1,1}, y; dist=StdT, meanspec=NoIntercept)
c = coef(Results_jl["tGARCH11"])
# ARCHModels coefficient order is (omega, beta, alpha)
@printf("omega: %.6f\nalpha: %.6f\nbeta: %.6f\nnu: %.6f\n", c[1], c[3], c[2], c[4])omega: 0.011564
alpha: 0.099801
beta: 0.894971
nu: 6.300256The degrees of freedom is the fourth parameter. A value in the single digits is the usual finding for daily equity returns and is what makes the Student-t worth fitting at all.
|Parameter | R| Python| Julia|
|:---------|--------:|--------:|--------:|
|omega | 0.011564| 0.020245| 0.020272|
|alpha | 0.099805| 0.006221| 0.006314|
|beta | 0.894969| 0.895360| 0.895225|
|shape/nu | 6.300055| 6.301738| 6.636328|18.3.3.3 Skewed Student-t GARCH(1,1)
The skewed Student-t distribution accommodates both heavy tails and asymmetry in the return distribution. This allows the model to capture the negative skewness often found in equity returns. The asymmetry here is in the innovation distribution and not in the variance dynamics, which stay symmetric — the leverage effect is the subject of the GJR-GARCH model below.
The three packages do not fit the same skewed Student-t distribution. R’s sstd is the skewed Student of Fernández and Steel, while Python’s skewt and Julia’s StdSkewT both implement the standardised skewed Student of Hansen (1994). The variance equation is the same in all three, but the distributions are different two-parameter families, so the tail and skewness parameters are not the same quantity and the log-likelihoods are not directly comparable across languages.
spec.3 = garch_modelspec(
y = y,
model = "garch",
order = c(1, 1),
distribution = "sstd"
)
Results$atGARCH11 = estimate(spec.3)
p = coef(Results$atGARCH11)
cat(sprintf("omega: %.6f\nalpha: %.6f\nbeta: %.6f\nshape: %.6f\nskew: %.6f\n", p["omega"], p["alpha1"], p["beta1"], p["shape"], p["skew"]))omega: 0.011870
alpha: 0.100437
beta: 0.893554
shape: 6.636310
skew: 0.901808model = arch_model(y, vol='GARCH', p=1, q=1, dist='skewt', mean='Zero')
Results['atGARCH11'] = model.fit(disp='off')
p = Results['atGARCH11'].params
print(f"omega: {p['omega']:.6f}")
print(f"alpha: {p['alpha[1]']:.6f}")
print(f"beta: {p['beta[1]']:.6f}")
print(f"eta: {p['eta']:.6f}")
print(f"lambda: {p['lambda']:.6f}")omega: 0.011866
alpha: 0.100441
beta: 0.893568
eta: 6.638628
lambda: -0.103133using Printf
Results_jl["atGARCH11"] = fit(GARCH{1,1}, y; dist=StdSkewT, meanspec=NoIntercept)
c = coef(Results_jl["atGARCH11"])
# ARCHModels coefficient order is (omega, beta, alpha)
@printf("omega: %.6f\nalpha: %.6f\nbeta: %.6f\nnu: %.6f\nlambda: %.6f\n", c[1], c[3], c[2], c[4], c[5])omega: 0.011870
alpha: 0.100437
beta: 0.893555
nu: 6.636328
lambda: -0.102989The variance parameters are comparable across the three languages and the table below shows them.
|Parameter | R| Python| Julia|
|:---------|--------:|--------:|--------:|
|omega | 0.011870| 0.020245| 0.020272|
|alpha | 0.100437| 0.006221| 0.006314|
|beta | 0.893554| 0.895360| 0.895225|The distribution parameters are not. Each language is reported against its own family, with the value the parameter takes under symmetry.
|Language |Skewed Student family | Tail parameter| Skew parameter| Value under symmetry|
|:--------|:---------------------|--------------:|--------------:|--------------------:|
|R |Fernández--Steel | 6.6363| 0.9018| 1|
|Python |Hansen | 6.6386| -0.1031| 0|
|Julia |Hansen | 6.6363| -0.1030| 0|R’s skew equals 1 under symmetry, while Python’s lambda and Julia’s λ equal 0. All three estimates sit on the left-skewed side of their own symmetry value, so the three fits agree on the direction of the asymmetry even though they do not estimate the same parameter.
18.3.3.4 Gaussian GJR-GARCH(1,1)
The GJR-GARCH specification introduces asymmetric volatility responses, where negative returns can have different impacts on future volatility than positive returns of the same magnitude. This captures the leverage effect common in equity markets.
The conditional variance becomes
\[\Vol_t^2 = \GARCHconst + \left(\ARCHcoeff + \GJRleverage I_{t-1}\right) \CompoundReturns_{t-1}^2 + \GARCHcoeff \Vol_{t-1}^2\]
where \(I_{t-1}\) equals one when \(\CompoundReturns_{t-1} < 0\) and zero otherwise. A positive return moves volatility by \(\ARCHcoeff\) and a negative return of the same size moves it by \(\ARCHcoeff + \GJRleverage\), so \(\GJRleverage\) is the extra response to bad news and the model reduces to GARCH(1,1) when \(\GJRleverage = 0\). Positivity now requires \(\ARCHcoeff + \GJRleverage \geq 0\) as well, and the stationarity condition becomes \(\ARCHcoeff + \GJRleverage / 2 + \GARCHcoeff < 1\) for symmetric innovations.
spec.4 = garch_modelspec(
y = y,
model = "gjrgarch",
order = c(1, 1),
distribution = "norm"
)
Results$GJRGARCH11 = estimate(spec.4)
p = coef(Results$GJRGARCH11)
cat(sprintf("omega: %.6f\nalpha: %.6f\ngamma: %.6f\nbeta: %.6f\n", p["omega"], p["alpha1"], p["gamma1"], p["beta1"]))omega: 0.020272
alpha: 0.006313
gamma: 0.154964
beta: 0.895228model = arch_model(y, vol='GARCH', p=1, o=1, q=1, dist='normal', mean='Zero')
Results['GJRGARCH11'] = model.fit(disp='off')
p = Results['GJRGARCH11'].params
print(f"omega: {p['omega']:.6f}")
print(f"alpha: {p['alpha[1]']:.6f}")
print(f"gamma: {p['gamma[1]']:.6f}")
print(f"beta: {p['beta[1]']:.6f}")omega: 0.020245
alpha: 0.006221
gamma: 0.154955
beta: 0.895360using Printf
Results_jl["GJRGARCH11"] = fit(TGARCH{1,1,1}, y; meanspec=NoIntercept)
c = coef(Results_jl["GJRGARCH11"])
@printf("omega: %.6f\ngamma: %.6f\nbeta: %.6f\nalpha: %.6f\n", c[1], c[2], c[3], c[4])omega: 0.020272
gamma: 0.154971
beta: 0.895225
alpha: 0.006314The leverage parameter is the extra term. A positive estimate means bad news raises tomorrow’s volatility by more than good news of the same size, which is the asymmetry the specification exists to capture.
|Parameter | R| Python| Julia|
|:---------|--------:|--------:|--------:|
|omega | 0.020272| 0.020245| 0.020272|
|alpha | 0.006313| 0.006221| 0.006314|
|gamma | 0.154964| 0.154955| 0.154971|
|beta | 0.895228| 0.895360| 0.895225|18.4 Comparing results
Now we compare the models. A single table of log-likelihoods shows how much each specification improves the fit.
|Model | R| Python| Julia|
|:-------------------|---------:|---------:|---------:|
|GARCH(1,1) | -11923.08| -11923.25| -11923.08|
|GARCH-t(1,1) | -11713.75| -11713.89| -11713.75|
|Skewed-t GARCH(1,1) | -11683.91| -11683.96| -11683.91|
|GJR-GARCH(1,1) | -11771.19| -11771.10| -11771.19|Read the first, second and fourth rows across the languages, since those three specifications are the same in all of them. The skewed-t row is not a cross-language comparison, because R maximises a Fernández–Steel likelihood while Python and Julia maximise Hansen’s. Within any one column all four log-likelihoods are computed on the same data and are comparable.
18.5 Diagnostics
A fitted model can look fine and still fail its own assumptions, so we test the residuals and compare nested specifications.
18.5.1 Likelihood ratio test
The likelihood ratio (LR) test compares two nested models by twice the difference in their log-likelihoods. Here we compare the Gaussian GARCH(1,1) with the GARCH-t(1,1), which restricts nothing in the variance equation and only changes the innovation distribution.
# Perform the LR test
LR_statistic = 2 * (as.numeric(logLik(Results$tGARCH11)) - as.numeric(logLik(Results$GARCH11)))
p_value = pchisq(LR_statistic, df = 1, lower.tail = FALSE)
cat(sprintf("Likelihood of Gaussian GARCH: %.1f\n", as.numeric(logLik(Results$GARCH11))))
cat(sprintf("Likelihood of GARCH-t: %.1f\n", as.numeric(logLik(Results$tGARCH11))))
cat(sprintf("2 * (Lu - Lr): %.2f\n", LR_statistic))
cat(sprintf("p-value: %.2e\n", p_value))Likelihood of Gaussian GARCH: -11923.1
Likelihood of GARCH-t: -11713.8
2 * (Lu - Lr): 418.65
p-value: 4.81e-93from scipy import stats
LR_statistic = 2 * (Results['tGARCH11'].loglikelihood - Results['GARCH11'].loglikelihood)
p_value = stats.chi2.sf(LR_statistic, df=1)
print(f"Likelihood of Gaussian GARCH: {Results['GARCH11'].loglikelihood:.1f}")
print(f"Likelihood of GARCH-t: {Results['tGARCH11'].loglikelihood:.1f}")
print(f"2 * (Lu - Lr): {LR_statistic:.2f}")
print(f"p-value: {p_value:.2e}")Likelihood of Gaussian GARCH: -11923.3
Likelihood of GARCH-t: -11713.9
2 * (Lu - Lr): 418.72
p-value: 4.64e-93using Printf, Distributions
LR_statistic = 2 * (loglikelihood(Results_jl["tGARCH11"]) - loglikelihood(Results_jl["GARCH11"]));
p_value = ccdf(Chisq(1), LR_statistic);
@printf("Likelihood of Gaussian GARCH: %.1f\n", loglikelihood(Results_jl["GARCH11"]))
@printf("Likelihood of GARCH-t: %.1f\n", loglikelihood(Results_jl["tGARCH11"]))
@printf("2 * (Lu - Lr): %.2f\n", LR_statistic)
@printf("p-value: %.2e\n", p_value)Likelihood of Gaussian GARCH: -11923.1
Likelihood of GARCH-t: -11713.8
2 * (Lu - Lr): 418.64
p-value: 4.81e-93The Gaussian model is a limiting case in the closure of the Student-t family as \(\DOF \to \infty\), a boundary point rather than an interior null. In the idealised boundary formulation the likelihood-ratio statistic’s asymptotic null distribution is the mixture \(\tfrac{1}{2}\chi^2_{(0)} + \tfrac{1}{2}\chi^2_{(1)}\), not \(\chi^2_{(1)}\).
For Julia, whose \(\DOF\) is unbounded, the \(\chi^2_{(1)}\) reference is conservative relative to this asymptotic boundary mixture. For the capped R and Python fits, the Gaussian limit is not even in the closure of the fitted family, so neither reference law formally applies and the \(\chi^2_{(1)}\) reference is only a heuristic approximation.
The p-value is extremely small but not zero, so we reject the null hypothesis that the Gaussian restriction is adequate. The GARCH-t model has a significantly larger likelihood than the Gaussian GARCH model, indicating that heavy tails improve the fit.
18.5.2 Residual analysis
Examining standardised residuals reveals whether our model successfully captures the volatility dynamics. Unlike price forecasting where we can directly compare predictions with observed outcomes, volatility is latent and unobservable. We can only validate our volatility models indirectly through residual behaviour. Well-specified models should produce residuals that are approximately independent and identically distributed, with no remaining volatility clustering.
18.5.2.1 GARCH
Although the likelihood ratio test above shows that the GARCH-t provides a significantly better fit, we use the Gaussian GARCH(1,1) here as a baseline example for residual analysis. The Gaussian model makes the strongest distributional assumption, so its residual diagnostics most clearly reveal departures from the assumed distribution — heavy tails will show up in the QQ plot, for instance. The same diagnostic procedures apply to any fitted specification.
We analyse the standardised residuals from our GARCH(1,1) model to check for remaining serial correlation and distributional assumptions. Patterns in these residuals show where the model fails and what we need to change.
The three tabs use the same conventions so that the figures can be read against each other. The ACF plots show 20 lags with 95% significance bands, and the QQ plots use the plotting positions \((i - 0.5) / \SampleSize\) with a 45 degree reference line.
std_resid = residuals(Results$GARCH11, standardize = TRUE)
df_resid = data.frame(t = 1:length(std_resid), resid = as.numeric(std_resid))
ggplot(df_resid, aes(x = t, y = resid)) +
geom_line() +
labs(x = "", y = "") +
theme_minimal()
myACF(as.numeric(std_resid), n = 21, main = "ACF of standardised residuals")
myACF(as.numeric(std_resid)^2, n = 21, main = "ACF of squared standardised residuals")
x = qqPlot(as.numeric(std_resid), distribution = "norm", envelope = FALSE, line = "none",
main = "Normal QQ of standardised residuals")
abline(0, 1, col = "red")
from plotnine import geom_bar, geom_hline, geom_point, geom_abline
from statsmodels.tsa.stattools import acf
from scipy.stats import norm
std_resid = Results['GARCH11'].std_resid
# Standardised residuals line plot
df_resid = pd.DataFrame({'t': range(len(std_resid)), 'resid': std_resid})
p = (ggplot(df_resid, aes(x='t', y='resid'))
+ geom_line()
+ labs(title='Standardised Residuals', x='', y='')
+ theme_minimal())
p.save("_figs/vol_resid_py.png", width=10, height=4, dpi=100)
# ACF of standardised residuals
acf_vals = acf(std_resid, nlags=20, fft=False)
sig = 1.96 / np.sqrt(len(std_resid))
df_acf = pd.DataFrame({'lag': range(1, 21), 'acf': acf_vals[1:21]})
p = (ggplot(df_acf, aes(x='lag', y='acf'))
+ geom_bar(stat='identity', fill='steelblue')
+ geom_hline(yintercept=sig, color='red')
+ geom_hline(yintercept=-sig, color='red')
+ labs(title='ACF of standardised residuals', x='Lag', y='ACF')
+ theme_minimal())
p.save("_figs/vol_acf_py.png", width=10, height=4, dpi=100)
# ACF of squared standardised residuals
acf_sq = acf(std_resid**2, nlags=20, fft=False)
df_acf_sq = pd.DataFrame({'lag': range(1, 21), 'acf': acf_sq[1:21]})
p = (ggplot(df_acf_sq, aes(x='lag', y='acf'))
+ geom_bar(stat='identity', fill='steelblue')
+ geom_hline(yintercept=sig, color='red')
+ geom_hline(yintercept=-sig, color='red')
+ labs(title='ACF of squared standardised residuals', x='Lag', y='ACF')
+ theme_minimal())
p.save("_figs/vol_acf_sq_py.png", width=10, height=4, dpi=100)
# QQ plot
sorted_resid = np.sort(std_resid)
n = len(sorted_resid)
theoretical = norm.ppf((np.arange(1, n+1) - 0.5) / n)
df_qq = pd.DataFrame({'theoretical': theoretical, 'sample': sorted_resid})
p = (ggplot(df_qq, aes(x='theoretical', y='sample'))
+ geom_point(alpha=0.5)
+ geom_abline(slope=1, intercept=0, color='red')
+ labs(title='Normal QQ of standardised residuals', x='Theoretical Quantiles', y='Sample Quantiles')
+ theme_minimal())
p.save("_figs/vol_qq_py.png", width=6, height=6, dpi=100)
std_resid_jl = residuals(Results_jl["GARCH11"]);
# Create residuals time series plot
df_resid = DataFrame(t = 1:length(std_resid_jl), resid = std_resid_jl);
p = ggplot(df_resid, @aes(x = t, y = resid)) +
geom_line(linewidth = 0.3) +
labs(title = "Standardised Residuals", x = "", y = "") +
theme_minimal();
ggsave("_figs/vol_resid_jl.png", p);
# Compute and plot ACF
sig = 1.96 / sqrt(length(std_resid_jl));
acf_vals = autocor(std_resid_jl, 1:20);
df_acf = DataFrame(lag = 1:20, acf = acf_vals);
p = ggplot(df_acf) +
geom_col(@aes(x = lag, y = acf)) +
geom_hline(yintercept = sig, color = "red") +
geom_hline(yintercept = -sig, color = "red") +
labs(title = "ACF of standardised residuals", x = "Lag", y = "ACF") +
theme_minimal();
ggsave("_figs/vol_acf_jl.png", p);
# ACF of squared residuals
acf_sq_vals = autocor(std_resid_jl.^2, 1:20);
df_acf_sq = DataFrame(lag = 1:20, acf = acf_sq_vals);
p = ggplot(df_acf_sq) +
geom_col(@aes(x = lag, y = acf)) +
geom_hline(yintercept = sig, color = "red") +
geom_hline(yintercept = -sig, color = "red") +
labs(title = "ACF of squared standardised residuals", x = "Lag", y = "ACF") +
theme_minimal();
ggsave("_figs/vol_acf_sq_jl.png", p);
# QQ plot
sorted_resid = sort(std_resid_jl);
theoretical_q = quantile.(Normal(), ((1:length(std_resid_jl)) .- 0.5) ./ length(std_resid_jl));
df_qq = DataFrame(theoretical = theoretical_q, sample = sorted_resid, line = theoretical_q);
p = ggplot(df_qq) +
geom_point(@aes(x = theoretical, y = sample), size = 0.5, alpha = 0.5) +
geom_line(@aes(x = theoretical, y = line), color = "red") +
labs(title = "Normal QQ of standardised residuals", x = "Theoretical", y = "Sample") +
theme_minimal();
ggsave("_figs/vol_qq_jl.png", p);
18.5.2.2 Formal residual tests
The plots suggest where a model fails and the tests put numbers on it. We run the Ljung–Box test for autocorrelation in the standardised residuals, the same test on their squares for remaining volatility clustering and the Jarque–Bera test for normality. All three languages use 20 lags and make no degrees of freedom adjustment for the estimated GARCH parameters.
lb = Box.test(as.numeric(std_resid), lag = 20, type = "Ljung-Box")
lb2 = Box.test(as.numeric(std_resid)^2, lag = 20, type = "Ljung-Box")
jb = jarque.bera.test(as.numeric(std_resid))
cat(sprintf("Ljung-Box residuals: Q = %.2f, df = %d, p = %.4g\n",
lb$statistic, lb$parameter, lb$p.value))
cat(sprintf("Ljung-Box squared residuals: Q = %.2f, df = %d, p = %.4g\n",
lb2$statistic, lb2$parameter, lb2$p.value))
cat(sprintf("Jarque-Bera: JB = %.2f, df = %d, p = %.4g\n",
jb$statistic, jb$parameter, jb$p.value))Ljung-Box residuals: Q = 32.33, df = 20, p = 0.03987
Ljung-Box squared residuals: Q = 23.17, df = 20, p = 0.2807
Jarque-Bera: JB = 1816.71, df = 2, p = 0from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import jarque_bera
lb = acorr_ljungbox(std_resid, lags=[20])
lb2 = acorr_ljungbox(std_resid**2, lags=[20])
jb_stat, jb_p = jarque_bera(std_resid)
print(f"Ljung-Box residuals: Q = {lb['lb_stat'].iloc[0]:.2f}, df = 20, p = {lb['lb_pvalue'].iloc[0]:.4g}")
print(f"Ljung-Box squared residuals: Q = {lb2['lb_stat'].iloc[0]:.2f}, df = 20, p = {lb2['lb_pvalue'].iloc[0]:.4g}")
print(f"Jarque-Bera: JB = {jb_stat:.2f}, df = 2, p = {jb_p:.4g}")Ljung-Box residuals: Q = 32.17, df = 20, p = 0.04148
Ljung-Box squared residuals: Q = 23.07, df = 20, p = 0.2854
Jarque-Bera: JB = 1816.47, df = 2, p = 0using Printf, HypothesisTests
lb = LjungBoxTest(std_resid_jl, 20);
lb2 = LjungBoxTest(std_resid_jl.^2, 20);
jb = JarqueBeraTest(std_resid_jl);
@printf("Ljung-Box residuals: Q = %.2f, df = %d, p = %.4g\n", lb.Q, lb.lag, pvalue(lb))
@printf("Ljung-Box squared residuals: Q = %.2f, df = %d, p = %.4g\n", lb2.Q, lb2.lag, pvalue(lb2))
@printf("Jarque-Bera: JB = %.2f, df = %d, p = %.4g\n", jb.JB, 2, pvalue(jb))Ljung-Box residuals: Q = 32.33, df = 20, p = 0.03986
Ljung-Box squared residuals: Q = 23.17, df = 20, p = 0.2807
Jarque-Bera: JB = 1816.68, df = 2, p = 0The squared standardised residuals show no autocorrelation at 20 lags that the test can detect, which suggests the variance equation has absorbed the volatility clustering — what we are asking of it. Absence of a rejection is not proof that no dependence remains. The residuals themselves are marginally autocorrelated, rejecting at 5% but not at 1%, which points at the mean equation rather than the volatility model, since we fitted a zero mean.
Treat that marginal result as marginal. The calibration here is weaker than the p-value suggests. These are fitted residuals, not observed data, and the test ignores the estimation of \(\GARCHconst\), \(\ARCHcoeff\) and \(\GARCHcoeff\) that produced them, while the residual distribution remains strongly non-normal. The nominal bands and critical values are the ones for independent, well-behaved observations. A result sitting between the 5% and 1% thresholds is exactly the case a modest calibration error can move either way, so it is a reason to look at the mean equation, not a finding about it.
The Jarque–Bera test rejects normality decisively. That is the same conclusion the QQ plot reaches by eye and the reason the GARCH-t fits better.
18.6 Which comparisons are legitimate
Before comparing the four fitted models, three qualifications decide which comparisons mean anything.
First, for the GJR-GARCH pair the \(\chi^2_{(1)}\) reference is standard only where zero leverage is an interior point of the parameter space. R and Python permit negative leverage, but the installed ARCHModels constrains \(\GJRleverage\) to \([0, 0.5]\), making zero a boundary there, so treat the \(\chi^2_{(1)}\) reference as a conservative approximation for the Julia fit specifically.
Second, the GARCH-t within skewed Student-t comparison is theoretical and local. R and Python cap \(\DOF\), at different values for the two variants, while Julia leaves it unbounded. The comparison is valid because the fitted optima lie inside the shared interior of these ranges, not because the bounds themselves match. It is also a within-language comparison only, since R restricts a Fernández–Steel likelihood while Python and Julia restrict Hansen’s.
Third, these three pairs are deliberately not an exhaustive list of nested relationships. Gaussian GARCH also sits in the closure of skewed Student-t GARCH via symmetry combined with \(\DOF \to \infty\), a nonstandard two-restriction boundary problem this exercise avoids, using AIC and BIC there instead.
A model that survives its diagnostics has been shown to fit the sample it was estimated on. That is not the same as describing the data well or forecasting it well. Chapter 19 takes the four fitted specifications and asks both questions.
18.7 Exercises
Model comparison: Not every pair of the four fitted models is nested, so a single pairwise LR sweep does not apply uniformly. Test the three one-parameter nested pairs with a likelihood ratio test — normal GARCH within GARCH-t (the boundary case discussed above), normal GARCH within GJR-GARCH (testing \(\GJRleverage = 0\)) and GARCH-t within the fitted skewed Student-t GARCH (testing the symmetry value). For every other direct comparison among the four models, use AIC and BIC instead, since those pairs are not nested. Comment on the results and discuss which models provide significant improvements.
Residual diagnostics: For the best-performing model, conduct residual analysis including normality tests, autocorrelation tests and ARCH effects tests.
Parameter stability: Investigate parameter stability by estimating your preferred model on rolling windows of 1000 observations. Plot the evolution of the main parameters over time.
Answer to exercise 1. Parameter counts must be hard-coded rather than read from the fitted objects. tsgarch’s stored parameter count runs one higher than its visible coefficients because it counts the mean specification internally even under a zero mean, while arch and ARCHModels count only the visible coefficients — mixing package-reported counts across languages would make AIC and BIC incomparable. The four models, in order, use \(k = (3, 4, 5, 4)\) visible parameters (GARCH, GARCH-t, skewed Student-t GARCH, GJR-GARCH), and \(\SampleSize\) is the sample size, which the code calls n.
k = c(3, 4, 5, 4)
n = length(y)
loglik = c(
as.numeric(logLik(Results$GARCH11)),
as.numeric(logLik(Results$tGARCH11)),
as.numeric(logLik(Results$atGARCH11)),
as.numeric(logLik(Results$GJRGARCH11))
)
model_names = c("GARCH", "GARCH-t", "skewed-t GARCH", "GJR-GARCH")
LR_garch_tgarch = 2 * (loglik[2] - loglik[1])
p_garch_tgarch = pchisq(LR_garch_tgarch, df = 1, lower.tail = FALSE)
LR_garch_gjr = 2 * (loglik[4] - loglik[1])
p_garch_gjr = pchisq(LR_garch_gjr, df = 1, lower.tail = FALSE)
LR_tgarch_atgarch = 2 * (loglik[3] - loglik[2])
p_tgarch_atgarch = pchisq(LR_tgarch_atgarch, df = 1, lower.tail = FALSE)
cat(sprintf("GARCH vs GARCH-t: LR = %.4f, p = %.4e\n", LR_garch_tgarch, p_garch_tgarch))
cat(sprintf("GARCH vs GJR-GARCH: LR = %.4f, p = %.4e\n", LR_garch_gjr, p_garch_gjr))
cat(sprintf("GARCH-t vs skewed-t: LR = %.4f, p = %.4e\n", LR_tgarch_atgarch, p_tgarch_atgarch))
AIC_vals = 2 * k - 2 * loglik
BIC_vals = k * log(n) - 2 * loglik
ic_table = data.frame(Model = model_names, LogLik = loglik, AIC = AIC_vals, BIC = BIC_vals)
print(ic_table)
cat("Best by AIC:", model_names[which.min(AIC_vals)], "\n")
cat("Best by BIC:", model_names[which.min(BIC_vals)], "\n")GARCH vs GARCH-t: LR = 418.6473, p = 4.8070e-93
GARCH vs GJR-GARCH: LR = 303.7730, p = 4.9635e-68
GARCH-t vs skewed-t: LR = 59.6914, p = 1.1096e-14
Model LogLik AIC BIC
1 GARCH -11923.08 23852.15 23873.45
2 GARCH-t -11713.75 23435.51 23463.90
3 skewed-t GARCH -11683.91 23377.82 23413.31
4 GJR-GARCH -11771.19 23550.38 23578.77
Best by AIC: skewed-t GARCH
Best by BIC: skewed-t GARCHk = [3, 4, 5, 4]
n = len(y)
loglik = [
Results['GARCH11'].loglikelihood,
Results['tGARCH11'].loglikelihood,
Results['atGARCH11'].loglikelihood,
Results['GJRGARCH11'].loglikelihood
]
model_names = ["GARCH", "GARCH-t", "skewed-t GARCH", "GJR-GARCH"]
LR_garch_tgarch = 2 * (loglik[1] - loglik[0])
p_garch_tgarch = stats.chi2.sf(LR_garch_tgarch, df=1)
LR_garch_gjr = 2 * (loglik[3] - loglik[0])
p_garch_gjr = stats.chi2.sf(LR_garch_gjr, df=1)
LR_tgarch_atgarch = 2 * (loglik[2] - loglik[1])
p_tgarch_atgarch = stats.chi2.sf(LR_tgarch_atgarch, df=1)
print(f"GARCH vs GARCH-t: LR = {LR_garch_tgarch:.4f}, p = {p_garch_tgarch:.4e}")
print(f"GARCH vs GJR-GARCH: LR = {LR_garch_gjr:.4f}, p = {p_garch_gjr:.4e}")
print(f"GARCH-t vs skewed-t: LR = {LR_tgarch_atgarch:.4f}, p = {p_tgarch_atgarch:.4e}")
AIC_vals = [2*k[i] - 2*loglik[i] for i in range(4)]
BIC_vals = [k[i]*np.log(n) - 2*loglik[i] for i in range(4)]
for i in range(4):
print(f"{model_names[i]:>15s}: LogLik={loglik[i]:.2f} AIC={AIC_vals[i]:.2f} BIC={BIC_vals[i]:.2f}")
print("Best by AIC:", model_names[int(np.argmin(AIC_vals))])
print("Best by BIC:", model_names[int(np.argmin(BIC_vals))])GARCH vs GARCH-t: LR = 418.7165, p = 4.6432e-93
GARCH vs GJR-GARCH: LR = 304.3105, p = 3.7904e-68
GARCH-t vs skewed-t: LR = 59.8583, p = 1.0194e-14
GARCH: LogLik=-11923.25 AIC=23852.50 BIC=23873.80
GARCH-t: LogLik=-11713.89 AIC=23435.79 BIC=23464.18
skewed-t GARCH: LogLik=-11683.96 AIC=23377.93 BIC=23413.42
GJR-GARCH: LogLik=-11771.10 AIC=23550.19 BIC=23578.59
Best by AIC: skewed-t GARCH
Best by BIC: skewed-t GARCHk = [3, 4, 5, 4]
n = length(y)
loglik = [
loglikelihood(Results_jl["GARCH11"]),
loglikelihood(Results_jl["tGARCH11"]),
loglikelihood(Results_jl["atGARCH11"]),
loglikelihood(Results_jl["GJRGARCH11"])
]
model_names = ["GARCH", "GARCH-t", "skewed-t GARCH", "GJR-GARCH"]
LR_garch_tgarch = 2 * (loglik[2] - loglik[1])
p_garch_tgarch = ccdf(Chisq(1), LR_garch_tgarch)
LR_garch_gjr = 2 * (loglik[4] - loglik[1])
p_garch_gjr = ccdf(Chisq(1), LR_garch_gjr)
LR_tgarch_atgarch = 2 * (loglik[3] - loglik[2])
p_tgarch_atgarch = ccdf(Chisq(1), LR_tgarch_atgarch)
@printf("GARCH vs GARCH-t: LR = %.4f, p = %.4e\n", LR_garch_tgarch, p_garch_tgarch)
@printf("GARCH vs GJR-GARCH: LR = %.4f, p = %.4e\n", LR_garch_gjr, p_garch_gjr)
@printf("GARCH-t vs skewed-t: LR = %.4f, p = %.4e\n", LR_tgarch_atgarch, p_tgarch_atgarch)
AIC_vals = 2 .* k .- 2 .* loglik
BIC_vals = k .* log(n) .- 2 .* loglik
for i in 1:4
@printf("%15s: LogLik=%.2f AIC=%.2f BIC=%.2f\n", model_names[i], loglik[i], AIC_vals[i], BIC_vals[i])
end
println("Best by AIC: ", model_names[argmin(AIC_vals)])
println("Best by BIC: ", model_names[argmin(BIC_vals)])GARCH vs GARCH-t: LR = 418.6445, p = 4.8137e-93
GARCH vs GJR-GARCH: LR = 303.7716, p = 4.9668e-68
GARCH-t vs skewed-t: LR = 59.6935, p = 1.1084e-14
GARCH: LogLik=-11923.08 AIC=23852.15 BIC=23873.45
GARCH-t: LogLik=-11713.75 AIC=23435.51 BIC=23463.90
skewed-t GARCH: LogLik=-11683.91 AIC=23377.82 BIC=23413.31
GJR-GARCH: LogLik=-11771.19 AIC=23550.38 BIC=23578.77
Best by AIC: skewed-t GARCH
Best by BIC: skewed-t GARCHThe three languages reach the same conclusions. Both normal-GARCH comparisons reject the restricted model overwhelmingly and the GARCH-t against skewed Student-t comparison also rejects symmetry. Skewed Student-t GARCH has the lowest AIC and BIC among the four models, so the extra skewness parameter over GARCH-t is justified by both criteria here. The last of these is three separate within-language results rather than one replicated result, because R restricts a Fernández–Steel likelihood while Python and Julia restrict Hansen’s.
Hansen, Bruce E. 1994. “Autoregressive Conditional Density Estimation.” International Economic Review 35, (3): 705–30.