Forecasting series temporales con python y scikitlearn

Forecasting series temporales con Python y Scikit-learn

Joaquín Amat Rodrigo
Febrero, 2021 (última actualización Mayo 2021)

Introducción


Una serie temporal (time series) es una sucesión de datos ordenados cronológicamente, espaciados a intervalos iguales o desiguales. El proceso de Forecasting consiste en predecir el valor futuro de una serie temporal, bien modelando la serie temporal únicamente en función de su comportamiento pasado (autorregresivo) o empleando otras variables externas a la serie temporal.

A lo largo de este documento, se describe cómo utilizar modelos de regresión de Scikit-learn para realizar forecasting sobre series temporales. En concreto, se hace uso de Skforecast, una librería sencilla (en construcción) que contiene las clases y funciones necesarias para adaptar cualquier modelo de regresón de Scikit-learn a problemas de forecasting.

Multi-Step Time Series Forecasting


Cuando se trabaja con series temporales, raramente se quiere predecir solo el siguiente elemento de la serie ($t_{+1}$), sino todo un intervalo futuro o un punto alejado en el tiempo ($t_{+n}$). A cada paso de predicción se le conoce como step.

Existen varias estrategias que permiten generar este tipo de predicciones múltiples.

Recursive multi-step forecasting

Dado que, para predecir el momento $t_{n}$ se necesita el valor de $t_{n-1}$, y $t_{n-1}$ se desconoce, es necesario hacer predicciones recursivas en las que, cada nueva predicción, se basa en la predicción anterior. A este proceso se le conoce como recursive forecasting o recursive multi-step forecasting.

La principal adaptación que se necesita hacer para aplicar modelos de scikit learn a problemas de recursive multi-step forecasting es transformar la serie temporal en un matriz en la que, cada valor, está asociado a la ventana temporal (lags) que le preceden. Esta estrategia de forecasting pueden generarse fácilmente con las clases ForecasterAutoreg y ForecasterCustom de la librería skforecast.

forecasting-python

Tranformación de una serie temporal en una matriz de 5 lags y un vector con el valor de la serie que sigue a cada fila de la matriz.

Direct multi-step forecasting

Esta estrategia consiste en entrenar un modelo distinto para cada step. Por ejemplo, si se quieren predecir los siguientes 5 valores de una serie temporal, se entrenan 5 modelos distintos, uno para cada step. Como resultado, las predicciones son independientes unas de otras. Esta estrategia de forecasting pueden generarse fácilmente con la clase ForecasterAutoregMultiOutput de la librería skforecast.

Multiple output forecasting

Determinados modelos son capaces de predecir de forma simultánea varios valores de una secuencia (one-shot). Un ejemplo de modelo con esta capacidad son las redes neuronales LSTM.

Forecasting autorregresivo recursivo


Se dispone de una serie temporal con el gasto mensual (millones de dólares) en fármacos con corticoides que tuvo el sistema de salud Australiano entre 1991 y 2008. Se pretende crear un modelo autoregresivo capaz de predecir el futuro gasto mensual. Para poder predecir múltiples valores a futuro, se emplea la estrategia recursiva (recursive multi-step forecasting).

Librerías


Las librerías utilizadas en este documento son:

In [17]:
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
%matplotlib inline

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')

Además de las anteriores, se utiliza skforecast, una librería sencilla (en construcción) que contiene las clases y funciones necesarias para adaptar cualquier modelo de regresón de scikit learn a problemas de forecasting.

Puede instalarse directamente desde Github con el comando:

!pip install git+https://github.com/JoaquinAmatRodrigo/skforecast@v0.1.8.1

Última versión (inestable):

!pip install git+https://github.com/JoaquinAmatRodrigo/skforecast#master

In [18]:
# Modelado y Forecasting
# ==============================================================================
#!pip install git+https://github.com/JoaquinAmatRodrigo/skforecast@v0.1.8.1

from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error

from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterCustom import ForecasterCustom
from skforecast.ForecasterAutoregMultiOutput import ForecasterAutoregMultiOutput
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import time_series_spliter
from skforecast.model_selection import cv_forecaster
from skforecast.model_selection import backtesting_forecaster_intervals

Datos


Los datos empleados en los ejemplos de este document se han obtenido del magnífico libro Forecasting: Principles and Practice by Rob J Hyndman and George Athanasopoulos.

In [19]:
# Descarga de datos
# ==============================================================================
url = 'https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/data/h2o.csv'
datos = pd.read_csv(url, sep=',')

La columna fecha se ha almacenado como string. Para convertirla en datetime, se emplea la función pd.to_datetime(). Una vez en formato datetime, y para hacer uso de las funcionalidades de pandas, se establece como índice. Además, dado que los datos son mensuales, se indica la frecuencia (Monthly Started 'MS').

In [20]:
# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y/%m/%d')
datos = datos.set_index('fecha')
datos = datos.rename(columns={'x': 'y'})
datos = datos.asfreq('MS')
datos = datos['y']
datos = datos.sort_index()

Se verifica que la serie temporal está completa.

