Generalised impulse response function in R.

This function is a replacement for vars:::.irf that estimates Pesaran and Shin (1998) orthogonalised and generalised IRFs and is intended to be slotted into the vars:::.boot function to bootstrap confidence intervals for IRFs.

I referred to the Matlab code posted by Ken Nyholm.

# orthog = FALSE gives GIRF.

psGIRF <- function(x, n.ahead = 20, cumulative = TRUE, orthog = FALSE){

  y.names <- colnames(x$y)
  impulse <- y.names
  response <- y.names
  
  # Ensure n.ahead is an integer
  n.ahead <- abs(as.integer(n.ahead))
  
  # Create arrays to hold calculations     
  # [1:nlags, 1:nvariables, shocked variable ] 
  
  IRF_o  = array(data = 0, dim = c(n.ahead,x$K,x$K),
                 dimnames = list(NULL,y.names,y.names))       
  IRF_g  = array(data = 0, dim = c(n.ahead,x$K,x$K),
                 dimnames = list(NULL,y.names,y.names))
  IRF_g1 = array(data = 0, dim = c(n.ahead,x$K,x$K))
  
  # Estimation of orthogonalised and generalised IRFs
  SpecMA <- Phi(x, n.ahead)
  params <- ncol(x$datamat[, -c(1:x$K)])
  sigma.u <- crossprod(resid(x))/(x$obs - params)
  P <- t(chol(sigma.u))
  sig_jj <- diag(sigma.u)
  
  for (jj in 1:x$K){
    indx_      <- matrix(0,x$K,1)
    indx_[jj,1] <- 1
    
    for (kk in 1:n.ahead){  #kk counts the lag
      
      IRF_o[kk, ,jj] <- SpecMA[, ,kk]%*%P%*%indx_  # Peseran-Shin eqn 7 (OIRF)
      
      IRF_g1[kk, ,jj] <- SpecMA[, ,kk]%*%sigma.u%*%indx_
      IRF_g[kk, ,jj] <- sig_jj[jj]^(-0.5)*IRF_g1[kk, ,jj]  # Peseran-Shin eqn 10 (GIRF)
      
    }
  }
  
  if(orthog==TRUE){
    irf <- IRF_o
  } else if(orthog==FALSE) {
    irf <- IRF_g
  } else {
    stop("\nError! Orthogonalised or generalised IRF?\n")
  }
  
  idx <- length(impulse)
  irs <- list()
  for (ii in 1:idx) {
    irs[[ii]] <- matrix(irf[1:(n.ahead), response, impulse[ii]], nrow = n.ahead)
    colnames(irs[[ii]]) <- response
    if (cumulative) {
      if (length(response) > 1) 
        irs[[ii]] <- apply(irs[[ii]], 2, cumsum)
      if (length(response) == 1) {
        tmp <- matrix(cumsum(irs[[ii]]))
        colnames(tmp) <- response
        irs[[ii]] <- tmp
      }
    }
  }
  names(irs) <- impulse
  result <- irs
  return(result)
 
}