R : VAR(1) Simulation

This post explains how to perform simulations of a VAR(1) model.

VAR(1) Simulation

R code

The following R code loads the data, estimates a VAR(1) model using the vars R package, and generates simulations with confidence bands.

#========================================================# # Quantitative Financial Econometrics & Deep Learning, # R, Python, Excel by Sang-Heon Lee  # # https://shleeai.blogspot.com #--------------------------------------------------------# # VAR(1) Simulation #========================================================# graphics.off(); rm(list = ls())   library(vars) library(MASS)   #---------------------------------------- # data #---------------------------------------- data(usmacro, package = "bvarsv") M <- usmacro[,c(1,3)]   thor <- 36+1 # with last one observation, VAR(1) nlag <- 1; nvar <- 2; nsim <- 5000   nobs <- nrow(M) nold <- nobs # historical data in graph   #---------------------------------------- # VAR(1) estimation #---------------------------------------- est <- vars::VAR(M, p = nlag, type = "const",                        season = NULL) ARsC <- do.call(rbind, lapply(est$varresult,                                \(x) x$coefficients)) a0   <- ARsC[,ncol(ARsC)]     # constant term a1   <- ARsC[,-ncol(ARsC)]    # AR(1) coef. matrix sig  <- summary(est)$covres   # covariance matrix   #---------------------------------------- # VAR(1) forecast and fanchart #---------------------------------------- fcst <- predict(est, n.ahead = thor, ci = 0.99)   nWd = 5*2; nHt = 6 x11(width=nWd, height=nHt); par(mai=rep(0.4,4)) plot(fcst) x11(width=nWd, height=nHt); par(mai=rep(0.4,4)) fanchart(fcst, cis=c(0.01,0.1,0.3,0.5,0.7,0.9,0.99))   #---------------------------------------- # VAR(1) simulation #---------------------------------------- Ysim <- array(NA, dim = c(nsim, thor, nvar)) for (s in 1:nsim) {     y1 <- matrix(0, nrow=thor, ncol=nvar)     y1[1,] <- t(tail(M, nlag))     for (i in 2:thor) {         mvr <- mvrnorm(1, c(0,0), sig)         y1[i,] <- a0 + a1 %*% y1[i-1,] + mvr     }     Ysim[s, , ] <- y1 } #---------------------------------------- # Percentiles #---------------------------------------- perc = c(0.995, 0.50, 0.005) Ypct <- array(NA, dim = c(length(perc), thor, nvar)) for (h in 1:thor) { for (k in 1:nvar) {     Ypct[, h, k] <- quantile(Ysim[, h, k], probs = perc,                               names = FALSE, type = 7) }}   #---------------------------------------- # Graph #---------------------------------------- draw.VAR1.sim <- function(nth) {     xfcst <- (nold):(nold+thor-1)     plot(1:nold, M[(nobs-nold+1):nobs,nth],           main = paste0("X",nth), xlim = c(1,nold+thor),           ylim = c(-2, max(Ysim[,,nth],M[,nth])),          type="l", lty=1, lwd=2, col="black",           ylab = paste0("X",nth), xlab="Date")          for (i in 1:nsim) {         col1 <- rgb(0.53,0.81,0.92,alpha=0.1)      # sky blue         if(nth == 2) col1 <- rgb(1,0.5,0.5,alpha=0.1)  # pink         lines(xfcst, Ysim[i,,nth], type="l",                lty=1, lwd=0.5, col=col1)     }     lines(xfcst, Ypct[1,,nth], type="l", lty=1, lwd=1, col="blue")     lines(xfcst, Ypct[2,,nth], type="l", lty=1, lwd=2, col="blue")     lines(xfcst, Ypct[3,,nth], type="l", lty=1, lwd=1, col="blue") }   x11(width=nWd, height=nHt) par(mfrow = c(2,1), mar = c(3,4,2,2), mgp = c(2,0.6,0)) draw.VAR1.sim(1) draw.VAR1.sim(2)   Colored by Color Scripter

cs

The simulation results with confidence bands are shown below.

For comparison, I also include prediction results and fan charts.

Since the focus is on the simulation functionality, only a VAR(1) model is considered. For cases with \( p > 1 \), only the coefficient matrix needs to be modified.