In [21]:
# Verificar que un índice temporal está completo
# ==============================================================================
(datos.index == pd.date_range(start=datos.index.min(),
                              end=datos.index.max(),
                              freq=datos.index.freq)).all()
Out[21]:
True
In [22]:
# Completar huecos en un índice temporal
# ==============================================================================
# datos.asfreq(freq=''30min'', fill_value=np.nan)

Se utilizan los últimos 36 meses como conjunto de test para evaluar la capacidad predictiva del modelo.

In [23]:
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]

fig, ax=plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
ax.legend();

ForecasterAutoreg


Se crea y entrena un modelo ForecasterAutoreg a partir de un regresor RandomForestRegressor y una ventana temporal de 6 lags. Esto último significa que, para cada predicción, se utilizan como predictores los 6 meses anteriores.

In [24]:
# Crear y entrenar forecaster
# ==============================================================================
forecaster_rf = ForecasterAutoreg(
                    regressor=RandomForestRegressor(random_state=123),
                    lags=6
                )

forecaster_rf.fit(y=datos_train)

forecaster_rf
Out[24]:
=======================ForecasterAutoreg=======================
Regressor: RandomForestRegressor(random_state=123)
Lags: [1 2 3 4 5 6]
Exogenous variable: False
Parameters: {'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'mse', 'max_depth': None, 'max_features': 'auto', 'max_leaf_nodes': None, 'max_samples': None, 'min_impurity_decrease': 0.0, 'min_impurity_split': None, 'min_samples_leaf': 1, 'min_samples_split': 2, 'min_weight_fraction_leaf': 0.0, 'n_estimators': 100, 'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0, 'warm_start': False}

Predicciones


Una vez entrenado el modelo, se predicen los datos de test (36 meses a futuro).

In [25]:
# Predicciones
# ==============================================================================
steps = 36
predicciones = forecaster_rf.predict(steps=steps)
# Se añade el índice temporal a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)
predicciones.head()
Out[25]:
fecha
2005-07-01    0.866263
2005-08-01    0.874688
2005-09-01    0.951697
2005-10-01    0.991223
2005-11-01    0.952589
Freq: MS, dtype: float64
In [26]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();

Error predicción


Se cuantifica el error que comete el modelo en sus predicciones. En este caso, se emplea como métrica el mean squared error (mse).

In [27]:
# Error
# ==============================================================================
error_mse = mean_squared_error(
                y_true = datos_test,
                y_pred = predicciones
            )
print(f"Error de test (mse): {error_mse}")
Error de test (mse): 0.0665147739321922

Ajuste de hiperparámetros (tuning)


El ForecasterAutoreg entrenado ha utilizado una ventana temporal de 6 lags y un modelo Random Forest con los hiperparámetros por defecto. Sin embargo, no hay ninguna razón por la que estos valores sean los más adecuados.

Para identificar la mejor combinación de lags e hiperparámetros, la librería skforecast dispone de la estrategia de validación cruzada temporal y de Backtesting.

In [28]:
# Grid search de hiperparámetros
# ==============================================================================
forecaster_rf = ForecasterAutoreg(
                    regressor = RandomForestRegressor(random_state=123),
                    lags      = 12 # Este valor será remplazado en el grid search
                 )

# Hiperparámetros del regresor
param_grid = {'n_estimators': [100, 500],
              'max_depth': [3, 5, 10]}

# Lags utilizados como predictores
lags_grid = [10, 20]

resultados_grid = grid_search_forecaster(
                        forecaster  = forecaster_rf,
                        y           = datos_train,
                        param_grid  = param_grid,
                        lags_grid   = lags_grid,
                        steps       = 10,
                        method      = 'cv',
                        metric      = 'neg_mean_squared_error',
                        initial_train_size    = int(len(datos_train)*0.5),
                        allow_incomplete_fold = False,
                        return_best = True,
                        verbose     = False
                   )
loop lags_grid:   0%|          | 0/2 [00:00<?, ?it/s]
loop param_grid:   0%|          | 0/6 [00:00<?, ?it/s]
loop param_grid:  17%|█▋        | 1/6 [00:01<00:06,  1.36s/it]
loop param_grid:  33%|███▎      | 2/6 [00:07<00:17,  4.46s/it]
loop param_grid:  50%|█████     | 3/6 [00:09<00:09,  3.07s/it]
loop param_grid:  67%|██████▋   | 4/6 [00:16<00:09,  4.61s/it]
loop param_grid:  83%|████████▎ | 5/6 [00:17<00:03,  3.49s/it]
loop param_grid: 100%|██████████| 6/6 [00:25<00:00,  4.79s/it]
loop lags_grid:  50%|█████     | 1/2 [00:25<00:25, 25.19s/it] 
loop param_grid:   0%|          | 0/6 [00:00<?, ?it/s]
loop param_grid:  17%|█▋        | 1/6 [00:01<00:07,  1.44s/it]
loop param_grid:  33%|███▎      | 2/6 [00:08<00:18,  4.71s/it]
loop param_grid:  50%|█████     | 3/6 [00:09<00:09,  3.26s/it]
loop param_grid:  67%|██████▋   | 4/6 [00:17<00:09,  4.93s/it]
loop param_grid:  83%|████████▎ | 5/6 [00:19<00:03,  3.73s/it]
loop param_grid: 100%|██████████| 6/6 [00:27<00:00,  5.28s/it]
loop lags_grid: 100%|██████████| 2/2 [00:52<00:00, 26.28s/it] 
2021-05-14 21:34:17,640 root       INFO  Refitting `forecaster` using the best found parameters and the whole data set: 
lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20] 
params: {'max_depth': 10, 'n_estimators': 500}

