When is the AIC a good model selection criterion for forecasting and when is it not?

I'm trying to wrap my head around why the AIC and other similar ICs work as proxies for out of sample error when trying to perform automated forecast generation.

So I performed an experiment on the Air Passengers data set, which is as forecastable as a real world time series can get. I used data through the end of 1958 as the training set and the data from 1959 and 1960 as the hold out set.

The results I got were surprising:

Using auto.arima() (with stepwise=FALSE and approximation=FALSE), the best fitting model selected was an $$ARIMA(0,1,3)(0,1,0)_{12}$$ model, with the AIC = 802.346 and a test RMSE = 69.235753.

I then fit an "over the top" ARIMA model with extreme parameters $$ARIMA(15,1,15)(4,1,4)_{12}$$. As expected, this lead to a higher AIC of 829.2997, but with a much lower RMSE = 27.871115.

Then I tried a more reasonable ARIMA model (i.e. one within the range that would have been considered by auto.arima()), $$ARIMA(2,1,2)(2,1,2)_{12}$$ and it gave me an AIC=804.56, and an RMSE = 42.051491.

From the plots, of the forecasts, it is even more clear that the model with the highest AIC is the one giving the best forecast. The AICc and the BIC gave a similar reverse behavior, with the model with the highest IC value giving the lowest RMSE and MAPE.

My questions:

• Isn't this the exact opposite of the intended behavior? I thought minimizing the AIC (or any of the other ICs) gives the best model, not maximizing it?
• Is there something wrong with my experiment, that's giving me the counter intuitive result?
• Does the fact that the Air Passengers time series is very regular and forecastable have something to do with this? Would the AIC work better for very noisy data?
• When are the AIC and other ICs appropriate for time series model selection if not in this case?

#Call the necessary libraries
library('ggplot2')
library('forecast')
library(zoo)
library(scales)
theme_set(theme_bw())

#load the airpassengers data
data("AirPassengers")

#Split the data into test and train
train <- window(AirPassengers, end = c(1958, 12))
test <- window(AirPassengers, start = c(1959, 1), end = c(1960,12))

#Fit the models
fit <- auto.arima(train, stepwise = FALSE, approximation = FALSE)
fit2 <- Arima(train, order=c(15, 1, 15), seasonal = list (order= c(4, 1, 4) , period = 12), method='ML')
fit3 <- Arima(train, order=c(2, 1, 2), seasonal = list (order= c(2, 1, 2) , period = 12), method='ML')

#Generate forecasts
#I am setting forecast intervals to 0 so that they are not displayed for better clarity
arima_fct <- forecast(fit,level = c (0,0), h=24)
arima_fct2 <- forecast(fit2,level = c (0,0), h=24)
arima_fct3 <- forecast(fit3,level = c (0,0), h=24)

fit$$aic fit2$$aic
fit3$aic accuracy(arima_fct,test) accuracy(arima_fct2,test) accuracy(arima_fct3,test) #Plot results autoplot(arima_fct , ylab = 'Passengers') + scale_x_yearmon() + autolayer(test, series="Test Data") + autolayer(arima_fct$$mean, series="ARIMA(0,1,3)(0,1,0)[12]: AIC = 802.3461, Test RMSE = 69.235753") + autolayer(arima_fct3$$mean, series="ARIMA(2,1,2)(2,1,2)[12]: AIC = 804.56, Test RMSE = 42.051491") + autolayer(arima_fct2$mean, series="ARIMA(15,1,15)(4,1,4)[12]: AIC = 829.2997, Test RMSE = 27.871115")

• For your first model, I get the same AIC as you report. The second model has been running a long time in fitting and shows no sign of finishing. For the third model, I get an AIC of 814.19, different from yours. For the first model, accuracy gives a test set RMSE of 259.8, for the third model 282.9, both very different from your numbers. Please edit your question to provide your code. In the meantime, I am voting to close as unclear. Aug 7, 2018 at 13:09
• @StephanKolassa I added my code. Aug 7, 2018 at 15:25
• @StephanKolassa note that I have switched stepwise and approximation off in auto.arima to avoid the issues mentioned in this post Aug 7, 2018 at 15:31
• Am I understanding correctly that you’re just comparing a single forecast? There is of course a lot of variability involved. I’d suggest doing a rolling forecasting exercise to reduce the impact of noise. Aug 7, 2018 at 18:37
• @hejseb understood - that's why I'm running a basic test on a very predictable series. The pattern is almost the same from one year to the next if you take into account the trend and the multiplicative seasonality. Aug 7, 2018 at 18:44

