# How to optimise an automatic ARIMA-model selection?

I've been using statsmodels.tsa.arima_model to fit the residual component of some data. I've written an algorithm to automatically select the ARIMA model. Results are not quite as good as I had hoped, so I am looking for suggestions on how could I improve things. Please find below a description of what I've tried thus far.

I am dividing the data into a training and a forecasting/testing set. Up to this point, I've based my choice of the ARIMA model only on the training set.

I allow the ARIMA parameters p and q to run from 0 to 7. This choice of range is arbitrary. The d parameter is allowed to be either 0 or 1. This function tries each of them and storages the results:

def iterative_ARIMA_fit(series):
""" Iterates within the allowed values of the p and q parameters

Returns a dictionary with the successful fits.
Keys correspond to models.
"""
ARIMA_fit_results = {}
for AR in ARrange :
for MA in MArange :
for Diff in Diffrange:
model = ARIMA(series, order = (AR,Diff,MA))
try:
results_ARIMA = model.fit(disp = -1, method = 'css')
RSS = sum((results_ARIMA.fittedvalues - series)**2)
if RSS > 0:
ARIMA_fit_results['%d-%d-%d' % (AR,Diff,MA)]=[RSS,results_ARIMA]
except:
continue
return ARIMA_fit_results


Next, I look for the model that miminises RSS (total squared residual) using:

def get_best_ARIMA_model_fit(series):
""" Returns a list with the best ARIMA model

The first element on the list contains the squared residual
The second element on the list contains the fit results
"""
if t.isstationary(series)[0]:
ARIMA_fit_results = iterative_ARIMA_fit(series)
best_ARIMA = min(ARIMA_fit_results, key = ARIMA_fit_results.get)
return ARIMA_fit_results[best_ARIMA]


I had initially tried to use as much data as possible for training. To my surprise, despite using up to x6 more data, the fitting quality worsens.

I made a scan of the training set length, plotting for each the total squared residual and the total squared residual per unit of training set. This is shown in 1. Here, clearly, the fit quickly deteriorates the more data I use, up to some point, where the ARIMA-fit-residual stabilizes. Also here, the best ARIMA model is consistently the highest-order available MA model. 0-0-7.

Summarising, my optimisation strategy yields an ARIMA fit that:

• (Fixed) Consistently (stubbornly) selects the MA model with of the highest available order
• Does not improve with the size of the dataset

Do any of you know whether this is somehow expected? Has any of you tried something different?

Update

The rejection of all the AR models was due to an error on the first function. These models return statsmodels.tsa.arima_model.ARMAResults.fittedvalues with less entries than the original series and, hence, the RSS returns nan. These new function does the trick:

def iterative_ARIMA_fit(series):
""" Iterates within the allowed values of the p and q parameters

Returns a dictionary with the successful fits.
Keys correspond to models.
"""
ARIMA_fit_results = {}
for AR in ARrange :
for MA in MArange :
for Diff in Diffrange:
model = ARIMA(series, order = (AR,Diff,MA))
fit_is_available = False
results_ARIMA = None
try:
results_ARIMA = model.fit(disp = -1, method = 'css')
fit_is_available = True
except:
continue
if fit_is_available:
safe_RSS = get_safe_RSS(series, results_ARIMA.fittedvalues)
ARIMA_fit_results['%d-%d-%d' % (AR,Diff,MA)]=[safe_RSS,results_ARIMA]

return ARIMA_fit_results


Plus this extra one:

def get_safe_RSS(series, fitted_values):
""" Checks for missing indices in the fitted values before calculating RSS

Missing indices are assigned as np.nan and then filled using neighboring points
"""
fitted_values_copy = fitted_values  # original fit is left untouched
missing_index = list(set(series.index).difference(set(fitted_values_copy.index)))
if missing_index:
nan_series = pd.Series(index = pd.to_datetime(missing_index))
fitted_values_copy = fitted_values_copy.append(nan_series)
fitted_values_copy.sort_index(inplace = True)
fitted_values_copy.fillna(method = 'bfill', inplace = True)  # fill holes
fitted_values_copy.fillna(method = 'ffill', inplace = True)
return sum((fitted_values_copy - series)**2)


The results are much better, albeit the overfitting. Here 2 is the updated plot showing the residuals and model choices.

## 1 Answer

You are . If you have the choice between an MA($q$) and an MA($q+1$) model, the larger model with more degrees of freedom will almost always fit the data better and yield smaller residual sums of squares. (I would have expected the same to happen for the AR orders, but that this does not happen may be due to the fact that you are modeling residuals.)

ARIMA models are typically selected based on information criteria, like , AICc, or , after deciding on whether to difference or not based on a statistical test. The documentation for the auto.arima() function in the forecast package for R may give you some inspiration as to what to look at.

Edit: Cagdas Ozgenc correctly notes that increasing the MA order will not necessarily always reduce the residual sums of squares, because the conditional sum of squares estimation is not convex. To illustrate this effect, I simulated 10,000 white noise time series of 100 realization each, fitted MA($q$) models for $q=0, \dots, 7$ and noted the RSS. Below are boxplots of $$\Delta(q) := \text{RSS}_{\text{MA}(q)}-\text{RSS}_{\text{MA}(q-1)}$$ against $q$. Out of the $10,000\times 7=70,000$ possible differences, $69,851 = 99.8%$ were negative, i.e., a larger model yielded smaller RSS - although there were zero moving average dynamics in the simulated series.

R code:

rm(list=ls())
library(forecast)

n.series <- 1e4
nn <- 100
ma.max <- 7

rss <- matrix(NA,nrow=n.series,ncol=ma.max+1)

pb <- winProgressBar(max=n.series)
for ( ii in 1:n.series ) {
setWinProgressBar(pb,ii,paste(ii,"of",n.series))
set.seed(ii)
xx <- ts(rnorm(nn))
for ( kk in 0:ma.max ) {
model <- Arima(xx,order=c(0,0,kk),method="CSS")
rss[ii,kk+1] <- sum(model\$residuals^2)
}
}
close(pb)

differences <- apply(rss,1,diff)
boxplot(t(differences),main="RSS differences between MA(q) and MA(q-1) models",xlab="q")
abline(h=0)

sum(differences<0)/prod(dim(differences))

• Thanks @StephanKolassa. Found the reason for any AR order being rejected. For future reference, I've added a second version of the example code and figure. – Nicolas Gutierrez Jun 7 '16 at 15:00
• Why do you think a higher order model will "always" fit the data better? Always is a very strong word. The fitting procedure is not convex, only for that reason a complex model may end up fitting worse. If it were convex, then I am more inclined towards your argument, but would be nice if you could provide a few pointers to such a proof. – Cagdas Ozgenc Jun 7 '16 at 15:30
• @CagdasOzgenc: thank you for your patience. I have edited the answer to include your point. – Stephan Kolassa Jun 11 '16 at 9:30
• Even in the convex case your argument is probably true only for a set of hierarchical linear models (which is the case here with ARMA without gaps in lags). If the more complex model is not subsuming the less complex model, it will not work. Also for non-linear models even if the fitting is convex I am not sure what will happen. Do you know any proofs or results? – Cagdas Ozgenc Jun 13 '16 at 10:57
• That may well be, but I'm only discussing nested models here, which is the situation the OP is in. No, I don't have any general proofs or results about time series overfitting "in general". – Stephan Kolassa Jun 14 '16 at 6:33