In [29]:
# Resultados Grid Search
# ==============================================================================
resultados_grid
Out[29]:
lags params metric
11 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 10, 'n_estimators': 500} 0.005271
10 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 10, 'n_estimators': 100} 0.005331
9 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 5, 'n_estimators': 500} 0.005354
8 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 5, 'n_estimators': 100} 0.005514
7 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 3, 'n_estimators': 500} 0.005744
6 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 3, 'n_estimators': 100} 0.005821
5 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] {'max_depth': 10, 'n_estimators': 500} 0.026603
4 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] {'max_depth': 10, 'n_estimators': 100} 0.028092
2 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] {'max_depth': 5, 'n_estimators': 100} 0.028693
3 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] {'max_depth': 5, 'n_estimators': 500} 0.028870
1 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] {'max_depth': 3, 'n_estimators': 500} 0.029387
0 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] {'max_depth': 3, 'n_estimators': 100} 0.030020

Los mejores resultados se obtienen si se utiliza una ventana temporal de 20 lags y una configuración de random forest {'max_depth': 10, 'n_estimators': 500}.

Modelo final


Finalmente, se entrena de nuevo un ForecasterAutoreg con la configuración óptima encontrada mediante validación. Este paso no es necesario si se indica return_best = True en la función grid_search_forecaster().

In [30]:
# Crear y entrenar forecaster con mejores hiperparámetros
# ==============================================================================
regressor = RandomForestRegressor(max_depth=10, n_estimators=500, random_state=123)

forecaster_rf = ForecasterAutoreg(
                    regressor = regressor,
                    lags      = 20
                )

forecaster_rf.fit(y=datos_train)
In [31]:
# Predicciones
# ==============================================================================
predicciones = forecaster_rf.predict(steps=steps)
# Se añade el índice a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)
In [32]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
In [33]:
# Error
# ==============================================================================
error_mse = mean_squared_error(
                y_true = datos_test,
                y_pred = predicciones
            )
print(f"Error de test (mse) {error_mse}")
Error de test (mse) 0.003977043111878048

Con la combinación óptima de hiperparámetros se consigue reducir notablemente el error de test.

Importancia predictores


Dado que el objeto ForecasterAutoreg utiliza modelos scikit-learn, una vez entrenado, se puede acceder a la importancia de los predictores. Cuando el regresor empleado es un LinearRegression(), Lasso() o Ridge(), la importancia queda reflejada en los coeficientes del modelo, que se obtienen con el método get_coef(). En regresores GradientBoostingRegressor() o RandomForestRegressor(), la importancia de los predictores está basada en la reducción de impureza y es accesible mediante el método get_feature_importances(). En ambos casos, el orden devuelto se corresponde con el de los lags.

In [34]:
# Importancia predictores
# ==============================================================================
impotancia = forecaster_rf.get_feature_importances()
dict(zip(forecaster_rf.lags, impotancia))
Out[34]:
{1: 0.012553886713487061,
 2: 0.08983951332807713,
 3: 0.010659102591406495,
 4: 0.002089457758242796,
 5: 0.0020659673443678165,
 6: 0.0025603607878369916,
 7: 0.002581234419075128,
 8: 0.00692666663137377,
 9: 0.011527319810863246,
 10: 0.02536813007770579,
 11: 0.0179467896103424,
 12: 0.7756201041018168,
 13: 0.0031870877064889024,
 14: 0.014718107972907108,
 15: 0.00787688226587692,
 16: 0.003480591965198101,
 17: 0.0027116142341500264,
 18: 0.0020417148550876244,
 19: 0.0018876654387458749,
 20: 0.004357802386950057}

Comparación con statsmodels


La librería statsmodels dispone de algunos de los principales modelos de forecasting (AutoReg, AR, ARMA, ARIMA). Se verifica que se obtienen los mismos resultados con AutoReg que con la clase ForecasterAutoreg si se utiliza un modelo de regresor lineal LinearRegression().

In [35]:
# Modelo autorregresivo lineal statsmodels
# ==============================================================================
from statsmodels.tsa.ar_model import AutoReg
lags = 15

modelo_ar = AutoReg(datos_train, lags=lags)
res = modelo_ar.fit()
predicciones_statsmodels = res.predict(start=datos_test.index[0], end=datos_test.index[-1])

# Modelo autorregresivo lineal Forecaster
# ==============================================================================
regressor = LinearRegression()
forecaster = ForecasterAutoreg(regressor=regressor, lags=lags)
forecaster.fit(y=datos_train)
predicciones_forecaster = forecaster.predict(steps=36)

