Deploying to gh-pages from @ 3ee5015 🚀

Author: gregorkastner

articles/Chapter06_files/figure-html/unnamed-chunk-7-1.png

@@ -206,10 +206,10 @@ Now we can compute their empirical quantiles, and those of the data.
 
 ``` r
 
-quants <- c(0.01, 0.05, 0.25, 0.4, 0.5, 0.6, 0.75, 0.95, 0.99)
-q_normal <- quantile(yf_normal, quants)
-q_t <- quantile(yf_t, quants)
-q_y <- quantile(y, quants)
+thesequants <- c(0.01, 0.05, 0.25, 0.4, 0.5, 0.6, 0.75, 0.95, 0.99)
+q_normal <- quantile(yf_normal, thesequants)
+q_t <- quantile(yf_t, thesequants)
+q_y <- quantile(y, thesequants)
 
 knitr::kable(round(t(cbind("Data" = q_y, "Student t" = q_t,
                            "Gaussian" = q_normal)), 3))
@@ -732,7 +732,7 @@ abline(h = 0, lty = 3)
 
 ![](Chapter09_files/figure-html/unnamed-chunk-34-1.png)
 
-### Section 2.5.3: Forecasting Volatility via GARCH(1,1) Models
+### Section 10.5.3: Forecasting Volatility via GARCH(1,1) Models
 
 We start by defining the parameter transformation.
 
@@ -757,7 +757,7 @@ loglik <- function(y, theta) {
   pars <- unrestricted2restricted(theta)
   sigma2 <- rep(NA_real_, length(y))
 
-  # initialize at unconditional variance
+  # Initialize at unconditional variance
   sigma2[1] <- pars[1] / (1 - pars[2] - pars[3])
 
   for (i in 2:length(y)) {
@@ -783,28 +783,32 @@ Finally, we define our random walk MH posterior sampler.
 ``` r
 
 sampler <- function(y, M = 20000L, nburn = 2000L, propsds = rep(0.15, 3)) {
-  # allocate storage
+  # Allocate storage
   thetas <- matrix(NA_real_, M, 3)
   sigma2s <- matrix(NA_real_, M, length(y))
 
-  # starting values
+  # Starting values
   theta <- rep(0, 3)
   currentlogpost <- loglik(y, theta) + logprior(theta)
 
   naccepts <- 0L
   for (m in seq_len(nburn + M)) {
-    if (m > nburn) sigma2s[m - nburn, ] <- attr(currentlogpost, "sigma2")
     prop <- rnorm(3, theta, propsds)
     proplogpost <- loglik(y, prop) + logprior(prop)
     logR <- proplogpost - currentlogpost
 
+    # Accept / reject
     if (log(runif(1)) < logR) {
       theta <- prop
       currentlogpost <- proplogpost
       if (m > nburn) naccepts <- naccepts + 1L
     }
 
-    if (m > nburn) thetas[m - nburn, ] <- theta
+    # Store
+    if (m > nburn) {
+      thetas[m - nburn, ] <- theta
+      sigma2s[m - nburn, ] <- attr(currentlogpost, "sigma2")
+    }
   }
 
   # Transform back to the usual parameterization
@@ -852,10 +856,10 @@ myhist <- function(dat, breaks = 30, freq = FALSE, ...) {
        main = "", ylab = "", ...)
 }
 
-layout(matrix(c(1, 1, 1, 2, 3, 4), nrow = 2, byrow = TRUE))
-ts.plot(y, col = "gray")
+ts.plot(y, col = "gray", ylab = "")
 lines(colMeans(sqrt(res$sigma2)))
-lines(-colMeans(sqrt(res$sigma2)))
+legend("topright", legend = c("Data", bquote(E(sigma[t] ~ "|" ~  bold(y)))),
+       lty = 1, col = c("gray", 1))
 myhist(res$para$alpha0, xlab = expression(alpha[0]))
 myhist(res$para$alpha1, xlab = expression(alpha[1]))
 myhist(res$para$gamma1, xlab = expression(gamma[1]))
