R/04GVARtests.R

In GVARX: Perform Global Vector Autoregression Estimation and Inference

.gvar.normality <- function(x, multivariate.only = TRUE){
  if(!inherits(x,c("varest","vec2var"))) {
    stop("\nPlease provide an object of class 'varest', generated by 'var()', or an object of class 'vec2var' generated by 'vec2var()'.\n")
  }
  obj.name <- deparse(substitute(x))
  K <- x$K
  obs <- x$obs
  resid <- resid(x)
  resids <- scale(resid, scale=FALSE)
  ## Jarque Bera Test (multivariate)
  jbm.resids <- .jb.multi(resids, obs = obs, K = K, obj.name = obj.name)
  if(multivariate.only){
    result <- list(resid = resid, jb.mul = jbm.resids)
  } else {
    ## Jarque Bera Test (univariate)
    jbu.resids <- apply(resids, 2, function(x) .jb.uni(x, obs = obs))
    for(i in 1 : K)
      jbu.resids[[i]][5] <- paste("Residual of", colnames(resids)[i], "equation")
    result <- list(resid = resid, jb.uni = jbu.resids, jb.mul = jbm.resids)
  }
  class(result) <- "varcheck"
  return(result)
}

.normality <- function(x, multivariate.only = TRUE){
  .Deprecated("normality.test", package = "vars", msg = "Function 'normality' is deprecated; use 'normality.test' instead.\nSee help(\"vars-deprecated\") and help(\"normality-deprecated\") for more information.")
  normality.test(x = x, multivariate.only = multivariate.only)
}

.gvar.serial <- function(x, lags.pt = 16, lags.bg = 5, type = c("PT.asymptotic", "PT.adjusted", "BG", "ES")){
  if(!inherits(x,c("varest","vec2var"))) {
    stop("\nPlease provide an object of class 'varest', generated by 'var()', or an object of class 'vec2var' generated by 'vec2var()'.\n")
  }
  obj.name <- deparse(substitute(x))
  type <- match.arg(type)
  K <- x$K
  obs <- x$obs
  resids <- resid(x)
  if((type == "PT.asymptotic") || (type == "PT.adjusted")){
    lags.pt <- abs(as.integer(lags.pt))
    ptm <- .pt.multi(x, K = K, obs = obs, lags.pt = lags.pt, obj.name = obj.name, resids = resids)
    ifelse(type == "PT.asymptotic", test <- ptm[[1]], test <- ptm[[2]])
  } else {
    lags.bg <- abs(as.integer(lags.bg))
    bgm <- .bgserial(x, K = K, obs = obs, lags.bg = lags.bg, obj.name = obj.name, resids = resids)
    ifelse(type == "BG", test <- bgm[[1]], test <- bgm[[2]])
  }
  result <- list(resid = resids, serial = test)
  class(result) <- "varcheck"
  return(result)
}

.serial <- function(x, lags.pt = 16, lags.bg = 5, type = c("PT.asymptotic", "PT.adjusted", "BG", "ES")){
  .Deprecated("serial.test", package = "vars", msg = "Function 'serial' is deprecated; use 'serial.test' instead.\nSee help(\"vars-deprecated\") and help(\"serial-deprecated\") for more information.")
  serial.test(x = x, lags.pt = lags.pt, lags.bg = lags.bg, type = type)
}

.gvar.stability <- function(x, type = c("OLS-CUSUM", "Rec-CUSUM", "Rec-MOSUM", "OLS-MOSUM", "RE", "ME", "Score-CUSUM", "Score-MOSUM", "fluctuation"), h = 0.15, dynamic = FALSE, rescale = TRUE){
  if(!inherits(x,c("varest"))){
    stop("\nPlease provide an object of class 'varest', generated by 'var()'.\n")
  }
  type <- match.arg(type)
  K <- x$K

  stability <- list()
  endog <- colnames(x$datamat)[1 : K]
  for(i in 1 : K){
    formula <- formula(x$varresult[[i]])
    data <- x$varresult[[i]]$model
    stability[[endog[i]]] <- strucchange::efp(formula = formula, data = data, type = type, h = h, dynamic = dynamic, rescale = rescale)
  }
  result <- list(stability = stability, names = endog, K = K)
  class(result) <- "varstabil"
  return(result)
}

