R/diagTests.R

In sstvars: Toolkit for Reduced Form and Structural Smooth Transition Vector Autoregressive Models

#' @title Perform adjusted Portmanteau test for a STVAR model
#'
#' @description \code{Portmanteau_test} performs adjusted Portmanteau test for remaining autocorrelation
#'   (or heteroskedasticity) in the residuals of a STVAR model.
#'
#' @param stvar an object of class \code{'stvar'} generated by \code{fitSTVAR} or \code{STVAR}.
#' @param nlags a strictly positive integer specifying the number of lags to be tested.
#' @param which_test should test for remaining autocorrelation or heteroskedasticity be calculated?
#' @details The implemented adjusted Portmanteau test is based on Lütkepohl (2005), Section 4.4.3.
#'   When testing for remaining heteroskedasticity, the Portmanteau test is applied to squared
#'   standardized residuals that are centered to have zero mean. Note that the validity of the
#'   heteroskedasticity test requires that the residuals are not autocorrelated.
#' @return A list with class "hypotest" containing the test results and arguments used to calculate the test.
#' @seealso \code{\link{LR_test}}, \code{\link{Rao_test}}, \code{\link{fitSTVAR}}, \code{\link{STVAR}},
#'   \code{\link{diagnostic_plot}}, \code{\link{profile_logliks}},
#' @references
#'   \itemize{
#'    \item Lütkepohl H. 2005. New Introduction to Multiple Time Series Analysis,
#'            \emph{Springer}.
#'   }
#' @examples
#' # Gaussian STVAR p=2, M=2, model with weighted relative stationary densities
#' # of the regimes as the transition weight function:
#' theta_222relg <- c(0.357, 0.107, 0.356, 0.086, 0.14, 0.035, -0.165, 0.387, 0.452,
#'  0.013, 0.228, 0.336, 0.239, 0.024, -0.021, 0.708, 0.063, 0.027, 0.009, 0.197,
#'  0.206, 0.005, 0.026, 1.092, -0.009, 0.116, 0.592)
#' mod222relg <- STVAR(data=gdpdef, p=2, M=2, d=2, params=theta_222relg,
#'  weight_function="relative_dens")
#'
#' # Test for remaining autocorrelation taking into account the first 20 lags:
#' Portmanteau_test(mod222relg, nlags=20)
#'
#' # Test for remaining heteroskedasticity taking into account the first 20 lags:
#' Portmanteau_test(mod222relg, nlags=20, which_test="het.sked")
#'
#' # Two-regime Student's t Threhold VAR p=3 model with the first lag of the second
#' # variable as the switching variable:
#' theta_322thres <- c(0.527, 0.039, 1.922, 0.154, 0.284, 0.053, 0.033, 0.453, 0.291,
#'  0.024, -0.108, 0.153, -0.108, 0.003, -0.128, 0.219, 0.195, -0.03, -0.893, 0.686,
#'  0.047, 0.016, 0.524, 0.068, -0.025, 0.044, -0.435, 0.119, 0.359, 0.002, 0.038,
#'  1.252, -0.041, 0.151, 1.196, 12.312)
#' mod322thres <- STVAR(data=gdpdef, p=3, M=2, d=2, params=theta_322thres,
#'  weight_function="threshold", weightfun_pars=c(2, 1), cond_dist="Student")
#'
#' # Test for remaining autocorrelation taking into account the first 25 lags:
#' Portmanteau_test(mod322thres, nlags=25)
#'
#' # Test for remaining heteroskedasticity taking into account the first 25 lags:
#' Portmanteau_test(mod322thres, nlags=25, which_test="het.sked")
#' @export

Portmanteau_test <- function(stvar, nlags=20, which_test=c("autocorr", "het.sked")) {
  which_test <- match.arg(which_test)
  if(is.null(stvar$residuals_raw) && which_test == "autocorr") {
    stop("The model needs to contain residuals_raw for the autocorrelation test")
  } else if(is.null(stvar$residuals_std) && which_test == "het.sked") {
    stop("The model needs to contain residuals_std for the heteroskedasticity test")
  }
  if(length(nlags) != 1 || nlags < 1 || nlags%%1 != 0) {
    stop("The argument nlags should be a strictly positive interger")
  }
  p <- stvar$model$p
  M <- stvar$model$M
  d <- stvar$model$d
  if(nlags <= p) {
    stop("The number of lags 'nlags' should be larger than the autoregressive order p")
  }

  if(which_test == "autocorr") {
    U <- t(stvar$residuals_raw)
  } else { # which_test == het.sked
    U <- t(stvar$residuals_std^2) # Test applied to squared standardized residuals
    U <- U - rowMeans(U) # Standardize zero mean
  }
  T_obs <- ncol(U) # U = (d x T_obs)

  # Function that calculates the matrix C_i given in Lütkepohl (2005), Equation (4.4.1).
  get_Ci <- function(i) { # T_obs and U taken from the parent environment
    U%*%tcrossprod(create_Fi_matrix(i=i, T_obs=T_obs), U)/T_obs
  }

  # Calculate the test statistic, Lütkepohl (2005), Equation (4.4.23)
  inv_C_0 <- solve(get_Ci(i=0))
  tmp_q <- numeric(nlags)
  for(i1 in 1:nlags) {
    C_i <- get_Ci(i=i1)
    tmp_q[i1] <- sum(diag(crossprod(C_i, inv_C_0%*%C_i%*%inv_C_0)))/(T_obs - i1)
  }
  test_stat <- T_obs^2*sum(tmp_q)

  # Calculate the p-value
  df <- d^2*(nlags - p)
  p_value <- pchisq(test_stat, df=df, lower.tail=FALSE)

  # Return
  structure(list(stvar=stvar,
                 df=df,
                 test_stat=test_stat,
                 p_value=p_value,
                 type="Portmanteau test",
                 which_test=which_test),
             class="hypotest")
}