Saltar al contenido

Modelado de series temporales con Scikit, Pandas y Numpy

septiembre 29, 2021

Modelado de series temporales con Scikit, Pandas y Numpy

Uso intuitivo de la estacionalidad para mejorar la precisión del modelo

Dr. Varshita Sher

14 de junio de 2020·9 min de lectura

¡Bienvenido a la Parte 2 del Análisis de series temporales! En esta publicación, trabajaremos en el modelado de datos de series de tiempo. Esta es una continuación de mi publicación anterior sobre datos de series temporales.

En nuestra publicación de blog anterior, hablamos sobre qué datos de series de tiempos, cómo formatear dichos datos para maximizar su utilidad y cómo manejar los datos faltantes. También aprendimos cómo volver a muestrear datos de series de tiempo por mes, semana, año, etc., y cómo calcular las medias continuas. Nos sumergimos en conceptos como tendencias, estacionalidad, diferenciación de primer orden y autocorrelación. Si está familiarizado con la mayoría de las cosas, está listo para comenzar. En caso de que necesite un repaso, puede hacer una búsqueda rápida en Google para estos temas o leer mi publicación anterior aquí.

Algunas palabras antes de comenzar:

Hay otros paquetes, indudablemente mejores, disponibles para los pronósticos de series de tiempo, como ARIMA o el software propietario de Facebook, Prophet. Sin embargo, este artículo se inspiró en la tarea para llevar a casa de una amiga que requería que usara solo Scikit, Numpy y Pandas (¡o se enfrentaría a la descalificación instantánea!).

Profundicemos en nuestro conjunto de datos

Trabajaremos con el conjunto de datos disponible públicamente Open Power System Data. Puede descargar los datos aquí. Contiene el consumo de electricidad, la producción de energía eólica y la producción de energía solar para 2006–2017.

url='https://raw.githubusercontent.com/jenfly/opsd/master/opsd_germany_daily.csv'
data = pd.read_csv(url,sep=",")

Después de configurar el Date columna como índice, así es como se ve nuestro conjunto de datos:

# to explicitly convert the date column to type DATETIME
data['Date'] = pd.to_datetime(data['Date'])
data = data.set_index('Date')
0*82pE15 BCQTYfMFe

Definición de la tarea de modelado

Objetivos de predicción

Nuestro objetivo es predecir Consumption (idealmente para futuras fechas no vistas) de este conjunto de datos de series de tiempo.

Conjunto de entrenamiento y prueba

Usaremos 10 años de datos para la capacitación, es decir, 2006–2016 y los datos del año pasado para las pruebas, es decir, 2017.

Medida de rendimiento

Para evaluar qué tan bueno es nuestro modelo, estaríamos usando R-cuadrado y Root Mean Squared Error (pero imprimiremos todas las métricas relevantes para que usted tome la última llamada).

Funciones auxiliares

Para imprimir todas las métricas de rendimiento relevantes para una tarea de regresión (como MAE y R-cuadrado), definiremos el regression_results función.

import sklearn.metrics as metrics
def regression_results(y_true, y_pred):
# Regression metrics
explained_variance=metrics.explained_variance_score(y_true, y_pred)
mean_absolute_error=metrics.mean_absolute_error(y_true, y_pred)
mse=metrics.mean_squared_error(y_true, y_pred)
mean_squared_log_error=metrics.mean_squared_log_error(y_true, y_pred)
median_absolute_error=metrics.median_absolute_error(y_true, y_pred)
r2=metrics.r2_score(y_true, y_pred)
print('explained_variance: ', round(explained_variance,4))
print('mean_squared_log_error: ', round(mean_squared_log_error,4))
print('r2: ', round(r2,4))
print('MAE: ', round(mean_absolute_error,4))
print('MSE: ', round(mse,4))
print('RMSE: ', round(np.sqrt(mse),4))

Ingeniería de funciones

Como línea de base, elegimos un modelo simplista, uno que predice el valor de consumo actual basado en

  • valor de consumo de ayer y;
  • diferencia entre el valor de consumo de ayer y anteayer.
# creating new dataframe from consumption column
data_consumption = data[['Consumption']]
# inserting new column with yesterday's consumption values
data_consumption.loc[:,'Yesterday'] =
data_consumption.loc[:,'Consumption'].shift()
# inserting another column with difference between yesterday and day before yesterday's consumption values.
data_consumption.loc[:,'Yesterday_Diff'] = data_consumption.loc[:,'Yesterday'].diff()
# dropping NAs
data_consumption = data_consumption.dropna()
1*RvqYTi5Gow5SheiPCMLMVQ

