22  GARCH(1,1) moments

Setup
library(tidyverse)
# **********************************************************************
#                     Simulate GARCH (1,1) 
# **********************************************************************
# Simulate GARCH process and variance 
simulate_GARCH11 <- function(B = 100, n = 1000, sigma2_0, omega, alpha, beta, mu, sd, p){
  scenarios <- list()
  for(b in 1:B){
    Bt <- rbinom(n, 1, p[1])
    ut <- rnorm(n, mu[1], sd[1])*Bt + rnorm(n, mu[2], sd[2])*(1-Bt)
    yt <- c(sigma2_0 * ut)
    sigma2_t <- c(sigma2_0)
    for(i in 2:n){
      sigma2_t[i] <- omega + alpha * (yt[i-1] - sqrt(sigma2_t[i-1])*sum(mu*p))^2 + beta * sigma2_t[i-1]
      yt[i] <- sqrt(sigma2_t[i]) * ut[i]
    }
    scenarios[[b]] <- dplyr::tibble(t = 1:n, j = b, x = yt, sigma2 = sigma2_t)
  }
  return(scenarios)
}
# **********************************************************************
#                     Simulate GARCH (1,1) Moments
# **********************************************************************
# Simulate GARCH and compute its moments
simulate_GARCH_moments <- function(B = 100, n = 1000, sigma2_0, omega, alpha, beta, mu, sd, p){
  m1 <- m2 <- m3 <- m4 <- m32 <- m12 <- v <- v12 <- c()
  for(b in 1:B){
    Bt <- rbinom(n, 1, p[1])
    ut <- rnorm(n, mu[1], sd[1])*Bt + rnorm(n, mu[2], sd[2])*(1-Bt)
    sigma2_t <- c(sigma2_0)
    for(i in 2:n){
      sigma2_t[i] <- omega + alpha * ut[i-1]^2 + beta * sigma2_t[i-1]
      ut[i] <- sqrt(sigma2_t[i]) * ut[i]
    }
    m1[b] <- mean(sigma2_t)
    m12[b] <- mean(sigma2_t^(1/2))
    m2[b] <- mean(sigma2_t^2)
    m32[b] <- mean(sigma2_t^(3/2))
    m3[b] <- mean(sigma2_t^3) 
    m4[b] <- mean(sigma2_t^4)
    v[b] <- var(sigma2_t)
    v12[b] <- var(sigma2_t^(1/2)) 
  }
  dplyr::tibble(m1 = m1, m2 = m2, m3 = m3, m4 = m4, m32 = m32, m12 = m12, v = v, v12 = v12) %>%
    dplyr::summarise_all(mean)
}
# **********************************************************************
#                              Inputs
# **********************************************************************
# Forecast horizon
h <- 30000
B <- 500
n <- h
# GARCH parameters
alpha <- 0.08
beta <- 0.35
omega <- 1.2*(1 - alpha - beta)
# Mixture moments
mu <- c(-0.6, 0.4)
sd <- c(0.8, 1.2)
p <- c(0.4, 0.6)
# Moments of a Gaussian Mixture 
e_x2 = sum((mu^2 + sd^2)*p)
e_x4 = sum((3 * sd^4 + 6 * mu^2 * sd^2 + mu^4) * p)
# Initial values for variance 
sigma2_t <- 1.2
sigma4_t <- sigma2_t^2 
# **********************************************************************
# Simulate GARCH and compute moments
mom <- simulate_GARCH_moments(B = B, n = n, sigma2_t, omega, alpha, beta, mu, sd, p)

Let’s consider a very general setup for GARCH(1,1) process where we do not assume that the GARCH residuals

  1. Has expected value equal to zero, i.e. \(\mathbb{E}\{u_t\} = 0\).
  2. Has NOT necessary a second moment that is constant and equal to one, i.e. \(\mathbb{E}\{u_t^2\}\) possibly time dependent (deterministic).
  3. Has NOT necessary a fourth moment that is constant and equal to 3, i.e. \(\mathbb{E}\{u_t^4\}\) possibly time dependent (deterministic).

The process we refer to has the form of a standard GARCH(1,1), i.e.  \[ \begin{aligned} {} & Y_t = \sigma_t u_t \\ & \sigma_t^2 = \omega + \alpha_1 u_{t-1}^2 + \beta_1 \sigma_{t-1}^2 \end{aligned} \] where \(u_t \sim \text{WN}(\mu_t, \sigma_t)\).

22.1 First moment \(\sigma_t^2\)

22.1.1 Short-term

Given the information at time \(t-1\), the expected value of the GARCH variance \(h\)-steps ahead, with \(h\ge1\), can be expanded as: \[ \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\} = \omega\left( 1 + \sum_{j=1}^{h-1} \prod_{i=1}^{j} \lambda_{t+h-i}\right) + \sigma_t^2 \prod_{j=1}^{h} \lambda_{t+h-j} \text{,} \tag{22.1}\] with the convention that when \(h=1\) the summation is equal to zero and where \[ \lambda_{t+h-i} = \alpha_1 \mathbb{E}\{u_{t+h-i}^2\} + \beta_1 \text{.} \tag{22.2}\] The iteration between two consecutive times \(t+s\) and \(t+s-1\) reads \[ \mathbb{E}\{\sigma_{t+s}^2 \mid \mathcal{F}_{t-1} \} = \omega + \lambda_{t+s-1} \mathbb{E}\{\sigma_{t+s-1}^2 \mid \mathcal{F}_{t-1}\} \text{.} \]

Proof. We start the proof by recalling that the conditional variance at time \(t\) is included in the information set \(\mathcal{F}_{t-1}\), i.e it is known, therefore at time \(t\) one obtain simply \[ \sigma_{t}^2 = \omega + \alpha_1 y_{t-1}^2 + \beta_1 \sigma_{t-1}^2 \text{.} \] Then, the GARCH variance at time \(t+1\) reads: \[ \sigma_{t+1}^2 = \omega + \alpha_1 y_{t}^2 + \beta_1 \sigma_{t}^2 \text{.} \] Taking the conditional expectation and recalling that \(y_{t} = \sigma_t u_t\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+1}^2 \mid \mathcal{F}_{t-1} \} {} & = \omega + \alpha_1 \mathbb{E}\{y_{t}^2\} + \beta_1 \sigma_{t}^2 = \\ & = \omega + \alpha_1 \sigma_t^2 \mathbb{E}\{u_{t}^2\} + \beta_1 \sigma_{t}^2 = \\ & = \omega + (\alpha_1 \mathbb{E}\{u_{t}^2\} + \beta_1) \sigma_{t}^2 \end{aligned} \] Therefore, denoting with \(\lambda_t = \alpha_1 \mathbb{E}\{u_t^2\} + \beta_1\) one obtain a recursive expression, i.e.  \[ \mathbb{E}\{\sigma_{t+1}^2 \mid \mathcal{F}_{t-1} \} = \omega + \lambda_t \sigma_{t}^2 \text{.} \] Iterating the expectation at time \(t+2\) gives \[ \begin{aligned} \mathbb{E}\{\sigma_{t+2}^2 \mid \mathcal{F}_{t-1} \} &= \omega + \lambda_{t+1} \mathbb{E}\{\sigma_{t+1}^2 \mid \mathcal{F}_{t-1}\} \\ &= \omega + \lambda_{t+1}\left(\omega + \lambda_{t} \sigma_t^2 \right) \\ &= \omega (1 + \lambda_{t+1}) + \lambda_{t+1}\lambda_{t} \sigma_t^2 \end{aligned} \] Then, at time \(t+3\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+3}^2 \mid \mathcal{F}_{t-1} \} &= \omega + \lambda_{t+2}\,\mathbb{E}\{\sigma_{t+2}^2 \mid \mathcal{F}_{t-1} \} \\ &= \omega + \lambda_{t+2}\left(\omega (1+\lambda_{t+1}) + \lambda_{t+1} \lambda_t \sigma_t^2 \right) \\ &= \omega\left(1 + \lambda_{t+2} + \lambda_{t+2} \lambda_{t+1}\right) + \lambda_{t+2} \lambda_{t+1} \lambda_t \sigma_t^2 \end{aligned} \]

