3. Forecasting with AR: intervals and walk-forward

Level 9: Time Series - ARIMA and GARCH

With the best model already selected we can generate a multi-step forecast. Here, 21 days of S&P 500 returns with 95% and 99% bands:

n_steps = 21

f = best_model.forecast(steps=n_steps)

# approximate standard error, assuming normality of the residuals
std_err = best_model.resid.std()

z_95 = 1.96
lower_ci_95 = f - z_95 * std_err
upper_ci_95 = f + z_95 * std_err

z_99 = 2.576
lower_ci_99 = f - z_99 * std_err
upper_ci_99 = f + z_99 * std_err

idx = pd.date_range(df.index[-1], periods=n_steps, freq='D')
fc_95 = pd.DataFrame({'forecast': f,
                      'lower_ci_95': lower_ci_95,
                      'upper_ci_95': upper_ci_95}, index=idx)

Two critical observations. First: these bands use a constant standard error, the same for day 1 as for day 21. It is a convenient but optimistic approximation, because the real uncertainty of a multi-step forecast grows with the horizon. Second: in a stationary AR the forecast converges rapidly to the mean of the series. If you zoom in, the 21-day return forecasts are an almost flat line near zero surrounded by an enormous band. That alone tells you how much directional signal there really is.

To evaluate on the test stretch we use predict with absolute indices:

predictions = best_model.predict(start=len(train),
                                 end=len(train) + len(test) - 1,
                                 dynamic=False)

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

dynamic=False means that at each step actual observations are used as lags; with dynamic=True the model feeds on its own forecasts and the error accumulates fast.

But there is an underlying problem: that model was estimated once, at the start of the test set, and never sees a new data point again. Obviously we want to rebuild the model each day as new data arrive.

Walk-forward validation. The statsmodels API does not make it easy to update the model with each observation. One option is to retrain every day, which is valid but computationally expensive. The alternative the professor proposes is to use the learned coefficients and compute the forecast by hand, keeping a history of previous observations:

# choose the optimal lag by BIC
best_bic = np.inf
best_lag = None
for lag in range(1, 51):
    result = AutoReg(train, lags=lag).fit()
    if result.bic < best_bic:
        best_bic = result.bic
        best_lag = lag

model_fit = AutoReg(train, lags=best_lag).fit()
window = len(model_fit.model._lags)
coef = model_fit.params

# walk forward over the test set
history = list(train[-window:])
predictions = []

for t in test:
    lag = history[-window:]
    yhat = coef[0]                       # intercept (trend='c' by default)
    for d in range(window):
        yhat += coef[d + 1] * lag[window - d - 1]
    predictions.append(yhat)
    history.append(t)                    # the ACTUAL value is appended, not the forecast

rmse = np.sqrt(mean_squared_error(test, predictions))
print('RMSE:', rmse)

The line history.append(t) is the heart of the method: the actual observation is appended after the forecast has been made, never before. That ordering is what prevents lookahead bias. Reverse it and your backtest will be predicting the past, and your results will be spectacular and useless.