.gvecm.fevd <-function(x, n.ahead=10, ...)
{ vars::fevd(x,n.ahead=10) }

.gvar.fevd <- function(x, n.ahead = 10, ...){
    UseMethod("fevd", x)
  }


.fevd.varest <- function(x, n.ahead=10, ...){
  if(!inherits(x,c("varest"))) {
    stop("\nPlease provide an object of class 'varest', generated by 'VAR()'.\n")
  }
  n.ahead <- abs(as.integer(n.ahead))
  K <- x$K
  p <- x$p
  ynames <- colnames(x$datamat[, 1 : K])
  msey <- .fecov(x, n.ahead = n.ahead)
  Psi <- Psi(x, nstep = n.ahead)
  mse <- matrix(NA, nrow = n.ahead, ncol = K)
  Omega <- array(0, dim = c(n.ahead, K, K))
  for(i in 1 : n.ahead){
    mse[i, ] <- diag(msey[, , i])
    temp <- matrix(0, K, K)
    for(l in 1 : K){
      for(m in 1 : K){
        for(j in 1 : i){
          temp[l, m] <- temp[l, m] + Psi[l , m, j]^2
        }
      }
    }
    temp <- temp / mse[i, ]
    for(j in 1 : K){
      Omega[i, ,j] <- temp[j, ]
    }
  }
  result <- list()
  for(i in 1 : K){
    result[[i]] <- matrix(Omega[, , i], nrow = n.ahead, ncol = K)
    colnames(result[[i]]) <- ynames
  }
  names(result) <- ynames
  class(result) <- "varfevd"
  return(result)
}





.gvar.arch <- function(x, lags.single = 16, lags.multi = 5, multivariate.only = TRUE){
  if(!inherits(x,c("varest","vec2var"))) {
    stop("\nPlease provide an object of class 'varest', generated by 'var()', or an object of class 'vec2var' generated by 'vec2var()'.\n")
  }
  obj.name <- deparse(substitute(x))
  lags.single <- abs(as.integer(lags.single))
  lags.multi <- abs(as.integer(lags.multi))
  K <- x$K
  obs <- x$obs
  resid <- resid(x)
  resids <- scale(resid)
  ## ARCH test (multivariate)
  archm.resids <- .arch.multi(resids, lags.multi = lags.multi, K = K, obs = obs, obj.name = obj.name)
  if(multivariate.only){
    result <- list(resid=resid, arch.mul = archm.resids)
  } else {
    ## ARCH test (univariate)
    archs.resids <- apply(resids, 2, function(x) .arch.uni(x, lags.single = lags.single))
    for(i in 1 : K)
      archs.resids[[i]][5] <- paste("Residual of", colnames(resids)[i], "equation")
    result <- list(resid=resid, arch.uni = archs.resids, arch.mul = archm.resids)
  }
  class(result) <- "varcheck"
  return(result)
}


.arch <- function(x, lags.single = 16, lags.multi = 5, multivariate.only = TRUE){
  .Deprecated("arch.test", package = "vars", msg = "Function 'arch' is deprecated; use 'arch.test' instead.\nSee help(\"vars-deprecated\") and help(\"arch-deprecated\") for more information.")
  arch.test(x = x, lags.single = lags.single, lags.multi = lags.multi, multivariate.only = multivariate.only)
}

