R/rolling-oos.R

In bivarhr: Bivariate Hurdle Regression with Bayesian Model Averaging

#' Rolling Out-of-Sample Forecast Evaluation
#'
#' Computes rolling out-of-sample (OOS) forecast accuracy for the
#' selected bivariate hurdle model by repeatedly truncating the sample
#' at different cut points \code{Tcut}, generating multi-step-ahead
#' predictive distributions, and summarizing them via RMSE for \code{I}
#' and \code{C}.
#'
#' @param best_fit A fitted model object as returned by \code{fit_one()},
#'   containing at least:
#'   \itemize{
#'     \item \code{$fit}: CmdStanR fit object with posterior draws.
#'     \item \code{$des}: design matrices used by the model.
#'     \item \code{$k}: lag order used in the fit.
#'   }
#'   This object is passed directly to \code{\link{predict_multistep}()}.
#' @param DT A \code{data.frame} or \code{data.table} containing the
#'   original time series and covariates used to fit the model, including
#'   at least columns \code{I} and \code{C}.
#' @param h Integer; maximum forecast horizon (number of steps ahead)
#'   requested at each cut. For a given \code{Tcut}, the effective horizon
#'   is \code{min(h, nrow(DT) - Tcut)}.
#' @param cuts Numeric vector of time indices (training end points) at
#'   which to perform the rolling evaluation. By default, a grid of five
#'   equally spaced cut points between 60\% and 90\% of the sample size is
#'   used:
#'   \code{seq(round(0.6 * nrow(DT)), round(0.9 * nrow(DT)), length.out = 5)}.
#' @param dir_csv Character scalar or \code{NULL}; directory where the
#'   rolling-OOS table (\code{"rolling_oos.csv"}) is written. If
#'   \code{NULL} (default), nothing is written to disk.
#'
#' @details
#' For each \code{Tcut} in \code{cuts}, the function:
#' \enumerate{
#'   \item Calls \code{\link{predict_multistep}()} with
#'         \code{fit_obj = best_fit}, the full \code{DT}, lag
#'         \code{k = best_fit$k}, and horizon
#'         \code{h_eff = min(h, nrow(DT) - Tcut)} to obtain posterior
#'         predictive paths \code{pred_I} and \code{pred_C}.
#'   \item Computes the posterior-mean forecast for each step
#'         (\code{mI}, \code{mC}) as the column means of \code{pred_I}
#'         and \code{pred_C}.
#'   \item Extracts the realized outcomes
#'         \code{yI = I[(Tcut + 1):(Tcut + h_eff)]} and analogously for
#'         \code{yC}.
#'   \item Computes RMSE for each series:
#'         \code{RMSE_I = sqrt(mean((yI - mI)^2))},
#'         \code{RMSE_C = sqrt(mean((yC - mC)^2))}.
#' }
#'
#' Progress is reported via \pkg{progressr}. When \code{dir_csv} is
#' supplied, the resulting table is written as \code{"rolling_oos.csv"} in
#' that directory.
#'
#' @return A \code{data.frame} with one row per \code{Tcut} and columns:
#' \itemize{
#'   \item \code{Tcut}: training end index.
#'   \item \code{RMSE_I}: rolling OOS RMSE for series \code{I}.
#'   \item \code{RMSE_C}: rolling OOS RMSE for series \code{C}.
#' }
#'
#' @examples
#' \donttest{
#' # Minimal synthetic example illustrating the expected data structure:
#' set.seed(123)
#' DT <- data.frame(
#'   id = rep(1:10, each = 2),
#'   t  = rep(1:2, times = 10),
#'   I  = rpois(20, lambda = 0.5),
#'   C  = rpois(20, lambda = 1.0)
#' )
#'
#' # Directory for CSV output (in practice, use a persistent path chosen
#' # by the user):
#' dir_csv <- file.path(tempdir(), "bivarhr_oos_csv")
#'
#' # Typical workflow (commented out to avoid heavy computation and
#' # external dependencies such as CmdStan during R CMD check):
#' #
#' # best_fit <- fit_one(
#' #   data = DT,
#' #   k    = 2,
#' #   spec = "C"
#' # )
#' #
#' # oos_res <- rolling_oos(
#' #   fit     = best_fit,
#' #   data    = DT,
#' #   h       = 6,
#' #   dir_csv = dir_csv
#' # )
#' # print(oos_res)
#' }
#'
#' @export

