####################
###  UNIT TESTS  ###
####################
#
# For exported "id" functions, 
# their auxiliary modules, 
# and residual diagnostics.


test_that("the residual diagnostics reproduce example from vars package", {
  tolerant = 0.0001  # for comparing p-values
  
  # VAR model #
  library("vars")
  data(Canada)
  var.2c <- vars::VAR(Canada, p = 2, type = "const") ### p=3 generates computationally singular matrices in test.serial(lag.h>2)
  
  # vars package #
  R.norm = vars::normality.test(var.2c, multivariate.only = FALSE)
  R.arch = vars::arch.test(var.2c, lags.multi = 5, lags.single = 5, multivariate.only = FALSE)
  R.seBG = vars::serial.test(var.2c, lags.bg = 2, type = "BG")
  R.seES = vars::serial.test(var.2c, lags.bg = 2, type = "ES")
  
  vars_norm = sapply(R.norm$jb.mul, function(x) x$p.value)
  vars_arch = R.arch$arch.mul$p.value
  vars_seri = c(R.seBG$serial$p.value, R.seES$serial$p.value)
  
  # pvars package: p-values from multivariate tests #
  resids = resid(var.2c)
  our_norm = test.normality(u=resids)
  our_arch = test.arch(u=resids, lag.h=5)
  our_seri = test.serial(u=resids, lag.h=2, x=var.2c$datamat[ ,-(1:ncol(var.2c$y))])
  
  # check #
  expect_equivalent(our_norm, vars_norm, tolerance=tolerant)
  expect_equivalent(our_arch, vars_arch, tolerance=tolerant)
  expect_equivalent(our_seri, vars_seri, tolerance=tolerant)
  
  # pvars package: p-values from univariate tests #
  # our_unorm = apply(resids, MARGIN=2, FUN=function(x) test.normality(u=x))
  # our_uarch = apply(resids, MARGIN=2, FUN=function(x) test.arch(u=x, lag.h=5))
  # our_useri = apply(resids, MARGIN=2, FUN=function(x) test.serial(u=x, lag.h=2, x=var.2c$datamat[ ,-(1:ncol(var.2c$y))]))
})


test_that("id.iv() returns the same B for all S2=c('MR', 'JL', 'NQ') if there is L=1 proxy", {
  # prepare data #
  data("PCIT")
  names_k = c("APITR", "ACITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")
  names_l = c("m_PI")  # name of single proxy
  
  # estimate and identify #
  R.vars = vars::VAR(PCIT[ , names_k], p=4, type="const")
  R.ivMR = id.iv(R.vars, iv=PCIT[-(1:4), names_l], cov_u="OMEGA", S2="MR")
  R.ivJL = id.iv(R.vars, iv=PCIT[-(1:4), names_l], cov_u="OMEGA", S2="JL")
  R.ivNQ = id.iv(R.vars, iv=PCIT[-(1:4), names_l], cov_u="OMEGA", S2="NQ")
  
  # check #
  expect_equal(R.ivJL$B, R.ivMR$B)
  expect_equal(R.ivJL$B, R.ivNQ$B[ , 1, drop=FALSE])
})


test_that("id.iv(S2='MR') can reproduce the point estimates of 'PCIT' in Jentsch,Lunsford 2019:2668, Ch.III", {
  ### Bootstrap results in the example of id.iv() are congruent with Fig.1 and 2. 
  tolerant = 0.1
  
  # prepare data #
  data("PCIT")
  R.norm = function(B) B / matrix(-diag(B), nrow(B), ncol(B), byrow=TRUE)
  
  # estimate and identify under ordering "blue" in Fig.1 and 2 #
  names_k1 = c("APITR", "ACITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")
  names_l1 = c("m_PI", "m_CI")  # proxy names
  R.var1 = vars::VAR(PCIT[ , names_k1], p=4, type="const")
  R.iv1  = id.iv(R.var1, iv=PCIT[-(1:4), names_l1], S2="MR", cov_u="OMEGA")
  R.irf1 = irf.varx(R.iv1, n.ahead=20, normf=R.norm)
  
  # check "blue" #
  our_B   = unlist(R.irf1$irf[1, -1])
  their_B = c(#1 | #2 Figure
    APITR = c(-1, 0.05),
    ACITR = c(0.6,  -1),
    PITB  = c(0.6,  NA), 
    CITB  = c(NA,  3.2), 
    GOV   = c(0.1, 0.6), 
    RGDP  = c(1.3, 0.4),
    DEBT  = c(NA,   NA))
  idx_na = is.na(their_B)
  expect_equivalent(our_B[!idx_na], their_B[!idx_na], tolerance=tolerant)
  ### cbind(unname(our_B[!idx_na]), unname(their_B[!idx_na]))
  
  # estimate and identify under ordering "red" in Fig.1 and 2 #
  names_k2 = c("ACITR", "APITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")
  names_l2 = c("m_CI", "m_PI")
  R.var2 = vars::VAR(PCIT[ , names_k2], p=4, type="const")
  R.iv2  = id.iv(R.var2, iv=PCIT[-(1:4), names_l2], S2="MR", cov_u="OMEGA")
  R.irf2 = irf.varx(R.iv2, n.ahead=20, normf=R.norm)
  
  # check "red" #
  our_B   = unlist(R.irf2$irf[1, -1])
  their_B = c(#2 | #1 Figure
    ACITR = c(-1,  0.03),
    APITR = c(0.07, -1),
    PITB  = c(NA,  0.7), 
    CITB  = c(3.2,  NA), 
    GOV   = c(0.6, 0.2), 
    RGDP  = c(0.3, 1.4),
    DEBT  = c(NA,   NA))
  idx_na = is.na(their_B)
  expect_equivalent(our_B[!idx_na], their_B[!idx_na], tolerance=tolerant)
  ### cbind(unname(our_B[!idx_na]), unname(their_B[!idx_na]))
})