1 Answer

I am not completely satisfied with my answer, but here goes.

1. To an extent, you are comparing apples and oranges. Your two calls to Arima() use method="ML", whereas your auto.arima() uses the default, which is method="CSS-ML". Then again, refitting everything with the default does not make a real difference.

2. Minimizing the AIC is asymptotically equivalent to minimizing the one-step ahead squared prediction error. (I don't have a reference at hand, sorry.) Note that this is an asymptotic result in a suitable statistical sense. It's quite possible for a handpicked model to outperform AIC on a limited length time series. And on a single one, at that.

3. Finally, as you write in a comment, the AirPassengers dataset exhibits strong multiplicative seasonality. ARIMA does not model multiplicative seasonality or trend; it can only deal with additive effects. Your overparameterized model gets the multiplicative trend and seasonality right, but it may also forecast this in a series that does not exhibit such effects. There are reasons why such large models are typically not considered.

To model multiplicative effects, allow auto.arima() to use Box-Cox transformations:

> (foo <- auto.arima(train,lambda="auto"))
Series: train
ARIMA(0,1,1)(0,1,1)[12]
Box Cox transformation: lambda= -0.3096628

Coefficients:
ma1     sma1
-0.3936  -0.5713
s.e.   0.1035   0.0863

> accuracy(forecast(foo,h=24,biasadj=TRUE),test)
ME      RMSE       MAE        MPE     MAPE      MASE       ACF1 Theil's U
Training set -0.7186038  8.915531  6.691014 -0.2079082 2.753580 0.2341638 0.04889565        NA
Test set     28.5600533 31.711896 28.884516  6.2710488 6.348486 1.0108644 0.17279165 0.6372069


I cut out the AIC, because that is not comparable to the AIC on nontransformed data. Note that we end up much closer to your large model in terms of the test RMSE, but the model is much more interpretable, and I personally would trust it a lot more than an ARIMA(15,1,15)(4,1,4)[12] one. Incidentally, searching through more possible ARIMA models yields the exact same model:

> (bar <- auto.arima(train,max.p=15,max.q=15,max.P=4,max.Q=4,
+ lambda="auto",stepwise=FALSE,approximation=FALSE))
Series: train
ARIMA(0,1,1)(0,1,1)[12]
Box Cox transformation: lambda= -0.3096628

Coefficients:
ma1     sma1
-0.3936  -0.5713
s.e.   0.1035   0.0863

• Dr Kolassa, this is a very informative answer, especially with regards to using the Box-Cox transformation and how to better use ARIMA models. However, I'm still hung up on the original question about the use of the AIC and the BIC: I get it that they are asymptotic approximations, but they are so widely used in practice (the AIC in the Forecast package, the BIC in the demand forecasting software that I use in my day job) that I assumed that they would work on a time series as generic as Air Passengers. The cases were they failed would be more unusual time series (with brakes or noise, etc...) Aug 9, 2018 at 6:17
• That is why I am not completely happy with my answer. Then again, it's a stretch to assume that data are actually generated by an ARIMA model, so even choosing the "best" ARIMA orders won't give you the data generating process. There is a reason why forecasting is still an active research field. Aug 9, 2018 at 6:29
• Maximum likelihood minimizes the error of fit to a set of observations, it is not an optimization with respect to least error of forecasting future tense withheld observations. So in general, polishing the ML apple with AIC/BIC, and the like, will not convert it into an orange. For that, you need to start with an orange, where in this case, the orange is an explanatory physical model of the cause and effect type.
– Carl
Oct 9, 2018 at 22:18
• You seem to contrast an estimation teachnique (MLE) with model definition/building (explanatory physical model of the cause and effect type). These are tangential, hence not comparable/contrastable. Oct 11, 2018 at 14:45