I have the following Structural VAR Model from the vars
package.
library(vars)
data(Canada)
var <- VAR(Canada, p = 2, type = "const")
svar <- BQ(var)
I'm interested in constructing counterfactual time series for the variables within my model, where I can evaluate hypothetical paths for each variable if shocks on other variables didn't occur. The svars
package has quite a nice implementation of this (bottom of this page), however I'm unsure of how to implement something similar within the vars
package.
Here is the source code to build a counterfactual series on an svars
, is there any way to do something similar on my model? My first guess is to change A_hat
to A
in the source code, unsure of what other parts of the code to amend in order to do this.
function (x, series = 1, transition = 0)
{
if (x$type == "const") {
A_hat <- x$A_hat[, -1]
}
else if (x$type == "trend") {
A_hat <- x$A_hat[, -1]
}
else if (x$type == "both") {
A_hat <- x$A_hat[, -c(1, 2)]
}
else {
A_hat <- x$A_hat
}
B_hat <- x$B
horizon <- x$n
IR <- array(unlist(IRF(A_hat, B_hat, horizon)), c(x$K, x$K,
horizon))
impulse <- matrix(0, ncol = dim(IR)[2]^2 + 1, nrow = dim(IR)[3])
colnames(impulse) <- rep("V1", ncol(impulse))
cc <- 1
impulse[, 1] <- seq(1, dim(IR)[3])
for (i in 1:dim(IR)[2]) {
for (j in 1:dim(IR)[2]) {
cc <- cc + 1
impulse[, cc] <- IR[i, j, ]
colnames(impulse)[cc] <- paste("epsilon[",
colnames(x$y)[j], "]", "%->%", colnames(x$y)[i])
}
}
y <- x$y
p <- x$p
obs <- x$n
k <- x$K
B <- x$B
A <- x$A_hat
Z <- t(YLagCr(y, p))
if (x$type == "const") {
Z <- rbind(rep(1, ncol(Z)), Z)
}
else if (x$type == "trend") {
Z <- rbind(seq(1, ncol(Z)), Z)
}
else if (x$type == "both") {
Z <- rbind(rep(1, ncol(Z)), seq(1, ncol(Z)), Z)
}
else {
Z <- Z
}
u <- t(y[-c(1:p), ]) - A %*% Z
s.errors <- solve(B_hat) %*% u
impulse <- impulse[, -1]
y_hat <- matrix(NA, nrow = obs, ncol = k)
if (transition == 0) {
for (i in 1:obs) {
for (j in 1:k) {
y_hat[i, j] <- impulse[1:i, j + series * k -
k] %*% t(s.errors)[i:1, j]
}
}
}
else {
for (i in (obs - floor(obs * (1 - transition))):obs) {
for (j in 1:k) {
y_hat[i, j] <- impulse[(obs - floor(obs * (1 -
transition))):i, j + series * k - k] %*% t(s.errors)[i:(obs -
floor(obs * (1 - transition))), j]
}
}
}
y_hat_a <- rowSums(y_hat)
yhat <- as.data.frame(cbind(seq(1, length(y_hat_a)), (y[-c(1:p),
series] - mean(y[-c(1:p), series])), y_hat_a, y_hat))
yhat_counter <- as.data.frame(matrix(NA, nrow = obs, ncol = k))
for (i in 1:k) {
yhat_counter[, i] <- (yhat[, 2] - yhat[, (3 + i)])
}
colnames(yhat)[3] <- paste("Constructed series ", colnames(y)[series])
colnames(yhat)[2] <- paste("Demeaned series ", colnames(y)[series])
for (i in 4:ncol(yhat)) {
colnames(yhat)[i] <- paste("Cumulative effect of ",
colnames(y)[i - 3], "shock on ", colnames(y)[series])
}
for (i in 1:ncol(yhat_counter)) {
colnames(yhat_counter)[i] <- paste(colnames(y)[series],
"with and without cumulative effect of ", colnames(y)[i],
"shock")
}
yhat_counter <- as.matrix(yhat_counter[, -series])
yhat <- yhat[, c(1, rep(2, ncol(yhat_counter)))]
if (inherits(x$y, "ts")) {
count <- ts(yhat[, -grep("V1", colnames(yhat))],
start = start(lag(x$y, k = -x$p)), end = end(x$y),
frequency = frequency(x$y))
counterfac <- list(actual = na.omit(count), counter = yhat_counter)
}
else {
counterfac <- list(actual = as.data.frame(na.omit(yhat)),
counter = yhat_counter)
}
class(counterfac) <- "cf"
return(counterfac)
}