R/ts_arimax.R

In tspredit: Time Series Prediction with Integrated Tuning

#'@title ARIMAX
#'@description Create a target-centered multivariate regressor based on ARIMA
#' with external regressors.
#'
#'@details
#' `ts_arimax()` is the singular multivariate counterpart of `ts_arima()`.
#'
#' The model keeps one target variable `y` as the main forecasting objective and
#' uses the aligned auxiliary variables `x1, ..., xn` as regressors through the
#' `xreg` mechanism of the `forecast` package.
#'
#' This is the natural choice when the user thinks in the following way:
#' - there is one main series `y`
#' - the remaining variables help explain or anticipate `y`
#' - the primary output is the future path of `y`
#'
#' In other words, `ts_arimax()` is a target-centered multivariate model, not a
#' symmetric system model like `ts_var()`.
#'
#' For multi-step forecasting, future auxiliary values can be supplied directly
#' or generated by the auxiliary univariate models stored in `models_x`.
#'
#' This makes `ts_arimax()` a natural member of the new `ts_reg_mv` branch:
#' - the target forecast remains central
#' - the aligned multivariate system is still available when `return_all = TRUE`
#' - auxiliary series can be modeled separately when their future path is not
#'   known beforehand
#'
#' The current implementation follows the same philosophy as `ts_arima()`:
#' - if `p`, `d`, and `q` are supplied, fit that specific order
#' - otherwise use `forecast::auto.arima()` on the target with `xreg`
#'
#' This keeps the singular multivariate branch consistent with the existing raw
#' univariate branch of the package.
#'
#'@param models_x Optional named list with one univariate model per auxiliary
#' variable.
#'@param p Optional integer autoregressive order. Leave `NULL` to let
#' `auto.arima()` choose it.
#'@param d Optional integer differencing order. Leave `NULL` to let
#' `auto.arima()` choose it.
#'@param q Optional integer moving-average order. Leave `NULL` to let
#' `auto.arima()` choose it.
#'@return A `ts_arimax` object inheriting from `ts_reg_mv`.
#'@references
#' - Box GEP, Jenkins GM, Reinsel GC, Ljung GM (2015). Time Series Analysis:
#'   Forecasting and Control. Wiley.
#' - Hyndman RJ, Athanasopoulos G (2021). Forecasting: Principles and Practice.
#'   Third Edition. OTexts. https://otexts.com/fpp3/
#'@examples
#'data(tsd)
#'x1 <- c(tsd$y[-1], tail(tsd$y, 1))
#'x2 <- stats::filter(tsd$y, rep(1/3, 3), sides = 1)
#'x2[is.na(x2)] <- tsd$y[is.na(x2)]
#'
#'mv <- ts_data_mv(data.frame(y = tsd$y, x1 = x1, x2 = as.numeric(x2)), y = "y")
#'samp <- ts_sample(mv, test_size = 5)
#'
#'model <- ts_arimax(models_x = list(x1 = ts_arima(), x2 = ts_arima()))
#'model <- daltoolbox::fit(model, samp$train)
#'predict(model, steps_ahead = 5)
#'@export
ts_arimax <- function(models_x = NULL, p = NULL, d = NULL, q = NULL) {
  obj <- ts_reg_mv(models_x = models_x)
  obj$p <- p
  obj$d <- d
  obj$q <- q
  obj$manual_order <- !is.null(p) && !is.null(d) && !is.null(q)
  obj$include_mean <- isTRUE(obj$manual_order && d == 0)
  obj$include_drift <- isTRUE(obj$manual_order && d == 1)
  class(obj) <- append("ts_arimax", class(obj))
  obj
}

#'@importFrom forecast auto.arima
#'@importFrom forecast Arima
#'@exportS3Method fit ts_arimax
#'@inheritParams do_fit
#'@return A fitted `ts_arimax` object.
#'@noRd
fit.ts_arimax <- function(obj, x, y = NULL, ...) {
  x <- reg_mv_validate_data(x)
  obj <- reg_mv_set_metadata(obj, x)
  obj <- reg_mv_validate_models_x(obj, x)
  obj <- reg_mv_fit_aux_models(obj, as.data.frame(x))

  data <- as.data.frame(x)
  y_series <- data[[obj$y_name]]
  xreg <- reg_mv_numeric_matrix(data, obj$x_names, context = "training auxiliary data")

  if (isTRUE(obj$manual_order)) {
    obj$model <- forecast::Arima(
      y_series,
      order = c(obj$p, obj$d, obj$q),
      xreg = xreg,
      include.mean = obj$include_mean,
      include.drift = obj$include_drift
    )
  } else {
    obj$model <- forecast::auto.arima(
      y_series,
      xreg = xreg,
      allowdrift = TRUE,
      allowmean = TRUE
    )
    order <- obj$model$arma[c(1, 6, 2, 3, 7, 4, 5)]
    obj$p <- order[1]
    obj$d <- order[2]
    obj$q <- order[3]
    obj$include_mean <- isTRUE(obj$d == 0 && !is.element("intercept", names(obj$model$coef)))
    obj$include_drift <- (NCOL(obj$model$xreg) == 1) && is.element("drift", names(obj$model$coef))
  }

  attr(obj, "params") <- list(p = obj$p, d = obj$d, q = obj$q)
  obj
}

#'@importFrom forecast forecast
#'@exportS3Method predict ts_arimax
#'@inheritParams do_predict
#'@param steps_ahead Integer. Forecast horizon.
#'@param return_all Logical. Ignored for compatibility. The method always
#' returns the target forecast as a vector with the full multivariate forecast
#' attached as attributes.
#'@return Numeric vector with target forecasts. The full multivariate forecast
#' is attached as attributes.
#'@noRd
predict.ts_arimax <- function(object, x = NULL, steps_ahead = 1, return_all = FALSE, ...) {
  steps_ahead <- as.integer(steps_ahead)
  future_x <- reg_mv_forecast_aux(object, x = x, steps_ahead = steps_ahead)
  forecast_h <- steps_ahead
  forecast_x <- future_x

  # forecast::forecast.Arima() can fail for h = 1 with xreg even when the
  # same call works for h > 1. Requesting one extra step avoids that internal
  # edge case; only the first step is used here.
  if (steps_ahead == 1L) {
    forecast_h <- 2L
    forecast_x <- rbind(future_x, future_x[1, , drop = FALSE])
    rownames(forecast_x) <- NULL
  }

  prediction_y <- as.vector(
    forecast::forecast(
      object$model,
      h = forecast_h,
      xreg = reg_mv_numeric_matrix(forecast_x, object$x_names, context = "future auxiliary data")
    )$mean
  )[seq_len(steps_ahead)]

  object$history <- reg_mv_update_history(object, future_x, prediction_y)

  prediction_x <- as.list(future_x)
  mv_compose_prediction(object, prediction_y, prediction_x)
}