6. Ex: Fit GARCH to SPY

See how close estimated sigma and rolling predicted sigma, using real data.

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

  getSymbols("SPY")     # SP500 download from Yahoo!
## [1] "SPY"
  SPY = Ad(SPY)["2022::"]
  is.ts(SPY)      # not ts object
  is.xts(SPY)     # its xts object

  plot( log(SPY) )
  plot( diff( log(SPY) ) )


6a. Fit GARCH and forecast

  library(fGarch)
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")
## 
## Attaching package: 'fGarch'
## The following object is masked from 'package:TTR':
## 
##     volatility
  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: 0x13808d200>
##  [data = Y]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1  
## 3.3724e-06  1.0003e-01  8.6956e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  3.372e-06   9.396e-07    3.589 0.000332 ***
## alpha1 1.000e-01   1.723e-02    5.804 6.48e-09 ***
## beta1  8.696e-01   2.129e-02   40.848  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  3782.322    normalized:  3.21626 
## 
## Description:
##  Sun Sep 13 16:21:47 2026 by user:
    sig.estim = xts(Fit01@sigma.t, order.by=index(Y))

    ##  Fit01@fit$par                  # estimated parameters
    ##  Fit01@residuals                # this is not GARCH residuals! This is same as Y.
    ##  Fit01@sigma.t                  # estimated sig_t
    ##  Fit01@residuals/Fit1@sigma.t   # this is the (standardized) GARCH residuals

    Fit01@fit$ics                      # AIC and BIC are here
##       AIC       BIC       SIC      HQIC 
## -6.427418 -6.414485 -6.427431 -6.422542
    res1  = Y/as.numeric(sig.estim)    # this is the GARCH residuals


6b. 10-step predicion of Sigma

  # Forecast sig.t
  sig.pred2 = predict(Fit01, 10)
  sig.pred2
##    meanForecast   meanError standardDeviation
## 1             0 0.007427117       0.007427117
## 2             0 0.007540356       0.007540356
## 3             0 0.007648550       0.007648550
## 4             0 0.007752013       0.007752013
## 5             0 0.007851027       0.007851027
## 6             0 0.007945852       0.007945852
## 7             0 0.008036725       0.008036725
## 8             0 0.008123864       0.008123864
## 9             0 0.008207469       0.008207469
## 10            0 0.008287727       0.008287727


6c. Rolling 1-step prediction of \(\sigma_t\)

  options(warn= -1)
  length(Y)
## [1] 1176
  # Rolling 1-step prediction of sig.t
  Y = Y                     # Original data
  window.size = 250         # Window size for estimation

    sig.pred1 = numeric(0)
    date.stp  = numeric(0)
    for (i in 1:(length(Y)-window.size)){
      T = Y[(1:window.size)+(i-1)]
      n = length(T)
      out1  =  garchFit(~ garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
      sig.pred1[i] = predict(out1)[1,3]   #- sig.t prediction
      date.stp[i]  = index(T[n])+1
    }
    sig.pred = xts(sig.pred1, order.by=as.Date(date.stp))


  # Plot and compare estimated vs rolling predicted
  plot(cbind(sig.estim['2022'],
             sig.pred[ '2022']),
     col=c("red", "blue"), lwd=c(1,1,1),
     main="sigma.t:  Estim (Red) vs Pred(Blue)")

  plot(cbind(sig.estim['2023'],
             sig.pred[ '2023']),
     col=c("red", "blue"), lwd=c(1,1,1),
     main="sigma.t:  Estim (Red) vs Pred(Blue)")

    sigma.upper.e = xts( 1.96*sig.estim, order.by=index(sig.estim))
    sigma.lower.e = xts(-1.96*sig.estim, order.by=index(sig.estim))
    sigma.upper.p = xts( 1.96*sig.pred, order.by=index(sig.pred))
    sigma.lower.p = xts(-1.96*sig.pred, order.by=index(sig.pred))

    plot( cbind(Y[            "2017::"],
                sigma.upper.e["2017::"], sigma.lower.e["2017::"],
                sigma.upper.p["2017::"], sigma.lower.p["2017::"]),
                main="Y with 95perc daily PI",
                col=c("black", "red", "red", "blue", "blue"), lwd=c(2,1,1,1,1) )

    time.period="2017::"
    plot( cbind(Y[            time.period],
                sigma.upper.e[time.period], sigma.lower.e[time.period],
                sigma.upper.p[time.period], sigma.lower.p[time.period]),
                main="Y with 95perc daily PI",
                col=c("black", "red", "red", "blue", "blue"), lwd=c(2,1,1,1,1) )

    time.period="2019::"
    plot( cbind(Y[            time.period],
                sigma.upper.e[time.period], sigma.lower.e[time.period],
                sigma.upper.p[time.period], sigma.lower.p[time.period]),
                main="Y with 95perc daily PI",
                col=c("black", "red", "red", "blue", "blue"), lwd=c(2,1,1,1,1) )

    time.period="2019::"
    plot(  as.numeric(            Y[time.period]), type="h", lwd=2, ylim=c(-.12, .12),
            xlab="diff( log(Y) )")
    lines( as.numeric(sigma.upper.e[time.period]), col="red")
    lines( as.numeric(sigma.lower.e[time.period]), col="red")
    lines( as.numeric(sigma.upper.p[time.period]), col="blue")
    lines( as.numeric(sigma.lower.p[time.period]), col="blue")


6d. More Examples in Books

Cryer : CREF stock values

Shumway : NYSE us GNP

Cowpertwait : SP500, Southern hem Temp.