test_coint.R
In pvars: VAR Modeling for Heterogeneous Panels

####################
###  UNIT TESTS  ###
####################
#
# For exported "coint" functions 
# and their auxiliary modules.



test_that("OLS of the transformed model and GLS are identical in aux_GLStrend.", {
  # 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
  L.data = sapply(names_i, FUN=function(i) ts(ERPT[ERPT$id_i==i, names_k], start=c(1995, 1), frequency=12), simplify=FALSE)
  
  # GLSE of the deterministic term #
  for(dim_p in 0:3){
    R.vecm = VECM(y=L.data$France, dim_r = 2, dim_p = 4, type = "Case4", t_D1=list(t_break=89))
    D_VAR  = aux_dummy(dim_T=R.vecm$dim_T+R.vecm$dim_p, t_break=R.vecm$t_D1$t_break, t_shift=R.vecm$t_D1$t_break, type="both")
    #alpha_oc = if(R.vecm$dim_r==0){ diag(dim_K) }else{ R.vecm$RRR$S00inv %*% R.vecm$RRR$S01 %*% R.vecm$RRR$V[ ,(R.vecm$dim_r+1):R.vecm$dim_K] }
    alpha_oc = MASS::Null(R.vecm$VECM$alpha)
    R.GLSd1 = aux_GLStrend(y=L.data$France, OMEGA=R.vecm$VECM$OMEGA, A=R.vecm$A, D=D_VAR, dim_p=R.vecm$dim_p)
    R.GLSd2 = aux_GLStrend(y=L.data$France, alpha=R.vecm$VECM$alpha, alpha_oc=alpha_oc,
                           OMEGA=R.vecm$VECM$OMEGA, A=R.vecm$A, D=D_VAR, dim_p=R.vecm$dim_p)
    
    # check #
    expect_equal(R.GLSd2$Q%*%t(R.GLSd2$Q), solve(R.vecm$VECM$OMEGA))
    expect_equal(R.GLSd1$MU, R.GLSd2$MU)
    
    # OLSE of the deterministic term #
    #R.GLSd1 = aux_GLStrend(y=L.data$France, OMEGA=diag(R.vecm$dim_K), A=R.vecm$A, D=D_VAR, dim_p=R.vecm$dim_p)
  }
})


test_that("aux_CointMoments can reproduce the examples of Doornik 1998:587, 
  Johansen et al. 2000:235, Trenkler et al. 2008:350, and  Kurita,Nielsen 2019:19", {
  # Johansen procedure with a weakly exogenous variable #
  our_JO   = aux_CointMoments(dim_K=3, dim_L=1, type="Case4")[ , 1:2]  # remove moments for ME
  their_JO = cbind(TR_EZ = c(36.52, 20.43, 8.265),
                   TR_VZ = c(59.55, 34.23, 14.40))  # Doornik 1998:587, Ch.10.2 for r_H0=0:2
  rownames(their_JO) = rownames(our_JO)
  expect_equal(our_JO, their_JO, tolerance = 0.005)
  
  # JMN (2000) procedure #
  moments = aux_CointMoments(dim_K=5, dim_T=91, t_D1=list(t_break=c(27,77)), type="JMN_Case4")
  m = moments[ , 1, drop=FALSE]  # means for TR
  v = moments[ , 2, drop=FALSE]  # variances for TR
  LR.stats  = c(256.46, 157.43, 71.96, 29.50, 10.42)  # Johansen et al. 2000:235, Tab.8
  our_JMN   = 1-pgamma(LR.stats, shape=m^2/v, rate=m/v)
  their_JMN = c(0, 0, 0.043, 0.621, 0.745)
  expect_equal(our_JMN, their_JMN, tolerance = 0.005)
  
  # TSL (2008) procedure #
  our_TSL   = aux_CointMoments(dim_K=3, dim_T=184, t_D1=list(t_break=125), r_H0=1, type="TSL_trend")[ , 1:2]  # remove moments for ME
  their_TSL = c(TR_EZ=11.3009, TR_VZ=17.5418)  # Trenkler et al. 2008:350
  expect_equal(our_TSL, their_TSL, tolerance = 0.0005)
  
  # KN (2019) procedure #
  # counted signs of coefficients in Kurita,Nielsen 2019:21/22, Tab.A1/A2
  expect_equal( colSums(coint_rscoef$KN_Case2 < 0), c(TR_lam=14, TR_del=17, TR_cov=18) )
  expect_equal( colSums(coint_rscoef$KN_Case2 > 0), c(TR_lam=17, TR_del=17, TR_cov=16) )
  expect_equal( colSums(coint_rscoef$KN_Case4 < 0), c(TR_lam=17, TR_del=16, TR_cov=13) )
  expect_equal( colSums(coint_rscoef$KN_Case4 > 0), c(TR_lam=18, TR_del=18, TR_cov=11) )
  
  # Kurita,Nielsen 2019:19, Ch.5, Tab.4 and supplementary spreadsheet
  moments = aux_CointMoments(dim_K=2, dim_L=3, dim_T=94, type="KN_Case4", t_D1=list(t_break=71))[ , 1:2]  # remove moments for ME
  m = moments[ , 1, drop=FALSE]  # means for TR
  v = moments[ , 2, drop=FALSE]  # variances for TR
  LR.stats = c(56.610, 21.964)
  KN_pvals = 1-pgamma(LR.stats, shape=m^2/v, rate=m/v)
  KN_cvals = c(qgamma(0.95, shape=m^2/v, rate=m/v))
  our_KN   = cbind(KN_pvals, KN_cvals)
  their_KN = cbind(c(0.014, 0.148), c(50.864, 26.334))
  dimnames(their_KN) = dimnames(our_KN)
  expect_equal(our_KN, their_KN, tolerance = 0.005)
  
  # approximate equivalence of JMN (2000) and KN (2019) 
  JMN_moments = aux_CointMoments(dim_K=4, dim_T=100, t_D1=list(t_break=c(25,75)), type="JMN_Case4")[ ,1:2]
  KN_moments  = aux_CointMoments(dim_K=4, dim_T=100, t_D1=list(t_break=c(25,75)), type="KN_Case4")[ ,1:2]
  zero_matrix = matrix(0, nrow=4, ncol=2, dimnames=dimnames(JMN_moments))
  expect_equal((JMN_moments-KN_moments)/JMN_moments, zero_matrix, tolerance = 0.09)
})


