In this markdown, we will fit different GARCH models to the daily log-returns of GameStop’s stock. Note that any data set can be used!
We initially load our data:

library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(ggplot2)
GME <- read.table("R_Diss\\GME(1).csv",header=T,sep=",")
attach(GME)
n=length(Adj.Close)
ret = log(Adj.Close[-1]/(Adj.Close[-n]-1))

We plot the daily close value over a period of 1 year and plot its respective log-returns:

myplot = ggplot(data = GME)+geom_line(mapping = aes(x=dmy(Date),y=Adj.Close)) + ggtitle("GME") + xlab("Time") + ylab("Adjusted Close ($USD)") 
myplot + theme_bw()

plot(ret)

hist(ret, breaks = 50, probability = T)

GARCH Models We will use the “rugarch” package to fit the GARCH models for the data. We create a 3x3 matrix to investigate all models from GARCH (1,1) up to and including GARCH (3,3) and tabulate their respective AIC values, observing which has the best relative fit.

library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
aic.garch <- matrix(0,3,3)
data <- log(Adj.Close[-1]/(Adj.Close[-n]-1))

for (i in 1:3) {
  for (j in 1:3) {
    garch.spec = ugarchspec(variance.model = list(garchOrder=c(i,j)), 
                            mean.model = list(armaOrder=c(0,0),include.mean = FALSE),
                            distribution.model = "std")
    garch.fit = ugarchfit(spec=garch.spec, data=data,
                          solver.control=list(trace = 1))
    aic.garch[i,j] <- infocriteria(garch.fit)[1]
  }
}
## 
## Iter: 1 fn: -23.8794  Pars:   0.009327  0.487846  0.410652 12.036651
## Iter: 2 fn: -23.8794  Pars:   0.009326  0.487847  0.410652 12.036373
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -24.0474  Pars:  9.310e-03 4.905e-01 4.077e-01 1.293e-10 1.191e+01
## Iter: 2 fn: -24.0474  Pars:  9.310e-03 4.905e-01 4.077e-01 5.442e-11 1.191e+01
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -23.9942  Pars:   0.00938285780  0.49012498616  0.40563712592  0.00000001789  0.00000002966 12.04058159116
## Iter: 2 fn: -23.9942  Pars:   0.009383009001  0.490121795874  0.405634069388  0.000000008152  0.000000010147 12.040947094846
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -25.8360  Pars:   0.015693  0.339536  0.499167  0.006144 15.282924
## Iter: 2 fn: -25.8360  Pars:   0.015693  0.339534  0.499166  0.006141 15.283282
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -25.8360  Pars:   0.015692958980  0.339532681577  0.499163970527  0.006140250751  0.000000004709 15.283538292251
## Iter: 2 fn: -25.8360  Pars:   0.015692850423  0.339534643718  0.499165617095  0.006141491220  0.000000002506 15.283195525393
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -25.8984  Pars:   0.01505135440  0.32317844543  0.50335307099  0.00000094479  0.00000003751  0.02292329427 15.08630637824
## Iter: 2 fn: -25.8984  Pars:   0.0150513308  0.3231812692  0.5033456260  0.0000006217  0.0000000148  0.0229269758 15.0861399766
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -25.8721  Pars:   0.015790495  0.332488563  0.505468619  0.002775972  0.000002948 15.364632233
## Iter: 2 fn: -25.8721  Pars:   0.015790563  0.332491331  0.505465426  0.002776368  0.000001715 15.364767977
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -25.8720  Pars:   0.01572795662  0.33255186642  0.50408986509  0.00000072323  0.00460029761  0.00000001201 15.38352545784
## Iter: 2 fn: -25.8720  Pars:   0.015728015555  0.332550801884  0.504089567806  0.000000766234  0.004599295385  0.000000007878 15.383823568599
## solnp--> Completed in 2 iterations
## 
## Iter: 1 fn: -25.8984  Pars:  1.505e-02 3.232e-01 5.033e-01 2.145e-07 4.371e-07 6.723e-10 2.293e-02 1.509e+01
## Iter: 2 fn: -25.8984  Pars:  1.505e-02 3.232e-01 5.033e-01 2.042e-07 4.133e-07 4.061e-12 2.293e-02 1.509e+01
## solnp--> Completed in 2 iterations
aic.garch
##            [,1]       [,2]       [,3]
## [1,] -0.1584014 -0.1517723 -0.1433797
## [2,] -0.1660239 -0.1580557 -0.1505847
## [3,] -0.1583437 -0.1503744 -0.1426166
min(aic.garch)
## [1] -0.1660239