Example 22.1 Expanding the expression in Equation 22.1 with \(h=1\) one obtain \[ \begin{aligned} \mathbb{E}\{\sigma_{t+1}^2 \mid \mathcal{F}_{t-1} \} & {} = \omega\left( 1 + \sum_{j=1}^{0} \prod_{i=1}^{j} \lambda_{t+1-i}\right) + \sigma_t^2 \prod_{j=1}^{1} \lambda_{t+1-j} = \\ & = \omega + \sigma_t^2 \lambda_{t} \end{aligned} \] with \(h=2\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+2}^2 \mid \mathcal{F}_{t-1} \} & {} = \omega\left( 1 + \sum_{j=1}^{1} \prod_{i=1}^{j} \lambda_{t+2-i}\right) + \sigma_t^2 \prod_{j=1}^{2} \lambda_{t+2-j} = \\ & = \omega (1 + \lambda_{t+1}) + \sigma_t^2 (\lambda_{t+1} \lambda_{t}) \end{aligned} \] with \(h=3\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+3}^2 \mid \mathcal{F}_{t-1} \} & {} = \omega\left( 1 + \sum_{j=1}^{2} \prod_{i=1}^{j} \lambda_{t+3-i}\right) + \sigma_t^2 \prod_{j=1}^{3} \lambda_{t+3-j} = \\ & = \omega (1 + \lambda_{t+2} + \lambda_{t+1} \lambda_{t+2}) + \sigma_t^2 (\lambda_{t+2} \lambda_{t+1} \lambda_{t}) \end{aligned} \] and so on.

Functions expectation of GARCH variance.
#' Formula for conditional first moment of GARCH variance
#' 
#' @param h integer, number of steps ahead. 
#' @param omega numeric, intercept of the model. 
#' @param alpha numeric, arch parameter of the model. 
#' @param beta numeric, garch parameter of the model. 
#' @param e_x2 numeric, second moment of ut. 
#' @param e_x4 numeric, fourth moment of ut. 
#' @param sigma2_t numeric, last value of the variance at time t-1. 
# Conditional first moment GARCH variance (exact)
e_sigma2_h_mix <- function(h, omega, alpha, beta, e_x2 = 1, sigma2_t){
  sigma2_h <- c(sigma2_t)
  # Second moment
  m2 <- e_x2
  if (length(e_x2) == 1 & h > 1){
    m2 <- rep(e_x2, h)
  }
  # Derived quantities
  lambda <- alpha * m2 + beta
  lambda_prod <- c(1, lambda[h:1][-1])
  #lambda_prod <- cumprod(lambda_prod)
  sigma2_h <- c(sigma2_t, omega * cumsum(cumprod(lambda_prod)) + sigma2_t * cumprod(lambda))
  names(sigma2_h) <- paste0("t+", 0:h)
  return(sigma2_h)
}
#' Iteration for conditional first moment of GARCH variance
#' 
#' @param h integer, number of steps ahead. 
#' @param omega numeric, intercept of the model. 
#' @param alpha numeric, arch parameter of the model. 
#' @param beta numeric, garch parameter of the model. 
#' @param e_x2 numeric, second moment of ut. 
#' @param e_x4 numeric, fourth moment of ut. 
#' @param sigma2_t numeric, last value of the variance at time t-1. 
e_sigma2_h_iter <- function(h, omega, alpha, beta, e_x2 = 1, sigma2_t){
  # Manual iteration 
  sigma2_h <- c(sigma2_t)
  # Second moment
  m2 <- e_x2
  if (length(e_x2) == 1 & h > 1){
    m2 <- rep(e_x2, h)
  }
  for(i in 1:h){
    sigma2_h[i+1] <- omega + (alpha * m2[i] + beta) * sigma2_h[i]
  }
  names(sigma2_h) <- paste0("t+", 0:h)
  return(sigma2_h)
}
Moment Step Formula Iteration Difference. (%)
\(\mathbb{E}\{\sigma_{t+0}^2\mid \mathcal{F}_{t-1}\}\) 0 1.200 1.200 0%
\(\mathbb{E}\{\sigma_{t+1}^2\mid \mathcal{F}_{t-1}\}\) 1 1.235 1.235 0%
\(\mathbb{E}\{\sigma_{t+2}^2\mid \mathcal{F}_{t-1}\}\) 2 1.250 1.250 0%
\(\mathbb{E}\{\sigma_{t+3}^2\mid \mathcal{F}_{t-1}\}\) 3 1.258 1.258 0%
\(\mathbb{E}\{\sigma_{t+4}^2\mid \mathcal{F}_{t-1}\}\) 4 1.261 1.261 0%
\(\mathbb{E}\{\sigma_{t+5}^2\mid \mathcal{F}_{t-1}\}\) 5 1.263 1.263 0%
\(\mathbb{E}\{\sigma_{t+6}^2\mid \mathcal{F}_{t-1}\}\) 6 1.263 1.263 0%
\(\mathbb{E}\{\sigma_{t+7}^2\mid \mathcal{F}_{t-1}\}\) 7 1.264 1.264 0%
\(\mathbb{E}\{\sigma_{t+8}^2\mid \mathcal{F}_{t-1}\}\) 8 1.264 1.264 0%
\(\mathbb{E}\{\sigma_{t+9}^2\mid \mathcal{F}_{t-1}\}\) 9 1.264 1.264 0%
\(\mathbb{E}\{\sigma_{t+10}^2\mid \mathcal{F}_{t-1}\}\) 10 1.264 1.264 0%
Table 22.1: Forecasted expectation of GARCH(1,1) variance with iteration and with formula.

22.1.2 Long-term

If \(\mathbb{E}\{u_t^2\}\) is constant for all \(t\) then, \(\lambda = \alpha_1 \mathbb{E}\{u_t^2\} + \beta_1\), became a constant and the formula simplifies to
\[ \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1} \} = \sigma_{\infty}^2 + \lambda^h (\sigma_t^2 - \sigma_{\infty}^2) \text{,} \tag{22.3}\] where \(\sigma_{\infty}^2\) denotes the long-term expected GARCH variance as \(h\to\infty\), i.e.  \[ \lim_{h\to\infty} \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1} \} = \frac{\omega}{1- \lambda} = \sigma_{\infty}^2 \text{.} \tag{22.4}\] It follows that, under the standard assumption that \(u_t \sim \text{WN}(0, 1)\), then one obtain the classic expression \[ \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\} = \sigma_{\infty}^2 + (\alpha_1 + \beta_1)^h (\sigma_t^2 - \sigma_{\infty}^2) \tag{22.5}\] where \(\sigma_{\infty}^2\) denotes the unconditional expectation as Equation 22.4 with \(\mathbb{E}\{u_t^2\} = 1\).

Proof. Let’s verify the formula in Equation 22.5 for constant \(\lambda_t\) (Equation 22.2) for all \(t\), i.e.  \[ \lambda = \alpha_1 \mathbb{E}\{u_t^2\} + \beta_1 \] In this case the iterative formula (Equation 22.1) simplifies, i.e.  \[ \begin{aligned} \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1} \} & {} = \omega\left(1 + \sum_{j=1}^{h-1} \prod_{i=1}^{j}\lambda \right) + \sigma_t^2 \prod_{j=1}^{h} \lambda = \\ & = \omega \sum_{j=0}^{h-1} \lambda^j + \sigma_t^2 \lambda^h = \\ & = \omega \left(\frac{1- \lambda^h}{1 - \lambda} \right) + \sigma_t^2 \lambda^h \\ & = \frac{\omega}{1- \lambda} + \lambda^h \left(\sigma_t^2 - \frac{\omega}{1-\lambda} \right) \end{aligned} \] Taking the limit as \(h \to \infty\) gives the long term stationary variance, i.e. \[ \lim_{h\to\infty} \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\} = \frac{\omega}{1- \lambda} = \sigma_{\infty}^2 \iff \lambda < 1 \] Hence, the general expression became \[ \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\} = \sigma_{\infty}^2 + \lambda^h (\sigma_t^2 - \sigma_{\infty}^2) \]