test_that("coint.JO() and VECM() can reproduce the basic examples with four seasons from urca package", {
  # prepare data #
  library("urca")
  data(denmark)
  sjd = denmark[, c("LRM", "LRY", "IBO", "IDE")]
  dim_K = ncol(sjd)  # number of endogenous variables
  
  # perform Johansen procedure on full VECM #
  sjd.vecm = urca::ca.jo(sjd, ecdet="const", type="eigen", K=2, spec="transitory", season=4)
  sjd.vars = vars::vec2var(sjd.vecm, r=1)
  R.JOrank = coint.JO(y=sjd,      dim_p=2, type="Case2", t_D2=list(n.season=4))
  R.JOvecm = VECM(y=sjd, dim_r=1, dim_p=2, type="Case2", t_D2=list(n.season=4))
  
  # perform Johansen procedure on partial VECM with structural shift #
  R.KNrank = coint.JO(y=sjd[ , c("LRM"), drop=FALSE],   dim_p=2, 
                      x=sjd[ , c("LRY", "IBO", "IDE")], dim_q=2, 
                      type="Case2", t_D1=list(t_shift=36), t_D2=list(n.season=4))
  R.KNvecm = VECM(y=sjd[ , c("LRM"), drop=FALSE],   dim_p=2, 
                  x=sjd[ , c("LRY", "IBO", "IDE")], dim_q=2, dim_r=1, 
                  type="Case2", t_D1=list(t_shift=36), t_D2=list(n.season=4))
  
  # check #
  our_statsME  = R.JOrank$stats_ME
  urca_statsME = sort(sjd.vecm@teststat, decreasing=TRUE)
  urca_resids  = t(resid(sjd.vars)); dimnames(urca_resids) = dimnames(R.JOvecm$resid)
  
  expect_equal(our_statsME,    urca_statsME)
  expect_equal(R.JOvecm$resid, urca_resids)
  expect_equal(R.JOrank$lambda[1:dim_K], R.JOvecm$RRR$lambda[1:dim_K])
  expect_equal(R.KNrank$lambda[1],       R.KNvecm$RRR$lambda[1])
})


### reproduce Oersal,Arsova 2016:22, Tab.7.5/8 ###
# R.JOrank  = coint.JO(y=sjd, dim_p=3, type="Case4")
# R.TSLrank = coint.SL(y=L.data$France, dim_p=3, type_SL="SL_trend", t_D=list(t_break=89))
# R.MSBrank = coint.MSB(eit=sjd, lag_max = 5, MIC="AIC", type_MSB = "MSB_trend")