Mejore el rendimiento del pronóstico de series temporales con una técnica de procesamiento de señales
Modelar datos de series de tiempo y pronosticarlos son temas complejos. Existen muchas técnicas que podrían usarse para mejorar el rendimiento del modelo para un trabajo de pronóstico. Aquí discutiremos una técnica que puede mejorar la forma en que los modelos de pronóstico absorben y utilizan características de tiempo y generalizan a partir de ellas. El enfoque principal será la creación de características estacionales que alimentan el modelo de pronóstico de series de tiempo en el entrenamiento; aquí se pueden obtener ganancias fáciles si se incluye la transformada de Fourier en el proceso de creación de características.
Este artículo supone que está familiarizado con los aspectos básicos del pronóstico de series de tiempo; no discutiremos este tema en general, sino solo un refinamiento de un aspecto del mismo. Para obtener una introducción al pronóstico de series de tiempo, consulte el curso de Kaggle sobre este tema: la técnica que se analiza aquí se basa en su lección sobre la estacionalidad.
Considerar el conjunto de datos trimestrales de producción de cemento Portland australianoque muestra la producción trimestral total de cemento Portland en Australia (en millones de toneladas) desde 1956:T1 hasta 2014:T1.
df = pd.read_csv('Quarterly_Australian_Portland_Cement_production_776_10.csv', usecols=['time', 'value'])
# convert time from year float to a proper datetime format
df['time'] = df['time'].apply(lambda x: str(int(x)) + '-' + str(int(1 + 12 * (x % 1))).rjust(2, '0'))
df['time'] = pd.to_datetime(df['time'])
df = df.set_index('time').to_period()
df.rename(columns={'value': 'production'}, inplace=True)
production
time
1956Q1 0.465
1956Q2 0.532
1956Q3 0.561
1956Q4 0.570
1957Q1 0.529
... ...
2013Q1 2.049
2013Q2 2.528
2013Q3 2.637
2013Q4 2.565
2014Q1 2.229[233 rows x 1 columns]
Se trata de datos de series de tiempo con una tendencia, componentes estacionales y otros atributos. Las observaciones se realizan trimestralmente y abarcan varias décadas. Primero echemos un vistazo a los componentes estacionales, utilizando el gráfico de periodograma (todo el código está incluido en el cuaderno complementario, vinculado al final).
El periodograma muestra la densidad de potencia de los componentes espectrales (componentes estacionales). El componente más fuerte es el que tiene un período igual a 4 trimestres o 1 año. Esto confirma la impresión visual de que las variaciones más fuertes hacia arriba y hacia abajo en el gráfico ocurren aproximadamente 10 veces por década. También hay un componente con un período de 2 trimestres; es el mismo fenómeno estacional, y simplemente significa que la periodicidad de 4 trimestres no es una onda sinusoidal simple, sino que tiene una forma más compleja. Ignoraremos los componentes con períodos superiores a 4, por simplicidad.
Si intenta ajustar un modelo a estos datos, tal vez para pronosticar la producción de cemento para momentos posteriores al final del conjunto de datos, sería una buena idea capturar esta estacionalidad anual en una característica o conjunto de características, y proporcionar esas al modelo en entrenamiento. Veamos un ejemplo.
Los componentes estacionales se pueden modelar de diferentes maneras. A menudo, se representan como variables ficticias del tiempo o como pares seno-coseno. Se trata de características sintéticas con un período igual a la estacionalidad que intenta modelar, o un submúltiplo de ella.
Maniquíes de tiempo anual:
seasonal_year = DeterministicProcess(index=df.index, constant=False, seasonal=True).in_sample()
print(seasonal_year)
s(1,4) s(2,4) s(3,4) s(4,4)
time
1956Q1 1.0 0.0 0.0 0.0
1956Q2 0.0 1.0 0.0 0.0
1956Q3 0.0 0.0 1.0 0.0
1956Q4 0.0 0.0 0.0 1.0
1957Q1 1.0 0.0 0.0 0.0
... ... ... ... ...
2013Q1 1.0 0.0 0.0 0.0
2013Q2 0.0 1.0 0.0 0.0
2013Q3 0.0 0.0 1.0 0.0
2013Q4 0.0 0.0 0.0 1.0
2014Q1 1.0 0.0 0.0 0.0[233 rows x 4 columns]
Pares seno-coseno anuales:
cfr = CalendarFourier(freq='Y', order=2)
seasonal_year_trig = DeterministicProcess(index=df.index, seasonal=False, additional_terms=[cfr]).in_sample()
with pd.option_context('display.max_columns', None, 'display.expand_frame_repr', False):
print(seasonal_year_trig)
sin(1,freq=A-DEC) cos(1,freq=A-DEC) sin(2,freq=A-DEC) cos(2,freq=A-DEC)
time
1956Q1 0.000000 1.000000 0.000000 1.000000
1956Q2 0.999963 0.008583 0.017166 -0.999853
1956Q3 0.017166 -0.999853 -0.034328 0.999411
1956Q4 -0.999963 -0.008583 0.017166 -0.999853
1957Q1 0.000000 1.000000 0.000000 1.000000
... ... ... ... ...
2013Q1 0.000000 1.000000 0.000000 1.000000
2013Q2 0.999769 0.021516 0.043022 -0.999074
2013Q3 0.025818 -0.999667 -0.051620 0.998667
2013Q4 -0.999917 -0.012910 0.025818 -0.999667
2014Q1 0.000000 1.000000 0.000000 1.000000[233 rows x 4 columns]
Los muñecos de tiempo pueden capturar formas de onda muy complejas del fenómeno estacional. Por otro lado, si el período es N, entonces necesita N variables ficticias de tiempo; por lo tanto, si el período es muy largo, necesitará muchas columnas ficticias, lo que puede no ser deseable. Para los modelos no penalizados, a menudo solo se utilizan dummies N-1; se elimina uno para evitar problemas de colinealidad (lo ignoraremos aquí).
Los pares seno/coseno pueden modelar períodos de tiempo arbitrariamente largos. Cada par modelará una forma de onda sinusoidal pura con alguna fase inicial. A medida que la forma de onda de la característica estacional se vuelve más compleja, necesitará agregar más pares (aumente el orden de la salida de CalendarFourier).
(Nota al margen: utilizamos un par seno/coseno para cada período que queremos modelar. Lo que realmente queremos es sólo una columna de A*sin(2*pi*t/T + phi) dónde A es el peso asignado por el modelo a la columna, t es el índice de tiempo de la serie, y T es el período del componente. Pero el modelo no puede ajustar la fase inicial. phi mientras se ajustan los datos, los valores de los senos se calculan previamente. Afortunadamente, la fórmula anterior equivale a: A1*sin(2*pi*t/T) + A2*cos(2*pi*t/T) y el modelo sólo necesita encontrar los pesos A1 y A2.)
TLDR: Lo que estas dos técnicas tienen en común es que representan la estacionalidad a través de un conjunto de características repetitivas, con valores que alternan hasta una vez por el período de tiempo que se modela (dummy de tiempo y el par base seno/coseno), o varias veces por período (seno/coseno de orden superior). Y cada una de estas características tiene valores que varían dentro de un intervalo fijo: de 0 a 1, o de -1 a 1. Estas son todas las características que tenemos para modelar la estacionalidad aquí.
Veamos qué sucede cuando ajustamos un modelo lineal a tales características estacionales. Pero primero, debemos agregar características de tendencia a la colección de características utilizadas para entrenar el modelo.
Creemos características de tendencia y luego unámoslas con las variables ficticias de tiempo estacional_año generadas anteriormente:
trend_order = 2
trend_year = DeterministicProcess(index=df.index, constant=True, order=trend_order).in_sample()
X = trend_year.copy()
X = X.join(seasonal_year)
const trend trend_squared s(1,4) s(2,4) s(3,4) s(4,4)
time
1956Q1 1.0 1.0 1.0 1.0 0.0 0.0 0.0
1956Q2 1.0 2.0 4.0 0.0 1.0 0.0 0.0
1956Q3 1.0 3.0 9.0 0.0 0.0 1.0 0.0
1956Q4 1.0 4.0 16.0 0.0 0.0 0.0 1.0
1957Q1 1.0 5.0 25.0 1.0 0.0 0.0 0.0
... ... ... ... ... ... ... ...
2013Q1 1.0 229.0 52441.0 1.0 0.0 0.0 0.0
2013Q2 1.0 230.0 52900.0 0.0 1.0 0.0 0.0
2013Q3 1.0 231.0 53361.0 0.0 0.0 1.0 0.0
2013Q4 1.0 232.0 53824.0 0.0 0.0 0.0 1.0
2014Q1 1.0 233.0 54289.0 1.0 0.0 0.0 0.0[233 rows x 7 columns]
Este es el marco de datos X (características) que se utilizará para entrenar/validar el modelo. Estamos modelando los datos con características de tendencia cuadrática, además de las cuatro variables ficticias necesarias para la estacionalidad anual. El marco de datos y (objetivo) serán solo las cifras de producción de cemento.
Hagamos un conjunto de validación a partir de los datos, que contenga las observaciones del año 2010. Ajustaremos un modelo lineal a las características X que se muestran arriba y al objetivo y representado por la producción de cemento (la parte de prueba), y luego evaluaremos el desempeño del modelo en la validación. También haremos todo lo anterior con un modelo de tendencia exclusiva que sólo se ajustará a las características de la tendencia pero ignorará la estacionalidad.
def do_forecast(X, index_train, index_test, trend_order):
X_train = X.loc[index_train]
X_test = X.loc[index_test]y_train = df['production'].loc[index_train]
y_test = df['production'].loc[index_test]
model = LinearRegression(fit_intercept=False)
_ = model.fit(X_train, y_train)
y_fore = pd.Series(model.predict(X_test), index=index_test)
y_past = pd.Series(model.predict(X_train), index=index_train)
trend_columns = X_train.columns.to_list()[0 : trend_order + 1]
model_trend = LinearRegression(fit_intercept=False)
_ = model_trend.fit(X_train[trend_columns], y_train)
y_trend_fore = pd.Series(model_trend.predict(X_test[trend_columns]), index=index_test)
y_trend_past = pd.Series(model_trend.predict(X_train[trend_columns]), index=index_train)
RMSLE = mean_squared_log_error(y_test, y_fore, squared=False)
print(f'RMSLE: {RMSLE}')
ax = df.plot(**plot_params, title='AUS Cement Production - Forecast')
ax = y_past.plot(color='C0', label='Backcast')
ax = y_fore.plot(color='C3', label='Forecast')
ax = y_trend_past.plot(ax=ax, color='C0', linewidth=3, alpha=0.333, label='Trend Past')
ax = y_trend_fore.plot(ax=ax, color='C3', linewidth=3, alpha=0.333, label='Trend Future')
_ = ax.legend()
do_forecast(X, index_train, index_test, trend_order)
RMSLE: 0.03846449744356434
El azul es tren, el rojo es validación. El modelo completo es la línea fina y nítida. El modelo de tendencia es la línea ancha y difusa.
Esto no está mal, pero hay un problema evidente: el modelo ha aprendido una estacionalidad anual de amplitud constante. Aunque la variación anual en realidad aumenta con el tiempo, el modelo sólo puede atenerse a variaciones de amplitud fija. Obviamente, esto se debe a que le dimos al modelo solo características estacionales de amplitud fija, y no hay nada más en las características (retrasos objetivo, etc.) para ayudarlo a superar este problema.
Profundicemos en la señal (los datos) para ver si hay algo que pueda ayudar con el problema de la amplitud fija.
El periodograma resaltará todos los componentes espectrales de la señal (todos los componentes estacionales de los datos) y proporcionará una descripción general de su intensidad general, pero es un agregado, una suma de todo el intervalo de tiempo. No dice nada sobre cómo la amplitud de cada componente estacional puede variar en el tiempo en todo el conjunto de datos.
Para capturar esa variación, hay que utilizar el espectrograma de Fourier. Es como el periodograma, pero se realiza repetidamente durante muchas ventanas de tiempo en todo el conjunto de datos. Tanto el periodograma como el espectrograma están disponibles como métodos en la biblioteca scipy.
Tracemos el espectrograma de los componentes estacionales con períodos de 2 y 4 trimestres, mencionados anteriormente. Como siempre, el código completo se encuentra en el cuaderno complementario vinculado al final.
spectrum = compute_spectrum(df['production'], 4, 0.1)
plot_spectrogram(spectrum, figsize_x=10)
Lo que este diagrama muestra es la amplitud, la “fuerza” de los componentes de 2 cuartos y 4 cuartos a lo largo del tiempo. Inicialmente son bastante débiles, pero se vuelven muy fuertes alrededor de 2010, lo que coincide con las variaciones que se pueden ver en el gráfico del conjunto de datos al principio de este artículo.
¿Qué pasa si, en lugar de alimentar el modelo con características estacionales de amplitud fija, ajustamos la amplitud de estas características en el tiempo, coincidiendo con lo que nos dice el espectrograma? ¿El modelo final funcionaría mejor en la validación?
Escojamos el componente estacional del cuarto trimestre. Podríamos ajustar un modelo lineal (llamado modelo envolvente) en la tendencia de orden=2 de este componente, en los datos del tren (model.fit()), y extender esa tendencia a la validación/prueba (model.predict()):
envelope_features = DeterministicProcess(index=X.index, constant=True, order=2).in_sample()spec4_train = compute_spectrum(df['production'].loc[index_train], max_period=4)
spec4_train
spec4_model = LinearRegression()
spec4_model.fit(envelope_features.loc[spec4_train.index], spec4_train['4.0'])
spec4_regress = pd.Series(spec4_model.predict(envelope_features), index=X.index)
ax = spec4_train['4.0'].plot(label='component envelope', color='gray')
spec4_regress.loc[spec4_train.index].plot(ax=ax, color='C0', label='envelope regression: past')
spec4_regress.loc[index_test].plot(ax=ax, color='C3', label='envelope regression: future')
_ = ax.legend()
La línea azul es la fuerza de la variación del componente de 4 trimestres en los datos del tren, ajustado como una tendencia cuadrática (orden=2). La línea roja es lo mismo, extendida (predicha) sobre los datos de validación.
Hemos modelado la variación en el tiempo del componente estacional de 4 trimestres. Podemos tomar el resultado del modelo envolvente y multiplicarlo por las variables ficticias de tiempo correspondientes al componente estacional de 4 trimestres:
spec4_regress = spec4_regress / spec4_regress.mean()season_columns = ['s(1,4)', 's(2,4)', 's(3,4)', 's(4,4)']
for c in season_columns:
X[c] = X[c] * spec4_regress
print(X)
const trend trend_squared s(1,4) s(2,4) s(3,4) s(4,4)
time
1956Q1 1.0 1.0 1.0 0.179989 0.000000 0.000000 0.000000
1956Q2 1.0 2.0 4.0 0.000000 0.181109 0.000000 0.000000
1956Q3 1.0 3.0 9.0 0.000000 0.000000 0.182306 0.000000
1956Q4 1.0 4.0 16.0 0.000000 0.000000 0.000000 0.183581
1957Q1 1.0 5.0 25.0 0.184932 0.000000 0.000000 0.000000
... ... ... ... ... ... ... ...
2013Q1 1.0 229.0 52441.0 2.434701 0.000000 0.000000 0.000000
2013Q2 1.0 230.0 52900.0 0.000000 2.453436 0.000000 0.000000
2013Q3 1.0 231.0 53361.0 0.000000 0.000000 2.472249 0.000000
2013Q4 1.0 232.0 53824.0 0.000000 0.000000 0.000000 2.491139
2014Q1 1.0 233.0 54289.0 2.510106 0.000000 0.000000 0.000000[233 rows x 7 columns]
Los 4 muñecos de tiempo ya no son 0 ni 1. Se han multiplicado por la envolvente del componente y ahora su amplitud varía en el tiempo, al igual que la envolvente.
Entrenemos el modelo principal nuevamente, ahora usando los muñecos de tiempo modificados.
do_forecast(X, index_train, index_test, trend_order)
RMSLE: 0.02546321729737165
No buscamos un ajuste perfecto aquí, ya que este es solo un modelo de juguete, pero es obvio que el modelo funciona mejor, sigue más de cerca las variaciones de 4 cuartos en el objetivo. La métrica de rendimiento también ha mejorado, en un 51%, lo que no está nada mal.
Mejorar el rendimiento del modelo es un tema muy amplio. El comportamiento de cualquier modelo no depende de una única característica, ni de una única técnica. Sin embargo, si está buscando extraer todo lo que pueda de un modelo determinado, probablemente debería proporcionarle características significativas. Las dummies de tiempo, o pares de Fourier seno/coseno, son más significativas cuando reflejan la variación en el tiempo de la estacionalidad que están modelando.
Ajustar la envolvente del componente estacional mediante la transformada de Fourier es particularmente efectivo para modelos lineales. Los árboles potenciados no se benefician tanto, pero aun así puedes ver mejoras casi sin importar el modelo que uses. Si está utilizando un modelo híbrido, donde un modelo lineal maneja características deterministas (calendario, etc.), mientras que un modelo de árboles potenciados maneja características más complejas, es importante que el modelo lineal “lo haga bien”, por lo que deja menos trabajo por hacer. Hazlo para el otro modelo.
También es importante que los modelos envolventes que utilice para ajustar las características estacionales solo se entrenen con los datos del tren y no vean ningún dato de prueba durante el entrenamiento, como cualquier otro modelo. Una vez que se ajusta la envolvente, las réplicas de tiempo contienen información del objetivo; ya no son características puramente deterministas, que pueden calcularse con anticipación en cualquier horizonte de pronóstico arbitrario. Entonces, en este punto, la información podría filtrarse de la validación/prueba a los datos de entrenamiento si no se tiene cuidado.
El conjunto de datos utilizado en este artículo está disponible aquí bajo la licencia de dominio público (CC0):
El código utilizado en este artículo se puede encontrar aquí:
Un cuaderno enviado al concurso Ventas en tienda: pronóstico de series temporales en Kaggle, utilizando las ideas descritas en este artículo:
Un repositorio de GitHub que contiene la versión de desarrollo del cuaderno enviado a Kaggle está aquí:
Todas las imágenes y el código utilizados en este artículo son creados por el autor.