Example 22.2  

Unconditional expected value of GARCH(1,1) variance
# Unconditional first moment GARCH variance (exact)
e_sigma2_inf <- function(omega, alpha, beta, e_x2 = 1){
  # Derived quantities
  lambda <- alpha * e_x2 + beta
  # Expectation 
  omega / (1 - lambda)
}
Moment Formula MonteCarlo Difference (%)
\(\sigma_{\infty}^2\) 1.264 1.264 -0.0037%
Table 22.2: Forecasted long-term expectation of GARCH(1,1) variance with formula and by Monte Carlo simulations.

22.2 Second moment \(\sigma_t^2\)

22.2.1 Short term

The second moment admits an iterative formula, i.e.  \[ \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1} \} = \sum_{i=1}^{h} b_{t+h-i} \prod_{j=1}^{i-1} \gamma_{t+h-j} + \sigma_t^4 \prod_{j=1}^{h} \gamma_{t+h-j} \tag{22.6}\] where
\[ \gamma_{t+h-j} = \alpha_1^2 \mathbb{E}\{u_{t+h-j}^4\} + \beta_1 (2\alpha_1 \mathbb{E}\{u_{t+h-j}^2\} + \beta_1) \tag{22.7}\] while \[ b_{t+h-i} = \omega (\omega + 2\lambda_{t+h-i} \mathbb{E}\{\sigma_{t+h-i}^2 \mid \mathcal{F}_{t-1}\}) \tag{22.8}\] with \(\lambda_{t+h-i}\) as in Equation 22.2.

Proof. Starting from \[ \begin{aligned} \sigma_{t+1}^4 & {} = (\omega + \alpha_1 y_{t}^2 + \beta_1 \sigma_t^2)^2 \end{aligned} \] and substitute the definition of \(y_t^2 = \sigma_t^2 u_t^2\), i.e. \[ \begin{aligned} \sigma_{t+1}^4 & {} = (\omega + (\alpha_1 u_{t}^2 + \beta_1) \sigma_t^2)^2 = \\ & = \omega^2 + 2 \omega(\alpha_1 u_{t}^2 + \beta_1) \sigma_t^2 + (\alpha_1 u_{t}^2 + \beta_1)^2 \sigma_t^4 \end{aligned} \] Then, let’s take the conditional expectation on both sides: \[ \begin{aligned} \mathbb{E}\{\sigma_{t+1}^4 \mid \mathcal{F}_{t-1}\} & {} = \omega^2 + \mathbb{E}\{(\alpha_1 u_{t}^2 + \beta_1)^2\} \mathbb{E}\{\sigma_{t}^4 \mid \mathcal{F}_{t-1}\} + 2 \omega (\alpha_1 \mathbb{E}\{u_t^2\} + \beta_1) \mathbb{E}\{\sigma_{t}^2 \mid \mathcal{F}_{t-1}\} = \\ & = \omega^2 + \mathbb{E}\{(\alpha_1 u_{t}^2 + \beta_1)^2\} \sigma_{t}^4 + 2 \omega \lambda_t \sigma_{t}^2 \end{aligned} \] where \(\lambda_t = \alpha_1 \mathbb{E}\{u_t^2\} + \beta_1\). Then note that: \[ \begin{aligned} \mathbb{E}\{(\alpha_1 u_{t}^2 + \beta_1)^2 \} & {} = \mathbb{E}\{\alpha_1^2 u_{t}^4 + \beta_1^2 + 2 \alpha_1 \beta_1 u_{t}^2\} = \\ & = \alpha_1^2 \mathbb{E}\{u_{t}^4\} + \beta_1 (2\alpha_1 \mathbb{E}\{u_{t}^2\} + \beta_1) = \gamma_t \end{aligned} \] Hence, we can write the expectation in terms of the previous expectation. \[ \begin{aligned} \mathbb{E}\{\sigma_{t+1}^4 \mid \mathcal{F}_{t-1}\} & {} = \omega^2 + \gamma_t \sigma_{t}^4 + 2 \omega \lambda_t \sigma_{t}^2 = \\ & = \omega (\omega + 2\lambda_t \sigma_{t}^2) + \gamma_t \sigma_{t}^4 = \\ & = b_t + \gamma_t \sigma_{t}^4 \end{aligned} \] where \(b_t = \omega (\omega + 2\lambda_t \sigma_{t}^2)\). Iterating at time \(t+2\): \[ \begin{aligned} \mathbb{E}\{\sigma_{t+2}^4 \mid \mathcal{F}_{t-1}\} & {} = b_{t+1} + \gamma_{t+1} \mathbb{E}\{\sigma_{t+1}^4 \mid \mathcal{F}_{t-1}\} = \\ & = b_{t+1} + \gamma_{t+1} b_t + \gamma_{t+1} \gamma_t \sigma_{t}^4 \end{aligned} \] where in this case \[ b_{t+1} = \omega (\omega + 2\lambda_t \mathbb{E}\{\sigma^2_{t+1} \mid \mathcal{F}_{t-1}\}) \] At time \(t+3\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+3}^4 \mid \mathcal{F}_{t-1}\} & {} = b_{t+2} + \gamma_{t+2} \mathbb{E}\{\sigma_{t+2}^4 \mid \mathcal{F}_{t-1}\} = \\ & = b_{t+2} + \gamma_{t+2} b_{t+1} + \gamma_{t+2} \gamma_{t+1} b_t + \gamma_{t+2} \gamma_{t+1} \gamma_t \sigma_{t}^4 \end{aligned} \]

Example 22.3 With \(h = 1\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+1}^4 \mid \mathcal{F}_{t-1} \} {} & = \sum_{i=1}^{1} b_{t+1-i} \prod_{j=1}^{i-1} \gamma_{t+1-j} + \sigma_t^4 \prod_{j=1}^{1} \gamma_{t+1-j} = \\ & = b_{t} + \gamma_{t} \sigma_t^4 \end{aligned} \] With \(h = 2\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+2}^4 \mid \mathcal{F}_{t-1} \} {} & = \sum_{i=1}^{2} b_{t+2-i} \prod_{j=1}^{i-1} \gamma_{t+2-j} + \sigma_t^4 \prod_{j=1}^{2} \gamma_{t+2-j} = \\ & = b_{t+1} \prod_{j=1}^{0} \gamma_{t+2-j} + b_{t} \prod_{j=1}^{1} \gamma_{t+2-j} + \sigma_t^4\gamma_{t+1} \gamma_{t} = \\ & = b_{t+1} + b_{t} \gamma_{t+1} + \sigma_t^4\gamma_{t+1} \gamma_{t} \end{aligned} \] With \(h = 3\) \[ \begin{aligned} \mathbb{E}\{\sigma_{t+3}^4 \mid \mathcal{F}_{t-1} \} {} & = \sum_{i=1}^{3} b_{t+3-i} \prod_{j=1}^{i-1} \gamma_{t+3-j} + \sigma_t^4 \prod_{j=1}^{3} \gamma_{t+3-j} = \\ & = b_{t+2} \prod_{j=1}^{0} \gamma_{t+3-j} + b_{t+1} \prod_{j=1}^{1} \gamma_{t+3-j} + b_{t} \prod_{j=1}^{2} \gamma_{t+3-j} + \sigma_t^4\gamma_{t+2} \gamma_{t+1} \gamma_{t} = \\ & = b_{t+2} + b_{t+1} \gamma_{t+2} + b_{t} \gamma_{t+2} \gamma_{t+1} + \sigma_t^4\gamma_{t+2} \gamma_{t+1} \gamma_{t} \end{aligned} \] and so on.

