Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multi-step fitted ARIMA faster solution #949

Merged
merged 2 commits into from
Dec 21, 2023

Conversation

danigiro
Copy link
Contributor

@danigiro danigiro commented Dec 20, 2023

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 robjhyndman merged commit f481632 into robjhyndman:master Dec 21, 2023
6 checks passed
@robjhyndman
Copy link
Owner

Thanks @danigiro

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants