#################### ### UNIT TESTS ### #################### # # For exported varx functions # and their auxiliary modules. test_that("the residuals from the rank-restricted VAR and VECM representation in VECM() are identical.", { # 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) dim_p = 3 for(type in c("Case1", "Case2", "Case3", "Case4", "Case5")){ # estimate VECM and VAR residuals # R.vecm = VECM(y=data_i, dim_r=2, dim_p=dim_p, type=type) type_A = switch(type, "Case1"="none", "Case2"="const", "Case3"="const", "Case4"="both", "Case5"="both") Y = t(data_i)[ ,-(1:dim_p)] D = aux_dummy(type=type_A, dim_T=nrow(data_i)) Z = aux_stack(data_i, dim_p=dim_p, D=D) if(type == "Case4"){ Z["trend", ] = Z["trend", ] - 1 } e = Y - R.vecm$A %*% Z # check identity of VECM and level-VAR representation # our_VAR = e # ... via residuals our_VECM = R.vecm$resid expect_equal(our_VAR, our_VECM) } }) test_that("VECM() can reproduce basic examples from urca package at different lag-orders and determinstic regressors D2.", { # prepare data # library("urca") data(denmark) sjd = denmark[, c("LRM", "LRY", "IBO", "IDE")] tolerant_u = 4e-08 # on residuals tolerant_a = 2e-06 # on LR-test stats and pvals for adjustment coefficients dim_K = ncol(sjd) dim_r = 2 for(dim_p in 2:9){ # lag-order of the endogenous variables in the level-VAR # estimate VECM # R.cajo = urca::ca.jo(sjd, ecdet="trend", type="trace", K=dim_p, spec="transitory") ### R.vecm_urca = urca::cajorls(R.cajo, r=dim_r) R.vars = vars::vec2var(R.cajo, r=dim_r) R.vecm = VECM(y=sjd, dim_r=dim_r, dim_p=dim_p, type="Case4") # check via residuals # our_VECM = R.vecm$resid urca_VECM = t(R.vars$resid) expect_equivalent(our_VECM, urca_VECM, tolerance = tolerant_u) # perform LR-test on weak exogeneity in conditional VECM # DA = diag(dim_K)[ ,-4, drop=FALSE] ### matrix(c(1,0,0,0), c(4,1)) R.alr = urca::alrtest(R.cajo, r=dim_r, A=DA) #### R.cvec_urca = urca::cajorls(R.alr, r=dim_r) R.cvec = VECM(y=sjd[ ,1:3], dim_p=dim_p, x=sjd[ , 4], dim_q=dim_p, dim_r=dim_r, type="Case4") our_LRalpha = aux_LRrest(lambda = R.vecm$RRR$lambda, lambda_rest = R.cvec$RRR$lambda, dim_r=dim_r, dim_T=R.cvec$dim_T, dim_L=R.cvec$dim_L) # check # urca_LRalpha = data.frame(stats=R.alr@teststat, pvals=R.alr@pval[1]) expect_equal(our_LRalpha, urca_LRalpha, tolerance = tolerant_a) } # estimate VECM with customized regressors in D2 # R.D2 = t(aux_dummy(dim_T=nrow(sjd), t_blip=40, t_impulse=c(10+0:3, 30), t_shift=10, n.season=4)) dim_p = 4 R.cajo_D2 = urca::ca.jo(sjd, ecdet="trend", type="trace", K=dim_p, spec="transitory", dumvar=R.D2) ### R.vecm_urca_D2 = urca::cajorls(R.cajo_D2, r=dim_r) R.vars_D2 = vars::vec2var(R.cajo_D2, r=dim_r) R.vecm_D2 = VECM(y=sjd, dim_r=dim_r, dim_p=dim_p, type="Case4", D2=R.D2) # check via residuals # our_VECM_D2 = R.vecm_D2$resid urca_VECM_D2 = t(R.vars_D2$resid) expect_equivalent(our_VECM_D2, urca_VECM_D2, tolerance = tolerant_u) }) test_that("auxVAR modules can reproduce OLS examples from vars package.", { # prepare data # library("vars") data(Canada) dim_T = nrow(Canada) # number of observations including presample to_resid = 2.8e-05 # VAR() is based on equation-wise lm(). tolerant = 3.8e-08 for(dim_p in 1:10){ # lag-order of the endogenous variables # stack operators # Y = t(Canada)[c(1,3,4),-(1:dim_p)] D = aux_dummy(dim_T=dim_T, center=TRUE, n.season=4, type="const") Z = aux_stack(y=Canada[ ,c(1,3,4)], dim_p=dim_p, x=Canada[ , 2, drop=FALSE], dim_q=0, D=D) B = tcrossprod(Y, Z) %*% solve(tcrossprod(Z)) # Y%*%t(Z) %*% solve(Z%*%t(Z)) R.def = aux_stackOLS(Canada[, c(1,3,4)], dim_p=dim_p, x=Canada[, 2, drop=FALSE], dim_q=0, D=D) R.vars = vars::VAR(y=Canada[,-2], p=dim_p, type="const", exogen=Canada[ , 2, drop=FALSE], season=4) names_slp = paste0(rownames(Y), ".l1") theirs = t(sapply(R.vars$varresult, FUN = function(k) k$coefficients))[ , c(names_slp, "prod")] expect_equal(B[ , names_slp], theirs[ , names_slp], tolerance = tolerant) # Block of seasonal dummies is shifted. expect_identical(Y, R.def$Y) expect_identical(Z, R.def$Z) # multivariate OLS estimation # R.sens = aux_dummy(dim_T=dim_T, center=TRUE, n.season=4) R.varx = aux_VARX(y=Canada[, -2], dim_p=dim_p, x=Canada[ , 2, drop=FALSE], dim_q=0, type="both", D=R.sens) R.vars = vars::VAR(y=Canada[,-2], p=dim_p, type="both", exogen=Canada[ , 2, drop=FALSE], season=4) our_resid = unname(R.varx$resid) varx_resid = unname(as.varx(R.vars)$resid) # coerce the VAR to varx object first vars_resid = unname(t(resid(R.vars))) expect_equal(our_resid, vars_resid, tolerance = to_resid) expect_equal(our_resid, varx_resid, tolerance = to_resid) our_slp = R.varx$A[ , c(names_slp, "prod.l0"), drop=FALSE] vars_slp = t(sapply(R.vars$varresult, FUN = function(k) k$coefficients))[ , c(names_slp, "prod")] expect_equivalent(our_slp, vars_slp, tolerance = tolerant) # Differing base season has no effect on slope coefficients! } # roots of the companion matrix # var.2c = vars::VAR(Canada, p = dim_p, type = "const") vars_root = vars::roots(var.2c) R.compani = aux_var2companion(A=do.call("cbind", Acoef(var.2c)), dim_p=var.2c$p) our_root = Mod(eigen( R.compani )$values) expect_equal(our_root, vars_root, tolerance = tolerant) # ordinary information criteria # R.criteria = NULL lag_max = 4 Y = t(Canada)[,-(1:lag_max)] dim_T = ncol(Y) for(p in 1:lag_max){ idx = (lag_max+1-p):nrow(Canada) Z = aux_stack(t(Canada)[ ,idx], dim_p=p, D=t(rep(1, dim_T))) B = Y%*%t(Z) %*% solve(Z%*%t(Z)) e = Y - B%*%Z Ucov = e%*%t(e)/dim_T IC_p = aux_MIC(Ucov, COEF=B, dim_T, lambda_H0=NULL) R.criteria = cbind(R.criteria, IC_p) rm(idx, Z, B, e, Ucov, IC_p) } R.selection = apply(R.criteria, MARGIN=1, function(x) which.min(x)) ours = list(selection=R.selection, criteria=R.criteria) theirs = vars::VARselect(y=Canada, lag.max=4, type="const") expect_equivalent(ours$selection, theirs$selection, tolerance = tolerant) expect_equivalent(ours$criteria, theirs$criteria, tolerance = tolerant) }) test_that("aux_vec2vma() provides the correct long-run multiplier matrix XI at different cointegration ranks", { # 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) x2 = urca::ca.jo(data_i, ecdet="trend", type="trace", K=2, spec="transitory") for(dim_r in 1:3){ # define # x = VECM(y=data_i, dim_r=dim_r, dim_p=2, type="Case4") dim_K = x$dim_K # number of endogenous variables dim_p = x$dim_p # number of lags dim_T = x$dim_T # number of observations without presample # calculate long-run multiplier matrix # alpha_oc = MASS::Null(x$VECM$alpha) beta_oc = MASS::Null(x$beta[1:dim_K, ]) our_XI = aux_vec2vma(GAMMA=x$VECM$GAMMA, alpha_oc=alpha_oc, beta_oc=beta_oc, dim_p=dim_p, n.ahead="XI") # define # r = dim_r K <- x2@P P <- x2@lag - 1 IK <- diag(K) vecr <- urca::cajorls(z = x2, r = r) # calculate long-run multiplier matrix, code chunk from vars::SVEC() # Coef.vecr <- coef(vecr$rlm) alpha <- t(Coef.vecr[1:r, ]) ifelse(r == 1, alpha.orth <- Null(t(alpha)), alpha.orth <- Null(alpha)) beta <- vecr$beta[1:K, ] beta.orth <- Null(beta) Gamma.array <- array(c(t(tail(coef(vecr$rlm), K*P))), c(K, K, P)) Gamma <- apply(Gamma.array, c(1, 2), sum) vars_Xi <- beta.orth %*% solve(t(alpha.orth) %*% (IK - Gamma) %*% beta.orth) %*% t(alpha.orth) # check # expect_equivalent(our_XI, vars_Xi) } })