robjhyndman / forecast

Forecasting Functions for Time Series and Linear Models
http://pkg.robjhyndman.com/forecast
1.1k stars 340 forks source link

Multi-step fitted ARIMA faster solution #949

Closed danigiro closed 6 months ago

danigiro commented 6 months ago

Multi-step fitted ARIMA faster solution

Hi @robjhyndman. Working with the forecast function hfitted I found a solution to make it faster in the ARIMA case by also looking at Fable's implementation. So, the hfitted function is splitted into two: hfitted.default and hfitted.Arima. The first is the general implementation, and the second is a re-implementation of the Fable's hfitted.ARIMA function with the KalmanRun and KalmanForecast functions.

If you're interested, here I did a comparison in terms of computational time (the fitted returns are equal):

library(forecast)
library(microbenchmark)
mod <- auto.arima(AirPassengers)
benchtime <- microbenchmark("old_version" = forecast:::hfitted(mod, h = 2), 
                            "new_version" = hfitted.Arima(mod, h = 2), 
                            check='equal', times = 10)
print(benchtime, digits =2)
#> Unit: milliseconds
#>        expr    min     lq   mean median     uq    max neval cld
#> old_version 821.38 841.44 859.82 850.38 861.93 935.31    10  a 
#> new_version  66.74  68.47  69.24  69.09  69.48  73.55    10   b

The hfitted.Arima function (already implemented in this pull request) is

hfitted.Arima <- function(object, h, ...) {
  if(h == 1){
    return(object$fitted)
  }
  y <- object$fitted+residuals(object, "innovation")
  yx <- residuals(object, "regression")
  # Get fitted model
  mod <- object$model
  # Reset model to initial state
  mod <-  stats::makeARIMA(mod$phi, mod$theta, mod$Delta)
  # Calculate regression component
  xm <- y-yx
  fits <- rep_len(NA_real_, length(y))

  start <- length(mod$Delta) + 1
  end <- length(yx) - h
  idx <- if(start > end) integer(0L) else start:end
  for(i in idx) {
    fc_mod <- attr(stats::KalmanRun(yx[seq_len(i)], mod, update = TRUE), "mod")
    fits[i + h] <- stats::KalmanForecast(h, fc_mod)$pred[h] + xm[i+h]
  }
  fits <- ts(fits)
  tsp(fits) <- tsp(object$x)
  fits
}
robjhyndman commented 6 months ago

Thanks @danigiro