R: VAR(1) Simulation

The article “R: VAR(1) Simulation” was originally published on SHLee AI Financial Model blog.

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

VAR(1) Simulation

#========================================================#
# 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)

Visit SHLee AI Financial Model blog for additional insights on this topic.