3.1 GARCH modelling in R

GARCH modelling in R involves following steps:

Step 1 – specification of GARCH model you want to use (along with selection of time lags \(p\) and \(q\) as well as returns distribution and possibly some external regressors in the mean or variance equations)

Step 2 – estimation of GARCH model by maximum likelihood method (numerical quasi–Newton algorithm is applied, such as BFGS or BHHH)

Step 3 – diagnostic checking if GARCH model meets the assumptions and parameters constraints

Step 4 – if initially specified model in step \(1\) does not meet the assumptions or parameters constraints, try with another GARCH specification and repeat steps \(1-3\)

Step 5 – after estimating several GARCH models choose the most appropriate one and use it for volatility predictions

The most appropriate GARCH model should: (a) fit data better than other models in terms of AIC or BIC infromation criteria, (b) meet parameter constraints, (c) meet all the assumptions (no ARCH effects are left in the residuals or assumed ditribution should be followed), and (d) estimated parameters are statistically significant (with or without robust standard errors)

TABLE 3.1: Commands for univariate GARCH modelling within package `rugarch`
Command	Description
`ugarchspec()`	univariate GARCH specification prior to estimation
`ugarchfit()`	estimating specified GARCH model
`print()`	displays the results of estimated GARCH model in console window
`fitted()`	produces predicted mean values of returns
`sigma()`	produces predicted volatilites in-the-sample and out-of-sample
`uncvariance()`	reports unconditional variance in the long-run
`ugarchforecast()`	forecasts the future volatility for h-steps ahead
`infocriteria()`	exctrcts information criteria, such as AIC, BIC,…
`residuals()`	extracts residuals from estimated GARCH model
`likelihood()`	extracts the maximum value of the likelihhod function
`plot()`	displys \(12\) possible plots by selecting numbers from \(1\) do \(12\)
`ugarchroll()`	re-estimates GARCH model within rolling window
`persistence()`	extracts the volatility persistence

Example 9. Using Tesla returns estimate a standard GARCH(\(1,1\)) model assuming normal distribution of innovations, and only constant term in the mean equation. Prior to estimation specify a GARCH model as model.spec object by ugarchspec() command. Afterwards, estimate the same model using garchfit() command as garch1 object. Display the results of object garch1 by employing summary() command. Write both equations analytically. Are parameters constraints meet? Are parameters statistically significant? What is the volatility persistence? Obtain the unconditional (long-run) volatility. Plot etimated volatilites in-the-sample. Obtain future volatilites \(5-\)steps ahead.

Solution

Copy the code lines below to the clipboard, paste them into an R Script file, and run them.
Standard GARCH(\(1,1\)) model in this example has \(1\) parameter in the mean equation and \(3\) parameters in the variance equation, thus \[\begin{equation} \begin{aligned} r_t=0.000887+u_t,~~~~~~~~~~\dfrac{u_t}{\sqrt{\sigma^2_t}}\sim N(0,1) \\ \sigma^2_t=0.000035+0.033023 u^2_{t-1}+0.942876 \sigma^2_{t-1} \end{aligned} \end{equation}\] A persistence close to \(1\) (here \(\alpha_1+\beta_1=0.9759\)) means that volatility shocks decay very slowly, i.e. volatility stays high or low for long periods. The long–run (unconditional) variance of returns is \(\bar{\sigma}^2=\dfrac{\omega}{\alpha_1+\beta_1}=0.00144\), and therefore the long–run daily volatility is around \(3.8\%\). It is evident that future volatilities at \(T+1\), \(T+2\), \(T+3\), \(\dots\) are converging to \(0.00144\).

# Installing and loading the package for univariate GARCH models
install.packages("rugarch")
library(rugarch)

# GARCH model specification
model.spec <- ugarchspec(mean.model = list(armaOrder = c(0,0), include.mean = TRUE, # conditional mean equation
                                           external.regressors = NULL), # additional variables in the mean equation
                         variance.model = list(model = "sGARCH", garchOrder = c(1,1), # conditional variance equation
                                               submodel = NULL, external.regressors = NULL), # submodel and additional variables
                         distribution.model = "norm") # standard normal distribution of innovations

# ML estimation of specified GARCH model
garch1 <- ugarchfit(spec = model.spec, data = returns)

# Displaying the results of estimated GARCH model
print(garch1)

# Individual information of the estimated GARCH model can also be extracted
garch1@fit$matcoef # displays coefficients with standard errors, t-statistics, and p-values
garch1@fit$robust.matcoef # displays coefficients with robust standard errors

persistence(garch1) # volatility persistence
uncvariance(garch1) # unconditinal volatility

plot(sigma(garch1)) # volatilites in-the-sample
plot(garch1, which =2) # conditional standard deviation vs absolute returns

# Forecasting volatility
forecast <- ugarchforecast(garch1, n.ahead = 5) # 5-days ahead
forecast_volatility <- sigma(forecast)
print(forecast_volatility)

\(~~~\)

Example 10. Using Tesla returns estimate two additional GARCH(\(1,1\)) models by changing theoretical distribution of innovations. Specify model.spec2 to account for standard t–distribution (argument distribution.model=“std”) and model.spec3 to account for skewed t–distribution (argument distribution.model=“sstd”). Afterwards, estimate two models as to separate objects garch2 and garch3 using garchfit() command. Display the results of both objects by employing print() command. Crate two functions to “tidy” and “glance” the GARCH objects so that you can easily compare the results of several GARCH models in a single well–formated table within modelsummary() command.

Solution

# Two additional GARCH specifications with respect to theoretical distributioin of innovations
model.spec2 <- ugarchspec(mean.model = list(armaOrder = c(0,0), include.mean = TRUE, 
                                           external.regressors = NULL), 
                         variance.model = list(model = "sGARCH", garchOrder = c(1,1), 
                                               submodel = NULL, external.regressors = NULL), # 
                         distribution.model = "std") # standard t-distribution

model.spec3 <- ugarchspec(mean.model = list(armaOrder = c(0,0), include.mean = TRUE, 
                                           external.regressors = NULL), 
                         variance.model = list(model = "sGARCH", garchOrder = c(1,1), 
                                               submodel = NULL, external.regressors = NULL), # 
                         distribution.model = "sstd") # skewed t-distribution

# ML estimation of both new models
garch2 <- ugarchfit(spec = model.spec2, data = returns)
garch3 <- ugarchfit(spec = model.spec3, data = returns)

# Displaying the results of estimated GARCH models
print(garch2)
print(garch3)

# Tabular display and comparison of multiple models requires installing the "broom" package
install.packages("broom")
library(broom)
# Two functions are defined that return two data frames: tidy and glance
# This is required as modesummary() command does not support GARCH model
tidy.uGARCHfit<- function(x, ...) {
  s <- x@fit$matcoef
  ret <- data.frame(
    term      = row.names(s),
    estimate  = s[, " Estimate"],
    std.error = s[, " Std. Error"],
    statistic = s[, " t value"],
    p.value =  s[,"Pr(>|t|)"])
  ret
}

glance.uGARCHfit <- function(x, ...) {
  s <- x@fit
  ret <- data.frame(
    Likelihood = round(s$LLH, digits=2),
    Persistence = round(s$persistence, digits=4),
    Akaike=round(infocriteria(x)[1], digits=4),
    Bayes=round(infocriteria(x)[2], digits=4))
  ret
}

modelsummary(list("GARCH(1,1) with normal dist."=garch1,
                  "GARCH(1,1) with t-dist."=garch2,
                  "GARCH(1,1) with skewed t-dist."=garch3),
             stars=TRUE,fmt=4)

\(~~~\)