Iterative formula for second moment of GARCH variance.
# GARCH second moment iteration
e_sigma4_h_iter <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_0){
  sigma4_iter <- numeric(h + 1)
  sigma4_iter[1] <- sigma4_0
  # Iterative formula for E[sigma^2]
  sigma2_iter <- numeric(h + 1)
  sigma2_iter[1] <- sqrt(sigma4_0)
  for(i in 1:h){
    lambda <- alpha * e_x2 + beta
    gamma <- alpha^2 * e_x4 + beta * (2 * alpha * e_x2 + beta)
    sigma2_iter[i+1] <- omega + lambda * sigma2_iter[i]
    b_t <- omega^2 + 2 * omega * lambda * sigma2_iter[i]
    sigma4_iter[i+1] <- b_t + gamma * sigma4_iter[i]
  }
  sigma4_iter
}
# GARCH second moment with formula
e_sigma4_h_mix <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_0){
  # Step 1: initialize vectors
  sigma2_iter <- numeric(h + 1)
  sigma4_mix  <- numeric(h + 1)
  sigma2_iter[1] <- sqrt(sigma4_0)        # E[sigma_t^2]
  sigma4_mix[1]  <- sigma4_0              # E[sigma_t^4]

  # Step 2: Compute lambda and gamma constants
  lambda <- alpha * e_x2 + beta
  gamma  <- alpha^2 * e_x4 + beta * (2 * alpha * e_x2 + beta)

  # Step 3: Forward recursion to get all E[sigma^2_{t+h}]
  for(i in 1:h){
    sigma2_iter[i + 1] <- omega + lambda * sigma2_iter[i]
  }
  # Step 4: Compute each E[sigma^4_{t+h}]
  for(k in 1:h){
    # Prepare sum for b_t terms
    b_sum <- 0
    for(i in 1:k){
      idx <- k - i + 1  # corresponds to time t + (k - i)
      b_t <- omega^2 + 2 * omega * lambda * sigma2_iter[idx]
      b_sum <- b_sum + b_t * gamma^(i - 1)
    }
    sigma4_mix[k + 1] <- b_sum + sigma4_0 * gamma^k
  }

  names(sigma4_mix) <- paste0("t+", 0:h)
  return(sigma4_mix)
}
Moment step Formula Iteration Difference
\(\mathbb{E}\{\sigma^4_{t+0} \mid \mathcal{F}_{t-1}\}\) 0 1.440000 1.440000 0%
\(\mathbb{E}\{\sigma^4_{t+1} \mid \mathcal{F}_{t-1}\}\) 1 1.559380 1.559380 0%
\(\mathbb{E}\{\sigma^4_{t+2} \mid \mathcal{F}_{t-1}\}\) 2 1.609123 1.609123 0%
\(\mathbb{E}\{\sigma^4_{t+3} \mid \mathcal{F}_{t-1}\}\) 3 1.630762 1.630762 0%
\(\mathbb{E}\{\sigma^4_{t+4} \mid \mathcal{F}_{t-1}\}\) 4 1.640413 1.640413 0%
\(\mathbb{E}\{\sigma^4_{t+5} \mid \mathcal{F}_{t-1}\}\) 5 1.644775 1.644775 0%
\(\mathbb{E}\{\sigma^4_{t+6} \mid \mathcal{F}_{t-1}\}\) 6 1.646762 1.646762 0%
\(\mathbb{E}\{\sigma^4_{t+7} \mid \mathcal{F}_{t-1}\}\) 7 1.647669 1.647669 0%
\(\mathbb{E}\{\sigma^4_{t+8} \mid \mathcal{F}_{t-1}\}\) 8 1.648085 1.648085 0%
\(\mathbb{E}\{\sigma^4_{t+9} \mid \mathcal{F}_{t-1}\}\) 9 1.648275 1.648275 0%
\(\mathbb{E}\{\sigma^4_{t+10} \mid \mathcal{F}_{t-1}\}\) 10 1.648363 1.648363 0%
Table 22.3: Forecasted second moment of GARCH(1,1) variance with iteration and with formula.

22.2.2 Long-term

If \(\mathbb{E}\{u_t^2\}\) and \(\mathbb{E}\{u_t^3\}\) are constant for all \(t\) then the formula simplifies to
\[ \begin{aligned} \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1} \} {} & = (\omega^2 + 2 \omega \lambda \sigma_{\infty}^2) \left(\frac{1-\gamma^h}{1-\gamma}\right) + 2 \omega \sigma_{\infty}^2 \lambda^{h} \left(\frac{\lambda(\lambda^h - (1 + \gamma)^h)}{\lambda^h (\lambda - 1 - \gamma)} \right) + \sigma_t^4 \gamma^h \end{aligned} \tag{22.9}\] where \(\sigma_{\infty}^2\) denotes the long-term expected GARCH variance as \(h\to\infty\), i.e.  \[ \sigma_{\infty}^4 = \lim_{h\to\infty} \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1} \} = \frac{\omega^2 (1 + \alpha_1 \mathbb{E}\{u_t^2\} + \beta_1)}{(1 - \alpha_1 \mathbb{E}\{u_t^2\} - \beta_1) (1-\alpha_1^2 \mathbb{E}\{u_{t}^4\} - 2\alpha_1 \beta_1 \mathbb{E}\{u_{t}^2\} - \beta_1^2)} \text{.} \tag{22.10}\] It follows that, under normality \(u_t \sim \mathcal{N}(0,1)\) we have that \(\mathbb{E}\{u_t^2\} = 1\) and \(\mathbb{E}\{u_t^4\} = 3\). Substituting, one obtain the same result as in Bollerslev (1986), i.e. \[ \sigma_{\infty}^4 = \frac{\omega^2 (1 + \alpha_1 + \beta_1)}{(1 - \alpha_1 - \beta_1) (1-3 \alpha_1^2 - 2 \alpha_1 \beta_1 - \beta_1^2)} \text{.} \]