# Verificación de que las predicciones de ambos modelos son iguales
# ==============================================================================
print(np.allclose(predicciones_statsmodels.values, predicciones_forecaster))

# Verificación de que los coeficients de ambos modelos son iguales
# ==============================================================================
print(np.allclose(res.params.values[1:], forecaster.get_coef()))
True
True

Forecasting autorregresivo recursivo con variables exógenas


En el ejemplo anterior, se han utilizado como predictores únicamente lags de la propia variable predicha. En ciertos escenarios, es posible disponer de información sobre otras variables, cuyo valor a futuro se conoce, y pueden servir como predictoreres adicionales en el modelo.

Siguiendo con el ejemplo anterior, se simula una nueva variable cuyo comportamiento está correlacionado con la serie temporal modelada y que, por lo tanto, se quiere incorporar como predictor.

Datos

In [36]:
# Descarga de datos
# ==============================================================================
url = 'https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/data/h2o.csv'
datos = pd.read_csv(url, sep=',')

# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y/%m/%d')
datos = datos.set_index('fecha')
datos = datos.rename(columns={'x': 'y'})
datos = datos.asfreq('MS')
datos = datos['y']
datos = datos.sort_index()
In [37]:
# Datos
# ==============================================================================
datos_exog = datos.rolling(window=10, closed='right').mean() + 0.5
datos_exog = datos_exog[10:]
datos = datos[10:]

fig, ax = plt.subplots(figsize=(9, 4))
datos.plot(ax=ax, label='y')
datos_exog.plot(ax=ax, label='variable exógena')
ax.legend();
In [38]:
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]

datos_exog_train = datos_exog[:-steps]
datos_exog_test  = datos_exog[-steps:]

ForecasterAutoreg

In [39]:
# Crear y entrenar forecaster
# ==============================================================================
forecaster_rf = ForecasterAutoreg(
                    regressor = RandomForestRegressor(random_state=123),
                    lags      = 8
                )

forecaster_rf.fit(y=datos_train, exog=datos_exog_train)

Predicciones


Si el ForecasterAutoreg se entrena con una variable exógena, hay que pasarle el valor de esta variable al predict(). Por lo tanto, solo es aplicable a escenarios en los que se dispone de información a futuro de la variable exógena.

In [40]:
# Predicciones
# ==============================================================================
predicciones = forecaster_rf.predict(steps=steps, exog=datos_exog_test)
# Se añade el índice a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)
In [41]:
# Gráfico
# ==============================================================================
fig, ax=plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();

Error predicciones

In [42]:
# Error
# ==============================================================================
error_mse = mean_squared_error(
                y_true = datos_test,
                y_pred = predicciones
            )
print(f"Error de test (mse): {error_mse}")
Error de test (mse): 0.03948116514043873

Tuning del modelo

In [43]:
# Grid search de hiperparámetros
# ==============================================================================
forecaster_rf = ForecasterAutoreg(
                    regressor = RandomForestRegressor(random_state=123),
                    lags      = 12 # Este valor será remplazado en el grid search
                 )

param_grid = {'n_estimators': [50, 100, 500],
              'max_depth': [3, 5, 10]}

lags_grid = [5, 12, 20]

resultados_grid = grid_search_forecaster(
                        forecaster  = forecaster_rf,
                        y           = datos_train,
                        exog        = datos_exog_train,
                        param_grid  = param_grid,
                        lags_grid = lags_grid,
                        steps       = 10,
                        method      = 'cv',
                        metric      = 'neg_mean_squared_error',
                        initial_train_size    = int(len(datos_train)*0.5),
                        allow_incomplete_fold = False,
                        return_best = True,
                        verbose     = False
                    )