Definición de conjuntos de entrenamiento y prueba

X_train = data_consumption[:'2016'].drop(['Consumption'], axis = 1)
y_train = data_consumption.loc[:'2016', 'Consumption']
X_test = data_consumption['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption.loc['2017', 'Consumption']

Validación cruzada de datos de series temporales

Una pregunta que surge a menudo durante las entrevistas de ciencia de datos es: ¿Qué técnica de validación cruzada usaría en datos de series de tiempo?

Es posible que sienta la tentación de gravitar hacia la validación cruzada de K-Fold favorita de todos los tiempos (créanme, hasta hace poco, ¡no pregunten cuándo! No sabía que existían otras técnicas de CV además de K-fold). Desafortunadamente, esa no sería la respuesta correcta. La razón es que no tiene en cuenta que los datos de series de tiempo tienen algún orden natural y la aleatorización en la validación cruzada estándar de k veces no conserva ese orden.

Una mejor alternativa para la validación cruzada en datos de series de tiempo (que K-fold CV) es la estrategia Forward Chaining.

En el encadenamiento hacia adelante, digamos con 3 pliegues, el tren y los conjuntos de validación se ven así:

  • pliegue 1: entrenamiento [1], validación [2]
  • doblez 2: entrenamiento [1 2], validación [3]
  • pliegue 3: entrenamiento [1 2 3], validación [4]

donde 1, 2, 3, 4 representan el año. De esta manera, los conjuntos de entrenamiento sucesivos son superconjuntos de los anteriores.

Afortunadamente para nosotros, sklearn tiene una disposición para implementar dicha división de prueba de tren usando TimeSeriesSplit.

from sklearn.model_selection import TimeSeriesSplit

los TimeSerieSplit La función toma como entrada el número de divisiones. Dado que nuestros datos de entrenamiento tienen 11 años únicos (2006-2016), estaríamos estableciendo n_splits = 10. De esta manera, tenemos conjuntos de validación y entrenamiento prolijos:

  • pliegue 1: entrenamiento [2006], validación [2007]
  • doblez 2: entrenamiento [2006 2007], validación [2008]
  • pliegue 3: entrenamiento [2006 2007 2008], validación [2009]
  • pliegue 4: entrenamiento [2006 2007 2008 2009], validación [2010]
  • pliegue 5: entrenamiento [2006 2007 2008 2009 2010], validación [2011]
  • pliegue 6: entrenamiento [2006 2007 2008 2009 2010 2011], validación [2012]
  • pliegue 7: entrenamiento [2006 2007 2008 2009 2010 2011 2012], validación [2013]
  • pliegue 8: entrenamiento [2006 2007 2008 2009 2010 2011 2012 2013], validación [2014]
  • pliegue 9: entrenamiento [2006 2007 2008 2009 2010 2011 2012 2013 2014], validación [2015]
  • pliegue 10: entrenamiento [2006 2007 2008 2009 2010 2011 2012 2013 2014 2015], validación [2016]

Algoritmos de verificación puntual

# Spot Check Algorithmsmodels = []
models.append(('LR', LinearRegression()))
models.append(('NN', MLPRegressor(solver = 'lbfgs'))) #neural network
models.append(('KNN', KNeighborsRegressor()))
models.append(('RF', RandomForestRegressor(n_estimators = 10))) # Ensemble method - collection of many decision trees
models.append(('SVR', SVR(gamma='auto'))) # kernel = linear
# Evaluate each model in turn
results = []
names = []
for name, model in models:
# TimeSeries Cross validation
tscv = TimeSeriesSplit(n_splits=10)

cv_results = cross_val_score(model, X_train, y_train, cv=tscv, scoring='r2')
results.append(cv_results)
names.append(name)
print('%s: %f (%f)' % (name, cv_results.mean(), cv_results.std()))

# Compare Algorithms
plt.boxplot(results, labels=names)
plt.title('Algorithm Comparison')
plt.show()
1*8ywSH5 NZpPI7gQbHUKB6Q

Tanto KNN como RF funcionan igualmente bien. Pero personalmente prefiero RF ya que este modelo de conjunto (combina varios modelos ‘individuales’ (diversos) juntos y ofrece un poder de predicción superior) casi puede funcionar de inmediato y esa es una de las razones por las que son muy populares.

Hiperparámetros de búsqueda de cuadrícula

Hablé de la necesidad de hiperparámetros de búsqueda de cuadrícula en uno de mis artículos anteriores.

Una combinación óptima de hiperparámetros maximiza el rendimiento de un modelo sin generar un problema de alta varianza (sobreajuste).

El código de Python para realizar la búsqueda en la cuadrícula es el siguiente:

from sklearn.model_selection import GridSearchCVmodel = RandomForestRegressor()
param_search = {
'n_estimators': [20, 50, 100],
'max_features': ['auto', 'sqrt', 'log2'],
'max_depth' : [i for i in range(5,15)]
}
tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)
gsearch.fit(X_train, y_train)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_

