test_Supp.R
In pvars: VAR Modeling for Heterogeneous Panels

####################
###  UNIT TESTS  ###
####################
#
# For exported "svarfevd", "svarirf",  
# "sboot" functions, and supporting modules.


test_that("irf.*() and fevd.*() provide the same results for 'svars', 'L.varx', and 'pvarx'", {
  # prepare data #
  data("PCAP")
  names_k = c("g", "k", "l", "y")  # variable names
  names_i = levels(PCAP$id_i)      # country names 
  L.data  = sapply(names_i[1], FUN=function(i) 
    ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), simplify=FALSE)
  
  # calculate single SVAR #
  ### Note that svars (version 1.3.11) does not count the number of deterministic regressors 
  ### ... when calculating the OLS covariance matrix of residuals. Only with a 
  ### ... single deterministic regressor, the results are identical to vars and pvars!
  L.vars  = lapply(L.data, FUN=function(x) vars::VAR(x, p=2, type="const"))
  R.svars = svars::id.chol(L.vars[[1]], order_k=names_k)
  R.pid   = pid.chol(L.vars, order_k=names_k)
  colnames(R.svars$B) = colnames(R.pid$L.varx[[1]]$B)  ### svars::id.*() do not define names_s for shocks.
  L.svars = list(AUS=R.svars)
  
  # calculate FEVD #
  n.ahead = 3
  R.fevd_svars = svars:::fevd.svars(R.svars, n.ahead=n.ahead)
  R.fevd_varx  = fevd.id(R.svars, n.ahead=n.ahead)
  
  # calculate IRF #
  R.irf_svars = svars:::irf.svars(R.svars, n.ahead=n.ahead)
  R.irf_varx  = irf.varx(R.svars,  n.ahead=n.ahead)
  R.irf_pid   = irf.pvarx(R.pid,   n.ahead=n.ahead)
  R.irf_pid_L = irf.pvarx(L.svars, n.ahead=n.ahead)
  R.irf_pid_i = irf.pvarx(R.pid,   n.ahead=n.ahead, idx_i=1)
  R.irf_id_i  = irf.varx(R.pid$L.varx[[1]], n.ahead=n.ahead)
  
  # check #
  ### Note that irf.svars (version 1.3.8) produces IRF over the interval [0, n.ahead-1]!
  idx_row = 1:n.ahead
  expect_equal(R.svars$B, R.pid$L.varx[[1]]$B)
  expect_equivalent(R.fevd_svars, R.fevd_varx)
  expect_equal(R.irf_svars$irf[idx_row, ], R.irf_pid$irf[idx_row, ])
  expect_equal(R.irf_varx, R.irf_pid)
  expect_equal(R.irf_pid,  R.irf_pid_L)
  expect_equal(R.irf_varx, R.irf_pid_i)
  expect_equal(R.irf_id_i, R.irf_pid_i)
})


test_that("PP.*() provide the same persistence profiles for objects of class 'varx' and 'vec2var'", {
  tolerant = 25e-08  # vec2var() is based on equation-wise lm().
  
  # prepare data #
  data("PCAP")
  names_k = c("g", "k", "l", "y")  # variable names
  names_i = levels(PCAP$id_i)      # country names 
  L.data  = sapply(names_i, FUN=function(i) ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), simplify=FALSE)
  
  # joint specifications #
  dim_r = 2
  dim_p = 2
  shock = cbind(c(1, 0, 0, 0), c(0, 0, 1, 0), c(0, 0, 1, 1))
  n.ahead = 10
  n.season = 4
  
  # calculate single reduced-rank VAR #
  R.vecm = VECM(y=L.data$DNK, dim_p=dim_p, dim_r=dim_r, type="Case4", t_D2=list(n.season=n.season))
  R.cajo = urca::ca.jo(L.data$DNK, ecdet="trend", type="eigen", K=dim_p, spec="transitory", season=n.season)
  R.vars = vars::vec2var(R.cajo, r=dim_r)
  
  # calculate persistence profiles #
  R.ppv_varx = PP.variable(R.vecm, n.ahead = n.ahead, shock = shock)
  R.pps_varx = PP.system(R.vecm,   n.ahead = n.ahead)
  R.ppv_vars = PP.variable(R.vars, n.ahead = n.ahead, shock = shock)
  R.pps_vars = PP.system(R.vars,   n.ahead = n.ahead)
  
  # check #
  expect_equivalent(R.vecm$resid, t(R.vars$resid))
  expect_equal(R.ppv_varx, R.ppv_vars, tolerance=tolerant)
  expect_equal(R.pps_varx, R.pps_vars, tolerance=tolerant)
})