Proof. Under the assumption of constant second and fourth moments of \(u_t\) one can simplify the expressions, i.e.
\[ \begin{aligned} {} & \gamma = \alpha_1^2 \mathbb{E}\{u_{t}^4\} + \beta_1 (2\alpha_1 \mathbb{E}\{u_{t}^2\} + \beta_1) \\ & b_{t+h-i} = \omega (\omega + 2 \lambda \mathbb{E}\{\sigma_{t+h-i}^2 \mid \mathcal{F}_{t-1}\}) \\ & \lambda = \alpha_1 \mathbb{E}\{u_{t}^2\} + \beta_1 \end{aligned} \] Recalling the expression of the expectation of the GARCH variance with constant moments in Equation 22.3 one can write \[ b_{t+h-i} = \omega^2 + 2 \omega \lambda \sigma_{\infty}^2 + 2 \omega \lambda^{h-i-1} \sigma_{\infty}^2 \] Substituting the above expression into Equation 22.6 one obtain \[ \begin{aligned} \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1} \} {} & = \sum_{i=1}^{h} \left( \omega^2 + 2 \omega \lambda \sigma_{\infty}^2 + 2 \omega \lambda^{h-i-1} \sigma_{\infty}^2 \right) \gamma^{i-1} + \sigma_t^4 \gamma^h \\ & = (\omega^2 + 2 \omega \lambda \sigma_{\infty}^2) \sum_{i=1}^{h} \gamma^{i-1} + 2 \omega \sigma_{\infty}^2 \lambda^{h} \sum_{i=1}^{h} \left(\frac{1}{\lambda}\right)^{i-1} \gamma^{i-1} + \sigma_t^4 \gamma^h = \\ & = (\omega^2 + 2 \omega \lambda \sigma_{\infty}^2) \sum_{i=0}^{h-1} \gamma^{i} + 2 \omega \sigma_{\infty}^2 \lambda^{h} \sum_{i=0}^{h-1} \left(\frac{1 + \gamma}{\lambda}\right)^{i} + \sigma_t^4 \gamma^h \end{aligned} \] Notably \[ \sum_{i=0}^{h-1} \gamma^{i} = \frac{1-\gamma^h}{1-\gamma} \] and \[ \sum_{i=0}^{h-1} \left(\frac{1 + \gamma}{\lambda}\right)^{i} = \frac{1 - \frac{(1 + \gamma)^h}{\lambda^h}}{1 - \frac{1 + \gamma}{\lambda}} = \frac{\lambda(\lambda^h - (1 + \gamma)^h)}{\lambda^h (\lambda - 1 - \gamma)} \] Hence, \[ \begin{aligned} \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1} \} {} & = (\omega^2 + 2 \omega \lambda \sigma_{\infty}^2) \left(\frac{1-\gamma^h}{1-\gamma}\right) + 2 \omega \sigma_{\infty}^2 \lambda^{h} \left(\frac{\lambda(\lambda^h - (1 + \gamma)^h)}{\lambda^h (\lambda - 1 - \gamma)} \right) + \sigma_t^4 \gamma^h \end{aligned} \] Taking the limit as \(h \to infty\), the second and third terms converges to zero if \(\lambda < 1\), therefore \[ \begin{aligned} \lim_{h\to\infty} \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1} \} {} & = \frac{\omega^2 + 2 \omega \lambda \sigma_{\infty}^2}{1-\gamma} \end{aligned} \] More explicitly, \[ \begin{aligned} \lim_{h\to\infty} \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1} \} {} & = \frac{\omega^2 (1 + 2 \frac{\lambda}{1-\lambda})}{1-\alpha_1^2 \mathbb{E}\{u_{t}^4\} - 2\alpha_1 \beta_1 \mathbb{E}\{u_{t}^2\} - \beta_1^2} = \\ & = \frac{\omega^2 (1 + \lambda)}{(1 - \lambda) (1-\alpha_1^2 \mathbb{E}\{u_{t}^4\} - 2\alpha_1 \beta_1 \mathbb{E}\{u_{t}^2\} - \beta_1^2)} = \\ & = \frac{\omega^2 (1 + \alpha_1 \mathbb{E}\{u_t^2\} + \beta_1)}{(1 - \alpha_1 \mathbb{E}\{u_t^2\} - \beta_1) (1-\alpha_1^2 \mathbb{E}\{u_{t}^4\} - 2\alpha_1 \beta_1 \mathbb{E}\{u_{t}^2\} - \beta_1^2)} \end{aligned} \]

Long-term formula for second moment of GARCH variance.
# Unconditional second moment GARCH variance (exact)
e_sigma4_inf <- function(omega, alpha, beta, e_x2 = 1, e_x4 = 3){
  # Expectation 
  e_sigma2 <- e_sigma2_inf(omega, alpha, beta, e_x2)
  # Derived quantities
  lambda <- alpha * e_x2 + beta
  gamma <- alpha^2 * e_x4 + 2 * alpha * beta + beta^2
  # Second moment  
  e_sigma2 * omega * (1 + lambda) / (1 - gamma)
}

# Unconditional second moment GARCH variance (exact)
e_sigma4_inf <- function(omega, alpha, beta, e_x2 = 1, e_x4 = 3){
  #1.352* omega^2 * (1 + alpha + beta) / ((1-alpha-beta) * (1 - beta^2 - 2*alpha*beta - e_x4 * alpha^2)) 
   #e_sigma2 * omega * (1 + alpha*e_x2 + beta) / ((1 - beta^2 - 2*alpha*beta*e_x2 - e_x4 * alpha^2)) 
  # Expectation 
  e_sigma2 <- e_sigma2_inf(omega, alpha, beta, e_x2)
  # Derived quantities
  lambda <- alpha * e_x2 + beta
  gamma <- alpha^2 * e_x4 + 2 * e_x2 * alpha * beta + beta^2
  # Second moment  
  e_sigma2 * omega * (1 + lambda) / (1 - gamma)
}
Moment Formula MonteCarlo Difference
\(\sigma_{\infty}^4\) 1.648437 1.648606 -0.0103%
Table 22.4: Forecasted long-term second moment of GARCH(1,1) variance with formula and by Monte Carlo simulations.

22.3 Variance \(\sigma_t^2\)

The conditional variance of \(\sigma_t^2\) with \(h \ge 1\) reads \[ \mathbb{V}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\} = \mathbb{E}\{\sigma_{t+h}^4 \mid \mathcal{F}_{t-1}\} - (\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\})^2 \text{.} \] Instead the long-term variance of \(\sigma_t^2\) reads explicitly \[ \mathbb{V}\{\sigma_t^2\} = \sigma_{\infty}^4 - (\sigma_{\infty}^2)^2 \text{.} \]

Example 22.4  

Iterative formula for variance of GARCH variance.
# Iterative GARCH variance formula
v_sigma2_h_mix <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_t){
  # Second moment GARCH variance
  e_sigma4 <- e_sigma4_h_mix(h, omega, alpha, beta, e_x2, e_x4, sigma4_t)
  # Expectation GARCH variance
  e_sigma2 <- e_sigma2_h_mix(h, omega, alpha, beta, e_x2, sqrt(sigma4_t))
  # Variance
  e_sigma4 - e_sigma2^2
}
# Iterative GARCH variance formula
v_sigma2_h_iter <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_t){
  # Second moment GARCH variance
  e_sigma4 <- e_sigma4_h_iter(h, omega, alpha, beta, e_x2, e_x4, sigma4_t)
  # Expectation GARCH variance
  e_sigma2 <- e_sigma2_h_iter(h, omega, alpha, beta, e_x2, sqrt(sigma4_t))
  # Variance
  e_sigma4 - e_sigma2^2
}
Long-term variance of GARCH variance.
# Unconditional variance GARCH std.dev (approximated)
v_sigma2_inf <- function(omega, alpha, beta, e_x2 = 1, e_x4 = 3){
  # Expectation GARCH variance 
  e_sigma2 <- e_sigma4_inf(omega, alpha, beta, e_x2, e_x4)
  # Expectation GARCH std. dev 
  e_sigma <- e_sigma2_inf(omega, alpha, beta, e_x2)
  # Variance 
  e_sigma2 - e_sigma^2
}
Moment step Formula Iteration Difference
\(\mathbb{V}\{\sigma_{t+1}^2\mid \mathcal{F}_{t-1}\}\) 1 0.0352420 0.0352420 0%
\(\mathbb{V}\{\sigma_{t+2}^2\mid \mathcal{F}_{t-1}\}\) 2 0.0455820 0.0455820 0%
\(\mathbb{V}\{\sigma_{t+3}^2\mid \mathcal{F}_{t-1}\}\) 3 0.0489759 0.0489759 0%
\(\mathbb{V}\{\sigma_{t+4}^2\mid \mathcal{F}_{t-1}\}\) 4 0.0502199 0.0502199 0%
\(\mathbb{V}\{\sigma_{t+5}^2\mid \mathcal{F}_{t-1}\}\) 5 0.0507180 0.0507180 0%
\(\mathbb{V}\{\sigma_{t+6}^2\mid \mathcal{F}_{t-1}\}\) 6 0.0509296 0.0509296 0%
\(\mathbb{V}\{\sigma_{t+7}^2\mid \mathcal{F}_{t-1}\}\) 7 0.0510227 0.0510227 0%
\(\mathbb{V}\{\sigma_{t+8}^2\mid \mathcal{F}_{t-1}\}\) 8 0.0510646 0.0510646 0%
\(\mathbb{V}\{\sigma_{t+9}^2\mid \mathcal{F}_{t-1}\}\) 9 0.0510835 0.0510835 0%
\(\mathbb{V}\{\sigma_{t+10}^2\mid \mathcal{F}_{t-1}\}\) 10 0.0510922 0.0510922 0%
Table 22.5: Forecasted variance of GARCH(1,1) variance with iteration and with formula.
Moment Formula MonteCarlo Difference
\(\mathbb{V}\{\sigma_{\infty}^2\}\) 0.0510995 0.0511493 -0.0975%
Table 22.6: Forecasted long-term variance of GARCH(1,1) variance with formula and by Monte Carlo simulations.