Si observa el código anterior, hemos definido un personalizado goleador configurando scoring = rmse_scoreen lugar de utilizar una de las métricas de puntuación comunes definidas en sklearn. Definimos nuestro marcador personalizado de la siguiente manera:

from sklearn.metrics import make_scorerdef rmse(actual, predict):predict = np.array(predict)
actual = np.array(actual)
distance = predict - actualsquare_distance = distance ** 2mean_square_distance = square_distance.mean()score = np.sqrt(mean_square_distance)return scorermse_score = make_scorer(rmse, greater_is_better = False)

Comprobación del mejor rendimiento del modelo en los datos de prueba

y_true = y_test.values
y_pred = best_model.predict(X_test)
regression_results(y_true, y_pred)
1*HhpaHNvhmOke00Ief UHbg

Esto no está mal para empezar. Veamos si podemos mejorar aún más nuestro modelo.

Devoluciones de ingeniería de funciones

Hasta ahora, hemos estado usando valores en (t-1)th día para predecir valores en t fecha. Ahora, usemos también valores de (t-2)días para predecir el consumo:

# creating copy of original dataframe
data_consumption_2o = data_consumption.copy()
# inserting column with yesterday-1 values
data_consumption_2o['Yesterday-1'] = data_consumption_2o['Yesterday'].shift()
# inserting column with difference in yesterday-1 and yesterday-2 values.
data_consumption_2o['Yesterday-1_Diff'] = data_consumption_2o['Yesterday-1'].diff()
# dropping NAs
data_consumption_2o = data_consumption_2o.dropna()
1*othDvMvqoehEE i

Restablecimiento del tren y el equipo de prueba

X_train_2o = data_consumption_2o[:'2016'].drop(['Consumption'], axis = 1)
y_train_2o = data_consumption_2o.loc[:'2016', 'Consumption']
X_test = data_consumption_2o['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption_2o.loc['2017', 'Consumption']

Comprobación para ver si el «mejor» bosque aleatorio que utiliza predictores «nuevos» funciona mejor

model = RandomForestRegressor()
param_search = {
'n_estimators': [20, 50, 100],
'max_features': ['auto', 'sqrt', 'log2'],
'max_depth' : [i for i in range(5,15)]
}
tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)
gsearch.fit(X_train_2o, y_train_2o)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_
y_true = y_test.values
y_pred = best_model.predict(X_test)
regression_results(y_true, y_pred)

¡¡Una gran noticia!! Hemos reducido significativamente los valores de RMSE y MAE, mientras que el valor de R cuadrado también ha aumentado.

La ingeniería de funciones contraataca

Veamos si agregar el valor de la producción solar es beneficioso de alguna manera para predecir el consumo de electricidad.

data_consumption_2o_solar = data_consumption_2o.join(data[['Solar']])data_consumption_2o_solar = data_consumption_2o_solar.dropna()
1*dd1t4Jc0HmkC6uP BDdPNg

Restablecimiento del tren / prueba + GridSearch + Comprobación del rendimiento

X_train_2o_solar = data_consumption_2o_solar[:'2016'].drop(['Consumption'], axis = 1)
y_train_2o_solar = data_consumption_2o_solar.loc[:'2016', 'Consumption']
X_test = data_consumption_2o_solar['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption_2o_solar.loc['2017', 'Consumption']
model = RandomForestRegressor()
param_search = {
'n_estimators': [20, 50, 100],
'max_features': ['auto', 'sqrt', 'log2'],
'max_depth' : [i for i in range(5,15)]
}
tscv = TimeSeriesSplit(n_splits=5)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)
gsearch.fit(X_train_2o_solar, y_train_2o_solar)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_
y_true = y_test.values
y_pred = best_model.predict(X_test)
close