loop lags_grid:   0%|          | 0/3 [00:00<?, ?it/s]
loop param_grid:   0%|          | 0/9 [00:00<?, ?it/s]
loop param_grid:  11%|█         | 1/9 [00:00<00:04,  1.66it/s]
loop param_grid:  22%|██▏       | 2/9 [00:01<00:06,  1.05it/s]
loop param_grid:  33%|███▎      | 3/9 [00:07<00:18,  3.15s/it]
loop param_grid:  44%|████▍     | 4/9 [00:08<00:10,  2.15s/it]
loop param_grid:  56%|█████▌    | 5/9 [00:09<00:07,  1.81s/it]
loop param_grid:  67%|██████▋   | 6/9 [00:15<00:09,  3.18s/it]
loop param_grid:  78%|███████▊  | 7/9 [00:15<00:04,  2.35s/it]
loop param_grid:  89%|████████▉ | 8/9 [00:17<00:01,  2.00s/it]
loop param_grid: 100%|██████████| 9/9 [00:23<00:00,  3.32s/it]
loop lags_grid:  33%|███▎      | 1/3 [00:23<00:46, 23.32s/it] 
loop param_grid:   0%|          | 0/9 [00:00<?, ?it/s]
loop param_grid:  11%|█         | 1/9 [00:00<00:05,  1.59it/s]
loop param_grid:  22%|██▏       | 2/9 [00:01<00:06,  1.03it/s]
loop param_grid:  33%|███▎      | 3/9 [00:07<00:19,  3.21s/it]
loop param_grid:  44%|████▍     | 4/9 [00:08<00:11,  2.21s/it]
loop param_grid:  56%|█████▌    | 5/9 [00:09<00:07,  1.87s/it]
loop param_grid:  67%|██████▋   | 6/9 [00:15<00:10,  3.34s/it]
loop param_grid:  78%|███████▊  | 7/9 [00:16<00:04,  2.48s/it]
loop param_grid:  89%|████████▉ | 8/9 [00:17<00:02,  2.11s/it]
loop param_grid: 100%|██████████| 9/9 [00:24<00:00,  3.47s/it]
loop lags_grid:  67%|██████▋   | 2/3 [00:47<00:23, 23.93s/it] 
loop param_grid:   0%|          | 0/9 [00:00<?, ?it/s]
loop param_grid:  11%|█         | 1/9 [00:00<00:05,  1.54it/s]
loop param_grid:  22%|██▏       | 2/9 [00:01<00:07,  1.01s/it]
loop param_grid:  33%|███▎      | 3/9 [00:08<00:19,  3.33s/it]
loop param_grid:  44%|████▍     | 4/9 [00:08<00:11,  2.29s/it]
loop param_grid:  56%|█████▌    | 5/9 [00:10<00:07,  1.94s/it]
loop param_grid:  67%|██████▋   | 6/9 [00:16<00:10,  3.54s/it]
loop param_grid:  78%|███████▊  | 7/9 [00:17<00:05,  2.62s/it]
loop param_grid:  89%|████████▉ | 8/9 [00:18<00:02,  2.23s/it]
loop param_grid: 100%|██████████| 9/9 [00:25<00:00,  3.66s/it]
loop lags_grid: 100%|██████████| 3/3 [01:13<00:00, 24.42s/it] 
2021-05-14 21:35:35,081 root       INFO  Refitting `forecaster` using the best found parameters and the whole data set: 
lags: [ 1  2  3  4  5  6  7  8  9 10 11 12] 
params: {'max_depth': 10, 'n_estimators': 500}

In [44]:
# Resultados Grid Search
# ==============================================================================
resultados_grid.head()
Out[44]:
lags params metric
17 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] {'max_depth': 10, 'n_estimators': 500} 0.005345
23 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 5, 'n_estimators': 500} 0.005403
26 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'max_depth': 10, 'n_estimators': 500} 0.005420
14 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] {'max_depth': 5, 'n_estimators': 500} 0.005423
16 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] {'max_depth': 10, 'n_estimators': 100} 0.005507

Los mejores resultados se obtienen utilizando una ventana temporal de 12 lags y una configuración de random forest {'max_depth': 10, 'n_estimators': 500}.

Modelo final


Como se ha indicado return_best = True en el grid_search_forecaster(), tras la búsqueda, el objeto ForecasterAutoreg ha sido modificado y entrenado con la mejor combinación encontrada.

In [45]:
# Predicciones
# ==============================================================================
predicciones = forecaster_rf.predict(steps=steps, exog=datos_exog_test)
# Se añade el índice a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)

# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
In [46]:
# Error
# ==============================================================================
error_mse = mean_squared_error(y_true = datos_test, y_pred = predicciones)
print(f"Error de test (mse) {error_mse}")
Error de test (mse) 0.004209073816980198

Forecasting autorregresivo recursivo con predictores custom


En determinados escenarios, puede ser interesante incorporar otras características de la serie temporal además de los lags, por ejemplo, la media movil de los últimos n valores puede servir para capturar la tendencia de la serie.

La clase ForecasterCustom se comporta de forma muy similar a la clase ForecasterAutoreg vista en los apartados anteriores pero con la diferencia de que, es el usuario, quien define la función empleada para crear los predictores.

Se repite el primer ejemplo del documento, predecir los últimos 36 meses de la serie temporal, pero esta vez, utilizando como predictores los 10 primeros lags y la media móvil de los últimos 20 meses.

Datos

In [47]:
# Descarga de datos
# ==============================================================================
url = 'https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/data/h2o.csv'
datos = pd.read_csv(url, sep=',')
In [48]:
# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y/%m/%d')
datos = datos.set_index('fecha')
datos = datos.rename(columns={'x': 'y'})
datos = datos.asfreq('MS')
datos = datos['y']
datos = datos.sort_index()
In [49]:
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]

ForecasterCustom


Se crea y entrena ForecasterCustom a partir de un regresor RandomForestRegressor. Para crear los predictores, se emplea la función create_predictors() que calcula los primeros 10 lags y la media móvil de los últimos 20 valores.

In [50]:
# Función para calcular los predictores a partir de la serie temporal
# ==============================================================================
def create_predictors(y):
    '''
    Crear los primeros 10 lags.
    Calcular la media móvil de los últimos 20 valores.
    '''
    
    X_train = pd.DataFrame({'y':y.copy()})
    for i in range(0, 10):
        X_train[f'lag_{i+1}'] = X_train['y'].shift(i)
        
    X_train['moving_avg'] = X_train['y'].rolling(20).mean()
    
    X_train = X_train.drop(columns='y').tail(1).to_numpy()  
    
    return X_train  

