5. ARIMA: integration, fitting and rolling-window forecasting

Level 9: Time Series - ARIMA and GARCH

The autoregressive integrated moving average model, or ARIMA, is the standard statistical model for analyzing and forecasting time series. Its acronym describes its three pieces exactly:

• AR: autoregression. It uses the dependent relationship between an observation and some number of lagged observations.
• I: integrated. It applies differencing to the raw observations (subtracting from each observation the one at the previous step) to make the series stationary.
• MA: moving average. It uses the dependency between an observation and the residual errors of a moving average model applied to lagged observations.

Its three parameters:

• p: number of lagged observations included, also called the lag order.
• d: number of times the raw observations are differenced, or degree of differencing.
• q: the size of the moving average window, or moving average order.

When d is zero we have an ARMA(p,q), which combines both effects:

# ARMA(p,q):
# x[t] = alpha_1*x[t-1] + ... + alpha_p*x[t-p]
#        + w[t] + beta_1*w[t-1] + ... + beta_q*w[t-q]

An ARMA usually needs fewer parameters than an AR(p) or an MA(q) separately for the same fit. And a warning that will bridge to the last two modules: an ARMA does not account for volatility clustering, a key empirical feature of many financial series.

Basic fitting (use the modern class, not the deprecated one):

from statsmodels.tsa.arima.model import ARIMA

model = ARIMA(Re, order=(5, 1, 0))
model_fit = model.fit()
print(model_fit.summary())

residuals = pd.DataFrame(model_fit.resid)
residuals.plot()
residuals.plot(kind='kde')
print(residuals.describe())

A residual mean close to zero indicates no bias; a positive or negative mean indicates a bias that can be corrected directly before using the forecasts.

Watch out for a detail the course code lets slide: Re is already log returns, that is, the series has already been differenced once. Adding d = 1 on top means over-differencing, which introduces an artificial MA structure into the residuals and inflates the variance. On returns, d = 0 is almost always the right choice.

Rolling-window ARIMA forecasting. Just as with the AR, we want to rebuild the model every day. Here it is refitted from scratch at every step:

from tqdm import tqdm

X = Re.values
size = int(len(X) * 0.66)
train, test = X[0:size], X[size:len(X)]
history = [x for x in train]
predictions = []

for t in tqdm(range(len(test)), desc='Processing'):
    model = ARIMA(history, order=(5, 1, 0))
    model_fit = model.fit()
    output = model_fit.forecast()
    yhat = output[0]
    predictions.append(yhat)
    history.append(test[t])       # only after forecasting

rmse = np.sqrt(mean_squared_error(test, predictions))
print('Test RMSE: %.3f' % rmse)

This is honest but costly: as many maximum likelihood fits as there are test days.

Hyperparameter grid search. Instead of fixing (5,1,0) by eye, we evaluate combinations and keep the one with the lowest out-of-sample RMSE:

def evaluate_arima_model(X, arima_order):
    train_size = int(len(X) * 0.66)
    train, test = X[0:train_size], X[train_size:]
    history = [x for x in train]
    predictions = list()
    for t in range(len(test)):
        model_fit = ARIMA(history, order=arima_order).fit()
        predictions.append(model_fit.forecast()[0])
        history.append(test[t])
    return np.sqrt(mean_squared_error(test, predictions))

def evaluate_models(dataset, p_values, d_values, q_values):
    dataset = dataset.astype('float32')
    best_score, best_cfg = float('inf'), None
    for p in p_values:
        for d in d_values:
            for q in q_values:
                order = (p, d, q)
                try:
                    rmse = evaluate_arima_model(dataset, order)
                    if rmse < best_score:
                        best_score, best_cfg = rmse, order
                    print('ARIMA%s RMSE=%.3f' % (order, rmse))
                except:
                    continue
    print('Best ARIMA%s RMSE=%.3f' % (best_cfg, best_score))

evaluate_models(Re.values, [0, 1, 2], [0, 1], [0, 1])

The bare except: continue is convenient (many combinations fail to converge) but dangerous: it also silences programming errors. In production, catch the exception and log it.