1. GARCH model

GARCH(1,1) model \[\begin{align*} Y_t &= \sigma_t \epsilon_t \hspace{10mm} \epsilon_t \sim_{iid} N(0,1) \\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align*}\]

Conditional Mean and Var of GARCH

Conditional Mean: \(0\)

Conditional Variance : \(\sigma_t^2\)


Mean and Var of GARCH

(unconditional) Mean: \(0\)

(unconditional) Variance : \(\sigma_t^2\)


1a. Ex: Daily SPY

Live donload Daily Price of SP500 ETF (SPY) from Yahoo!

  library(quantmod)
  source('https://nmimoto.github.io/R/TS-00.txt')

  getSymbols("SPY")    #- Download SP500 now from Yahoo!
## [1] "SPY"
  SPY = Ad(SPY)['2010::']

  is.ts(SPY)      # not ts object
## [1] FALSE
  is.xts(SPY)     # its xts object
## [1] TRUE
  plot( SPY )

  plot( log(SPY) )

  plot( diff( log(SPY) ) )


1b. Fit GARCH(1,1)

  library(fGarch)   # for garchFit()

  Y = diff( log(SPY) )[-1]     # remove the first diff for NA

    Fit01 =  garchFit(~ garch(1,1), data=Y, cond.dist="norm", include.mean = FALSE, trace = FALSE)
    Fit01
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = Y, cond.dist = "norm", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x11aa29b80>
##  [data = Y]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1  
## 3.7798e-06  1.5098e-01  8.1387e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  3.780e-06   4.755e-07    7.949 1.78e-15 ***
## alpha1 1.510e-01   1.290e-02   11.703  < 2e-16 ***
## beta1  8.139e-01   1.396e-02   58.291  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  13917.88    normalized:  3.31615 
## 
## Description:
##  Sun Sep 13 15:21:21 2026 by user:
    ##  Fit01@fit$par                  # estimated parameters
    ##  Fit01@residuals                # this is not GARCH residuals! This is same as Y.
    ##  Fit01@sigma.t                  # estimated sig_t
    ##  Fit01@fit$ics                   # AIC and BIC are here
    ##  Fit01@residuals/Fit1@sigma.t   # this is the (standardized) GARCH residuals

    Fit01@fit$ics                   # AIC and BIC are here
##       AIC       BIC       SIC      HQIC 
## -6.630870 -6.626336 -6.630871 -6.629267
    res1  = Y/Fit01@sigma.t        #- this is GARCH residuals

1c. plotting resutls

    plot(Y) # log difference

    plot(Fit01@sigma.t, type="l")   # estimated Sigma_t

    sigma.upper = xts(1.96*Fit01@sigma.t,  order.by=index(Y))
    sigma.lower = xts(-1.96*Fit01@sigma.t, order.by=index(Y))

    # Plot the two together
    plot( cbind(Y, sigma.upper, sigma.lower),
        col=c("black", "red", "red"), lwd=c(2,1,1) )

    # Zoom in
    plot(  as.numeric(          Y["2018::2021"]), type="h", lwd=2)
    lines( as.numeric(sigma.upper["2018::2021"]), col="red")
    lines( as.numeric(sigma.lower["2018::2021"]), col="red")

# Test standardized log difference for randomness
Randomness.tests(res1)

##   B-L test H0: the series is uncorrelated
##   M-L test H0: the square of the series is uncorrelated
##   J-B test H0: the series came from Normal distribution
##   SD         : Standard Deviation of the series
##       BL15  BL20 BL25  ML15  ML20 JB    SD
## [1,] 0.365 0.221 0.03 0.369 0.663  0 0.998


2. Parameter Estimation

2a. Quasi-MLE

Quasi-Maximum Likelihood Estimator, means we will use MLE for normal errors, even in the case of non-normality. Theoretically, for GARCH parameters, it is shown that Normal MLE is still consistent under non-normal errors.


GARCH(1,1)

\[\begin{align*} Y_t &= \sigma_t \epsilon_t \hspace{10mm} \epsilon_t \sim_{iid} N(0,1)\\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align*}\]

Even thorugh we can use quasi-MLE based on Normal distribution, using correct conditional distribution is importatnt when estimating parameters of GARCH.


2b. Importatnce of Conditional Distribution

\[\begin{align*} Y_t &= \sigma_t \epsilon_t \hspace{10mm} \epsilon_t \sim_{iid} N(0,1)\\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align*}\]

[Conditonal distribution] \(=\) [Distribution of \(\epsilon_t\)]

Guessing correct distribution for \(\epsilon_t\) is very important in GARCH parameter estimation.

2c. Four Built-in distributions

  x = seq(-10,10, .1)

  #- Standardized t-distribution
  plot(x, dstd(x, mean = 0, sd = 1, nu = 5), type="l", ylim=c(0,.5))
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )

  #- Generalized Error Distribution
  plot(x, dged(x, mean = 0, sd = 1, nu = 4), type="l", ylim=c(0,.5))
  lines(x, dnorm(x, mean = 0, sd=1) , col="red" )

  #- Skewed-Standardized t-distribution
  plot( x, dsstd(x, mean = 0, sd = 1, nu = 5, xi = 1.5), type="l", ylim=c(0,.5))
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )

  #- Skewed-Generalized Error Distribution
  plot(x, dsged(x, mean = 0, sd = 1, nu = 5, xi = 1.5) , type="l", ylim=c(0,.5))
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )


2d. Simulation Demo of Parameter Estimation Efficiency

True parameter (.024, .1, .8)

True conditional distribution is Normal. Estimated using Normal.

True conditional distribution is std(5), but estimated using Normal.