@@ -868,25 +872,88 @@ a function that predicts a number of steps ahead.
 
 ``` r
 
-prediction <- function(obj, nahead = 1L) {
-  sigma2f <- yf <- matrix(NA_real_, obj$ndraws, nahead)
-  sigma2f[, 1] <- obj$para$alpha0 + obj$para$alpha1 * tail(obj$y, 1)^2 +
-    obj$para$gamma1 * obj$sigma2[, ncol(obj$sigma2)]
-  yf[, 1] <- rnorm(obj$ndraws, 0, sqrt(sigma2f[, 1]))
+prediction <- function(obj, nahead = 1L, drawsmult = 1L) {
+  # Allocate storage space
+  sigma2f <- yf <- matrix(NA_real_, obj$ndraws * drawsmult, nahead)
+  
+  # Predict one step ahead
+  sigma2f[, 1] <- rep(obj$para$alpha0 + obj$para$alpha1 * tail(obj$y, 1)^2 +
+    obj$para$gamma1 * obj$sigma2[, ncol(obj$sigma2)], drawsmult)
+  yf[, 1] <- rnorm(obj$ndraws * drawsmult, 0, sqrt(sigma2f[, 1]))
+  
+  # Predict two and more steps ahead
   if (nahead > 1L) {
     for (i in 2:nahead) {
       sigma2f[, i] <- obj$para$alpha0 + obj$para$alpha1 * yf[, i - 1]^2 +
         obj$para$gamma1 * sigma2f[, i - 1]
-      yf[, i] <- rnorm(obj$ndraws, 0, sqrt(sigma2f[, i]))
+      yf[, i] <- rnorm(obj$ndraws * drawsmult, 0, sqrt(sigma2f[, i]))
     }
   }
+  
+  # Return
   list(y = yf, sigma2 = sigma2f)
 }
 ```
 
-Let’s try it out.
+Let’s try it out. We cut the data a some point, re-estimate and predict.
+
+``` r
+
+cutoff <- 2990
+nahead <- 50
+ytrain <- head(y, cutoff)
+restrain <- sampler(ytrain)
+predstatic <- prediction(restrain, nahead, drawsmult = 10)
+predstaticquants <- apply(predstatic$y, 2, quantile, thesequants)
+```
+
+Alternatively, we might want to do sequential one-step-ahead
+predictions. Note that this amounts to re-estimating the model very
+often, thus we reduce the number of posterior draws.
 
 ``` r
 
-pred <- prediction(res, 2)
+predquants <- matrix(NA_real_, nrow = length(thesequants), nahead)
+predsigmamean <- rep(NA_real_, nahead)
+longrunsdmean <- rep(NA_real_, nahead)
+for (i in seq_len(nahead)) {
+  y1 <- head(y, cutoff - 1 + i)
+  res1 <- sampler(y1, M = 2000)
+  longrunsdmean[i] <- mean(sqrt(res1$para$alpha0 /
+                                   (1 - res1$para$alpha1 - res1$para$gamma1)))
+  pred1 <- prediction(res1, 1, drawsmult = 100)
+  predsigmamean[i] <- mean(sqrt(pred1$sigma2[, 1]))
+  predquants[, i] <- quantile(pred1$y[, 1], thesequants)
+}
 ```
+
+And visualize.
+
+``` r
+
+plotwindow <- (cutoff - nahead):(cutoff + nahead)
+plot(plotwindow, y[plotwindow], xlab = "", ylab = "", col = "gray", type = "l",
+     main = "Static predictions", ylim = range(predstaticquants, y[plotwindow]))
+abline(h = 0, lty = 3, col = "gray")
+abline(v = cutoff, lty = 3)
+lines((cutoff - nahead):cutoff, colMeans(sqrt(restrain$sigma2[, (cutoff - nahead):cutoff])))
+lines(cutoff + 1:nahead, colMeans(sqrt(predstatic$sigma2)))
+abline(h = mean(sqrt(restrain$para$alpha0 / (1 - restrain$para$alpha1 - restrain$para$gamma1))),
+       lty = 2)
+for (i in seq_along(thesequants)) {
+  lines(cutoff + 1:nahead, predstaticquants[i, ], col = 2)
+}
+
+plot(plotwindow, y[plotwindow], xlab = "", ylab = "", col = "gray", type = "l",
+     main = "Dynamic predictions", ylim = range(predstaticquants, y[plotwindow]))
+abline(h = 0, lty = 3, col = "gray")
+abline(v = cutoff, lty = 3)
+lines((cutoff - nahead):cutoff, colMeans(sqrt(restrain$sigma2[, (cutoff - nahead):cutoff])))
+lines(cutoff + 1:nahead, predsigmamean)
+lines(cutoff + 1:nahead, longrunsdmean, lty = 2)
+for (i in seq_along(thesequants)) {
+  lines(cutoff + 1:nahead, predquants[i, ], col = 2)
+}
+```
+
+![](Chapter09_files/figure-html/unnamed-chunk-43-1.png)

articles/Chapter07_files/figure-html/unnamed-chunk-27-1.png

@@ -13,7 +13,7 @@ articles:
   Chapter09: Chapter09.html
   Chapter10: Chapter10.html
   Chapter11: Chapter11.html
-last_built: 2026-05-22T02:54Z
+last_built: 2026-05-22T11:55Z
 urls:
   reference: https://gregorkastner.github.io/BayesianLearningCode/reference
   article: https://gregorkastner.github.io/BayesianLearningCode/articles