R/cf.R

In svars: Data-Driven Identification of SVAR Models

#' Counterfactuals for SVAR Models
#'
#' Calculation of Counterfactuals for an identified SVAR object 'svars' derived by function id.st( ), id.cvm( ),id.cv( ),id.dc( ) or id.ngml( ).
#'
#' @param x SVAR object of class "svars"
#' @param series Integer. indicating the series for which the counterfactuals should be calculated.
#' @param transition Numeric. Value from [0, 1] indicating how many initial values should be discarded, i.e., 0.1 means that the first 10 per cent observations of the sample are considered as transient.
#'
#' @return A list with class attribute "hd" holding the Counterfactuals as data frame.
#'
#' @references Kilian, L., Luetkepohl, H., 2017. Structural Vector Autoregressive Analysis, Cambridge University Press.
#'
#' @seealso \code{\link{id.cvm}}, \code{\link{id.dc}}, \code{\link{id.ngml}}, \code{\link{id.cv}}, \code{\link{id.garch}} or \code{\link{id.st}}
#'
#' @examples
#' \donttest{
#' v1 <- vars::VAR(USA, lag.max = 10, ic = "AIC" )
#' x1 <- id.dc(v1)
#' x2 <- cf(x1, series = 2)
#' plot(x2)
#' }
#'
#' @export

cf <- function(x, series = 1, transition = 0){

  ## Step 1: Calculate MA coefficients

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

  # Step 2: Calculate structural errors

  # gathering informations from vars object
  y <- x$y
  p <- x$p
  obs <- x$n
  k <- x$K
  B <- x$B

  # calculating covariance from actual VAR
  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

  # Step 3: Match up structural shocks with appropriate impuslse response
  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))

  # Calculating counterfactuals
  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)
}