##
## univariate ARCH test
##
.arch.uni <- function(x, lags.single){
  lags.single <- lags.single + 1
  mat <- embed(scale(x)^2, lags.single)
  arch.lm <- summary(lm(mat[, 1] ~ mat[, -1]))
  STATISTIC <- arch.lm$r.squared*length(resid(arch.lm))
  names(STATISTIC) <- "Chi-squared"
  PARAMETER <- lags.single - 1
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "ARCH test (univariate)"
  result <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = deparse(substitute(x)))
  class(result) <- "htest"
  return(result)
}
##
## multivariate ARCH test
##
.arch.multi <- function(x, lags.multi, obj.name, K, obs){
  col.arch.df <- 0.5 * K * (K + 1)
  arch.df <- matrix(NA, nrow = obs, ncol = col.arch.df)
  for( i in 1 : obs){
    temp <- outer(x[i,], x[i,])
    arch.df[i,] <- temp[lower.tri(temp, diag=TRUE)]
  }
  lags.multi <- lags.multi + 1
  arch.df <- embed(arch.df, lags.multi)
  archm.lm0 <- lm(arch.df[ , 1:col.arch.df] ~ 1)
  archm.lm0.resids <- resid(archm.lm0)
  omega0 <- cov(archm.lm0.resids)
  archm.lm1 <- lm(arch.df[ , 1 : col.arch.df] ~ arch.df[ , -(1 : col.arch.df)])
  archm.lm1.resids <- resid(archm.lm1)
  omega1 <- cov(archm.lm1.resids)
  R2m <- 1 - (2 / (K * (K + 1))) * sum(diag(omega1 %*% solve(omega0)))
  n <- nrow(archm.lm1.resids)
  STATISTIC <- 0.5 * n * K * (K+1) * R2m
  names(STATISTIC) <- "Chi-squared"
  lags.multi <- lags.multi - 1
  PARAMETER <- lags.multi * K^2 * (K + 1)^2 / 4
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "ARCH (multivariate)"
  result <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(result) <- "htest"
  return(result)
}
##
## univariate normality test
##
.jb.uni <- function(x, obs){
  x <- as.vector(x)
  m1 <- sum(x) / obs
  m2 <- sum((x - m1)^2) / obs
  m3 <- sum((x - m1)^3) / obs
  m4 <- sum((x - m1)^4) / obs
  b1 <- (m3 / m2^(3 / 2))^2
  b2 <- (m4/m2^2)
  STATISTIC <- obs * b1 / 6 + obs * (b2 - 3)^2 / 24
  names(STATISTIC) <- "Chi-squared"
  PARAMETER <- 2
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = 2)
  METHOD <- "JB-Test (univariate)"
  result <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = deparse(substitute(x)))
  class(result) <- "htest"
  return(result)
}
##
## multivariate normality test
##
.jb.multi <- function(x, obs, K, obj.name){
  P <- chol(crossprod(x) / obs)
  resids.std <- x %*% solve(P)
  b1 <- apply(resids.std, 2, function(x) sum(x^3) / obs)
  b2 <- apply(resids.std, 2, function(x) sum(x^4) / obs)
  s3 <- obs * t(b1) %*% b1 / 6
  s4 <- obs * t(b2 - rep(3, K)) %*% (b2 - rep(3, K)) / 24
  STATISTIC <- s3 + s4
  names(STATISTIC) <- "Chi-squared"
  PARAMETER <- 2 * K
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "JB-Test (multivariate)"
  result1 <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(result1) <- "htest"
  STATISTIC <- s3
  names(STATISTIC) <- "Chi-squared"
  PARAMETER <- K
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "Skewness only (multivariate)"
  result2 <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(result2) <- "htest"
  STATISTIC <- s4
  names(STATISTIC) <- "Chi-squared"
  PARAMETER <- K
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "Kurtosis only (multivariate)"
  result3 <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(result3) <- "htest"
  result <- list(JB = result1, Skewness = result2, Kurtosis = result3)
  return(result)
}
##
## Convenience function for computing lagged x
##