22.4 First moment \(\sigma_t\)

22.4.1 Short term

The expected value of the GARCH std. deviation can be approximated as \(\sigma_{t+h} = (\sigma_{t+h}^2)^{1/2}\) with a Taylor expansion \[ \mathbb{E}\{\sigma_{t+h} \mid \mathcal{F}_{t-1} \} \approx \sqrt{\mathbb{E}\{\sigma_{t+h}^{2} \mid \mathcal{F}_{t-1} \}} - \frac{1}{8} \frac{\mathbb{V}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1} \}}{\mathbb{E}\{\sigma_{t+h}^{2} \mid \mathcal{F}_{t-1} \}^{\frac{3}{2}}} \text{.} \]

Proof. Let’s \(X\) be a non-negative random variable, and let’s say one want to approximate: \[ \mathbb{E}\{\sqrt{X}\} \text{.} \] Let \(f(x) = \sqrt{x}\) and expand \(f(x)\) around the point \(\mu = \mathbb{E}\{X\}\), i.e. \[ \sqrt{X} \approx \sqrt{\mu} +\frac{1}{2} \mu^{-1/2} (X - \mu) - \frac{1}{8} \mu^{-3/2} (X - \mu)^2 \text{,} \] and take the expectation \[ \mathbb{E}\{\sqrt{X}\} \approx \sqrt{\mu} + \frac{1}{2} \mu^{-1/2} \left(\mathbb{E}\{X\} - \mu \right) - \frac{1}{8} \mu^{-3/2} \mathbb{E}\{(X - \mu)^2\} \text{,} \] where \(\mathbb{E}\{X\} - \mu = 0\) and \(\mathbb{E}\{(X - \mu)^2\} = \mathbb{V}\{X\}\). Applying this result to the random variable \(\sigma_{t+h}^2\) with \(1 < h\) and let’s approximate around the expected value with a Taylor expansion, i.e.  \[ \sigma_{t+h} = \sqrt{\sigma_{t+h}^2} \approx \sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}} + \frac{1}{2} \frac{(\sigma_{t+h}^2 -\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}) }{\sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}}} - \frac{1}{8} \frac{(\sigma_{t+h}^2 -\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\})^2 }{\sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}^3}} \text{,} \] Taking the expectation gives the result. \[ \mathbb{E}\{\sigma_{t+h} \mid \mathcal{F}_{t-1}\} \approx \sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}} - \frac{1}{8} \frac{\mathbb{V}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}}{\sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}^3}} \text{.} \]

Example 22.5  

Iterative formula for expectation of GARCH std. deviation.
# Conditional first moment GARCH std. dev (approximated)
e_sigma12_h_mix <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_t){
  # Expectation GARCH variance
  e_sigma2 <- e_sigma2_h_mix(h, omega, alpha, beta, e_x2, sqrt(sigma4_t))
  # Variance GARCH variance
  v_sigma2 <- v_sigma2_h_mix(h, omega, alpha, beta, e_x2, e_x4, sigma4_t)
  # Moment to power 1/2 (Approximated)
  e_sigma2^(1/2) - (1/8) * v_sigma2 / sqrt(e_sigma2)^3
}
# Conditional first moment GARCH std. dev (approximated)
e_sigma12_h_iter <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_t){
  # Expectation GARCH variance
  e_sigma2 <- e_sigma2_h_mix(h, omega, alpha, beta, e_x2, sqrt(sigma4_t))
  # Variance GARCH variance
  v_sigma2 <- v_sigma2_h_iter(h, omega, alpha, beta, e_x2, e_x4, sigma4_t)
  # Moment to power 1/2 (Approximated)
  e_sigma2^(1/2) - (1/8) * v_sigma2 / sqrt(e_sigma2)^3
}
Moment step Formula Iteration Difference
\(\mathbb{E}\{\sigma_{t+0}\mid \mathcal{F}_{t-1}\}\) 0 1.095445 1.095445 0%
\(\mathbb{E}\{\sigma_{t+1}\mid \mathcal{F}_{t-1}\}\) 1 1.107896 1.107896 0%
\(\mathbb{E}\{\sigma_{t+2}\mid \mathcal{F}_{t-1}\}\) 2 1.114145 1.114145 0%
\(\mathbb{E}\{\sigma_{t+3}\mid \mathcal{F}_{t-1}\}\) 3 1.117128 1.117128 0%
\(\mathbb{E}\{\sigma_{t+4}\mid \mathcal{F}_{t-1}\}\) 4 1.118522 1.118522 0%
\(\mathbb{E}\{\sigma_{t+5}\mid \mathcal{F}_{t-1}\}\) 5 1.119168 1.119168 0%
\(\mathbb{E}\{\sigma_{t+6}\mid \mathcal{F}_{t-1}\}\) 6 1.119466 1.119466 0%
\(\mathbb{E}\{\sigma_{t+7}\mid \mathcal{F}_{t-1}\}\) 7 1.119603 1.119603 0%
\(\mathbb{E}\{\sigma_{t+8}\mid \mathcal{F}_{t-1}\}\) 8 1.119666 1.119666 0%
\(\mathbb{E}\{\sigma_{t+9}\mid \mathcal{F}_{t-1}\}\) 9 1.119694 1.119694 0%
\(\mathbb{E}\{\sigma_{t+10}\mid \mathcal{F}_{t-1}\}\) 10 1.119708 1.119708 0%
Table 22.7: Forecasted expectation of GARCH(1,1) std. deviation with iteration and with formula.

22.4.2 Long term

The unconditional expected GARCH(1,1) std. deviation \[ \sigma_{\infty} \approx \sqrt{\sigma_{\infty}^2} + \frac{1}{8} \frac{\mathbb{V}\{\sigma_{\infty}^2\}}{(\sigma_{\infty}^2)^{\frac{3}{2}}} \text{.} \]

Example 22.6  

Long-term expectation of GARCH std. deviation.
# Unconditional first moment GARCH std. dev (approximated)
e_sigma12_inf <- function(omega, alpha, beta, e_x2 = 1, e_x4 = 3){
  # Expectation GARCH variance 
  e_sigma2 <- e_sigma2_inf(omega, alpha, beta, e_x2)
  # Variance GARCH variance 
  v_sigma2 <- v_sigma2_inf(omega, alpha, beta, e_x2, e_x4)
  # Moment to power 1/2 (Approximated) 
  e_sigma2^(1/2) - (1/8) * v_sigma2 / sqrt(e_sigma2)^3
}
Moment Formula MonteCarlo Difference
\(\sigma_{\infty}\) 1.119719 1.120488 -0.0687%
Table 22.8: Long-term expectation of GARCH(1,1) std. deviation with approximated formula and by Monte Carlo simulations.

22.5 Variance \(\sigma_t\)

22.5.1 Short term

The variance of the GARCH std. deviation can be approximated \[ \mathbb{V}^{\mathbb{P}}\{\sigma_{t+h} \mid \mathcal{F}_{t-1} \} \approx \frac{\mathbb{V}^{\mathbb{P}}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}}{4\mathbb{E}^{\mathbb{P}}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}} \text{,} \]

Example 22.7  

Iterative formula for variance of GARCH std. deviation.
# Iterative GARCH variance formula (approximated)
v_sigma_h_mix <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_t){
  # Expectation GARCH variance
  v_sigma2 <- v_sigma2_h_mix(h, omega, alpha, beta, e_x2, e_x4, sqrt(sigma4_t))
  # Expectation GARCH std. dev
  e_sigma2 <- e_sigma2_h_mix(h, omega, alpha, beta, e_x2, sqrt(sigma4_t))
  # Variance
  v_sigma2/(4*e_sigma2)
}

