3. Predecir con AR: intervalos y walk-forward

Nivel 9: Series Temporales - ARIMA y GARCH

Con el mejor modelo ya seleccionado podemos generar una previsión a varios pasos. Aquí, 21 días de retornos del S&P 500 con bandas al 95% y al 99%:

n_steps = 21

f = best_model.forecast(steps=n_steps)

# error estandar aproximado, asumiendo normalidad de los residuos
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)

Dos observaciones críticas. Primera: estas bandas usan un error estándar constante, el mismo para el día 1 que para el día 21. Es una aproximación comoda pero optimista, porque la incertidumbre real de una previsión a varios pasos crece con el horizonte. Segunda: en un AR estacionario la previsión converge rápidamente a la media de la serie. Si haces zoom, las previsiones a 21 días de retornos son una línea casi plana cerca de cero rodeada por una banda enorme. Eso ya te dice cuánta señal direccional hay realmente.

Para evaluar sobre el tramo de test se usa predict con índices absolutos:

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 significa que en cada paso se usan observaciones reales como retardos; con dynamic=True el modelo se alimenta de sus propias predicciones y el error se acumula rápido. Ojo: esa distinción solo opera dentro de la muestra de ajuste. Como aquí el modelo se estimó únicamente con train, todo lo que va de len(train) en adelante es fuera de muestra y statsmodels lo calcula siempre de forma recursiva, así que ese RMSE no es el de una predicción a un paso, sino el de una previsión a len(test) pasos.

Pero hay un problema de fondo: ese modelo se estimó una sola vez, al principio del test, y no vuelve a ver un dato nuevo nunca más. Evidentemente queremos reconstruir el modelo cada día según llegan datos nuevos.

Validación walk-forward. La API de statsmodels no facilita actualizar el modelo con cada observación. Una opción es reentrenar cada día, válida pero cara en cómputo. La alternativa que propone el profesor es usar los coeficientes aprendidos y calcular la predicción a mano, manteniendo el historial de observaciones previas:

# elegir el retardo optimo por 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

# recorrido hacia delante sobre el test
history = list(train[-window:])
predictions = []

for t in test:
    lag = history[-window:]
    yhat = coef[0]                       # intercepto (trend='c' por defecto)
    for d in range(window):
        yhat += coef[d + 1] * lag[window - d - 1]
    predictions.append(yhat)
    history.append(t)                    # se incorpora el valor REAL, no la prediccion

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

La línea history.append(t) es el corazón del método: se añade la observación real después de haber hecho la predicción, nunca antes. Ese orden es lo que impide el sesgo de anticipación. Si lo inviertes, tu backtest predecira el pasado y tus resultados serán fantásticos e inútiles.