Are models identified by auto.arima() parsimonious?

I have been trying to learn and apply ARIMA models. I have been reading an excellent text on ARIMA by Pankratz - Forecasting with Univariate Box - Jenkins Models: Concepts and Cases. In the text the author especially emphasizes the priciple of parsimony in choosing ARIMA models.

I started playing with auto.arima() function in R package forecast. Here is what I did, I simulated ARIMA and then applied auto.arima(). Below are 2 examples. As you can see in both example auto.arima() clearly identified a model that many would consider non-parsimonious. Especially in example 2, where auto.arima() identified ARIMA(3,0,3) when actually ARIMA(1,0,1) would be sufficient and parsimonious.

Below are my questions. I would appreciate any suggestions and recommendations.

1. Are there any guidance on when to use/modify the models identified using automatic algorithms such as auto.arima()?
2. Are there any pit falls in just using AIC (which is what I think auto.arima() uses) to identify models?
3. Can an automatic algorithm built that is parsimonious?

By the way I used auto.arima() just as an example. This would apply to any automatic algorithm.

Below is Example #1:

set.seed(182)
y <- arima.sim(n=500,list(ar=0.2,ma=0.6),mean = 10)

auto.arima(y)

qa <- arima(y,order=c(1,0,1))
qa

Below are the results from auto.arima(). Please note that all the coefficients are insignificant. i.e., $t$ value < 2.

ARIMA(1,0,2) with non-zero mean

Coefficients:
ar1     ma1      ma2  intercept
0.5395  0.2109  -0.3385    19.9850
s.e.  0.4062  0.4160   0.3049     0.0878

sigma^2 estimated as 1.076:  log likelihood=-728.14
AIC=1466.28   AICc=1466.41   BIC=1487.36

Below are the results from running regular arima() with order ARIMA(1,0,1)

Series: y
ARIMA(1,0,1) with non-zero mean

Coefficients:
ar1     ma1  intercept
0.2398  0.6478    20.0323
s.e.  0.0531  0.0376     0.1002

sigma^2 estimated as 1.071:  log likelihood=-727.1
AIC=1462.2   AICc=1462.28   BIC=1479.06

Example 2:

set.seed(453)
y <- arima.sim(n=500,list(ar=0.2,ma=0.6),mean = 10)

auto.arima(y)

qa <- arima(y,order=c(1,0,1))
qa

Below are the results from auto.arima():

ARIMA(3,0,3) with non-zero mean

Coefficients:
ar1      ar2     ar3     ma1     ma2     ma3  intercept
0.7541  -1.0606  0.2072  0.1391  0.5912  0.5491    20.0326
s.e.  0.0811   0.0666  0.0647  0.0725  0.0598  0.0636     0.0939

sigma^2 estimated as 1.027:  log likelihood=-716.84
AIC=1449.67   AICc=1449.97   BIC=1483.39

Below are the results running regular arima() with order ARIMA(1,0,1)

Series: y
ARIMA(1,0,1) with non-zero mean

Coefficients:
ar1     ma1  intercept
0.2398  0.6478    20.0323
s.e.  0.0531  0.0376     0.1002

sigma^2 estimated as 1.071:  log likelihood=-727.1
AIC=1462.2   AICc=1462.28   BIC=1479.06
• You appear to be on to something here as the AIC procedure might be mis-identifying the simplest model. Additionally the AIC criteria premises no pulses/no level shifts/no seasonal pulses/no local time trends/constancy of parameters/constancy of error variance. – IrishStat Sep 26 '13 at 14:52
• The AIC and BIC are proportional to the variance of the errors from a guessed model. This approach in my experience is just to simple due to some of the issues I have already raised and the invertibility issues raise here. There is no replacement for an intelligent identification system that constructs iterative models concluding when no uneeded parameters remain and no information/structure is evident in the residuals. In summary automatic model identification is an iterative process not a one-step process just as it is in all statistical analysis. – IrishStat Sep 26 '13 at 17:57
• @Irishstat what is the sign of ma coefficient should it be interpretted as -0.1391 -0.5912 -0.5491 ? – forecaster Sep 26 '13 at 19:36
• in terms of the Pankratz Text .. yes ! . You can check the roots of the ma polynomial to see if they meet the invertibility requirements. – IrishStat Sep 26 '13 at 20:41
• so the above auto.arima is not invertible i.e, -0.1391-0.5912-0.5491 is <1 therefore the model is fine. – forecaster Sep 26 '13 at 20:57