test_that("pid.iv(c('pool','indiv')) can reproduce the id.iv(S2='NQ') example for 'PCIT' under duplicated N=1", {
  # prepare data #
  data("PCIT")
  names_k = c("APITR", "ACITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")
  names_l = c("m_PI", "m_CI")  # proxy names
  L.iv = list(USA=PCIT[-(1:4), names_l], USB=PCIT[-(1:4), names_l])
  
  # estimate and identify #
  R.vars = vars::VAR(PCIT[ , names_k], p=4, type="const")
  L.vars = list(USA=R.vars, USB=R.vars )  # list of N=2 duplication of a single individual
  R.iv1  = id.iv(L.vars[[1]], iv=L.iv[[1]], S2="NQ")
  R.iv2  = pid.iv(L.vars, S2="NQ", cov_u="OMEGA", combine="indiv", iv=L.iv[[1]])  # duplicates the common proxy
  R.iv3  = pid.iv(L.vars, S2="NQ", cov_u="OMEGA", combine="pool", iv=L.iv)
  R.iv4  = as.pvarx(list(USA=R.iv1, USB=R.iv1))
  
  # normalized IRF #
  R.norm = function(B) B / matrix(-diag(B), nrow(B), ncol(B), byrow=TRUE)
  R.irf1 = irf(R.iv1, normf=R.norm)
  R.irf2 = irf(R.iv2, normf=R.norm)
  R.irf3 = irf(R.iv3, normf=R.norm)
  R.irf4 = irf(R.iv4, normf=R.norm)
  
  # check #
  expect_equal(R.iv2$args_pid$L.iv, R.iv3$args_pid$L.iv)
  expect_equal(R.iv2$args_pid$L.iv, R.iv4$args_pid$L.iv)
  
  expect_equal(R.iv1$B, R.iv2$B)
  expect_equal(R.iv1$B, R.iv3$B)
  expect_equal(R.iv1$B, R.iv4$B)
  
  expect_equal(R.irf1, R.irf2)
  expect_equal(R.irf1, R.irf3)
  expect_equal(R.irf1, R.irf4)
})


test_that("the impact matrix from id.grt() fulfills BB'=OMEGA in a just-identified SVECM with period-pecific deterministic term", {
  # prepare data #
  data("PCAP")
  names_k = c("g", "k", "l", "y")  # variable names
  data_i  = ts(PCAP[PCAP$id_i=="DEU", names_k], start=1960, end=2019, frequency=1)
  
  # specifications #
  LR = matrix(NA, nrow=4, ncol=4); LR[  , 1:2] = 0
  SR = matrix(NA, nrow=4, ncol=4); SR[ 1, 2] = 0; SR[ 3, 4] = 0
  startv = c(-0.363, -1.448, -0.365, 0.696, -0.090, 0.957, 1.044, 0.083, -0.169, 0.044)
  
  # estimate and identify #
  R.vecm = VECM(y=data_i, dim_p=2, dim_r=2, type="Case4", t_D1=list(t_break=c(23, 49)))
  R.svec = id.grt(R.vecm, LR=LR, SR=SR, start=startv)
  
  # check #
  B = unname(R.svec$B)
  vecm_OMEGA = unname(R.vecm$VECM$OMEGA)
  svec_OMEGA = unname(R.svec$OMEGA)
  expect_equivalent(B%*%t(B), vecm_OMEGA)
  expect_equivalent(B%*%t(B), svec_OMEGA)
})


