R/predict.R

In garma: Fitting and Forecasting Gegenbauer ARMA Time Series Models

#' Predict future values.
#'
#' Predict ahead using algorithm of Godet (2009).
#'
#' @param object (garma_model) The garma_model from which to predict the values. This should have been generated by the [garma()]
#'   function.
#' @param n.ahead (int) The number of time periods to predict ahead. Default: 1
#' @param newdata (real vector or matrix) If the original model was fitted with the 'xreg=' option then this will provide the xreg
#'   values for predictions. If this is a vector then its length should be 'n.ahead'; if it is a matrix then it should have
#'   'n.ahead' rows.
#'
#'   It should have columns with the same names as the original xreg matrix.
#' @param ... Other parameters. Ignored.
#' @return A "ts" object containing the requested forecasts.
#' @references
#' Godet, F. Linear prediction of long-range dependent time series, ESAIM: PS (2009) 13 115-134.
#' DOI: https://doi.org/10.1051/ps:2008015
#' @examples
#' data(AirPassengers)
#' ap <- as.numeric(diff(AirPassengers, 12))
#' mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE)
#' predict(mdl, n.ahead = 12)
#' @export
predict.garma_model <- function(object, n.ahead = 1, newdata=NULL, ...) {
  internal_build_weights <- function(len) {
    wgts <- 1
    for (gf in object$model$ggbr_factors) {
      u <- cos(2*pi*gf$freq)
      fctr <- internal_ggbr_coef(len + 2, -gf$fd, u)
      wgts <- pracma::conv(wgts, fctr)
    }
    # Next line multiplies and divides the various polynomials to get psi = theta * delta * ggbr / phi
    # pracma::conv gives polynomial multiplication, and pracma::deconv gives polynomial division.
    # we don't bother with the remainder. For non-ggbr models this may be a mistake.
    wgts <- pracma::conv(phi_vec, wgts)
    if (length(theta_vec) > 1) wgts <- pracma::deconv(wgts, theta_vec)$q

    return(wgts[(len + 1):2])
  }

  ## Start of Function "predict"

  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  # Validations first

  if (n.ahead <= 0 | !is.numeric(n.ahead) | is.na(n.ahead)) {
    rlang::abort("The parameter 'n.ahead' must be an integer g.t. 0.")
    return(invisible(NULL))
  }

  if (is.null(newdata)) {
    if ("lm_reg" %in% names(object$lm_xreg)) {
      # if intercept and drift only then no need for newdata...
      other_xreg_factors <- setdiff(names(object$xreg), c("intercept", "drift"))
      if (length(other_xreg_factors) > 0) {
        rlang::abort(
          paste0("To predict you will need to provide some xreg values for the following factors: ",
                 paste(other_xreg_factors, collapse = ", "), ".")
        )
      }
    }
  } else {
    # newdata is not null
    if (is.null(object$xreg)) {
      rlang::abort("You have provided the 'newdata' parameter but the original 'fit' did not include an xreg.")
    }

    if (is.vector(newdata)) {
      if (length(newdata)  != n.ahead) {
        rlang::abort(paste0("The length of newdata must be ", n.ahead))
      }
    } else if (is.matrix(newdata) || is.data.frame(newdata)) {
      if (nrow(newdata) != n.ahead) {
        rlang::abort(paste0("The length of newdata must be ", n.ahead))
      }
    } else {
      rlang::abort("The 'newdata' parameter must be a vector, matrix or dataframe.")
    }
  }

  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  #

  # next code will transform newdata to ensure it is a dataframe - this is a requirement for the
  # predict.lm() function. It will add on columns for the intercept and drift components.
  # First we check for a vector and if so, then we convert it to a 1-column dataframe, and we
  # add a default column name.
  if (is.vector(newdata)) {
    newdata_name <- internal_make_column_name(as.character(deparse(substitute(newdata))))
    newdata <- data.frame(X1 = newdata)
    colnames(newdata) <- newdata_name
  } else if (is.matrix(newdata))  {
    newdata <- as.data.frame(newdata)
  }

  if (object$include.mean) {
    if (is.null(newdata)) {
      newdata <- data.frame(intercept = rep(1, n.ahead))
    } else {
      newdata$intercept <- 1
    }
  }
  if (object$include.drift) {
    if (is.null(newdata)) {
      newdata <- data.frame(drift = 1:n.ahead)
    } else {
      newdata$drift <- 1:n.ahead
    }
  }

  coef <- unname(object$coef[1, ])
  p <- object$order[1]
  id <- object$order[2]
  q <- object$order[3]
  k <- object$k
  y <- as.numeric(object$diff_y)
  orig_y <- as.numeric(object$y)
  n <- length(orig_y)
  resid <- as.numeric(object$resid)

  mean_y <- mean(y)
  phi_vec <- c(1, -object$model$phi)
  theta_vec <- c(1, -object$model$theta)
  if (any(Mod(polyroot(phi_vec)) < 1) | any(Mod(polyroot(theta_vec)) < 1)) {
    rlang::warn("model estimates are not Stationary! Forecasts may become unbounded.\n")
  }


  if (length(object$model$ggbr_factors) > 0) {
    wgts <- internal_build_weights(length(y) + n.ahead + id)
    y_dash <- y
    for (h in 1:(n.ahead)) {
      yy <- y_dash
      wgts1 <- tail(wgts, length(yy))
      next_forecast <- -(sum(yy * wgts1))
      y_dash <- c(y_dash, next_forecast)
    }

    pred <- y_dash[(length(y) + 1):length(y_dash)]
    if (id > 0) {
      pred <- diffinv(pred + mean_y, differences = id, xi = tail(orig_y, id))
      if (length(pred) > n.ahead) pred <- tail(pred, n.ahead)
    } else {
      if (length(pred) > n.ahead) pred <- tail(pred, n.ahead)
      n <- length(orig_y)
    }
  } else { # ARIMA forecasting only
    y_dash <- y
    phi_vec <- rev(-phi_vec[2:length(phi_vec)])
    if (length(theta_vec) > 1) {
      theta_vec <- rev(-theta_vec[2:length(theta_vec)])
    } else {
      theta_vec <- numeric(0)
    }
    pp <- length(phi_vec)
    qq <- length(theta_vec)
    pred <- rep(0, n.ahead)

    for (i in 1:n.ahead) {
      if (p > 0) {
        if (i > 1) {
          ar_vec <- tail(c(rep(0, pp), y_dash, pred[1:(i - 1)]), pp)
        } else {
          ar_vec <- tail(c(rep(0, pp), y_dash), pp)
        }
        pred[i] <- pred[i] + sum(phi_vec * ar_vec)
      }
      if (q > 0) {
        if (i > qq) {
          ma_vec <- rep(0, qq)
        } else if (i > 1) {
          ma_vec <- tail(c(rep(0, qq), resid, rep(0, i - 1)), qq)
        } else {
          ma_vec <- rep(0, qq)
        }
        pred[i] <- pred[i] + sum(theta_vec * ma_vec)
      }
    }
    if (id > 0) {
      pred <- diffinv(pred, differences = id, xi = tail(orig_y, id))
      if (length(pred) > n.ahead) pred <- tail(pred, n.ahead)
    } else {
      n <- length(orig_y)
      if (object$include.drift) pred <- pred + ((n + 1):(n + length(pred))) * object$model$drift
    }
  }

  # Now we have the forecasts, we set these up as a "ts" object
  y_end <- object$y_end
  if (length(y_end) > 1) {
    if (object$y_freq >= y_end[2]) {
      y_end[1] <- y_end[1] + 1
      y_end[2] <- y_end[2] - object$y_freq + 1
    } else {
      y_end[2] <- y_end[2] + 1
    }
  } else {
    y_end <- y_end + 1
  }

  # Now add on the regression predictions
  if ("lm_xreg" %in% names(object) & !is.null(newdata)) {
    pred <- pred + predict(object$lm_xreg, newdata)
  }

  res <- stats::ts(pred, start = y_end, frequency = object$y_freq)
  return(list(pred = res))
}