Simple GARCH Models In R

This page shows you how to fit a GARCH(1,1) model using the rugarch package in R. We will then retrieve the fitted standard deviations, which can be used as a (potentially not very good) estimate of volatility for whatever external exercise you happen to be doing.

Load the package and load the DJIA daily return data:

library(rugarch)
data(dji30ret)

Create the default GARCH specification, which is an ARMA(1,1) model of the conditional mean and a GARCH(1,1) model of the variance of the error term. You can read the documentation for ugarchspec-methods to see the many other options available for the model specification.

spec <- ugarchspec()

Estimate the model:

fit <- ugarchfit(spec=spec, data=dji30ret[,1])
fit

This prints a lot of output. We are interested in the conditional variance equation, for which the parameters are omega, alpha1, and beta1, in the “Optimal Parameters” section of the output. The model is \[\sigma_{t}^{2} = \omega + \alpha_{1} \varepsilon_{t-1}^{2} + \beta_{1} \sigma_{t-1}^{2}\] If the variance of the residual did not change through time (at least not in a way that is predictable), \(\sigma_{t}^{2} = \omega\), which is constant. We have the output:

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000646    0.000254   2.5474 0.010852
ar1    -0.540238    0.107065  -5.0459 0.000000
ma1     0.599170    0.101543   5.9006 0.000000
omega   0.000005    0.000002   2.0265 0.042718
alpha1  0.051440    0.007274   7.0719 0.000000
beta1   0.939775    0.009838  95.5293 0.000000

Which implies the following model of the variance: \[\sigma_{t}^{2} = 0.000005 + 0.051440 \varepsilon_{t-1}^{2} + 0.939775 \sigma_{t-1}^{2}\] Although this is more significant digits than I’d report in a paper, we are going to calculate fitted values of the volatility, so we want to keep all six digits for accuracy.

We can also get the fitted variances and the residual terms. We’ll include only the first few observations here.

> s <- sigma(fit)
> s[1:5]
                          [,1]
1970-01-01 18:00:00 0.02300389
1970-01-02 18:00:00 0.02311680
1970-01-03 18:00:00 0.02323291
1970-01-04 18:00:00 0.02263498
1970-01-05 18:00:00 0.02259049


> r <- residuals(fit)
> r[1:5]
                            [,1]
1970-01-01 18:00:00 -0.025111960
1970-01-02 18:00:00  0.025300044
1970-01-03 18:00:00 -0.002936391
1970-01-04 18:00:00  0.021692065
1970-01-05 18:00:00 -0.005649212

The series s is the estimate of \(\sigma^{2}\) at each point in time. We can calculate it ourselves. For the second observation, we have lagged standard deviation \(\sigma_{t-1} = 0.02300389\) and lagged residual \(\varepsilon_{t-1} = -0.025111960\). The prediction for the second observation’s residual variance is therefore

> 0.000005 + 0.051440*(-0.025111960)^2 + 0.939775*0.02300389^2
[1] 0.0005347478

Taking the square root, we get something close (accounting for rounding error) to the 0.02311680 for the second observation:

> sqrt(0.0005347478)
[1] 0.02312461

We can plot the fitted \(\sigma\) values through time, keeping in mind that the dates on the plot will not be correct, because the original dataset dji30ret does not have any time series properties.

plot(sigma(fit))

The big spike at the beginning of the dataset captures the effect of “Black Monday” on October 1987 and the spike at the end captures the financial crisis in 2008 and 2009.

Last updated: July 28, 2017