.matlag1 <- function(x, lag = 1){
  totcols <- ncol(x)
  nas <- matrix(NA, nrow = lag, ncol = totcols)
  x <- rbind(nas, x)
  totrows <- nrow(x)
  x <- x[-c((totrows - lag + 1) : totrows), ]
  return(x)
}
##
## Multivariate Portmanteau Statistic
##
.pt.multi <- function(x, K, obs, lags.pt, obj.name, resids){
  C0 <- crossprod(resids) / obs
  C0inv <- solve(C0)
  tracesum <- rep(NA, lags.pt)
  for(i in 1 : lags.pt){
    Ut.minus.i <- .matlag1(resids, lag = i)[-c(1 : i), ]
    Ut <- resids[-c(1 : i), ]
    Ci <- crossprod(Ut, Ut.minus.i) / obs
    tracesum[i] <- sum(diag(t(Ci) %*% C0inv %*% Ci %*% C0inv))
  }
  vec.adj <- obs - (1 : lags.pt)
  Qh <- obs * sum(tracesum)
  Qh.star <- obs^2 * sum(tracesum / vec.adj)
  nstar <- K^2 * x$p
  ## htest objects for Qh and Qh.star
  STATISTIC <- Qh
  names(STATISTIC) <- "Chi-squared"
  if(identical(class(x), "varest")){
    PARAMETER <- (K^2 * lags.pt - nstar)
  } else {
    PARAMETER <- (K^2 * lags.pt - nstar + x$K)
  }
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "Portmanteau Test (asymptotic)"
  PT1 <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(PT1) <- "htest"
  STATISTIC <- Qh.star
  names(STATISTIC) <- "Chi-squared"
  if(identical(class(x), "varest")){
    PARAMETER <- (K^2 * lags.pt - nstar)
  } else {
    PARAMETER <- (K^2 * lags.pt - nstar + x$K)
  }
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "Portmanteau Test (adjusted)"
  PT2 <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(PT2) <- "htest"
  result <- list(PT1 = PT1, PT2 = PT2)
  return(result)
}
##
## Breusch-Godfrey and Edgerton-Shukur Test
##

.bgserial <- function(x, K, obs, lags.bg, obj.name, resids){
  ylagged <- as.matrix(x$datamat[, -c(1 : K)])
  resids.l <- .matlag2(resids, lag = lags.bg)
  if(is.null(x$restrictions)){
    regressors <- as.matrix(cbind(ylagged, resids.l))
    lm0 <- lm(resids ~ -1 + regressors)
    lm1 <- lm(resids ~ -1 + ylagged)
    sigma.1 <- crossprod(resid(lm1)) / obs
    sigma.0 <- crossprod(resid(lm0)) / obs
  } else {
    resid0 <- matrix(NA, ncol = K, nrow = obs)
    resid1 <- matrix(NA, ncol = K, nrow = obs)
    for(i in 1 : K){
      datares <- data.frame(ylagged[, which(x$restrictions[i, ] == 1)])
      regressors <- data.frame(cbind(datares, resids.l))
      lm0 <- lm(resids[, i] ~ -1 + ., data=regressors)
      lm1 <- lm(resids[, i] ~ -1 + ., data=datares)
      resid0[, i] <- resid(lm0)
      resid1[, i] <- resid(lm1)
      sigma.0 <- crossprod(resid0) / obs
      sigma.1 <- crossprod(resid1) / obs
    }
  }
  LMh.stat <- obs * (K - sum(diag(crossprod(solve(sigma.1), sigma.0))))
  STATISTIC <- LMh.stat
  names(STATISTIC) <- "Chi-squared"
  PARAMETER <- lags.bg * K^2
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
  METHOD <- "Breusch-Godfrey LM test"
  LMh <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(LMh) <- "htest"
  ## small sample correction of Edgerton Shukur
  R2r <- 1 - det(sigma.0) / det(sigma.1)
  m <- K * lags.bg
  q <- 0.5 * K * m - 1
  n <- ncol(x$datamat) - K
  N <- obs - n - m - 0.5 * (K - m + 1)
  r <- sqrt((K^2 * m^2 - 4)/(K^2 + m^2 - 5))
  LMFh.stat <- (1 - (1 - R2r)^(1 / r))/(1 - R2r)^(1 / r) * (N * r - q) / (K * m)
  STATISTIC <- LMFh.stat
  names(STATISTIC) <- "F statistic"
  PARAMETER1 <- lags.bg * K^2
  names(PARAMETER1) <- "df1"
  PARAMETER2 <- floor(N * r - q)
  names(PARAMETER2) <- "df2"
  PVAL <-   1 - pf(LMFh.stat, PARAMETER1, PARAMETER2)
  METHOD <- "Edgerton-Shukur F test"
  LMFh <- list(statistic = STATISTIC, parameter = c(PARAMETER1, PARAMETER2), p.value = PVAL, method = METHOD, data.name = paste("Residuals of VAR object", obj.name))
  class(LMFh) <- "htest"
  return(list(LMh = LMh, LMFh = LMFh))
}
##
## Duplication matrix
##
.duplicate <- function(n){
  D <- matrix(0, nrow = n^2, ncol = n * (n + 1) / 2)
  count <- 0
  for(j in 1 : n){
    D[(j - 1) * n + j, count + j] <- 1
    if((j + 1) <= n){
      for(i in (j + 1):n){
        D[(j - 1) * n + i, count + i] <- 1
        D[(i - 1) * n + j, count + i] <- 1
      }
    }
    count <- count + n - j
  }
  return(D)
}
##
## Convenience function for computing lagged residuals
##
.matlag2 <- function(x, lag = 1){
  K <- ncol(x)
  obs <- nrow(x)
  zeromat <- matrix(0, nrow = obs, ncol = K * lag)
  idx1 <- seq(1, K * lag, K)
  idx2 <- seq(K, K * lag, K)
  for(i in 1:lag){
    lag <- i + 1
    res.tmp <- embed(x, lag)[, -c(1 : (K * i))]
    zeromat[-c(1 : i), idx1[i] : idx2[i]] <- res.tmp
  }
  resids.l <- zeromat
  return(resids.l)
}



