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.

 A: You are overfitting. 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 aic, AICc, or bic, 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))

