fable::ARIMA() compared to forecast::auto.arima()

[robjhyndman]

Different numbers of differences chosen, and excessive memory used.

library(fable)
#> Loading required package: fabletools
#> Registered S3 method overwritten by 'tsibble':
#>   method               from 
#>   as_tibble.grouped_df dplyr
library(forecast)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

# Different order of differences
auto.arima(AirPassengers, lambda = 0, biasadj = TRUE)
#> Series: AirPassengers 
#> ARIMA(0,1,1)(0,1,1)[12] 
#> Box Cox transformation: lambda= 0 
#> 
#> Coefficients:
#>           ma1     sma1
#>       -0.4018  -0.5569
#> s.e.   0.0896   0.0731
#> 
#> sigma^2 = 0.001371:  log likelihood = 244.7
#> AIC=-483.4   AICc=-483.21   BIC=-474.77
as_tsibble(AirPassengers) |>
  model(ARIMA(log(value))) |>
  report()
#> Series: value 
#> Model: ARIMA(2,0,0)(0,1,1)[12] w/ drift 
#> Transformation: log(value) 
#> 
#> Coefficients:
#>          ar1     ar2     sma1  constant
#>       0.5754  0.2614  -0.5553    0.0193
#> s.e.  0.0843  0.0842   0.0771    0.0015
#> 
#> sigma^2 estimated as 0.001323:  log likelihood=249.65
#> AIC=-489.29   AICc=-488.82   BIC=-474.88

# Benchmarks
bench::mark(
  forecast = auto.arima(AirPassengers, lambda = 0, biasadj = TRUE) |>
    forecast(h = 12),
  fable = as_tsibble(AirPassengers) |>
    model(ARIMA(log(value))) |>
    forecast(h = 12),
  check = FALSE
) |>
  select(expression, min, median, `itr/sec`, mem_alloc)
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> # A tibble: 2 × 5
#>   expression      min   median `itr/sec` mem_alloc
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>
#> 1 forecast      1.59s    1.59s     0.628  398.74MB
#> 2 fable          5.5s     5.5s     0.182    1.27GB

^{Created on 2025-05-30 with reprex v2.1.1}

[mitchelloharawild]

The difference with the number of differences is that the forecast package uses trunc(3 * sqrt(n) / 13), while the feasts package uses the urca default of trunc(4 * (n/100)^0.25). For the AirPassengers time series, the forecast package uses 2 lags while feasts uses 4.

The default lag for the test in fable is consistent with urca::ur.kpss() and tseries::kpss.test(). Was there any motivating reason for the default lag in the forecast package being different?

---

The performance differences are still being investigated.

[mitchelloharawild]

The performance difference is due to fable::ARIMA() considering more models in the stepwise search. If I recall correctly, the automatic search algorithm in fable more accurately matches the procedure in the 2008 paper. I think the search procedure in the forecast package misses some of the diagonal steps.

library(forecast)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
auto.arima(AirPassengers, d = 1L, lambda = 0, biasadj = TRUE, stepwise = TRUE, trace = TRUE)
#> 
#>  ARIMA(2,1,2)(1,1,1)[12]                    : Inf
#>  ARIMA(0,1,0)(0,1,0)[12]                    : -434.799
#>  ARIMA(1,1,0)(1,1,0)[12]                    : -474.6299
#>  ARIMA(0,1,1)(0,1,1)[12]                    : -483.2101
#>  ARIMA(0,1,1)(0,1,0)[12]                    : -449.8857
#>  ARIMA(0,1,1)(1,1,1)[12]                    : -481.5957
#>  ARIMA(0,1,1)(0,1,2)[12]                    : -481.6451
#>  ARIMA(0,1,1)(1,1,0)[12]                    : -477.2164
#>  ARIMA(0,1,1)(1,1,2)[12]                    : Inf
#>  ARIMA(0,1,0)(0,1,1)[12]                    : -467.4644
#>  ARIMA(1,1,1)(0,1,1)[12]                    : -481.582
#>  ARIMA(0,1,2)(0,1,1)[12]                    : -481.2991
#>  ARIMA(1,1,0)(0,1,1)[12]                    : -481.3006
#>  ARIMA(1,1,2)(0,1,1)[12]                    : -481.5633
#> 
#>  Best model: ARIMA(0,1,1)(0,1,1)[12]
#> Series: AirPassengers 
#> ARIMA(0,1,1)(0,1,1)[12] 
#> Box Cox transformation: lambda= 0 
#> 
#> Coefficients:
#>           ma1     sma1
#>       -0.4018  -0.5569
#> s.e.   0.0896   0.0731
#> 
#> sigma^2 = 0.001371:  log likelihood = 244.7
#> AIC=-483.4   AICc=-483.21   BIC=-474.77

