3.2 GARCH modelling in R

GARCH modelling in R involves following steps:

Step 1 – specification of GARCH model you want to estimate (along with selection of lags $p$ and $q$ as well as distribution of innovations and possibly some external regressors in the mean and/or variance equation)

Step 2 – estimation of GARCH model by maximum likelihood method (numerical quasi–Newton iterative algorithm is applied, such as SQP, BFGS or BHHH). If one algorithm does not converge you can change initial parameter values (close to the true values; zeros are not good starting points) or change stopping criteria (commonly used tolerance is $0.0001$ ) or consider another algorithm

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 log-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 11. 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 print() 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 fitted (estimated) volatilites in–the–sample. Obtain future volatilites

$5-$ steps ahead. Comment diagnostic tests.

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 a long period. The long–run (unconditional) variance of returns is

$\bar{\sigma}^2=\dfrac{\omega}{1-(\alpha_1+\beta_1)}=0.00144$ , and therefore the long–run daily standard deviation is around

$3.8\%$ . It is evident that future volatilities at

$T+1$ ,

$T+2$ ,

$T+3$ ,

$\dots$ are converging to

$0.00144$ .
The weighted Ljung–Box tests on both standardized residuals and their squares show high

$p-$ values at various lags, indicating no significant serial correlation, which supports model adequacy in capturing volatility clustering and heteroscedasticity, also confirmed by ARCH LM tests. The Nyblom stability test shows that both the joint and individual statistics are below critical values, implying parameter stability over time. However, the adjusted Pearson goodness–of–fit test reveals very low

$p-$ values across all groupings, indicating that assumed distribution of innovations does not fit well to the empirical distribution.
Note

$1$ : the Sign Bias Test is missing (not automatically included) in the standard GARCH specification unless asymmetric GARCH model is specified.
Note

$2$ : the quantiles used for post–computing the

$1\%$ VaR are based on the distribution specified in argument distribution.model within ugarchspec() command.

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

# Loading already installed package
library(quantmod)
# Fetching Tesla stock data from Yahoo Finance source and calculating returns
getSymbols("TSLA", src="yahoo", from = "2021-01-01", to = "2024-12-31")
returns <- dailyReturn(Cl(TSLA), type="log")

# Standard 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)) # plotting volatilites in-the-sample
plot(garch1, which =3) # conditional standard deviation vs absolute returns
plot(garch1, which =2) # returns with VaR limits

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

$~~~$

Example 12. 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. To visualize the distribution fit use the QQ plot (Quantile–Quantile plot) of the standardized residuals for each estimated GARCH model. Which distribution is a better fit? Which GARCH model is appropriate according to AIC and BIC? Create 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

Copy the code lines below to the clipboard, paste them into an R Script file, and run them.
According to the QQ plots standard Student’s t–distribution with

$4.94706$ degrees of freedom (shape parameter) is more appropriate than standard normal distribution, while skewed Student’s t–distribution contributes very poorly in terms of the positive skew parameter (although statistically significant). The same conclusion emerges when comparing AIC and BIC values, i.e. model garch2 is the most optimal since it minimizes both criteria.
Note: functions tidy.uGARCHfit and glance.uGARCHfit are defined to create compatibility with modelsummary(), which does not support uGARCHfit class objects. For that purpose broom package is required.

# 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)

# QQ plots
plot(garch1, which = 9)
plot(garch2, which = 9)
plot(garch3, which = 9)

# Tabular display and comparison of multiple models requires installing the "broom" package
# for models which are not suported by modelsummary() command
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 distribution"=garch1,
                  "GARCH(1,1) with t-distribution"=garch2,
                  "GARCH(1,1) with skewed t-distribution"=garch3),
             stars = c("***"=0.01, "**"=0.05, "*"= 0.1),fmt=4)

$~~~$