R/inars_forecasting.R

In rdmatheus/inars: INARS(1) Model for Integer-Valued Time Series

#' @export
forecast <- function(fit, h, xh = NULL, B = 1000, level = 0.95){

  y <- fit$response
  n <- length(y)
  p <- ncol(fit$x)
  est <- fit$estimates

  # Predictor
  pred <- vector()
  if (p == 1){
    m.h <- rep(est[2], h)
  } else {
    m.h <- xh%*%est[2:(p + 1)]
  }

  for(k in 1:h){
    pred[k] <- est[1]^k * y[n] + sum(est[1]^(0:(k - 1)) * m.h[k:1])
  }

  ## Step 1
  e <- fit$residuals

  ## Step 2
  Fr <- round(e)

  ## Step 3

  # Signed binomial thinning operation
  sign_thinning <- function(alpha, y){
    sign(alpha) * sign(y) * stats::rbinom(1, abs(y), abs(alpha))
  }

  forecast <- matrix(NA, B, h)
  for(j in 1:B){

    y_star <- vector()
    e_star <- sample(Fr, n, replace = TRUE)
    y_star[1] <- e_star[1]

    for (t in 2:n){
      y_star[t] <- sign_thinning(est[1], y_star[t-1]) + e_star[t]
    }

    ## Step 4
    alpha_star <- inars_est(y_star, X = NULL, innovation = fit$innovation,
                            method = "MM")$par[1]

    ## Step 5
    yh <- vector()
    yh[1] <- y[n]
    for (t in 2:(h + 1)){
      yh[t] <-  sign_thinning(alpha_star, yh[t-1]) + sample(Fr, 1, replace = TRUE)
    }

    forecast[j, ] <- yh[-1]

  }

  level <- 1 - level
  MD <- apply(forecast, 2, quantile, probs = 0.5)
  LI <- apply(forecast, 2, quantile, probs = level/2)
  LS <- apply(forecast, 2, quantile, probs = 1 - level/2)
  list(forecast = pred,
       LI = LI,
       LS = LS)
}