Constructing a counterfactual series on an SVAR in the vars package

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)
}

1 Like

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.