R/irf.nlVar.R

In tsDyn: Nonlinear Time Series Models with Regime Switching

#' @importFrom vars irf
#' @export
vars::irf

#' Impulse response function
#' 
#' Compute the impulse response coefficients (IRF) of a VAR(p) (or transformed VECM to VAR(p)) for
#' \code{n.ahead} steps.
#' For TVECM and TVAR model, argument \code{regime} offers regime-specific IRF. 
#' 
#' @aliases irf.nlVar
#' @param x Object of class \sQuote{\code{VAR}}; generated by
#' \command{lineVar()}, or object of class \sQuote{\code{VECM}}; generated
#' by\command{VECM()}.
#' @param impulse,response Not used!
#' @param n.ahead Integer specifying the steps.
#' @param ortho Logical, if \code{TRUE} (the default) the orthogonalised
#' impulse response coefficients are computed .
#' @param cumulative Logical, if \code{TRUE} the cumulated impulse response
#' coefficients are computed. The default value is false.
#' @param boot,ci,runs,seed Arguments for the bootstrap, see  \code{\link[vars]{irf.varest}}
#' @param ... Currently not used.
#' @return A list of class \sQuote{\code{varirf}} see  \code{\link[vars]{irf.varest}}
#' @author Matthieu Stigler
#' @seealso \code{\link{plot}} for the plot method. \code{\link{lineVar}},
#' \code{\link{VECM}} for the models.
#' Hamilton, J. (1994), \emph{Time Series Analysis}, Princeton University
#' Press, Princeton.
#' 
#' Lutkepohl, H. (2006), \emph{New Introduction to Multiple Time Series
#' Analysis}, Springer, New York.
#' @keywords regression
#' @examples
#' 
#' data(barry)
#' 
#' ## For VAR
#' mod_var <- lineVar(barry, lag = 2)
#' irf(mod_var, impulse = "dolcan", response = c("dolcan", "cpiUSA", "cpiCAN"), boot =
#' FALSE)
#' 
#' ## For VECM
#' mod_VECM <- VECM(barry, lag = 2, estim="ML", r=2)
#' irf(mod_VECM, impulse = "dolcan", response = c("dolcan", "cpiUSA", "cpiCAN"), boot =
#' FALSE)
#' 

#'@rdname irf.nlVar
#'@export
irf.VAR <-  function(x, impulse = NULL, response = NULL, n.ahead = 10, 
                     ortho = TRUE, cumulative = FALSE, boot = TRUE, ci = 0.95, 
                     runs = 100, seed = NULL, ...) {
  irf_any(x=x, n.ahead = n.ahead, cumulative = cumulative, boot = boot, ci = ci, runs = runs,
          regime = "all", seed = seed, ortho = ortho, ...) 
}

#'@rdname irf.nlVar
#'@export
irf.VECM <-  function(x, impulse = NULL, response = NULL, n.ahead = 10, 
                      ortho = TRUE, cumulative = FALSE, boot = TRUE, ci = 0.95, 
                      runs = 100, seed = NULL, ...) {
  if(x$inputArgs$LRinclude !="none") stop("irf not available for VECM with LRinclude != 'none'")
  if(x$inputArgs$I !="level") stop("irf not available for VECM with I != 'level'")
  
  irf_any(x=x, n.ahead = n.ahead, cumulative = cumulative, boot = boot, ci = ci, runs = runs,
          regime = "all", seed = seed, ortho = ortho, ...) 
}

#'@rdname irf.nlVar
#'@export
irf.TVAR <-  function(x, impulse = NULL, response = NULL, n.ahead = 10, 
                      ortho = TRUE, cumulative = FALSE, boot = TRUE, ci = 0.95, 
                      runs = 100, seed = NULL, 
                      regime = c("L", "M", "H"), ...) {
  regime <- match.arg(regime)
  irf_any(x=x, n.ahead = n.ahead, cumulative = cumulative, boot = boot, ci = ci, runs = runs,
          regime = regime, seed = seed, ortho = ortho, ...) 
}