library(fable)
#> Loading required package: fabletools
as_tsibble(AirPassengers) |>
  model(ARIMA(log(value) ~ pdq(d = 1), stepwise = TRUE, trace = TRUE))
#> Model specification      Selection metric
#> ARIMA(2,1,2)(1,1,1)[12]      Inf
#> ARIMA(0,1,0)(0,1,0)[12]      -434.798957
#> ARIMA(1,1,0)(1,1,0)[12]      -474.629862
#> ARIMA(0,1,1)(0,1,1)[12]      -483.210085
#> ARIMA(0,1,1)(0,1,0)[12]      -449.885657
#> ARIMA(0,1,0)(0,1,1)[12]      -467.464383
#> ARIMA(0,1,2)(0,1,1)[12]      -481.299086
#> ARIMA(1,1,1)(0,1,1)[12]      -481.582022
#> ARIMA(0,1,1)(0,1,2)[12]      -481.645052
#> ARIMA(0,1,1)(1,1,1)[12]      -481.595677
#> ARIMA(0,1,2)(0,1,0)[12]      -447.896470
#> ARIMA(1,1,1)(0,1,0)[12]      -448.062767
#> ARIMA(1,1,0)(0,1,1)[12]      -481.300623
#> ARIMA(1,1,2)(0,1,1)[12]      -481.563300
#> ARIMA(0,1,0)(0,1,2)[12]      -465.375718
#> ARIMA(0,1,2)(0,1,2)[12]      -479.675871
#> ARIMA(1,1,1)(0,1,2)[12]      -479.877452
#> ARIMA(0,1,1)(1,1,0)[12]      -477.216363
#> ARIMA(0,1,0)(1,1,1)[12]      -465.374845
#> ARIMA(0,1,2)(1,1,1)[12]      -479.625161
#> ARIMA(1,1,1)(1,1,1)[12]      -479.830857
#> ARIMA(0,1,1)(1,1,2)[12]      Inf
#> # A mable: 1 x 1
#>   `ARIMA(log(value) ~ pdq(d = 1), stepwise = TRUE, trace = TRUE)`
#>                                                           <model>
#> 1                                       <ARIMA(0,1,1)(0,1,1)[12]>

^{Created on 2025-06-03 with reprex v2.1.1}

In total, fable considers 22 models while forecast considers 14.

library(forecast)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
library(fable)
#> Loading required package: fabletools
bench::mark(
  auto.arima(AirPassengers, d=1, lambda = 0, biasadj = TRUE),
  as_tsibble(AirPassengers) |>
    model(ARIMA(log(value) ~ pdq(d=1))),
  check = FALSE
)
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> # A tibble: 2 × 6
#>   expression                             min median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>                           <bch> <bch:>     <dbl> <bch:byt>    <dbl>
#> 1 auto.arima(AirPassengers, d = 1, la… 1.76s  1.76s     0.567     401MB     5.67
#> 2 model(as_tsibble(AirPassengers), AR… 2.13s  2.13s     0.470     604MB     6.58

^{Created on 2025-06-03 with reprex v2.1.1}

So 0.126s/model for forecast, and 0.09s/model for fable. The internal Arima code is more efficient in fable, but there is some overhead in the generality of fable for handling data frames, guessing seasonality, symbolic inversion of the transformation, etc. So the single model performance is slightly worse for fable.

[robjhyndman]

Thanks. Now I recall playing with the urca::ur.kpss arguments -- see this post. I remember that one of the reasons I didn't want to change it was to keep the d=D=1 result for the AirPassengers data.

[mitchelloharawild]

I'm open to changing the default test for d, but I think more evidence than it giving better results on a single classic example is needed. Are there other newer unit root tests that behave better?

[robjhyndman]

I agree. I think sticking with KPSS is ok, but we could investigate the various suggestions for lags that have been proposed, and test the resulting forecasting performance on a large data set.