True conditional distribution is sged(skew=.7, shape=1.45) but estimated using Normal.

  # Simulation showing how error distribution affect parameter estimation in GARCH(1,1) ---
  # (Takes about 3min)
  n   = 1000
  itt = 1000
  theta = c(.024, .1, .8)

  set.seed(62192)  #36621)
  Param11 = matrix(rep(0,3*itt), itt, 3)
  for (i in 1:itt) {

    spec = garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]),
                                                               cond.dist="norm")
    x    = garchSim(spec, n = n, extended=FALSE)

    est2 =  garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    Param11[i,] = coef(est2)
  }

  Param12 = matrix(rep(0,3*itt), itt, 3)
  for (i in 1:itt) {

    spec = garchSpec(model = list(shape=5, omega=theta[1], alpha=theta[2], beta=theta[3]),
                                                                cond.dist="std")
    x     = garchSim(spec, n = n, extended=FALSE)

    est2 =  garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    Param12[i,] = coef(est2)
  }

  Param13 = matrix(rep(0,3*itt), itt, 3)
  for (i in 1:itt) {

    spec = garchSpec(model = list(skew=.9,  shape=1.45, omega=theta[1], alpha=theta[2], beta=theta[3]),
                                                                cond.dist="sged")
    x    = garchSim(spec, n = n, extended=FALSE)

    est2 =  garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    Param13[i,] = coef(est2)
  }
## Warning in sqrt(diag(fit$cvar)): NaNs produced
  layout(matrix(1:9, 3, 3, byrow=T))
  hist(Param11[,1], 10, xlim=c(0,.5)); hist(Param11[,2], 10, xlim=c(0,.5));
    hist(Param11[,3], xlim=c(0.4,1), breaks=10)
  hist(Param12[,1], 10, xlim=c(0,.5)); hist(Param12[,2], 10, xlim=c(0,.5));
    hist(Param12[,3], xlim=c(0.4,1), breaks=10)
  hist(Param13[,1], 10, xlim=c(0,.5)); hist(Param13[,2], 10, xlim=c(0,.5));
    hist(Param13[,3], xlim=c(0.4,1), breaks=10)

  P11 = (Param11 - t(matrix(theta, 3, itt) ))^2
  P12 = (Param12 - t(matrix(theta, 3, itt) ))^2
  P13 = (Param13 - t(matrix(theta, 3, itt) ))^2

  # Raw MSE
  MSE1 = apply(P11, 2, mean)
  MSE2 = apply(P12, 2, mean)
  MSE3 = apply(P13, 2, mean)

  # Relative MSE
  rbind(MSE1 / MSE1,
        MSE2 / MSE1,
        MSE3 / MSE1)
##          [,1]     [,2]     [,3]
## [1,] 1.000000 1.000000 1.000000
## [2,] 1.770278 2.513438 1.858293
## [3,] 1.490498 1.324999 1.431167


3. Stationarity Cond’n for GARCH(1,1)

GARCH(1,1) is weakly stationary if \[ E( \log(\beta + \alpha \epsilon_t^2)) < 0. \]

This is satisfied if \[ \alpha + \beta < 1. \]

4. Predicting \(\sigma_t\)

4a. Residuals of GARCH

GARCH(1,1) \[\begin{align*} Y_t &= \sigma_t \epsilon_t \hspace{10mm} \epsilon_t \sim_{iid} N(0,1)\\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align*}\]

Using observation \(\{Y_1, \ldots, Y_n\}\), the residuals are \[ \hat \epsilon_t = Y_t/ \hat \sigma_t \]

how can we get \(\hat \sigma_t\)?

4b. Write \(\sigma^2_t\) using \(Y_t^2\)

Starting with the definition of GARCH(1,1),

\[\begin{align*} \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \\\\ &= \omega + \alpha Y_{t-1}^2 + \beta \Big(\omega + \alpha Y_{t-2}^2 + \beta \sigma^2_{t-2} \Big)\\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \sigma^2_{t-2} \\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \Big(\omega + \alpha Y_{t-3}^2 + \beta \sigma^2_{t-3} \Big) \\ \\ &= \omega + \beta \omega + \beta^2 \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \alpha Y_{t-3}^2 + \beta^3 \sigma^2_{t-3} \end{align*}\]

Continuing, we get \[\begin{align*} &= \omega (1 + \beta + \beta^2 + \cdots ) + \alpha \sum_{i=0}^k \beta^i Y_{t-1-i}^2 + \beta^{k+1} \sigma^2_{t-1-k} \\\\ &= \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \end{align*}\]

That means if we keep going, we can write \[ \sigma_t^2 \hspace{3mm} = \hspace{3mm} \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \]

We will use the truncated and estimated version of this, \[ \hat \sigma_t^2 \hspace{3mm} = \hspace{3mm} \frac{\hat \omega}{1- \hat \beta} \hspace{3mm} + \hspace{3mm} \hat \alpha \sum_{i=0}^{t-1} \hat \beta^i Y_{t-1-i}^2 \]

4c. Now We Can Calculate Residuals

GARCH(1,1) says, \[ Y_t \hspace{3mm} = \hspace{3mm} \sigma_t \epsilon_t \]

Using observation \(\{Y_1, \ldots, Y_n\}\), the residuals are \[ \hat \epsilon_t = Y_t/ \hat \sigma_t \]


4d. Forecast \(\sigma_{t+1}\)

GARCH(1,1) \[\begin{align*} Y_t &= \sigma_t \epsilon_t \hspace{10mm} \epsilon_t \sim_{iid} N(0,1)\\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align*}\]

Now have 1-step forecast of tomorrow’s volatility \(\sigma_{t+1}\), \[\begin{align*} \sigma_{t}^2(1) &= \hat \omega + \hat \alpha Y_{t}^2 + \hat \beta \hat \sigma^2_{t} \end{align*}\]

How well can GARCH predict \(\sigma_t\) of the future?