#'@rdname irf.nlVar
#'@export
irf.TVECM <-  function(x, impulse = NULL, response = NULL, n.ahead = 10, 
                       ortho = TRUE, cumulative = FALSE, boot = TRUE, ci = 0.95, 
                       runs = 100, seed = NULL, 
                       regime = c("L", "M", "H"), ...) {
  regime <- match.arg(regime)
  irf_any(x=x, n.ahead = n.ahead, cumulative = cumulative, boot = boot, ci = ci, runs = runs,
          regime = regime, seed = seed, ortho = ortho, ...) 
}

#' @importFrom vars irf
irf_old <- function(x, impulse=NULL, response=NULL, n.ahead=10, ortho=TRUE, cumulative=FALSE, boot=TRUE, ci=0.95, runs=100, seed=NULL, ...){
  model <- attr(x, "model")
  if(model=="VECM"){
    LRinc <- x$model.specific$LRinclude
    inc <- x$include
    if(LRinc=="both"|inc=="none"&LRinc=="none") stop("Sorry, irf() is not available for this specification of deterministic terms.")
  }
  irf(vec2var.tsDyn(x), impulse=impulse, response=response, n.ahead = n.ahead, ortho=ortho, cumulative=cumulative, boot=boot, ci=ci, runs=runs, seed=seed, ...)
}


#' @importFrom stats df.residual

irf_1.VAR <-  function(x, n.ahead=10, cumulative=FALSE, ortho =FALSE, regime = "all",  ...) {
  irf_1.nlVar(x, n.ahead = n.ahead,
              cumulative = cumulative,
              regime = regime,
              ortho = ortho, ...)

}

irf_1.VECM <-  function(x, n.ahead=10, cumulative=FALSE, ortho =FALSE, regime = "all",  ...) {
  irf_1.nlVar(x, n.ahead = n.ahead,
              cumulative = cumulative,
              regime = regime,
              ortho = ortho, ...)
  
}

irf_1.TVAR <-  function(x, n.ahead=10, cumulative=FALSE, ortho =FALSE,  
                        regime = c("L", "M", "H"), ...) {
  regime <-  match.arg(regime)
  irf_1.nlVar(x, n.ahead = n.ahead, 
              cumulative = cumulative, 
              regime = regime, 
              ortho = ortho, ...)
  
}

irf_1.TVECM <-  function(x, n.ahead=10, cumulative=FALSE, ortho =FALSE,  
                        regime = c("L", "M", "H"), ...) {
  regime <-  match.arg(regime)
  irf_1.nlVar(x, n.ahead = n.ahead, 
              cumulative = cumulative, 
              regime = regime, 
              ortho = ortho, ...)
  
}


irf_1.nlVar <-  function(x, n.ahead=10, cumulative=FALSE, regime = c("all", "L", "M", "H"), ortho =FALSE,  ...) {
  regime <-  match.arg(regime)
  series <- get_series(x)
  
  ##
  model <- class(x)[1]
  mod_type <- switch(model, 
                     "VAR" = "VAR", "TVAR" = "VAR",
                     "VECM" = "VECM", "TVECM" = "VECM",
                     "ERROR")
  lag <-  x$lag
  k <-  x$k
  thDelay <-  x$thDelay# unused so far!?!
  beta <- x$model.specific$beta

  coefs <- coef(x, regime = regime)
  
  ## 
  n_start <- ifelse(mod_type == "VAR", lag, lag+1)
  start <-  matrix(0, nrow= n_start, ncol = k)
  innovs <-  matrix(0, nrow= n.ahead +1, ncol = k)
  
  if(any(grepl("Const|Intercept|Trend", colnames(coefs)))) coefs <-  coefs[,-grep("Const|Intercept|Trend", colnames(coefs))]
  init_vals_M <- diag(k)
  
  ## orthogonoalise
  if(ortho) {
    sigma.u <- crossprod(residuals(x)) / df.residual(x)
    P <- t(chol(sigma.u))
    init_vals_M <- init_vals_M %*% P
  } 

  ## shock first
  gen_1_VAR <- function(k) {
    innov_k <- innovs
    innov_k[1, ] <-  init_vals_M[, k]
    out <- VAR.gen(B = coefs, n=n.ahead + 1, lag=lag, include = "none",  
                      starting = start,
                      innov = innov_k, 
                      returnStarting = FALSE, ...)
    colnames(out) <- series
    out
  }
  gen_1_VECM <- function(k) {
    innov_k <- innovs
    innov_k[1, ] <-  init_vals_M[, k]
    out <- VECM.gen(B = coefs, n=n.ahead + 1, lag=lag, include = "none",  
                    starting = start,
                    beta = beta,
                    innov = innov_k, 
                    returnStarting = FALSE, ...)
    colnames(out) <- series
    out
  }
  if(mod_type == "VAR") {
    gen_1 <- gen_1_VAR
  } else if (mod_type == "VECM") {
    gen_1 <- gen_1_VECM
  } else {
    stop("Problem: not a recognized model?")
  }
  res_M <- lapply(seq_len(k), gen_1 )
  names(res_M) <- series
  
  if(cumulative) stop("TODO") #res <-  apply(res, 2, cumsum)
  
  ### results
  res_df <- as.data.frame(do.call("rbind", res_M))
  res_df$impulse <-  rep(names(res_M), each = n.ahead +1)
  res_df
}