rolling_oos <- function(best_fit, DT, h = 5,
                        cuts = seq(round(0.6 * nrow(DT)),
                                   round(0.9 * nrow(DT)), length.out = 5),
                        dir_csv = NULL) {
  res <- list()
  
  progressr::with_progress({
    p <- progressr::progressor(steps = length(cuts))
    for (Tcut in cuts) {
      pr <- predict_multistep(
        fit_obj = best_fit,   
        DT = DT,
        k = best_fit$k,
        Tcut = Tcut,
        h = min(h, nrow(DT) - Tcut),
        ndraws = 400
      )
      
      mI <- colMeans(pr$pred_I, na.rm = TRUE)
      mC <- colMeans(pr$pred_C, na.rm = TRUE)
      yI <- DT$I[(Tcut + 1):(Tcut + length(mI))]
      yC <- DT$C[(Tcut + 1):(Tcut + length(mC))]
      rmseI <- sqrt(mean((yI - mI)^2, na.rm = TRUE))
      rmseC <- sqrt(mean((yC - mC)^2, na.rm = TRUE))
      
      res[[as.character(Tcut)]] <- data.frame(
        Tcut = Tcut,
        RMSE_I = rmseI,
        RMSE_C = rmseC
      )
      
      p(message = sprintf("Rolling OOS: cut T=%d", Tcut))
    }
  })
  
  tab <- dplyr::bind_rows(res)
  if (!is.null(dir_csv)) {
    if (!dir.exists(dir_csv)) dir.create(dir_csv, recursive = TRUE)
    readr::write_csv(tab, file.path(dir_csv, "rolling_oos.csv"))
  }
  tab
}

#' Print summary of FLOOR smoke test (ELPD ranking invariance)
#'
#' Nicely prints a summary of the FLOOR smoke test produced by
#' \code{\link{smoketest_floor_elpd_invariance}}, indicating whether
#' the ELPD-based ranking of models is invariant across different
#' \code{FLOOR} constants and listing the combined results.
#'
#' @param st A list returned by \code{\link{smoketest_floor_elpd_invariance}},
#'   containing at least:
#'   \itemize{
#'     \item \code{same_order}: logical flag indicating whether the ELPD
#'           ranking is identical for all \code{FLOOR} values.
#'     \item \code{combined}: data frame or tibble with columns
#'           \code{FLOOR}, \code{fit_id}, \code{elpd}, \code{elpd_se},
#'           and \code{rank_elpd}, among others.
#'   }
#'
#' @details
#' The function uses \pkg{cli} to print a section header and an info
#' message stating whether the ELPD ranking is invariant across values
#' of \code{FLOOR}. It then arranges the combined table by \code{FLOOR}
#' and decreasing \code{elpd}, selects a subset of columns, and prints
#' it to the console.
#'
#' This is a convenience/reporting helper and does not modify \code{st}.
#'
#' @return Invisibly returns the input object \code{st}, so it can be
#'   used in pipes if desired.
#'
#' @examples
#' \donttest{
#' # 1. Define dummy data inside the example so it runs on CRAN checks
#' st_dummy <- list(
#'   same_order = TRUE,
#'   combined = data.frame(
#'     FLOOR     = rep(c(-1e6, -1e4), each = 2),
#'     fit_id    = rep(c("model_1", "model_2"), 2),
#'     elpd      = c(-100.1, -101.3, -100.1, -101.3),
#'     elpd_se   = c(1.2, 1.3, 1.2, 1.3),
#'     rank_elpd = c(1L, 2L, 1L, 2L)
#'   )
#' )
#'
#' # 2. Run the function
#' print_floor_smoketest(st_dummy)
#' }
#'
#' @export
print_floor_smoketest <- function(st) {
  cli::cli_h2("Smoke-test FLOOR -> ELPD ranking invariance")
  
  cli::cli_alert_info(sprintf(
    "Identical ranking across all FLOORS? %s",
    if (isTRUE(st$same_order)) "YES" else "NO"
  ))
  
  # Ensure your DESCRIPTION has Depends: R (>= 4.1.0) for the pipe |>
  tab <- st$combined |>
    dplyr::arrange(FLOOR, dplyr::desc(elpd)) |>
    dplyr::select(FLOOR, fit_id, elpd, elpd_se, rank_elpd)
  
  print(tab, row.names = FALSE)
  
  invisible(st)
}