There are a couple of issues here. Firstly, don't presume that the simulated ARIMA is truly of the order you specify; you are taking a sample from the specified model and due to randomness, the best fitting model for the particular sample drawn may not be the one from which the simulations were drawn.

I mention this because of the second and more important issue: the auto.arima() function can estimate models via a more efficient fitting algorithm, using conditional sums of squares, to avoid excessive computational time for long series or for complex seasonal models. When this estimation process is in use, auto.arima() approximates the information criteria for a model (because the log likelihood of the model has not been computed). A simple heuristic is used to determine whether the conditional sums of squares estimation is active, if the user does not indicate which approach should be used.

The behaviour is controlled via argument approximation and the simple heuristic is (length(x)>100 | frequency(x)>12), hence approximation takes a value TRUE if the length of the series is greater than $n = 100$, or there are more than 12 observations within each year. As you simulated series with $n = 500$ but did not specify a value for the approximation argument, you ran auto.arima() with approximation = TRUE. This explains the apparently erroneous selection of a model with larger AIC, AICc, and BIC than the simpler model you fitted with arima().

For your example 1, we should have

> auto.arima(y, approximation = FALSE)
Series: y
ARIMA(0,0,1) with non-zero mean

Coefficients:
ma1  intercept
0.7166    19.9844
s.e.  0.0301     0.0797

sigma^2 estimated as 1.079:  log likelihood=-728.94
AIC=1463.87   AICc=1463.92   BIC=1476.52
> qa
Series: y
ARIMA(1,0,1) with non-zero mean

Coefficients:
ar1     ma1  intercept
0.0565  0.6890    19.9846
s.e.  0.0626  0.0456     0.0830

sigma^2 estimated as 1.078:  log likelihood=-728.53
AIC=1465.06   AICc=1465.14   BIC=1481.92

Hence auto.arima() has selected a more parsimonious model than the true model; an ARIMA(0, 0, 1) is chosen. But this is based on the information criteria and now they are in accordance; the selected model has lower AIC, AICc, and BIC, although the differences for AIC and AICc are small. At least now the selection is consistent with the norms for choosing models based on information criteria.

The reason for the MA(1) being chosen, I believe, relates to the first issue I mentioned; namely that the best fitting model to a sample drawn from a stated ARIMA(p, d, q) may not be of the same order as the true model. This is due to random sampling. Taking a longer series or a longer burn in period may help increase the chance that the true model is selected, but don't bank on it.

Regardless, the moral here is that when something looks obviously wrong, like in your question, do read the associated man page or documentation to assure yourself that you understand how the software works.

• thanks for the detailed response. I did use your approach for the second example: set.seed(453) y <- arima.sim(n=500,list(ar=0.2,ma=0.6),mean = 10) auto.arima(y,approximation=FALSE) and this is what i get, it is clearly overfitting the data -- ARIMA(2,0,4) with non-zero mean Coefficients: ar1 ar2 ma1 ma2 ma3 ma4 intercept 0.5369 -0.9627 0.3681 0.6799 0.7065 0.1701 20.0329 s.e. 0.0278 0.0499 0.0533 0.0630 0.0793 0.0574 0.0927 sigma^2 estimated as 1.024: log likelihood=-716.17 AIC=1448.33 AICc=1448.63 BIC=1482.05 – forecaster Sep 26 '13 at 15:48
• I have not yet read the chapter that has something called "invertiblity". Does the auto.arima in the second case violate the "invertiblity" rule in model diagnosis?, I may be off, can you please correct if I'm wrong? – forecaster Sep 26 '13 at 15:59
• @forecaster Assume you don't know the truth, all you have is the sample to hand. In that circumstance, AIC etc are doing what you are asking them to do. As I said above, the simulated series may not be the one you asked for, or rather it may not be possible to successfully identify the true series from the small sample you took. There is nothing implied here or elsewhere in stats that you'll get the one, TRUE model. As the statistician, if you think the series is overfitted (say from prior knowledge), place a restriction on the order of the terms searched. Or use BIC as the stopping criterion. – Gavin Simpson Sep 26 '13 at 16:06
• I don't know sufficiently well what invertibility is to be able to answer that question. Perhaps ask this as a new question here? – Gavin Simpson Sep 26 '13 at 16:07
• @forecaster (Reason I say BIC, is that it places an extra penalty on complexity. In this case, the BIC of the ARIMA(1, 0, 1) has lower BIC than the model auto.arima() settled on, so if those were the only two candidate models, the simpler one would have been selected...) – Gavin Simpson Sep 26 '13 at 16:08