Al crear el forecaster, el argumento window_size debe ser un valor como mínimo tan grande como la ventana que utiliza la función que crea los predictores. En este caso, 20.

In [51]:
# Crear y entrenar forecaster
# ==============================================================================
forecaster_rf = ForecasterCustom(
                    regressor      = RandomForestRegressor(random_state=123),
                    fun_predictors = create_predictors,
                    window_size    = 20
                )

forecaster_rf.fit(y=datos_train)
forecaster_rf
Out[51]:
=======================ForecasterCustom=======================
Regressor: RandomForestRegressor(random_state=123)
Predictors created with: create_predictors
Window size: 20
Exogenous variable: False
Parameters: {'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'mse', 'max_depth': None, 'max_features': 'auto', 'max_leaf_nodes': None, 'max_samples': None, 'min_impurity_decrease': 0.0, 'min_impurity_split': None, 'min_samples_leaf': 1, 'min_samples_split': 2, 'min_weight_fraction_leaf': 0.0, 'n_estimators': 100, 'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0, 'warm_start': False}

Predicciones

In [52]:
# Predicciones
# ==============================================================================
steps = 36
predicciones = forecaster_rf.predict(steps=steps)
# Se añade el índice temporal a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)
In [53]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();

Error predicción

In [54]:
# Error
# ==============================================================================
error_mse = mean_squared_error(
                y_true = datos_test,
                y_pred = predicciones
            )
print(f"Error de test (mse): {error_mse}")
Error de test (mse): 0.04487765885818191

Ajuste de hiperparámetros (tuning)


Al utilizar la función grid_search_forecaster() con un ForecasterCustom, no se indica el argumento lags_grid.

In [55]:
# Grid search de hiperparámetros
# ==============================================================================
forecaster_rf = ForecasterCustom(
                    regressor      = RandomForestRegressor(random_state=123),
                    fun_predictors = create_predictors,
                    window_size    = 20
                )

# Hiperparámetros del regresor
param_grid = {'n_estimators': [100, 500],
              'max_depth': [3, 5, 10]}

resultados_grid = grid_search_forecaster(
                        forecaster  = forecaster_rf,
                        y           = datos_train,
                        param_grid  = param_grid,
                        steps       = 10,
                        method      = 'cv',
                        metric      = 'neg_mean_squared_error',
                        initial_train_size    = int(len(datos_train)*0.5),
                        allow_incomplete_fold = True,
                        return_best = True,
                        verbose     = False
                    )
loop lags_grid:   0%|          | 0/1 [00:00<?, ?it/s]
loop param_grid:   0%|          | 0/6 [00:00<?, ?it/s]
loop param_grid:  17%|█▋        | 1/6 [00:06<00:34,  6.97s/it]
loop param_grid:  33%|███▎      | 2/6 [00:19<00:41, 10.38s/it]
loop param_grid:  50%|█████     | 3/6 [00:26<00:26,  8.83s/it]
loop param_grid:  67%|██████▋   | 4/6 [00:40<00:21, 10.60s/it]
loop param_grid:  83%|████████▎ | 5/6 [00:47<00:09,  9.34s/it]
loop param_grid: 100%|██████████| 6/6 [01:00<00:00, 10.78s/it]
loop lags_grid: 100%|██████████| 1/1 [01:00<00:00, 60.73s/it] 
2021-05-14 21:36:39,509 root       INFO  Refitting `forecaster` using the best found parameters and the whole data set: 
lags: custom predictors 
params: {'max_depth': 10, 'n_estimators': 500}

In [56]:
# Resultados Grid Search
# ==============================================================================
resultados_grid
Out[56]:
lags params metric
5 custom predictors {'max_depth': 10, 'n_estimators': 500} 0.022736
3 custom predictors {'max_depth': 5, 'n_estimators': 500} 0.022742
4 custom predictors {'max_depth': 10, 'n_estimators': 100} 0.023564
2 custom predictors {'max_depth': 5, 'n_estimators': 100} 0.024030
1 custom predictors {'max_depth': 3, 'n_estimators': 500} 0.025694
0 custom predictors {'max_depth': 3, 'n_estimators': 100} 0.026545

Modelo final

In [57]:
# Predicciones
# ==============================================================================
predicciones = forecaster_rf.predict(steps=steps)
# Se añade el índice a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)

# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
In [58]:
# Error
# ==============================================================================
error_mse = mean_squared_error(y_true = datos_test, y_pred = predicciones)
print(f"Error de test (mse) {error_mse}")
Error de test (mse) 0.044590618568342955

Direct multi-step forecasting


Los modelos ForecasterAutoreg y ForecasterCustom siguen una estrategia de predicción recursiva en la que, cada nueva predicción, se basa en la predicción anterior. Una alternativa es entrenar un modelo para cada uno de los steps que se desea predecir. Esta estrategia, normalmente conocida como direct multi-step forecasting, es computacionalmente más costosa que la recursiva puesto que requiere entrenar múltiples modelos. Sin embargo, en algunos escenarios, consigue mejores resultados. Este tipo de modelos pueden obtenerse con la clase ForecasterAutoregMultiOutput y pueden incluir también una o múltiples variables exógenas.

