Notes on GARCH

Simulating GARCH

Let’s simulate a GARCH(1,1) process over time

library(fGarch)
rm(list = ls())
seed <- 17
# set the number simulated returns
n <- 10^2
# set the model parameters
con <- 0
w <- 0.1
a <- 0.2
b <- 0.5
par_0 <- c(con,w,a,b)
sigma_t <- sqrt(w/(1-a-b)) # start with initial volatility 
ht_seq <- numeric()
R_seq <- numeric()
set.seed(seed)
for(i in 1:n) {
  R_t1 <- rnorm(1,con,sigma_t)
  R_seq <- c(R_seq,R_t1)
  ht1 <- w + a*(R_t1^2) + b*sigma_t^2
  ht_seq <- c(ht_seq,ht1)
  sigma_t <- sqrt(ht1)
  }
plot(sqrt(ht_seq), type = "l")

To make sure the simulation is correct, let’s try to estimate the coefficients using the data given by R_seq:

g <- garchFit(formula = ~ garch(1,1), data = R_seq,trace = F)
coef_hat <- g@fit$matcoef
coef_hat[,1]

        mu      omega     alpha1      beta1 
0.05969198 0.09070289 0.37858626 0.44038373 

The above is based on 100 observations. One of the properties of MLE is consistency, for this reason, we expect the estimated parameters to converge to the true ones as the sample increases. To demonstrate this, consider the following:

sim_garch_f <- function(n) {
  set.seed(seed)
  sigma_t <- sqrt(w/(1-a-b)) 
  ht_seq <- numeric()
  R_seq <- numeric()
  for(i in 1:n) {
    R_t1 <- rnorm(1,con,sigma_t)
    R_seq <- c(R_seq,R_t1)
    ht1 <- w + a*(R_t1^2) + b*sigma_t^2
    ht_seq <- c(ht_seq,ht1)
    sigma_t <- sqrt(ht1)
  }
  
return(R_seq)
}
n_seq <- 10^(2:4)
R_list <- lapply(n_seq, sim_garch_f)

The above simulates the R_seq for different n values. For each sample size, we can check the difference between the estimated parameters and the true ones

garch_list <- lapply(R_list, function(R) garchFit(formula = ~ garch(1,1), data = R,trace = F))
coef_hat <- sapply(garch_list, function(g) g@fit$matcoef[,1])
colnames(coef_hat) <- paste("n",n_seq,sep = "_")
data.frame(coef_hat)

As we can see, the larger the sample is the more the closer the estimates approximates the true parameters. The mu parameter denotes the average mean return, which by construction we set to zero. As a result, its estimate should converge to zero as \(n \rightarrow \infty\).

Note that the above example is rather a reverse engineering problem. We ``cook’’ the data that imposing it is generated from a GARCH model. After doing so, we estimate a GARCH model using the cooked data. Note that if the simulation is well executed, consistency should hold true. Otherwise, you need to double check your simulation procedure.

Forecasting Volatility

Let’s stick with our cooked data to demonstrate volatility forecasting. As discussed in the slides, the volatility forecast is given by

w <- 0.1
gamma <- 0.7
k <- 20 # number of days ahead
h_for <-  last(ht_seq)
h_seq <- ht_seq
for(i in 1:k){ 
  h_for <- w + gamma*h_for
  h_seq <- c(h_seq,h_for)
}
in_index <- 1:n
out_index <- (n+1):(n+k)
plot(h_seq,col = 1,type = "l")
lines(h_seq[out_index] ~ out_index,col = 2,lwd = 2)

abline(h = w/(1-gamma),lty = 2)
# now add forecast for a different gamma
gamma <- 0.9
w <- 0.1/3
h_for <-  last(ht_seq)
h_seq <- ht_seq
for(i in 1:k){ 
  h_for <- w + gamma*h_for
  h_seq <- c(h_seq,h_for)
}
lines(h_seq,col = 1,type = "l")

lines(h_seq[out_index] ~ out_index,col = 3,lwd = 2)
# finally, add legend
legend("topleft",c("Fitted","Forecast Gamma = 0.7","Forecast Gamma = 0.9","Long Term Variance"),col = 1:3,lwd = c(1,2,2,1),lty = c(1,1,1,2))

In fact, the red and green curves are commonly referred to as the volatility term structure. This indicates how far the current volatility is deviating from its long-term mean and how long it will take to revert back to. As demonstrated above, the green line takes longer to revert since its \(\gamma\) is larger. This means that the persistence of recent shocks (large squared returns) will last longer. On the other hand, the red line has \(\gamma = 0.7\) and, as a result, recent shocks will be less persistence, i.e. it will take the process shorter time period to revert back to its long-term mean.

Let’s repeat the above for a different simulation, however, in which the recent volatility is above the long-term average volatility

seed <- 1713134
# set the number simulated returns
n <- 10^2
# set the model parameters
con <- 0
w <- 0.1
a <- 0.2
b <- 0.5
par_0 <- c(con,w,a,b)
sigma_t <- sqrt(w/(1-a-b)) # start with initial volatility 
ht_seq <- numeric()
R_seq <- numeric()
set.seed(seed)
for(i in 1:n) {
  R_t1 <- rnorm(1,con,sigma_t)
  R_seq <- c(R_seq,R_t1)
  ht1 <- w + a*(R_t1^2) + b*sigma_t^2
  ht_seq <- c(ht_seq,ht1)
  sigma_t <- sqrt(ht1)
  }
w <- 0.1
gamma <- 0.7
k <- 20 # number of days ahead
h_for <-  last(ht_seq)
h_seq <- ht_seq
for(i in 1:k){ 
  h_for <- w + gamma*h_for
  h_seq <- c(h_seq,h_for)
}
in_index <- 1:n
out_index <- (n+1):(n+k)
plot(h_seq,col = 1,type = "l")
lines(h_seq[out_index] ~ out_index,col = 2,lwd = 2)

abline(h = w/(1-gamma),lty = 2)
# now add forecast for a different gamma
gamma <- 0.9
w <- 0.1/3
h_for <-  last(ht_seq)
h_seq <- ht_seq
for(i in 1:k){ 
  h_for <- w + gamma*h_for
  h_seq <- c(h_seq,h_for)
}
lines(h_seq,col = 1,type = "l")

lines(h_seq[out_index] ~ out_index,col = 3,lwd = 2)
# finally, add legend
legend("topleft",c("Fitted","Forecast Gamma = 0.7","Forecast Gamma = 0.9","Long Term Variance"),col = 1:3,lwd = c(1,2,2,1),lty = c(1,1,1,2))

We observe a similar trend, however, now the volatility term structure is down sloping. The reason is the recent volatility is above the long-term mean average and, hence, the forecast should eventually converge to the long-term level, which is lower than the current level.