test_that("VECM(), id.grt(), and pid.grt() can reproduce the basic example from vars package under r=1 and N=1", {
  # prepare data #
  library(vars)
  data(Canada)
  names_k = c("prod", "e", "U", "rw")
  CAN = Canada[ , names_k]
  
  # joint specifications #
  dim_p = 3  # lag-order
  dim_r = 1  # cointegration rank
  SR = matrix(NA, nrow=4, ncol=4, dimnames=list(names_k, names_k))
  SR[4, 2] = 0
  LR = matrix(NA, nrow=4, ncol=4, dimnames=list(names_k, names_k))
  LR[1, 2:4] = 0
  LR[2:4, 4] = 0
  startv = c(-0.150, 1.011, -1.457, 0.044, 0.085, 0.467, 0.053, -1.073, 2.073, -1.888)
  
  # pvars: estimate and identify #
  set.seed(4679)
  R.vecm = VECM(y=CAN, dim_p=dim_p, dim_r=dim_r, type="Case4")
  R.grt  = id.grt(R.vecm, LR=LR, SR=SR, start=startv)
  
  set.seed(4679)
  L.vecm = list(A1=R.vecm)  # list of N=1 individual
  R.pgrt = pid.grt(L.vecm, LR=LR, SR=SR, start=startv)
  
  # urca and vars: estimate and identify #
  set.seed(4679)
  R.cajo = urca::ca.jo(CAN, type="trace", ecdet="trend", K=dim_p, spec="transitory")
  R.svec = vars::SVEC(R.cajo, LR=LR, SR=SR, r=dim_r, start=startv, lrtest=FALSE, boot=FALSE)
  
  # check #
  expect_equal(R.pgrt$B, as.varx(R.svec)$B)
  expect_equal(R.grt$B,  as.varx(R.svec)$B)
  expect_equal(R.grt$SR, R.svec$SR)
  expect_equal(R.grt$LR, R.svec$LR)
})


test_that("aux_idGRT() is correctly grafted under r=2", {
  # prepare data #
  data("PCAP")
  names_k = c("g", "k", "l", "y")  # variable names
  data_i  = ts(PCAP[PCAP$id_i=="AUT", names_k], start=1960, end=2019, frequency=1)
  
  # joint specifications #
  dim_p = 2  # lag-order
  dim_r = 2  # cointegration rank
  LR = matrix(NA, nrow=4, ncol=4, dimnames=list(names_k, names_k)); LR[  , 1:2] = 0
  SR = matrix(NA, nrow=4, ncol=4, dimnames=list(names_k, names_k)); SR[ 1, 2] = 0; SR[ 3, 4] = 0
  startv  = c(-0.363, -1.448, -0.365, 0.696, -0.090, 0.957, 1.044, 0.083, -0.169, 0.044)
  n.ahead = 10
  n.boot  = 10
  
  # pvars: estimate, identify, and bootstrap #
  set.seed(8349)
  R.vecm = VECM(y=data_i, dim_p=dim_p, dim_r=dim_r, type="Case4")
  R.grt  = id.grt(R.vecm, LR=LR, SR=SR, start=startv)
  R.boot = sboot.mb(R.grt, b.length=1, n.ahead=n.ahead, n.boot=n.boot, n.cores=1)
  
  # urca and vars: estimate, identify, and bootstrap #
  set.seed(8349)
  R.cajo = urca::ca.jo(data_i, type="eigen", ecdet="trend", K=dim_p, spec="transitory")
  R.svec = vars::SVEC(R.cajo, LR=LR, SR=SR, r=dim_r, start=startv, lrtest=FALSE, boot=FALSE)
  R.irf  = vars::irf(R.svec, n.ahead=n.ahead, boot=TRUE, ci=0.9, runs=n.boot)
  
  # align confidence bands, from svars::plot.sboot #
  kk <- ncol(R.boot$true$irf)
  intervals <- array(0, dim = c(n.ahead+1, kk, n.boot))
  for(i in 1:n.boot){
    intervals[,,i] <- as.matrix(R.boot$bootstrap[[i]]$irf)
  }
  R.bootLW <- array(0, dim = c(n.ahead, kk, 1))
  R.bootUP <- array(0, dim = c(n.ahead, kk, 1))
  for(i in 1:n.ahead){
    for(j in 1:kk){
      R.bootLW[i,j, ] <- quantile(intervals[i,j, ], probs = 0.05)
      R.bootUP[i,j, ] <- quantile(intervals[i,j, ], probs = 0.95)
    }
  }
  
  # check #
  expect_equal(R.grt$SR, R.svec$SR)
  expect_equal(R.grt$LR, R.svec$LR)
  #expect_equal(R.bootLW, R.irf$Lower)
  #expect_equal(R.bootUP, R.irf$Upper)
})


test_that("aux_sico() returns the correct permutation with names", {
  ### This is rather a reminder not to change the output. ###
  # define #
  names_s = paste0("eps_", 1:10)
  L.mat = list(A=cbind(2:3, -c(2:3)),
               B=cbind(-c(3:5), c(5, 0, 1), c(1:3), -c(4, 2, 1)),
               C=diag(5))
  
  for(mat_B in L.mat){
    dim_K = nrow(mat_B)
    dim_S = ncol(mat_B)
    
    # apply permutation matrix and its alternative #
    R.sico = aux_sico(mat_B, names_s=names_s[1:dim_S])
    mat_si = matrix(1, nrow=dim_K, ncol=dim_S)
    mat_si[ , R.sico$idx_si] = mat_si[ , R.sico$idx_si, drop=FALSE] * -1
    B.sico = mat_B[ , R.sico$idx_co, drop=FALSE] * mat_si
    
    # check #
    expect_equivalent(R.sico$B, B.sico)
    expect_equal(colnames(R.sico$B), names_s[1:dim_S])
    expect_equal(colnames(R.sico$mat_perm), names_s[1:dim_S])
  }
})