# Iterative GARCH variance formula (approximated)
v_sigma_h_iter <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_t){
  # Expectation GARCH variance
  v_sigma2 <- v_sigma2_h_iter(h, omega, alpha, beta, e_x2, e_x4, sqrt(sigma4_t))
  # Expectation GARCH std. dev
  e_sigma2 <- e_sigma2_h_iter(h, omega, alpha, beta, e_x2, sqrt(sigma4_t))
  # Variance
  v_sigma2/(4*e_sigma2)
}
Moment step Formula Iteration Difference
\(\mathbb{V}\{\sigma_{t+1}\mid \mathcal{F}_{t-1}\}\) 1 0.0059471 0.0059471 0%
\(\mathbb{V}\{\sigma_{t+2}\mid \mathcal{F}_{t-1}\}\) 2 0.0082691 0.0082691 0%
\(\mathbb{V}\{\sigma_{t+3}\mid \mathcal{F}_{t-1}\}\) 3 0.0092727 0.0092727 0%
\(\mathbb{V}\{\sigma_{t+4}\mid \mathcal{F}_{t-1}\}\) 4 0.0097250 0.0097250 0%
\(\mathbb{V}\{\sigma_{t+5}\mid \mathcal{F}_{t-1}\}\) 5 0.0099319 0.0099319 0%
\(\mathbb{V}\{\sigma_{t+6}\mid \mathcal{F}_{t-1}\}\) 6 0.0100270 0.0100270 0%
\(\mathbb{V}\{\sigma_{t+7}\mid \mathcal{F}_{t-1}\}\) 7 0.0100707 0.0100707 0%
\(\mathbb{V}\{\sigma_{t+8}\mid \mathcal{F}_{t-1}\}\) 8 0.0100908 0.0100908 0%
\(\mathbb{V}\{\sigma_{t+9}\mid \mathcal{F}_{t-1}\}\) 9 0.0101000 0.0101000 0%
\(\mathbb{V}\{\sigma_{t+10}\mid \mathcal{F}_{t-1}\}\) 10 0.0101042 0.0101042 0%
Table 22.9: Forecasted variance of GARCH(1,1) std. deviation with iteration and with formula.

22.5.2 Long term

The long-term variance of the GARCH std. deviation can be approximated as \[ \mathbb{V}\{\sigma_{\infty}\} \approx \frac{\mathbb{V}\{\sigma_{\infty}^2\}}{\sigma_{\infty}^2} \text{.} \]

Example 22.8  

Long-term variance of GARCH std. deviation.
# Unconditional variance GARCH std.dev (approximated)
v_sigma_inf <- function(omega, alpha, beta, e_x2 = 1, e_x4 = 3){
  # Expectation GARCH variance 
  v_sigma2 <- v_sigma2_inf(omega, alpha, beta, e_x2, e_x4)
  # Expectation GARCH std. dev 
  e_sigma2 <- e_sigma2_inf(omega, alpha, beta, e_x2)
  # Variance 
  v_sigma2 / (4 * e_sigma2^2)
}
Moment Formula MonteCarlo Difference
\(\mathbb{V}\{\sigma_{\infty}\}\) 0.0079976 0.0084112 -5.1714%
Table 22.10: Long-term variance of GARCH(1,1) std. deviation with approximated formula and by Monte Carlo simulations.

22.6 Third moment \(\sigma_t\)

22.6.1 Short term

The expected value of \(\sigma_{t+h}^3\) can be approximated with a Taylor expansion, i.e. \(\sigma_{t+h}^3 = (\sigma_{t+h}^2)^{3/2}\). Then, we can approximate
\[ \mathbb{E}\{\sigma_t^3 \mid \mathcal{F}_{t-1} \} \approx \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1} \}^{\frac{3}{2}} + \frac{3}{8} \frac{\mathbb{V}\{\sigma_{t+h} ^2 \mid \mathcal{F}_{t-1}\}}{\sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\}}} \text{.} \]

Example 22.9  

Iterative formula for third moment of GARCH std. deviation.
# Conditional third moment GARCH std. dev (approximated)
e_sigma32_h_mix <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_0){
  # Expectation variance
  e_sigma2 <- e_sigma2_h_mix(h, omega, alpha, beta, e_x2, sqrt(sigma4_0))
  # Variance variance
  v_sigma2 <- v_sigma2_h_mix(h, omega, alpha, beta, e_x2, e_x4, sigma4_0)
  # Moment to power 3/2 (Approximated)
  e_sigma2^(3/2) + (3/8) * v_sigma2 / sqrt(e_sigma2)
}