Datos

In [59]:
# Descarga de datos
# ==============================================================================
url = 'https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/data/h2o.csv'
datos = pd.read_csv(url, sep=',')

# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y/%m/%d')
datos = datos.set_index('fecha')
datos = datos.rename(columns={'x': 'y'})
datos = datos.asfreq('MS')
datos = datos['y']
datos = datos.sort_index()
In [60]:
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]

datos_exog_train = datos_exog[:-steps]
datos_exog_test  = datos_exog[-steps:]

ForecasterAutoregMultiOutput


A diferencia de cuando se utiliza ForecasterAutoreg o ForecasterCustom, en los modelos de tipo ForecasterAutoregMultiOutput hay que indicar, en el momento de su creación, el número de steps que se quieren predecir. Esto significa que, el número de predicciones obtenidas al ejecutar el método predict() es siempre el mismo.

In [61]:
# Grid search de hiperparámetros
# ==============================================================================
forecaster_rf = ForecasterAutoregMultiOutput(
                    regressor = Lasso(random_state=123),
                    steps     = 36,
                    lags      = 8 # Este valor será remplazado en el grid search
                )

param_grid = {'alpha': np.logspace(-5, 5, 10)}

lags_grid = [5, 12, 20]

resultados_grid = grid_search_forecaster(
                        forecaster  = forecaster_rf,
                        y           = datos_train,
                        param_grid  = param_grid,
                        lags_grid = lags_grid,
                        steps       = 36,
                        method      = 'cv',
                        metric      = 'neg_mean_squared_error',
                        initial_train_size    = int(len(datos_train)*0.5),
                        allow_incomplete_fold = False,
                        return_best = True,
                        verbose     = False
                    )
loop lags_grid:   0%|          | 0/3 [00:00<?, ?it/s]
loop param_grid:   0%|          | 0/10 [00:00<?, ?it/s]
loop param_grid:  20%|██        | 2/10 [00:00<00:00, 10.82it/s]
loop param_grid:  40%|████      | 4/10 [00:00<00:00, 10.93it/s]
loop param_grid:  60%|██████    | 6/10 [00:00<00:00, 11.01it/s]
loop param_grid:  80%|████████  | 8/10 [00:00<00:00, 11.03it/s]
loop param_grid: 100%|██████████| 10/10 [00:00<00:00, 11.12it/s]
loop lags_grid:  33%|███▎      | 1/3 [00:00<00:01,  1.10it/s]   
loop param_grid:   0%|          | 0/10 [00:00<?, ?it/s]
loop param_grid:  20%|██        | 2/10 [00:00<00:00, 10.97it/s]
loop param_grid:  40%|████      | 4/10 [00:00<00:00, 11.16it/s]
loop param_grid:  60%|██████    | 6/10 [00:00<00:00, 11.14it/s]
loop param_grid:  80%|████████  | 8/10 [00:00<00:00, 11.09it/s]
loop param_grid: 100%|██████████| 10/10 [00:00<00:00, 11.13it/s]
loop lags_grid:  67%|██████▋   | 2/3 [00:01<00:00,  1.10it/s]   
loop param_grid:   0%|          | 0/10 [00:00<?, ?it/s]
loop param_grid:  20%|██        | 2/10 [00:00<00:00, 10.75it/s]
loop param_grid:  40%|████      | 4/10 [00:00<00:00, 10.92it/s]
loop param_grid:  60%|██████    | 6/10 [00:00<00:00, 11.03it/s]
loop param_grid:  80%|████████  | 8/10 [00:00<00:00, 11.08it/s]
loop param_grid: 100%|██████████| 10/10 [00:00<00:00, 11.13it/s]
loop lags_grid: 100%|██████████| 3/3 [00:02<00:00,  1.10it/s]   
2021-05-14 21:36:45,300 root       INFO  Refitting `forecaster` using the best found parameters and the whole data set: 
lags: 5 
params: {'estimator__alpha': 0.0001291549665014884}

In [62]:
# Resultados Grid Search
# ==============================================================================
resultados_grid.head()
Out[62]:
lags params metric
1 5 {'estimator__alpha': 0.0001291549665014884} 0.015163
21 20 {'estimator__alpha': 0.0001291549665014884} 0.015163
11 12 {'estimator__alpha': 0.0001291549665014884} 0.015163
0 5 {'estimator__alpha': 1e-05} 0.015340
20 20 {'estimator__alpha': 1e-05} 0.015340

Los mejores resultados se obtienen utilizando una ventana temporal de 5 lags y una configuración de Lasso {'alpha': 0.000129}.

In [63]:
# Predicciones
# ==============================================================================
predicciones = forecaster_rf.predict()
# Se añade el índice a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)

# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
In [64]:
# Error
# ==============================================================================
error_mse = mean_squared_error(y_true = datos_test, y_pred = predicciones)
print(f"Error de test (mse) {error_mse}")
Error de test (mse) 0.0073864936455776565

