8c: Simulated GARCH

5 Ex: Simulated Data

See how closely true sigma is estimated and predicted, using simulated GARCH

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

  theta = c(.5, .7, .2)
  spec1 = garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]),  cond.dist="norm")

  # Generate GARCH
  set.seed(66343)
  X = garchSim(spec1, n=1000, extended=TRUE)   #extended=T:  now  Y has 3 columns: garch, sigma, eps
  X = as.xts(X)
  X = xts(X, order.by=as.Date(index(X)))

5a. Fit GARCH and forecast

  Y        = X[,1]    # simulated GARCH
  sig.true = X[,2]    # this is true sigma

  plot(Y, main="Simulated GARCH")

  plot(sig.true, main="True Sigma_t")

  # Estimate GARCH
  Fit01  =  garchFit(~ garch(1,1), data=Y, cond.dist="norm", include.mean = FALSE, trace = FALSE)
  sig.estim = xts(Fit01@sigma.t, order.by=index(Y))

  # compare true sigma vs estimated sigma
  head( cbind(sig.true, sig.estim) )

## Warning: object timezone ('UTC') is different from system timezone ('')
##   NOTE: set 'options(xts_check_TZ = FALSE)' to disable this warning
##     This note is displayed once per session

##                sigma sig.estim
## 2023-12-18 0.9844680 1.7773189
## 2023-12-19 0.8719402 1.0198346
## 2023-12-20 0.8528739 0.9322072
## 2023-12-21 0.9465185 1.0167695
## 2023-12-22 0.9430606 1.0033031
## 2023-12-23 1.1075998 1.1624516

  plot(cbind(sig.estim, sig.true), col=c("red", "black"),
                    main="True Sigma (Blk) vs Estimated (Red)")

  # Zoom in
  plot(cbind(sig.estim["2024"], sig.true["2024"]), col=c("red", "black"),
                    main="True Sigma (Blk) vs Estimated (Red)")

  plot(sig.estim/as.numeric(sig.true),
       main="Estimated Sigma / True Sigma", ylim=c(.8, 1.2) )

5b. Predict Sigma

  # Forecast sig.t for tomorrow (forecast sigma_t) by hand
  w = Fit01@fit$par[1]
  a = Fit01@fit$par[2]
  b = Fit01@fit$par[3]
  sig.pred = xts( sqrt( w + a*Y[1000]^2 + b*Fit01@sigma.t[1000]^2 ), order.by=index(Y[1000])+1)

  # Forecast sig.t (same as above)
  sig.pred2 = predict(Fit01, 1)
  sig.pred2

##   meanForecast meanError standardDeviation
## 1            0  1.061806          1.061806

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

  # Rolling 1-step prediction of sig.t
  sig.pred1 = numeric(0)
  date.stp = numeric(0)
  for (i in 1:750){
    T = Y[(1:250)+(i-1)]
    n = length(T)
    Fit01  =  garchFit(~ garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
    sig.pred1[i] = predict(Fit01)[1,3]   #- sig.t prediction
    date.stp[i]  = as.Date(index(Y[(250+(i-1)+1)]), format="%Y-%m-%d")
  }
  sig.pred = xts(sig.pred1, order.by=as.Date(date.stp))

  # Plot and compare estimated vs rolling predicted
  plot(cbind(sig.pred, sig.estim, sig.true),
    col=c("blue", "red", "black"),
    main="sigma.t:  True (Blk) vs Estim (Red) vs Pred(Blue)")

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

Under Construction

    # has problem with Rmd
    sigma.upper.t = xts( 1.96*sig.true,  order.by=index(Y))
    sigma.lower.t = xts(-1.96*sig.true,  order.by=index(Y))
    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, sigma.upper.t, sigma.lower.t),
        col=c("black", "black", "black"), lwd=c(2,1,1) )

    plot( cbind(Y, sigma.upper.t, sigma.lower.t,
                   sigma.upper.e, sigma.lower.e,
                   sigma.upper.p, sigma.lower.p),
        col=c("black", "black", "black", "red", "red", "blue", "blue"),
                                                      lwd=c(2,1,1,1,1) )

    plot( cbind(Y["2025"],
                sigma.upper.t["2025"], sigma.lower.t["2025"],
                sigma.upper.e["2025"], sigma.lower.e["2025"],
                sigma.upper.p["2025"], sigma.lower.p["2025"]),
        col=c("black", "black", "black", "red", "red", "blue", "blue"),
                                                      lwd=c(2,1,1,1,1) )

    plot(  as.numeric(            Y["2025::"]), type="h", lwd=2,
            xlab="diff( log(Y) )")
    lines( as.numeric(sigma.upper.e["2025::"]), col="red")
    lines( as.numeric(sigma.lower.e["2025::"]), col="red")
    lines( as.numeric(sigma.upper.p["2025::"]), col="blue")
    lines( as.numeric(sigma.lower.p["2025::"]), col="blue")