Thanks so much@Gavin, @Irishstat and @Rob for responding to my question. It is clear that if I need a parsimonious model from an automatic algorithms like auto.arima BIC information criterion should be used as opposed to AIC especially after looking this post and @Gavin's post above.

I also very much agree with @Irishstat that choosing a model based on IC criterion has limitations in that it does not choose a better model to fit a data with outliers and level shifts. In my opinion, outliers + level shifts + messy data = real word business data, anything else is textbook datasets. Any automatic model that doesn't consider outliers + level shifts, again in my opinion should be used with caution.

Coming to the code - auto.arima has an option to choose between the AIC or BIC. See below the code has been modified from teh above questions.

Many thanks Cross-validated community. I learn new and intresting things every day.

###############
set.seed(453)
y <- arima.sim(n=500,list(ar=0.2,ma=0.6),mean = 10)

## Adequetly describes the unknown data
fit.aic <- auto.arima(y,ic = c("aic"))
fit.aic

## Selects the model that is parsimonious
fit.bic <- auto.arima(y,ic = c("bic"))
fit.bic

BIC IC chooses an MA(2) model.

> fit.bic
Series: y
ARIMA(0,0,2) with non-zero mean

Coefficients:
ma1     ma2  intercept
0.9256  0.2335    20.0326
s.e.  0.0453  0.0444     0.0992

sigma^2 estimated as 1.059:  log likelihood=-724.19
AIC=1456.39   AICc=1456.47   BIC=1473.24

I took the 500 values into AUTOBOX ( a piece of commercial software that I have helped develop) and received the following advisory based upon the Chow Test for constancy of parameters. A very basic mistake that is made in the study of a time series is to assume that the data is driven by a particular model with constant parameters. AUTOBOX detected a break point at period 246 which might reflect a simulation that had not been "warmed up". When simulating data good practice is to delete the first "n" values and to then study the remaining ones. I took the data and segmented it into two sections ; the first 245 and the remaining 255 . Here are the two very dissimilar acf plots .

Returning to the analysis: Here is the model that was identified for the last 246 values and here with the following statistics . The Actual/Fit and Forecast is here with residual plot here . The ACF of the residuals suggests sufficiency . Note that the 5 pulses that were identified had a very small effect and could easily be disregarded ( in this case ! ). In summary the lesson learned here is that sometimes we have too much data and we need to consider time changing coefficients. In this case we are identifying a change in parameters which (apparently) doesn’t have a major impact on the resultant model/parameters but it points out a generally needed process improvement in time series analysis. My experience with auto.arima suggests that since it explicitely doesn't treat/remedy gaussian violations it tends to over-model by leaning too much on historical values rather than extracting structure from the data. In this case since it was a tightly controlled simulation without gaussian violations it worked but I would be generally suspect of such a bandwith limited and single step approach to ARIMA model identification. Trust but Verify !

Akaike_information_criterion says that $\ \exp^{(AIC_{min} − AIC_i)/2}$
"can be interpreted as the relative probability that the ith model minimizes the (estimated) information loss".

IF this is so, it would help users to see these relative probabilities, together with the AICs (?) from auto.arima( ... trace=TRUE ). For example, the eggs data run as in this question gives

# relprob = exp( (AICmin - AIC) / 2 ) * 100
ARIMA(0,1,0) with drift : 784.5    100
ARIMA(0,1,1) with drift : 784.8     86
ARIMA(1,1,0) with drift : 784.9     82
ARIMA(0,1,0)            : 792.4      2
ARIMA(2,1,2) with drift : Inf    0
ARIMA(1,1,1) with drift : Inf    0
• This is an answer or a question? – Tim May 11 '15 at 10:35
• @Tim, I'm so tentative becuase I don't know if these relative probs are "real" or at least common. – denis May 11 '15 at 10:52