R/irf.R

In jetroant/sbetel: Estimation of moment condition models with smoothed BETEL

irf <- function(model, output, 
                horizon = 40, N = 10000,
                cumulate = c(), secondary = NULL) {
  
  m <- ncol(model$y)
  p <- model$args$p
  
  A_indices <- 1:(m*(1+m*p))
  B_indices <- (m*(1+m*p)+1):(m*(1+m*p)+m^2) 
  
  ret <- list()
  rows <- sample.int(nrow(output$sample), N, replace = TRUE)
  A_post <- output$sample[rows, A_indices]
  B_post <- output$sample[rows, B_indices]
  
  for(shock_index in 1:m) {
    
    e <- rep(0, m)
    e[shock_index] <- 1
    irfs <- array(NA, dim = c(m, horizon+1, nrow(A_post)))
    
    print(paste0("Computing impulse responses... (", shock_index, "/", m,")"))
    pb <- txtProgressBar(min = 0, max = nrow(A_post), style = 3)
    for(row_index in 1:nrow(A_post)) {
      
      B <- matrix(B_post[row_index,], ncol = m)
      A <- matrix(A_post[row_index,], ncol = m)
      AA <- sbetel:::stackA(A)
      
      for(h in 1:(horizon+1)) {
        if(h == 1) {
          zero <- B %*% e
          zero_long <- c(zero, rep(0, (ncol(AA)-length(zero))))
          irfs[,h,row_index] <- zero
        } else {
          irfs[,h,row_index] <- (expm::`%^%`(AA, (h-1)) %*% zero_long)[1:length(zero)]
        }
      }
      
      setTxtProgressBar(pb, row_index)
    }
    close(pb)
    
    if(length(cumulate) > 0) {
      for(i in 1:length(cumulate)) {
        for(j in 1:N) {
          irfs[cumulate[i],,j] <- cumsum(irfs[cumulate[i],,j])
        }
      }
    }
    
    ret[[shock_index]] <- irfs
  }
  
  if(!is.null(secondary)) {
    cat("Secondary impulse responses: \n")
    secondary <- sbetel:::irf(model = model, 
                              output = secondary,
                              horizon = horizon,
                              N = N,
                              cumulate = cumulate)
    ret$secondary <- secondary
  }
  
  ret
}

#Stacks VAR(p) coefficient matrix to VAR(1) coefficient matrix
stackA <- function(A) {
  A <- t(A)[,-1]
  m <- nrow(A)
  lags <- ncol(A)/m
  eye <- diag(m*lags-m)
  A <- rbind(A, cbind(eye, matrix(0, ncol = m, nrow= nrow(eye))))
  A
}

#Plots impulse response functions
irf_plot <- function(irf_obj, varnames, probs = NULL) {
  
  m <- nrow(irf_obj[[1]]) 
  par(mar = c(2,4,2,1))
  par(mfrow = c(m, m))
  indexmat <- matrix(1:m^2, ncol = m)
  row <- 0
  col <- 0
  for(fig_index in 1:m^2) {
    
    if((fig_index-1) %% m == 0) {
      row <- row + 1
      col <- 1
    } else {
      col <- col + 1
    }
    
    sub_irfs <- t(irf_obj[[row]][col,,])
    mean_sub_irfs <- ts(apply(sub_irfs, 2, mean), start = 0)
    
    if(is.null(probs)) {
      p <- c(0.0249, 0.025, seq(0.1, 0.9, 0.1), 0.975, 0.9751)
    } else {
      p <- probs
    }
    quant <- function(column, probs = p) quantile(column, probs = probs)
    quantiles_sub_irfs <- apply(sub_irfs, 2, quant)
    
    color <- "tomato"
    ylims <- c(min(quantiles_sub_irfs), max(quantiles_sub_irfs))
    
    if(!is.null(irf_obj$secondary)) {
      sub_irfs_second <- t(irf_obj$secondary[[row]][col,,])
      mean_sub_irfs_second <- ts(apply(sub_irfs_second, 2, mean), start = 0)
      if(is.null(probs)) {
        p_second <- c(0.025, 0.975)
      } else {
        p_second <- c(probs[2], probs[(length(probs)-1)])
      }
      quantiles_sub_irfs_second <- apply(sub_irfs_second, 2, quant, probs = p_second)
      ylims <- c(min(c(quantiles_sub_irfs, quantiles_sub_irfs_second)), 
                 max(c(quantiles_sub_irfs, quantiles_sub_irfs_second)))
    }
      
    plot(mean_sub_irfs, lwd = 2, lty = 2, col = color, ylab = "", xlab = "", 
         main = paste0("Shock ", row, " on ", varnames[col]), 
         ylim = ylims)
    grid()
    fanplot::fan(data = quantiles_sub_irfs, data.type = "values", probs = p,
                 start = 0, fan.col = colorRampPalette(c(color, "white")),
                 rlab = NULL, ln = NULL)
    abline(h = 0, lwd = 2, lty = 2)
    
    if(!is.null(irf_obj$secondary)) {
      lines(ts(quantiles_sub_irfs_second[1,]), lwd = 1, lty = 2)
      lines(ts(quantiles_sub_irfs_second[2,]), lwd = 1, lty = 2)
    }
    
    post_mass <- (max(p[-length(p)]) - min(p[-1]))*100
    if(col == 1 & row == 1) legend("topright", c(paste0(post_mass,"% of post. prob. mass")), lwd = 0, bty = "n", col = "tomato")
  }
  par(mfrow = c(1,1))
}