test_that("aux_check() warns or stops in case of incorrect arguments", {
  # prepare data #
  data("ERPT")
  names_k = c("lpm5", "lfp5", "llcusd")  # variable names for "Chemicals and related products"
  names_i = levels(ERPT$id_i)[c(1,6,2,5,4,3,7)]  # ordered country names
  dim_N   = length(names_i)
  
  # arguments #
  L.data  = sapply(names_i, FUN=function(i) ts(ERPT[ERPT$id_i==i, names_k], start=c(1995, 1), frequency=12), simplify=FALSE)
  L.dataF = L.data; L.dataF[[2]] = L.data[[2]][ , -2, drop=FALSE]
  R.lags  = c(3, 3, 3, 4, 3, 3, 3); names(R.lags)=names_i  # lags of VAR model by MAIC
  R.lagsF = R.lags; names(R.lagsF) = names_i[dim_N:1]  # lags with false names
  t_DF    = list(t_impulse=20, t_error=10)
  DF_sq   = matrix((1:100)^2, nrow=1)
  
  # test the stops and warnings: pvarx #
  expect_error(pvarx.VAR(L.dataF, lags=R.lags, type="both"),
               "The number of variables in 'L.data' must be the same for all individuals.")
  expect_error(pvarx.VAR(DF_sq, lags=R.lags, type="both"),
               "Argument 'L.data' must be a list of N individual data.frame objects.")
  expect_error(pvarx.VAR(L.data, lags=R.lags, type="both", D=list(A=DF_sq, B=DF_sq)),
               "Argument 'D' must be either a single data.frame object or a list of N individual data.frame objects.")
  
  # test the stops and warnings: pcoint # 
  expect_error(pcoint.BR(L.dataF, lags=R.lags, type="Case4"),
               "The number of variables in 'L.data' must be the same for all individuals.")
  expect_error(pcoint.JO(L.data, lags=1:2, type="Case4"),
               "Argument 'lags' must be either a single integer or a vector of N integers specific to each individual.")
  expect_warning(pcoint.SL(L.data, lags=R.lagsF, type="SL_trend"),
                 "Arguments 'L.data' and lags have mismatching names for individuals.")
  expect_warning(pcoint.SL(L.data, lags=2, type="SL_trend", t_D=list(t_impulse=10), n.factors=2),
                 "The factor estimation by PCA ignores period-specific deterministic regressors.")
  expect_warning(pcoint.CAIN(L.data, lags=2, type="SL_trend", t_D=t_DF),
                 "Unrecognized elements in 't_D' have been removed: t_error")
})


test_that("speci.VAR() can reproduce basic examples from vars and urca package.", {
  # prepare data #
  library("urca")
  library("vars")
  data("denmark")
  sjd = denmark[, c("LRM", "LRY", "IBO", "IDE")]
  
  # use a single lag-order to determine only the break period #
  R.vars  = vars::VAR(sjd, type="both", p=1, season=4)
  R.speci = speci.VAR(R.vars, lag_set=2, dim_m=1, trim=3, add_dummy=FALSE)
  
  # perform cointegration test procedure #
  R.t_D   = list(t_shift=8, n.season=4)
  R.coint = coint.SL(sjd, dim_p=2, type_SL="SL_trend", t_D=R.t_D)
  R.urca  = urca::cajolst(sjd, trend=TRUE, K=2, season=4)
  
  # use m=0 to determine only the lag-order #
  data(Canada)
  R.vars_ic = vars::VARselect(Canada, lag.max=5, type="const")
  R.spec_ic = speci.VAR(vars::VAR(Canada, p=1, type="const"), lag_set=1:5, dim_m=0)
  
  # joint determination of lag-order and break period #
  R.spec_jnt = speci.VAR(R.vars, lag_set=1:2, dim_m=2, trim=0.2)
  ### has no original results to replicate and test against.
  
  # check #
  expect_equivalent(R.speci$selection[2,2], R.urca@bp)
  ### expect_equal(rev(R.coint$stats_TR),       R.urca@teststat)
  ### Trace statistics are not equal because cajolst() uses just OLS-detrending,
  ### ... while coint.SL() uses GLS-detrending ('urca' version 1.3-4).
  expect_equivalent(t(R.spec_ic$df[ ,-1]),  R.vars_ic$criteria)
  expect_equivalent(R.spec_ic$selection,    R.vars_ic$selection)
  
  # test the stops and warnings #
  R.vars_none = VAR(sjd, type="none", p=1, season=4)
  expect_warning(speci.VAR(R.vars_none, lag_set=2, dim_m=1, trim=3, type_break="const"),
                 "'type_break' pertains the constant, which is missing in the provided VAR model 'x'.")
  expect_warning(speci.VAR(R.vars_none, lag_set=2, dim_m=1, trim=4, type_break="trend"),
                 "'type_break' pertains the linear trend, which is missing in the provided VAR model 'x'.")
})