## small function for internal checks
irf_1_check <-  function(x) {
  res <-  list()
  res$irf <-  split(x[, -ncol(x)], x$impulse)
  to_mat <-  function(x){
    x2 <- as.matrix(x)
    rownames(x2) <- NULL
    x2
  }
  res$irf <- lapply(res$irf, to_mat)
  res$boot <- FALSE
  class(res) <- "varirf"
  res
}

#' @importFrom vars vec2var
irf_comp_tsD_vars <-  function(tsD, vars,r=1,  ortho = TRUE) {
  irf_tsD <- irf_any(tsD, ortho = ortho, boot = FALSE)
  if(inherits(vars, "ca.jo")) vars <-  vec2var(vars, r = r)
  irf_vars <- irf(vars , ortho = ortho, boot = FALSE)
  test <- all.equal(irf_tsD$irf, irf_vars$irf)
  if(!isTRUE(test)) {
    print(irf_tsD$irf[[1]])
    print( irf_vars$irf[[1]])
    print(irf_tsD$irf[[1]] - irf_vars$irf[[1]])
  }
  test
}


## #' @rdname irf.nlVar
### #' @export
# irf.VAR <-  function(x, impulse=NULL, response=NULL, n.ahead=10, ortho=TRUE, cumulative=FALSE, 
#                      boot=TRUE, ci=0.95, runs=100, seed=NULL, ...) {
#   irf_any(x=x, n.ahead = n.ahead, cumulative = cumulative, 
#           boot = boot, ci = ci, runs = runs, seed = seed, ...)
# }