##
## Forecast variance-covariance matrix (VAR)
##
.fecov <- function(x, n.ahead) {
  n.par<-sapply(x$varresult, function(x) summary(x)$df[2])
  sigma.u <- crossprod(resid(x))/n.par
  Sigma.yh <- array(NA, dim = c(x$K, x$K, n.ahead))
  Sigma.yh[, , 1] <- sigma.u
  Phi <- Phi(x, nstep = n.ahead)
  if (n.ahead > 1) {
    for (i in 2:n.ahead) {
      temp <- matrix(0, nrow = x$K, ncol = x$K)
      for (j in 2:i) {
        temp <- temp + Phi[, , j] %*% sigma.u %*% t(Phi[, , j])
      }
      Sigma.yh[, , i] <- temp + Sigma.yh[, , 1]
    }
  }
  return(Sigma.yh)
}
##
## Forecast variance-covariance matrix (SVAR)
##
.fecovsvar <- function(x, n.ahead) {
  Sigma.yh <- array(NA, dim = c(x$var$K, x$var$K, n.ahead))
  Phi <- Phi(x, nstep = n.ahead)
  Sigma.yh[, , 1] <- Phi[, , 1]%*%t(Phi[, , 1])
  if (n.ahead > 1) {
    for (i in 2:n.ahead) {
      temp <- matrix(0, nrow = x$var$K, ncol = x$var$K)
      for (j in 2:i) {
        temp <- temp + Phi[, , j]%*%t(Phi[, , j])
      }
      Sigma.yh[, , i] <- temp + Sigma.yh[, , 1]
    }
  }
  return(Sigma.yh)
}
##
## Forecast variance-covariance matrix (vec2var)
##
.fecovvec2var <- function(x, n.ahead) {
  sigma.u <- crossprod(resid(x))/x$obs
  Sigma.yh <- array(NA, dim = c(x$K, x$K, n.ahead))
  Sigma.yh[, , 1] <- sigma.u
  Phi <- Phi(x, nstep = n.ahead)
  if (n.ahead > 1) {
    for (i in 2:n.ahead) {
      temp <- matrix(0, nrow = x$K, ncol = x$K)
      for (j in 2:i) {
        temp <- temp + Phi[, , j] %*% sigma.u %*% t(Phi[, , j])
      }
      Sigma.yh[, , i] <- temp + Sigma.yh[, , 1]
    }
  }
  return(Sigma.yh)
}
##
## Forecast variance-covariance matrix (SVEC)
##
.fecovsvec <- function(x, n.ahead, K) {
  Sigma.yh <- array(NA, dim = c(K, K, n.ahead))
  Phi <- Phi(x, nstep = n.ahead)
  Sigma.yh[, , 1] <- Phi[, , 1]%*%t(Phi[, , 1])
  if (n.ahead > 1) {
    for (i in 2:n.ahead) {
      temp <- matrix(0, nrow = K, ncol = K)
      for (j in 2:i) {
        temp <- temp + Phi[, , j]%*%t(Phi[, , j])
      }
      Sigma.yh[, , i] <- temp + Sigma.yh[, , 1]
    }
  }
  return(Sigma.yh)
}