We see that the GARCH (2,1) model has the best relative fit, so we will begin to define it in order to use it in future steps:

garch21.spec = ugarchspec(variance.model = list(garchOrder=c(2,1)), 
                          mean.model = list(armaOrder=c(0,0),include.mean = FALSE),
                          distribution.model = "std")
garch21.fit = ugarchfit(spec=garch21.spec, data=data,
                        solver.control=list(trace = 1))
## 
## Iter: 1 fn: -25.8360  Pars:   0.015693  0.339536  0.499167  0.006144 15.282924
## Iter: 2 fn: -25.8360  Pars:   0.015693  0.339534  0.499166  0.006141 15.283282
## solnp--> Completed in 2 iterations
garch21.fit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## omega   0.015693    0.005364 2.925520 0.003439
## alpha1  0.339534    0.152346 2.228699 0.025834
## alpha2  0.499166    0.192702 2.590345 0.009588
## beta1   0.006141    0.127789 0.048059 0.961670
## shape  15.283282   10.696923 1.428755 0.153075
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## omega   0.015693    0.006372 2.462881 0.013783
## alpha1  0.339534    0.124029 2.737547 0.006190
## alpha2  0.499166    0.170535 2.927061 0.003422
## beta1   0.006141    0.072185 0.085078 0.932199
## shape  15.283282   14.355750 1.064611 0.287052
## 
## LogLikelihood : 25.836 
## 
## Information Criteria
## ------------------------------------
##                       
## Akaike       -0.166024
## Bayes        -0.095796
## Shibata      -0.166797
## Hannan-Quinn -0.137762
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      6.526 1.063e-02
## Lag[2*(p+q)+(p+q)-1][2]    11.538 7.499e-04
## Lag[4*(p+q)+(p+q)-1][5]    25.797 4.783e-07
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.5883  0.4431
## Lag[2*(p+q)+(p+q)-1][8]     2.1112  0.8414
## Lag[4*(p+q)+(p+q)-1][14]    4.5940  0.8179
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]    0.8076 0.500 2.000  0.3688
## ARCH Lag[6]    0.8442 1.461 1.711  0.7928
## ARCH Lag[8]    2.5089 2.368 1.583  0.6387
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.3948
## Individual Statistics:             
## omega  0.1229
## alpha1 0.3379
## alpha2 0.1599
## beta1  0.2789
## shape  2.3811
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.28 1.47 1.88
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value    prob sig
## Sign Bias           1.7679 0.07832   *
## Negative Sign Bias  1.1271 0.26079    
## Positive Sign Bias  0.1874 0.85152    
## Joint Effect        4.0216 0.25914    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     257.8    8.250e-44
## 2    30     272.6    1.985e-41
## 3    40     281.7    1.549e-38
## 4    50     303.8    1.874e-38
## 
## 
## Elapsed time : 0.1671369

Here we obtain the MLE parameters. Using these MLEs, we can obtain the conditional variance for the GARCH(2,1) model and use it to plot the time-varying value at risk:

data <- log(Adj.Close[-1]/(Adj.Close[-n]-1))
log.likt <- function(par){
  mu <- par[1]; sigma <- exp(par[2]); nu <- exp(par[3])
  out <- -sum( dt( (data -  mu)/sigma, df = nu , log=TRUE) - log(sigma) )
  return(out)
}

OPTt <- optim(par = c(0,0,0), fn = log.likt)
MLEt <- c(OPTt$par[1], exp(OPTt$par[2]), exp(OPTt$par[3]))

p <- qplot(y = ret , x = dmy(Date[-1]) , geom = 'point') + geom_point(colour = 'grey30' , size = 2) + 
  geom_line(aes(x = dmy(Date[-1]), y = garch21.fit@fit$sigma*qt(0.05, MLEt[3])+MLEt[1] ) , colour = 'red') +
  geom_hline(yintercept = MLEt[2]*qt(0.05, MLEt[3])+MLEt[1] , colour = 'darkgreen' , size = 1.2) + theme_bw() + xlab("Time") + ylab("Daily Logarithmic Returns") + ggtitle("GARCH(2,1) VaR")

p + theme(
  plot.title = element_text(color="Black", size=20, face="bold.italic"),
  axis.title.x = element_text(color="black", size=20, face="bold"),
  axis.title.y = element_text(color="black", size=20, face="bold")
)

For further reading, please refer to my undergraduate dissertation: