Fit a GARCH (1,1) - model with covariates in R

I have some experiences with time series modelling, in the form of simple ARIMA models and so on. Now I have some data that exhibits volatility clustering, and I would like to try to start with fitting a GARCH (1,1) model on the data.

I have a data series and a number of variables I think influence it. So in basic regression terms, it looks like:

$$y_t = \alpha + \beta_1 x_{t1} + \beta_2 x_{t2} + \epsilon_t .$$

But I am at a complete loss at how to implement this into a GARCH (1,1) - model? I've looked at the rugarch-package and the fGarch-package in R, but I haven't been able to do anything meaningful besides the examples one can find on the internet.

Here is an example of implementation using the rugarch package and with to some fake data. The function ugarchfit allows for the inclusion of external regressors in the mean equation (note the use of external.regressors in fit.spec in the code below).

To fix notations, the model is \begin{align*} y_t &= \lambda_0 + \lambda_1 x_{t,1} + \lambda_2 x_{t,2} + \epsilon_t, \\ \epsilon_t &= \sigma_t Z_t , \\ \sigma_t^2 &= \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma_{t-1}^2 , \end{align*} where $x_{t,1}$ and $x_{t,2}$ denote the covariate at time $t$, and with the "usual" assumptions/requirements on parameters and the innovation process $Z_t$.

The parameter values used in the example are as follows.

## Model parameters
nb.period  <- 1000
omega      <- 0.00001
alpha      <- 0.12
beta       <- 0.87
lambda     <- c(0.001, 0.4, 0.2)


The image below shows the series of covariate $x_{t,1}$ and $x_{t,2}$ as well as the series $y_t$. The R code used to generate them is provided below.

## Dependencies
library(rugarch)

## Generate some covariates
set.seed(234)
ext.reg.1 <- 0.01 * (sin(2*pi*(1:nb.period)/nb.period))/2 + rnorm(nb.period, 0, 0.0001)
ext.reg.2 <- 0.05 * (sin(6*pi*(1:nb.period)/nb.period))/2 + rnorm(nb.period, 0, 0.001)
ext.reg   <- cbind(ext.reg.1, ext.reg.2)

## Generate some GARCH innovations
sim.spec    <- ugarchspec(variance.model     = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model         = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "norm",
fixed.pars         = list(omega = omega, alpha1 = alpha, beta1 = beta))
path.sgarch <- ugarchpath(sim.spec, n.sim = nb.period, n.start = 1)
epsilon     <- as.vector(fitted(path.sgarch))

## Create the time series
y <- lambda[1] + lambda[2] * ext.reg[, 1] + lambda[3] * ext.reg[, 2] + epsilon

## Data visualization
par(mfrow = c(3,1))
plot(ext.reg[, 1], type = "l", xlab = "Time", ylab = "Covariate 1")
plot(ext.reg[, 2], type = "l", xlab = "Time", ylab = "Covariate 2")
plot(y, type = "h", xlab = "Time")
par(mfrow = c(1,1))


A fit is done with ugarchfit as follows.

## Fit
fit.spec <- ugarchspec(variance.model     = list(model = "sGARCH",
garchOrder = c(1, 1)),
mean.model         = list(armaOrder = c(0, 0),
include.mean = TRUE,
external.regressors = ext.reg),
distribution.model = "norm")
fit      <- ugarchfit(data = y, spec = fit.spec)


Parameter estimates are

## Results review
fit.val     <- coef(fit)
fit.sd      <- diag(vcov(fit))
true.val    <- c(lambda, omega, alpha, beta)
fit.conf.lb <- fit.val + qnorm(0.025) * fit.sd
fit.conf.ub <- fit.val + qnorm(0.975) * fit.sd
> print(fit.val)
#     mu       mxreg1       mxreg2        omega       alpha1        beta1
#1.724885e-03 3.942020e-01 7.342743e-02 1.451739e-05 1.022208e-01 8.769060e-01
> print(fit.sd)
#[1] 4.635344e-07 3.255819e-02 1.504019e-03 1.195897e-10 8.312088e-04 3.375684e-04


And corresponding true values are

> print(true.val)
#[1] 0.00100 0.40000 0.20000 0.00001 0.12000 0.87000


The following figure shows a parameter estimates with 95% confidence intervals, and the true values. The R code used to generate it is provided is below.

plot(c(lambda, omega, alpha, beta), pch = 1, col = "red",
ylim = range(c(fit.conf.lb, fit.conf.ub, true.val)),
xlab = "", ylab = "", axes = FALSE)
box(); axis(1, at = 1:length(fit.val), labels = names(fit.val)); axis(2)
points(coef(fit), col = "blue", pch = 4)
for (i in 1:length(fit.val)) {
lines(c(i,i), c(fit.conf.lb[i], fit.conf.ub[i]))
}
legend( "topleft", legend = c("true value", "estimate", "confidence interval"),
col = c("red", "blue", 1), pch = c(1, 4, NA), lty = c(NA, NA, 1), inset = 0.01)

• how did you estimate parameters( lambda,omega,alpha,beta)??
– chs
Nov 28, 2017 at 5:38
• @chs The parameter estimates were obtained with the ugarchfit function. Nov 30, 2017 at 20:08
• is there a similar guide for the fgarch package? Jun 16, 2020 at 16:03