R/tvar-star.R

In bivarhr: Bivariate Hurdle Regression with Bayesian Model Averaging

#' Fit VARX model with diagnostics for I and C
#'
#' Estimates a bivariate VAR model for \code{I} and \code{C} with
#' exogenous covariates (VARX), and computes a set of standard
#' diagnostics (stability, serial correlation, normality, ARCH). The
#' fitted model and diagnostics are saved to disk and also returned.
#'
#' @param DT A \code{data.table} (or \code{data.frame}) containing at
#'   least the following columns:
#'   \itemize{
#'     \item \code{I}, \code{C}: endogenous variables for the VAR.
#'     \item \code{EconCycle}, \code{PopDensity}, \code{Epidemics},
#'           \code{Climate}, \code{War}, \code{t_norm}: exogenous
#'           regressors included in the VARX.
#'   }
#' @param p Integer; lag order of the VAR part (number of lags for
#'   \code{I} and \code{C}).
#' @param dir_out Character scalar or \code{NULL}; directory where the
#'   fitted model and diagnostics (\code{"varx_fit.rds"}) are saved. If
#'   \code{NULL} (default), nothing is saved to disk.
#'
#' @details
#' This function requires the \pkg{vars} package (listed under
#' \code{Suggests}); an informative error is raised at call time if it is
#' not installed.
#'
#' The endogenous vector is \eqn{y_t = (I_t, C_t)'} and the exogenous
#' regressors are:
#' \code{EconCycle}, \code{PopDensity}, \code{Epidemics},
#' \code{Climate}, \code{War}, \code{t_norm}. The model is fit using
#' \code{vars::VAR()} with \code{type = "const"} and the exogenous
#' matrix passed via \code{exogen}.
#'
#' After estimation, the following diagnostics from \pkg{vars} are
#' (attempted to be) computed:
#' \itemize{
#'   \item \code{vars::stability(fit, type = "OLS-CUSUM")} for stability.
#'   \item \code{vars::serial.test(fit, lags.pt = 10, type = "PT.asymptotic")}
#'         for serial correlation.
#'   \item \code{vars::normality.test(fit)} for residual normality.
#'   \item \code{vars::arch.test(fit, lags.multi = 5)} for ARCH effects.
#' }
#' Each diagnostic call is wrapped in \code{try()}, so if a diagnostic
#' fails, the corresponding element in the output will contain a
#' \code{"try-error"} instead of stopping the function.
#'
#' When \code{dir_out} is supplied, the result is saved as an RDS file
#' named \code{"varx_fit.rds"} in that directory.
#'
#' @return A list with components:
#' \itemize{
#'   \item \code{fit}: the estimated VAR model (\code{vars} object).
#'   \item \code{stability}: result of \code{vars::stability()} (or
#'         \code{"try-error"} on failure).
#'   \item \code{serial}: result of \code{vars::serial.test()} (or
#'         \code{"try-error"} on failure).
#'   \item \code{normal}: result of \code{vars::normality.test()} (or
#'         \code{"try-error"} on failure).
#'   \item \code{arch}: result of \code{vars::arch.test()} (or
#'         \code{"try-error"} on failure).
#' }
#'
#' @examples
#' \donttest{
#' # This example runs only when 'vars' is installed.
#' if (requireNamespace("vars", quietly = TRUE)) {
#'   DT <- data.table::data.table(
#'     I = rpois(50, lambda = 10),
#'     C = rpois(50, lambda = 8),
#'     EconCycle = rnorm(50),
#'     PopDensity = rnorm(50),
#'     Epidemics = rnorm(50),
#'     Climate = rnorm(50),
#'     War = rnorm(50),
#'     t_norm = seq(-1, 1, length.out = 50)
#'   )
#'
#'   # p = 1 keeps the example fast and stable.
#'   res_varx <- run_varx(DT, p = 1)
#'   if (!inherits(res_varx$fit, "try-error")) {
#'     print(res_varx$fit)
#'   }
#' }
#' }
#'
#' @export

run_varx <- function(DT, p = 2, dir_out = NULL) {
  if (!requireNamespace("vars", quietly = TRUE)) {
    stop("Package 'vars' is required for run_varx(). Please install it.",
         call. = FALSE)
  }
  DT <- data.table::as.data.table(DT)
  y_cols <- c("I", "C")
  x_cols <- c("EconCycle", "PopDensity", "Epidemics", "Climate", "War", "t_norm")

  y <- as.matrix(DT[, ..y_cols])
  exog <- as.matrix(DT[, ..x_cols])
  colnames(exog) <- x_cols

  fit <- vars::VAR(y, p = p, type = "const", exogen = exog)
  st <- try(vars::stability(fit, type = "OLS-CUSUM"), silent = TRUE)
  serial <- try(vars::serial.test(fit, lags.pt = 10, type = "PT.asymptotic"), silent = TRUE)
  normal <- try(vars::normality.test(fit), silent = TRUE)
  arch <- try(vars::arch.test(fit, lags.multi = 5), silent = TRUE)
  res <- list(fit = fit, stability = st, serial = serial, normal = normal, arch = arch)
  if (!is.null(dir_out)) {
    if (!dir.exists(dir_out)) dir.create(dir_out, recursive = TRUE)
    saveRDS(res, file.path(dir_out, "varx_fit.rds"))
  }
  res
}