tests/testthat/test_pid.R

In pvars: VAR Modeling for Heterogeneous Panels

####################
###  UNIT TESTS  ###
####################
#
# For exported "pid" functions 
# and their subordinated modules.


test_that("pid.cvm() and pid.dc() return identical objects for group-ICA with one common factor", {
  # 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)
    
  # estimate #
  L.vars = lapply(L.data, FUN=function(x) vars::VAR(x, p=1, type="both"))
  R.pcvm = pid.cvm(L.vars, combine="group", n.factors=1)
  R.pdc  = pid.dc(L.vars,  combine="group", n.factors=1)
  
  # check #
  R.pcvm$args_pid = NULL
  R.pdc$args_pid  = NULL
  R.pcvm$L.varx = lapply(R.pcvm$L.varx, function(x_i) x_i$args_id = NULL)
  R.pdc$L.varx  = lapply(R.pdc$L.varx,  function(x_i) x_i$args_id = NULL)
  expect_equal(R.pcvm, R.pdc)
})


test_that("pid.dc(c('pool','indiv')) can reproduce the id.dc() example from svars package under N=1", {
  # estimate #
  library(svars)
  L.vars = list(USA = vars::VAR(USA, p=6, type="const"))
  
  # identify #
  set.seed(23211)
  R.pool = pid.dc(L.vars, combine="pool",  n.iterations=100, PIT=FALSE)
  set.seed(23211)
  R.indi = pid.dc(L.vars, combine="indiv", n.iterations=100, PIT=FALSE)
  set.seed(23211)
  R.svar = svars::id.dc(L.vars$USA, PIT=FALSE)
  
  # baseline shocks #
  R.varx = as.varx(L.vars$USA)
  R.eps  = t(solve(t(chol(R.varx$SIGMA))) %*% R.varx$resid)
  
  # check #
  R.svar$B = aux_sico(R.svar$B)$B
  expect_equivalent(R.svar$B, R.indi$B, tolerance=5e-05)  # identical on Linux Ubuntu (22.04)
  expect_equivalent(R.svar$B, R.pool$B, tolerance=5e-05)  # win-builder needs higher tolerance than other platforms
  expect_equivalent(R.eps, R.pool$eps)  
  ### View(R.eps - R.pool$eps)  # small rounding errors lead to slightly different rotation angles
})


test_that("pid.cvm(c('pool','indiv')) can reproduce the id.cvm() example from svars package under N=1", {
  # estimate #
  library(svars)
  L.vars = list(USA = vars::VAR(USA, p=6, type="const"))
  
  # identify #
  set.seed(23211)
  dd = copula::indepTestSim(L.vars$USA$obs, L.vars$USA$K, verbose=FALSE)
  set.seed(23211)
  R.pool = pid.cvm(L.vars, combine="pool",  dd=dd, itermax=25, steptol=100, iter2=5)
  set.seed(23211)
  R.indi = pid.cvm(L.vars, combine="indiv", dd=dd, itermax=25, steptol=100, iter2=5)
  set.seed(23211)
  R.svar = svars::id.cvm(L.vars$USA, dd=dd, itermax=25, steptol=100, iter2=5)
  
  # check #
  R.svar$B = aux_sico(R.svar$B)$B
  expect_equivalent(R.svar$B, R.indi$B)
  expect_equivalent(R.svar$B, R.pool$B)
})


test_that("aux_cvmICA() is correctly grafted", {
  # prepare data #
  data("USA")
  
  # joint arguments and model estimation #
  set.seed(23211)
  v1 = vars::VAR(USA, lag.max = 10, ic = "AIC" )
  dd = copula::indepTestSim(v1$obs, v1$K, verbose=FALSE)
  dim_p   = unname(v1$p)  # lag-order
  itermax = 500/10
  steptol = 100/10
  iter2   = 75/15
  n.ahead = 10
  n.boot  = 10
  n.cores = 1
  
  # pvars: identify and bootstrap #
  set.seed(23211)
  R.varx = as.varx.varest(v1)
  R.chol = t(chol(R.varx$SIGMA))
  R.eps  = solve(R.chol) %*% R.varx$resid  # pre-whitened shocks form baseline decomposition
  R.ICA  = aux_cvmICA(t(R.eps), dd=dd, itermax=itermax, steptol=steptol, iter2=iter2)
  R.varx$B = R.chol %*% R.ICA$W
  #R.varx$arg_id = list(dd=dd, itermax=30, steptol=20, iter2=5)
  #R.boot1 = sboot.mb(R.varx, b.length=1, n.ahead=n.ahead, n.boot=n.boot, n.cores=n.cores)
  
  # svars: identify and bootstrap #
  set.seed(23211)
  R.svars = svars::id.cvm(v1, dd=dd, itermax=itermax, steptol=steptol, iter2=iter2)
  #R.boot2 = svars::mb.boot(R.svars, b.length=1, n.ahead=n.ahead, nboot=n.boot, nc=n.cores,
  #                         dd=dd, itermax=30, steptol=20, iter2=5)
  
  # check #
  expect_equal(R.ICA$theta, R.svars$rotation_angles)
  expect_equal(R.varx$B, R.svars$B)
  #expect_equal(R.boot1, R.boot2)
})


test_that("aux_UGroup() and aux_UPool() provide whitened residuals", {
  # prepare data and estimate #
  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)
  L.cajo = lapply(L.data, FUN=function(x) urca::ca.jo(x, K=2, ecdet="trend", spec="transitory", type="trace"))
  L.vars = lapply(L.cajo, FUN=function(x) vars::vec2var(x, r=3))
  
  # apply #
  L.varx  = as.pvarx(L.vars)$L.varx
  L.Ucov  = lapply(L.varx, FUN=function(x) x$OMEGA) # use MLE covariance matrix
  L.resid = lapply(L.varx, FUN=function(x) x$resid)
  dim_K   = nrow(L.varx[[1]]$A)
  R.group = aux_UGroup(L.resid, n.factors=dim_K)
  R.pool  = aux_UPool(L.resid, L.Ucov)
  
  # check #
  our_UGroup = t(R.group$eps) %*% (R.group$eps) / nrow(R.group$eps)
  our_UPool  = t(R.pool$eps)  %*% (R.pool$eps)  / nrow(R.pool$eps) 
  ### Do not correct for degrees-of-freedom when using MLE covariance matrix.
  expect_equivalent(our_UGroup, diag(dim_K))
  expect_equivalent(our_UPool,  diag(dim_K))
})