test_that("sboot.pmb() under N=1 and sboot.mb() return identical Chol-VAR and BQ-VAR results", {
  # prepare data and arguments #
  data("PCAP")
  names_k = c("g", "k", "l", "y")  # variable names
  names_i = levels(PCAP$id_i)      # country names 
  L.data  = sapply(names_i, FUN=function(i) 
    ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), 
    simplify=FALSE)
  dim_p   = 2
  n.boot  = 3
  n.ahead = 3
  
  # estimate, identify, and bootstrap panel VAR #
  set.seed(8349)
  L.vars = lapply(L.data[1], FUN=function(x) vars::VAR(x, p=dim_p, type="none"))
  L.idBQ = lapply(L.vars[1], FUN=function(x) vars::BQ(x))
  R.pbBQ = sboot.pmb(L.idBQ, b.dim=c(5, FALSE), n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  
  set.seed(8349)
  R.pvar = pvarx.VAR(L.data[1], lags=dim_p, type="none")
  R.pid  = pid.chol(R.pvar)
  R.pmb  = sboot.pmb(R.pid, b.dim=c(5, FALSE), n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  R.pmb2 = sboot.pmb(R.pmb, b.dim=c(5, FALSE), n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  
  # estimate, identify, and bootstrap individual VAR #
  set.seed(8349)
  R.id  = R.pid$L.varx[[1]]
  R.mb  = sboot.mb(R.id, b.length=5, n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  R.mb2 = sboot.mb(R.mb, b.length=5, n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  
  set.seed(8349)
  R.idBQ = L.idBQ[[1]]
  R.mbBQ = sboot.mb(R.idBQ, b.length=5, n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  
  ###tolerant = 1.04e-04  
  ### Note that vars::VAR() is based on equation-wise lm() instead of multivariate LS. 
  ### ... The resulting computational differences are most notable in some bootstrapped IRF.
  ###R.svars = svars::id.chol(L.vars[[1]])  # dim_Kpn == Tob - 1 - k * p under type="const"
  ### Note that svars (version 1.3.11) does not count the number of deterministic regressors 
  ### ... when calculating the OLS covariance matrix of residuals. Only with a 
  ### ... single deterministic regressor, the results are identical to vars and pvars!
  ###colnames(R.svars$B) = colnames(R.pid$B)  # no automatic naming by id.chol()
  
  # check #
  expect_equal(R.pmb$true, R.mb$true)
  expect_equal(R.pmb$bootstrap, R.mb$bootstrap)
  expect_equal(R.pmb$A, R.mb$A)
  expect_equal(R.pmb$B, R.mb$B)
  
  expect_equal(R.pmb2$true, R.mb2$true)
  expect_equal(R.pmb2$bootstrap, R.mb2$bootstrap)
  expect_equal(R.pmb2$A, R.mb2$A)
  expect_equal(R.pmb2$B, R.mb2$B)
  
  expect_equal(R.pbBQ$true, R.mbBQ$true)
  expect_equal(R.pbBQ$bootstrap, R.mbBQ$bootstrap)
  expect_equal(R.pbBQ$A, R.mbBQ$A)
  expect_equal(R.pbBQ$B, R.mbBQ$B)
})


test_that("sboot.pmb() under N=1 and sboot.mb() return identical SVECM results", {
  # prepare data #
  data("PCAP")
  names_k = c("g", "k", "l", "y")  # variable names
  names_i = levels(PCAP$id_i)      # country names 
  L.data  = sapply(names_i, FUN=function(i) 
    ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), 
    simplify=FALSE)
  
  # arguments #
  type = "Case4"
  t_D1 = list(t_break=c(23, 49), t_shift=30)
  t_D2 = list(t_impulse=c(10), t_blip=12)
  ### A single element in t_D1 or t_D2 would be caught by aux_check() in pvarx.VEC(), ...
  ### ... but a panel model with single individual is not relevant anyway.
  LR = matrix(NA, nrow=4, ncol=4); LR[  , 1:2] = 0  # transitory shocks from g_t and k_t 
  SR = matrix(NA, nrow=4, ncol=4); SR[ 1, 2] = 0; SR[ 3, 4] = 0  # no instantaneous reaction of g_t to k_t resp. l_t to y_t 
  dim_r   = 2
  dim_p   = 2
  n.boot  = 3
  n.ahead = 3
  
  # estimate, identify, and bootstrap panel VAR #
  set.seed(8349)
  L.vecm = lapply(L.data["DNK"], FUN=function(x) VECM(x, dim_r=dim_r, dim_p=dim_p, type=type, t_D1=t_D1, t_D2=t_D2))
  L.id   = lapply(L.vecm["DNK"], FUN=function(x) id.grt(x, LR=LR, SR=SR))
  R.mb   = sboot.mb(L.id[["DNK"]], b.length=5, n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  
  set.seed(8349)
  R.pvec = pvarx.VEC(L.data["DNK"], dim_r=dim_r, lags=dim_p, type=type, t_D1=t_D1, t_D2=t_D2)
  R.pid  = pid.grt(R.pvec, LR=LR, SR=SR)
  R.pmb  = sboot.pmb(R.pid, b.dim=c(5, FALSE), n.boot=n.boot, n.ahead=n.ahead, n.cores=1)
  
  # extract a single individual from the panel application
  R.mg = sboot.mg(R.pid, n.ahead=n.ahead, idx_i="DNK")
  ### sboot.pmb could therewith show the result for a single ... 
  ### ... individual within the sampling of the complete panel ...
  ### ... (including potential selective pooled estimation).
  
  # check #
  expect_equal(R.mb$true, R.pmb$true)
  expect_equal(R.mb$bootstrap, R.pmb$bootstrap)
  expect_equal(R.mb$beta, R.pmb$beta)
  expect_equal(R.mb$A, R.pmb$A)
  expect_equal(R.mb$B, R.pmb$B)
})


test_that("sboot.pmb() return identical SVECM results under weighting and subsetting", {
  # prepare data #
  data("PCAP")
  names_k = c("g", "k", "l", "y")  # variable names
  names_i = levels(PCAP$id_i)      # country names 
  L.data  = sapply(names_i, FUN=function(i) 
    ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), 
    simplify=FALSE)
  
  # arguments #
  type = "Case4"
  t_D1 = list(t_break=c(23, 49), t_shift=30)
  t_D2 = list(t_impulse=c(10), t_blip=12)
  LR = matrix(NA, nrow=4, ncol=4); LR[  , 1:2] = 0  # transitory shocks from g_t and k_t 
  SR = matrix(NA, nrow=4, ncol=4); SR[ 1, 2] = 0; SR[ 3, 4] = 0  # no instantaneous reaction of g_t to k_t resp. l_t to y_t 
  dim_r   = 2
  dim_p   = 2
  n.boot  = 3
  n.ahead = 3
  b.dim   = c(FALSE, 2)
  w1 = c(TRUE, TRUE, FALSE, FALSE, FALSE)
  w2 = c(1, 1, 0, 0, 0)
  
  # estimate, identify, and bootstrap panel VAR #
  R.pvec = pvarx.VEC(L.data[1:5], dim_r=dim_r, lags=dim_p, type=type, t_D1=t_D1, t_D2=t_D2)
  R.pid  = pid.grt(R.pvec, LR=LR, SR=SR)
  set.seed(8349)
  R.pmb1 = sboot.pmb(R.pid, b.dim=b.dim, n.boot=n.boot, n.ahead=n.ahead, n.cores=1, fix_beta=FALSE, w=w1)
  set.seed(8349)
  R.pmb2 = sboot.pmb(R.pid, b.dim=b.dim, n.boot=n.boot, n.ahead=n.ahead, n.cores=1, fix_beta=FALSE, w=w2)
  
  # check #
  expect_equal(R.pmb1$true, R.pmb2$true)
  expect_equal(R.pmb1$bootstrap, R.pmb2$bootstrap)
  expect_equal(R.pmb1$beta, R.pmb2$beta)
  expect_equal(R.pmb1$A, R.pmb2$A)
  expect_equal(R.pmb1$B, R.pmb2$B)
})