# Conditional third moment GARCH std. dev (approximated)
e_sigma32_h_iter <- function(h, omega, alpha, beta, e_x2 = 1, e_x4 = 3, sigma4_0){
  # Expectation variance
  e_sigma2 <- e_sigma2_h_iter(h, omega, alpha, beta, e_x2, sqrt(sigma4_0))
  # Variance variance
  v_sigma2 <- v_sigma2_h_iter(h, omega, alpha, beta, e_x2, e_x4, sigma4_0)
  # Moment to power 3/2 (Approximated)
  e_sigma2^(3/2) + (3/8) * v_sigma2 / sqrt(e_sigma2)
}
Moment step Formula Iteration Difference
\(\mathbb{E}\{\sigma_{t+0}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 0 1.314534 1.314534 0%
\(\mathbb{E}\{\sigma_{t+1}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 1 1.383623 1.383623 0%
\(\mathbb{E}\{\sigma_{t+2}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 2 1.413526 1.413526 0%
\(\mathbb{E}\{\sigma_{t+3}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 3 1.426837 1.426837 0%
\(\mathbb{E}\{\sigma_{t+4}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 4 1.432849 1.432849 0%
\(\mathbb{E}\{\sigma_{t+5}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 5 1.435585 1.435585 0%
\(\mathbb{E}\{\sigma_{t+6}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 6 1.436836 1.436836 0%
\(\mathbb{E}\{\sigma_{t+7}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 7 1.437408 1.437408 0%
\(\mathbb{E}\{\sigma_{t+8}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 8 1.437670 1.437670 0%
\(\mathbb{E}\{\sigma_{t+9}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 9 1.437791 1.437791 0%
\(\mathbb{E}\{\sigma_{t+10}^{\frac{3}{2}}\mid \mathcal{F}_{t-1}\}\) 10 1.437846 1.437846 0%
Table 22.11: Forecasted third moment of GARCH(1,1) std. deviation with iteration and with formula.

22.6.2 Long term

With a Taylor approximation, the long term third moment of the GARCH std. deviation reads \[ \sigma^3_{\infty} \approx (\sigma_{\infty}^2)^{\frac{3}{2}} + \frac{3}{8} \frac{\mathbb{V}\{\sigma_{\infty}^2\}}{\sqrt{\sigma_{\infty}^2}} \text{.} \]

Example 22.10  

Long-term third moment of GARCH std. deviation.
# Unconditional third moment GARCH std. dev (approximated)
e_sigma32_inf <- function(omega, alpha, beta, e_x2 = 1, e_x4 = 3){
  # Expectation GARCH variance 
  e_sigma2 <- e_sigma2_inf(omega, alpha, beta, e_x2)
  # Variance GARCH variance 
  v_sigma2 <- v_sigma2_inf(omega, alpha, beta, e_x2, e_x4)
  # Moment to power 3/2 (Approximated)
  e_sigma2^(3/2) + (3/8) * v_sigma2 / sqrt(e_sigma2)
}
Moment Formula MonteCarlo Difference
\(\sigma_{\infty}^{\frac{3}{2}}\) 1.437893 1.436856 0.0721%
Table 22.12: Long-term third moment of GARCH(1,1) std. deviation with approximated formula and by Monte Carlo simulations.

22.7 Covariance

Given the information at time \(t-1\), the covariance between two GARCH variances at time \(t+h\) and \(t+s\) reads \[ \mathbb{C}v\{\sigma_{t+s}^2 \cdot \sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\} = \left( \prod_{i=1}^{|h-s|} \lambda_{t+\max(s,h)-i} \right) \mathbb{V}\{\sigma_{\min(s,h)}^2 \mid \mathcal{F}_{t-1}\} \] It follows that the covariance between any variance at time \(t+h\) with \(h>1\) and the variance at time \(t\), given the infomation at time \(t-1\), is zero, i.e.  \[ \begin{aligned} \mathbb{C}v\{\sigma_t^2 \cdot \sigma_{t+h}^2 \mid \mathcal{F}_{t-1} \} {} & = \left( \prod_{i=1}^{h} \lambda_{t+h-i} \right) \mathbb{V}\{\sigma_t^2 \mid \mathcal{F}_{t-1}\} = 0 \end{aligned} \] for all \(h \ge 1\) since the variance of \(\sigma_t^2\), given \(\mathcal{F}_{t-1}\) is zero.

Proof. Without loss of generality, consider \(s < h\). Then, leveraging the tower property of the conditional expectation we can write: \[ \mathbb{E}\{\sigma_{t+s}^2 \cdot \sigma_{t+h}^2 \mid \mathcal{F}_{t-1} \} = \mathbb{E}\{ \sigma_{t+s}^2 \cdot \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t+s-1}\} \mid \mathcal{F}_{t-1}\} \] From our previous result we have: \[ \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t+s-1}\} = \omega\left( 1 + \sum_{j=1}^{h-s-1} \prod_{i=1}^{j} \lambda_{t+h-i}\right) + \sigma_{t+s}^2 \prod_{j=1}^{h-s} \lambda_{t+h-j} \text{,} \] Therefore, denoting as \[ \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_{t+s-1}\} = A_{h|s-1} + \sigma_{t+s}^2 B_{h|s} \text{,} \] We obtain that the covariance by definition is: \[ \mathbb{C}v\{\sigma_{t+s}^2, \sigma_{t+h}^2 \mid \mathcal{F}_{t-1}\} = \mathbb{V}\{\sigma_{t+s}^2 \mid \mathcal{F}_{t-1}\} \prod_{j=1}^{h-s} \lambda_{t+h-j} \text{,} \]

Example 22.11  

GARCH covariances bewtween variance
# Covariance 
cov_sigma2_h_mix <- function(h, s, omega, alpha, beta, e_x2, e_x4, sigma4_t){
  # Compute the variances
  v_sigma2_s <- v_sigma2_h_mix(h-1, omega, alpha, beta, e_x2, e_x4, sigma4_t)[-c(1:(s))]
  # Extract the moments from the indexes
  m2 <- e_x2[s:(h-1)]
  # Compute lambda
  lambda <- (alpha * m2 + beta)
  # Compute cumulative product
  prod_lambda <- cumprod(lambda)[(h-s):1]
  # Compute covariances
  cv_sigma2_hs <- prod_lambda * v_sigma2_s
  # Assign correct names
  names(cv_sigma2_hs) <- paste0(names(v_sigma2_s),"|",paste0("t+",h))
  cv_sigma2_hs
}

# Simulate multiple ACF 
sim_acf_garch <- function(B, h){
  acf_sigma2 <- 0
  acf_sigma <- 0
  for(i in 1:B){
    sim <- simulate_GARCH11(B = 1, n = 10000, sigma2_t, omega, alpha, beta, mu, sd, p)
    sim <- bind_rows(sim)
    acf_sigma2 <- acf_sigma2 + acf(sim$sigma2, plot = F, lag.max = h+1, type = "cov")$acf[,,1][-1]
    acf_sigma <- acf_sigma + acf(sqrt(sim$sigma2), plot = F, lag.max = h+1, type = "cov")$acf[,,1][-1]
  }
  structure(
    list(
      acf_sigma = acf_sigma/B,
      acf_sigma2 = acf_sigma2/B
    )
  )
}
Moment Formula MonteCarlo Difference
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+10}^2\mid \mathcal{F}_{t-1}\}\) 0.0234411 0.0235074 -0.2827%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+9}^2\mid \mathcal{F}_{t-1}\}\) 0.0107530 0.0107950 -0.3913%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+8}^2\mid \mathcal{F}_{t-1}\}\) 0.0049316 0.0049715 -0.808%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+7}^2\mid \mathcal{F}_{t-1}\}\) 0.0022608 0.0022927 -1.4142%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+6}^2\mid \mathcal{F}_{t-1}\}\) 0.0010353 0.0010490 -1.317%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+5}^2\mid \mathcal{F}_{t-1}\}\) 0.0004730 0.0005244 -10.8517%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+4}^2\mid \mathcal{F}_{t-1}\}\) 0.0002149 0.0002327 -8.2978%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+3}^2\mid \mathcal{F}_{t-1}\}\) 0.0000962 0.0000829 13.7854%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+2}^2\mid \mathcal{F}_{t-1}\}\) 0.0000411 0.0000072 82.3895%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+1}^2\mid \mathcal{F}_{t-1}\}\) 0.0000146 0.0000138 5.2527%
Table 22.13: Forecasted covariance between GARCH(1,1) variances with formula and by monte Carlo simulations

For the covariance between the GARCH std. deviations, we apply a Taylor approximation around the mean of \(\sigma_{t+s}\) and \(\sigma_{t+h}\), using the delta method: \[ \mathbb{C}v\{\sigma_{t+s}, \sigma_{t+h} \mid \mathcal{F}_t\} \approx \frac{\mathbb{C}v\{\sigma_{t+s}^2, \sigma_{t+h}^2 \mid \mathcal{F}_t\}}{4 \sqrt{\mathbb{E}\{\sigma_{t+s}^2 \mid \mathcal{F}_t\} \mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_t\}}} \text{.} \]

Proof. Considering the product of \[ \sigma_{t+h} = \sqrt{\sigma_{t+h}^2} \approx \sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_t\}} + \frac{1}{2} \frac{(\sigma_{t+h}^2 -\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_t\}) }{\sqrt{\mathbb{E}\{\sigma_{t+h}^2 \mid \mathcal{F}_t\}}} \text{,} \] and applying the covariance between \(\sigma_{t+h}\) and \(\sigma_{t+s}\) with \(1 < s < h\) one obtain the result.

Example 22.12  

Moment Formula MonteCarlo Difference
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+10}^2\mid \mathcal{F}_{t-1}\}\) 0.0046915 0.0039642 15.5031%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+9}^2\mid \mathcal{F}_{t-1}\}\) 0.0021521 0.0018440 14.3165%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+8}^2\mid \mathcal{F}_{t-1}\}\) 0.0009870 0.0008559 13.2839%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+7}^2\mid \mathcal{F}_{t-1}\}\) 0.0004525 0.0003956 12.5632%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+6}^2\mid \mathcal{F}_{t-1}\}\) 0.0002072 0.0001809 12.6954%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+5}^2\mid \mathcal{F}_{t-1}\}\) 0.0000947 0.0000890 6.0413%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+4}^2\mid \mathcal{F}_{t-1}\}\) 0.0000430 0.0000390 9.3954%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+3}^2\mid \mathcal{F}_{t-1}\}\) 0.0000193 0.0000137 29.0796%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+2}^2\mid \mathcal{F}_{t-1}\}\) 0.0000082 0.0000014 83.2984%
\(\mathbb{C}v\{\sigma_{t+11}^2, \sigma_{t+1}^2\mid \mathcal{F}_{t-1}\}\) 0.0000029 0.0000021 28.4146%
Table 22.14: Forecasted covariance between GARCH(1,1) std. deviation with formula and by monte Carlo simulations