if(FALSE){
  library(tsDyn)
  library(tidyverse)
  
  library(devtools)
  load_all()
  
  irf_1 <-  tsDyn:::irf_1
  irf_1.nlVar <-  tsDyn:::irf_1.nlVar
  irf_1_check <-  tsDyn:::irf_1_check
  irf_any <-  tsDyn:::irf_any
  
  ## vars
  library(vars)
  # data(Canada)
  irf.ca.jo <- function(x, r=1, ...) irf(vec2var(x, r = r), ...)
  
  var_vars_l1 <- VAR(Canada, p = 1, type = "const")
  var_vars_l4 <- VAR(Canada, p = 4, type = "const")
  irf_l1_vars <- irf(var_vars_l1, runs = 10, ortho = FALSE)
  irf_l1_vars_ortho <- irf(var_vars_l1, boot = FALSE, ortho = TRUE)
  irf_l4_vars <- irf(var_vars_l4, boot = FALSE, ortho = FALSE)
  irf_l4_vars_ortho <- irf(var_vars_l4, boot = FALSE, ortho = TRUE)
  lapply(irf_l1_vars$irf, head, 2)
  
  # tsDyn
  library(tsDyn)
  # x <- lineVar(Canada, lag=4)
  var_tsD_l1 <- lineVar(Canada, lag=1)
  var_tsD_l4 <- lineVar(Canada, lag=4)
  coef(var_tsD_l1)[1,]
  coef(var_vars_l1)$e[, 1]
  
  ## compare
  
  ## deg freedom
  all.equal(var_vars_l1$obs,   var_tsD_l1$t)
  df.residual(var_tsD_l1)
  
  ##
  irf_l1_tsD <- irf_1(var_tsD_l1) %>% irf_1_check
  irf_l4_tsD <- irf_1(var_tsD_l4) %>% irf_1_check
  irf_l4_tsD$irf
  irf_l1_tsD
  # plot(irf_l1_tsD)
  
  ## irf tsD ortho
  irf_l1_tsD_ortho <- irf_1(x=var_tsD_l1, ortho = TRUE) %>% irf_1_check
  irf_l4_tsD_ortho <- irf_1(var_tsD_l4, ortho = TRUE) %>% irf_1_check
  
  
  ## compare
  irf_comp_tsD_vars(tsD = var_tsD_l1, var_vars_l1, ortho = TRUE)
  irf_comp_tsD_vars(tsD = var_tsD_l4, var_vars_l4, ortho = TRUE)
  
  irf_comp_tsD_vars(tsD = var_tsD_l1, var_vars_l1, ortho = FALSE)
  irf_comp_tsD_vars(tsD = var_tsD_l4, var_vars_l4, ortho = FALSE)
  
  
  ##
  a <- irf_1(x=var_tsD_l1, ortho = TRUE)
  
  #### Now irf_any
  irf_any_VAR <- irf_any(x=var_tsD_l1, boot = TRUE, ortho = FALSE)
  irf_any_VAR <- irf_any(x=var_tsD_l1, boot = TRUE, ortho = FALSE)
  irf_any_VAR
  
  ##
  r1 <- irf_l1_vars$irf[[1]]
  all.equal(irf_l1_vars$irf, irf_any_VAR$irf)
  irf_any_VAR$irf[[1]]
  irf_any_VAR$Lower[[1]]
  all.equal(irf_l1_vars$Lower[[1]], irf_any_VAR$Lower[[1]])
  
  
  ## TVAR
  tv_l1_th1 <- TVAR(zeroyld, lag =1, trace = FALSE)
  tv_l2_th1 <- TVAR(zeroyld, lag =2, trace = FALSE)
  tv_l1_th2 <- TVAR(zeroyld, lag =1, trace = FALSE, nthresh = 2)
  
  irf_1(x=tv_l1_th1, show.parMat = FALSE, n.ahead = 2, regime = "L")
  irf_1(x=tv_l1_th1, show.parMat = FALSE, n.ahead = 2, regime = "H")
  irf_1(tv_l2_th1, n.ahead = 2, regime = "L")
  irf_1(tv_l1_th2, n.ahead = 2, regime = "L")
  
  #
  irf_any(tv_l1_th1, show.parMat = FALSE, runs = 2, regime = "L")$irf
  irf_any(tv_l1_th1, show.parMat = FALSE, runs = 2, regime = "H")$irf
  
  ## VECM
  vec_l1_r1 <- VECM(barry, estim = "ML", lag = 1)
  vec_l1_r2 <- VECM(barry, estim = "ML", lag = 1, r = 2)
  
  vecm_l1_co_var <- ca.jo(barry, K=2, ecdet="none", spec="transitory")
  vecm_l3_co_var <- ca.jo(barry, K=4, ecdet="none", spec="transitory")
  
  irf_1(x=vec_l1_r1, show.parMat = FALSE, n.ahead = 2)
  irf_1(x=vec_l1_r2, show.parMat = FALSE, n.ahead = 2)
  
  
  irf_comp_tsD_vars(tsD=vec_l1_r1, vars=vecm_l1_co_var, ortho = FALSE)
  irf_comp_tsD_vars(tsD=vec_l1_r2, vars=vecm_l1_co_var, ortho = FALSE, r = 2)
  
  irf_comp_tsD_vars(tsD=vec_l1_r1, vars=vecm_l1_co_var, ortho = TRUE)
  irf_comp_tsD_vars(tsD=vec_l1_r2, vars=vecm_l1_co_var, ortho = TRUE, r = 2)
  
  ## TVECM
  tvec_l1_r1 <- TVECM(barry[, 1:2], lag = 1, trace = FALSE, plot = FALSE)
  irf_1(x=tvec_l1_r1, show.parMat = FALSE, n.ahead = 2, regime = "L")
  irf_1(x=tvec_l1_r1, show.parMat = FALSE, n.ahead = 2, regime = "H")
}