Intervalos de predicción


Un intervalo de predicción define el intervalo dentro del cual es de esperar que se encuentre el verdadero valor de $y$ con una determinada probabilidad.

Rob J Hyndman y George Athanasopoulos, listan en su libro Forecasting: Principles and Practice mútiples formas de estimar intervalos de predicción, la mayoría los cuales requieren que los resudios (errores) del modelo se distribuyan de forma normal. Cuando no se puede asumir esta propiedad, se puede recurrir a bootstrapping, que solo asume que los residuos no están correlacionados. Este es el método utilizado en la librería skforecast para los modelos de tipo ForecasterAutoreg y ForecasterCustom.

In [65]:
# Descarga de datos
# ==============================================================================
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/data/h2o.csv')
datos = pd.read_csv(url, sep=',')

# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y/%m/%d')
datos = datos.set_index('fecha')
datos = datos.rename(columns={'x': 'y'})
datos = datos.asfreq('MS')
datos = datos['y']
datos = datos.sort_index()

# Split train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]
In [66]:
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
                    regressor=LinearRegression(),
                    lags=15
                )

forecaster.fit(y=datos_train)

# Intervalos de predicción
# ==============================================================================
predictions = forecaster.predict_interval(
                    steps    = steps,
                    interval = [1, 99],
                    n_boot   = 1000
              )

# Se añade índice datetime
predictions = pd.DataFrame(data=predictions, index=datos_test.index)

# Error de predicción
# ==============================================================================
error_mse = mean_squared_error(
                y_true = datos_test,
                y_pred = predictions.iloc[:, 0]
            )
print(f"Test error (mse): {error_mse}")

# Gráfico
# ==============================================================================
fig, ax=plt.subplots(figsize=(9, 4))
#datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predictions.iloc[:, 0].plot(ax=ax, label='predictions')
ax.fill_between(predictions.index,
                predictions.iloc[:, 1],
                predictions.iloc[:, 2],
                alpha=0.5)
ax.legend();
Test error (mse): 0.011051937043503771
In [67]:
# Backtest con intervalos de predicción
# ==============================================================================
n_test = 36*3 + 1
datos_train = datos[:-n_test]
datos_test  = datos[-n_test:]

steps = 36
regressor = LinearRegression()
forecaster = ForecasterAutoreg(regressor=regressor, lags=15)

metric, predictions = backtesting_forecaster_intervals(
                            forecaster = forecaster,
                            y          = datos,
                            initial_train_size = len(datos_train),
                            steps      = steps,
                            metric     = 'neg_mean_squared_error',
                            interval            = [1, 99],
                            n_boot              = 100,
                            in_sample_residuals = True,
                            verbose             = True
                       )

print(metric)

# Se añade índice datetime
predictions = pd.DataFrame(data=predictions, index=datos_test.index)

# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
#datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predictions.iloc[:, 0].plot(ax=ax, label='predictions')
ax.fill_between(predictions.index,
                predictions.iloc[:, 1],
                predictions.iloc[:, 2],
                alpha=0.5)
ax.legend();
Number of observations used for training: 95
Number of observations used for testing: 109
    Number of folds: 4
    Number of steps per fold: 36
    Last fold only includes 1 observations
[0.02150972]

Información de sesión

In [68]:
from sinfo import sinfo
sinfo()
-----
PIL         8.0.1
gif         NA
matplotlib  3.3.2
numpy       1.19.5
pandas      1.2.3
sinfo       0.3.1
skforecast  0.1.8.1
sklearn     0.24.1
statsmodels 0.11.1
-----
IPython             7.23.1
jupyter_client      6.1.11
jupyter_core        4.7.1
notebook            6.2.0
-----
Python 3.7.9 (default, Aug 31 2020, 12:42:55) [GCC 7.3.0]
Linux-5.4.0-1048-aws-x86_64-with-debian-buster-sid
4 logical CPU cores, x86_64
-----
Session information updated at 2021-05-14 21:36

Bibliografía


Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. libro

Time Series Analysis and Forecasting with ADAM Ivan Svetunkov libro

Python Data Science Handbook by Jake VanderPlas libro

Python for Finance: Mastering Data-Driven Finance libro

Markus Löning, Anthony Bagnall, Sajaysurya Ganesh, Viktor Kazakov, Jason Lines, Franz Király (2019): “sktime: A Unified Interface for Machine Learning with Time Series”

Markus Löning, Tony Bagnall, Sajaysurya Ganesh, George Oastler, Jason Lines, ViktorKaz, …, Aadesh Deshmukh (2020). alan-turing-institute/sktime. Zenodo. http://doi.org/10.5281/zenodo.3749000

skforecast

¿Cómo citar este documento?

Forecasting series temporales con Python y Scikitlearn by Joaquín Amat Rodrigo, available under a Attribution 4.0 International (CC BY 4.0) at https://www.cienciadedatos.net/py27-forecasting-series-temporales-python-scikitlearn.html


¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊


Creative Commons Licence
This work by Joaquín Amat Rodrigo is licensed under a Creative Commons Attribution 4.0 International License.