Time Series Python Web

By Prof. Angel Colmenares, 05-2025

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 únicamente en función de su comportamiento pasado (autorregresivo) o empleando otras variables externas.

Ingeniería de características para series temporales

El objetivo de la ingeniería de características es proporcionar relaciones sólidas e idealmente sencillas entre las nuevas características de entrada y las características de salida para que el algoritmo las modele. En efecto, estamos desplazando la complejidad. La complejidad reside en las relaciones entre los datos de entrada y de salida.

En el caso de las series temporales, no existe el concepto de variables de entrada y salida. Podemos apoyarnos en la capacidad de modelos sofisticados para descifrar la complejidad del problema. Podemos facilitar el trabajo a estos modelos (e incluso utilizar modelos más sencillos) si podemos exponer mejor la relación inherente entre entradas y los resultados de los datos. La dificultad estriba en que no conocemos la relación funcional inherente subyacente entre los datos de entrada y de salida que intentamos exponer. Si lo supiéramos, probablemente no necesitaríamos aprendizaje automático. En su lugar, la única información que tenemos es el rendimiento de los modelos desarrollados sobre los conjuntos de datos de aprendizaje supervisado o las visiones del problema que creamos. En efecto, la mejor estrategia por defecto es utilizar todos los conocimientos disponibles para crear muchos conjuntos de datos y utilizar el rendimiento del modelo.

Datos

En esta seccion, utilizaremos como ejemplo el conjunto de datos Temperaturas mínimas diarias. Este conjunto de datos describe las temperaturas mínimas diarias a lo largo de 10 años (1981-1990) en la ciudad de Melbourne Australia.

https://github.com/upul/WhiteBoard/blob/master/data/daily-minimum-temperatures-in-me.csv

Funciones de fecha y hora

Empecemos por algunas de las funciones más sencillas que podemos utilizar. Se trata de características de la fecha/hora de cada observación. De hecho, estas pueden ser bastante complejas. Dos características con las que podemos empezar son el mes y el día enteros de cada observación.

from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-minimum-temperatures.csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
dataframe = DataFrame()
dataframe['month'] = [series.index[i].month for i in range(len(series))]
dataframe['day'] = [series.index[i].day for i in range(len(series))]
dataframe['temperature'] = [series[i] for i in range(len(series))]
print(dataframe.head(5))

   month  day  temperature
0      1    1         20.7
1      1    2         17.9
2      1    3         18.8
3      1    4         14.6
4      1    5         15.8

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\3700144450.py:9: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

temps = DataFrame(series.values)
dataframe = concat([temps.shift(3), temps.shift(2), temps.shift(1), temps], axis=1)
dataframe.columns = ['t-2', 't-1', 't', 't+1']
print(dataframe.head(5))

    t-2   t-1     t   t+1
0   NaN   NaN   NaN  20.7
1   NaN   NaN  20.7  17.9
2   NaN  20.7  17.9  18.8
3  20.7  17.9  18.8  14.6
4  17.9  18.8  14.6  15.8

Estadísticas de ventana móvil

Podemos calcular estadísticas resumidas de los valores de la ventana móvil e incluirlas como características en nuestro conjunto de datos. Quizá la más útil sea la media de los valores anteriores, también llamada media móvil.

Podemos calcular la media de los valores actuales y anteriores y utilizarla para predecir el valor siguiente. Para los datos de temperatura, tendríamos que esperar 3 pasos de tiempo antes de tener 2 valores de los que sacar la media antes de poder utilizar ese valor para predecir un tercer valor.

# create a rolling mean feature
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-minimum-temperatures.csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
temps = DataFrame(series.values)
shifted = temps.shift(1)
window = shifted.rolling(window=2)
means = window.mean()
dataframe = concat([means, temps], axis=1)
dataframe.columns = ['mean(t-1,t)', 't+1']
print(dataframe.head(5))
temps = DataFrame(series.values)

   mean(t-1,t)   t+1
0          NaN  20.7
1          NaN  17.9
2        19.30  18.8
3        18.35  14.6
4        16.70  15.8

Expansión de estadísticas de ventanas

Otro tipo de ventana que puede resultar útil incluye todos los datos anteriores de la serie. Esto se denomina ventana expansiva y puede ayudar a seguir los límites de los datos observables. Al igual que la función rolling() en DataFrame, Pandas proporciona una función expanding() que recoge conjuntos de todos los valores anteriores para cada paso de tiempo.

# create expanding window features
window = temps.expanding()
dataframe = concat([window.min(), window.mean(), window.max(), temps.shift(-1)], axis=1)
dataframe.columns = ['min', 'mean', 'max', 't+1']
print(dataframe.head(5))

    min       mean   max   t+1
0  20.7  20.700000  20.7  17.9
1  17.9  19.300000  20.7  18.8
2  17.9  19.133333  20.7  14.6
3  14.6  18.000000  20.7  15.8
4  14.6  17.560000  20.7  15.8

Visualización de series temporales

La visualización desempeña un papel importante en el análisis y la previsión de series temporales. Los gráficos de los datos pueden proporcionar valiosos diagnósticos para identificar estructuras temporales como tendencias, ciclos y estacionalidad, que pueden influir en la elección del modelo, y estacionalidad que pueden influir en la elección del modelo. El problema es que muchos principiantes de la previsión de series temporales se detienen en los gráficos de líneas. En este tutorial, echaremos un vistazo a 6 tipos de visualizaciones que puede utilizar con sus propios datos de series temporales. Son los siguientes

• Gráficos de líneas.

• Histogramas y gráficos de densidad.

• Gráficos de cajas y bigotes.

• Mapas de calor.

• Lag Plots o gráficos de dispersión.

• Gráficos de autocorrelación.

La atención se centra en las series temporales univariantes, pero las técnicas son igualmente aplicables a las series temporales multivariantes, cuando se dispone de más de una observación en cada paso temporal. multivariadas, cuando se tiene más de una observación en cada paso temporal. A continuación, echemos un vistazo al conjunto de datos que utilizaremos para demostrar la visualización de series temporales en este material.

Gráfico de líneas

La primera visualización de series temporales, y quizás la más popular, es el gráfico de líneas. En este gráfico, el tiempo se muestra en el eje x con los valores de observación a lo largo del eje y. A continuación se muestra un ejemplo de visualización de la Serie Pandas del conjunto de datos de Temperaturas Mínimas Diarias directamente como un gráfico de líneas.

# create a line plot
from pandas import read_csv
from matplotlib import pyplot
from skforecast.plot import set_dark_theme
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
set_dark_theme()
series.plot()
pyplot.show()

# create a dot plot
series.plot(style='k.')
pyplot.show()

# create stacked line plots
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
groups = series.groupby(Grouper(freq='Y'))
years = DataFrame()
for name, group in groups:
   years[name.year] = group.values
years.plot(subplots=True, legend=False)
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\2150578903.py:8: FutureWarning:

'Y' is deprecated and will be removed in a future version, please use 'YE' instead.

Histograma y gráficos de densidad

Otra visualización importante es la de la distribución de las propias observaciones. Es decir, un gráfico de los valores sin la ordenación temporal. Algunos métodos de previsión lineal de series temporales suponen que la distribución de las observaciones es correcta (es decir, una curva de campana o una distribución normal). Esto puede comprobarse explícitamente con herramientas como las pruebas de hipótesis estadísticas. Pero los gráficos pueden proporcionar una primera comprobación útil de la distribución de las observaciones, tanto en las observaciones brutas como después de realizar cualquier tipo de transformación de los datos.

# create a histogram plot
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
series.hist()
pyplot.show()
pyplot.close()

# create a density plot
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
series.plot(kind='kde')
pyplot.show()
pyplot.close()

Gráficos de caja y bigotes por intervalo

Los histogramas y los gráficos de densidad proporcionan una visión de la distribución de todas las observaciones, pero podemos estar interesados en la distribución de los valores por intervalo de tiempo. Otro tipo de gráfico útil para resumir la distribución de las observaciones es el gráfico de cajas y bigotes. Este gráfico dibuja una caja alrededor de los percentiles 25 y 75 de los datos que captura el 50% medio de las observaciones. Se traza una línea en el percentil 50 (la mediana) y los bigotes se dibujan por encima y por debajo de la caja para resumir la extensión general de las observaciones. Se dibujan puntos para los valores atípicos fuera de los bigotes o de la extensión de los datos.

# create a boxplot of yearly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
groups = series.groupby(Grouper(freq='Y'))
years = DataFrame()
for name, group in groups:
  years[name.year] = group.values
years.boxplot()
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\1192400507.py:8: FutureWarning:

'Y' is deprecated and will be removed in a future version, please use 'YE' instead.

# create a boxplot of monthly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
from pandas import concat
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
one_year = series['1990']
groups = one_year.groupby(Grouper(freq='M'))
months = concat([DataFrame(x[1].values) for x in groups], axis=1)
months = DataFrame(months)
months.columns = range(1,13)
months.boxplot()
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\842634312.py:10: FutureWarning:

'M' is deprecated and will be removed in a future version, please use 'ME' instead.

Mapas de calor

Una matriz de números puede trazarse como una superficie, en la que a los valores de cada celda de la matriz se les asigna un color único. Esto se denomina mapa de calor, ya que los valores más grandes se pueden dibujar con colores más cálidos (amarillos y rojos) y los valores más pequeños se pueden dibujar con colores más fríos (azules y verdes). Al igual que los gráficos de caja y bigotes, podemos comparar observaciones entre intervalos utilizando un mapa de calor.

# create a heat map of yearly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
groups = series.groupby(Grouper(freq='Y'))
years = DataFrame()
for name, group in groups:
  years[name.year] = group.values
years = years.T
pyplot.matshow(years, interpolation=None, aspect='auto')
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\3935000226.py:8: FutureWarning:

'Y' is deprecated and will be removed in a future version, please use 'YE' instead.

# create a heat map of monthly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
from pandas import concat
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
one_year = series['1990']
groups = one_year.groupby(Grouper(freq='M'))
months = concat([DataFrame(x[1].values) for x in groups], axis=1)
months = DataFrame(months)
months.columns = range(1,13)
pyplot.matshow(months, interpolation=None, aspect='auto')
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\192252240.py:10: FutureWarning:

'M' is deprecated and will be removed in a future version, please use 'ME' instead.

Gráficos de dispersión de retrasos

La modelización de series temporales presupone una relación entre una observación y la observación anterior. Las observaciones anteriores en una serie temporal se denominan retardos, con la observación en el paso de tiempo anterior denominada lag1, la observación en dos pasos de tiempo hace lag=2, y así sucesivamente. Un tipo de gráfico útil para explorar la relación entre cada observación y un retardo de esa observación se llama gráfico de dispersión. Pandas tiene una función incorporada para exactamente esto llamado el gráfico de retardo. Se traza la observación en el tiempo t en el eje x y la observación en el siguiente paso de tiempo (t + 1) en el eje y.

# create a scatter plot
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import lag_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
lag_plot(series)
pyplot.show()

# create multiple scatter plots
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
values = DataFrame(series.values)
lags = 7
columns = [values]
for i in range(1,(lags + 1)):
  columns.append(values.shift(i))
dataframe = concat(columns, axis=1)
columns = ['t']
for i in range(1,(lags + 1)):
  columns.append('t-' + str(i))
dataframe.columns = columns
pyplot.figure(1)
for i in range(1,(lags + 1)):
  ax = pyplot.subplot(240 + i)
  ax.set_title('t vs t-' + str(i))
  pyplot.scatter(x=dataframe['t'].values, y=dataframe['t-'+str(i)].values)
pyplot.show()

Gráficos de autocorrelación

Podemos cuantificar la fuerza y el tipo de relación entre las observaciones y sus retardos. En estadística, esto se denomina correlación y, cuando se calcula en función de los valores de retardo de las series temporales, se denomina autocorrelación. Un valor de correlación calculado entre dos grupos de números, como las observaciones y sus valores de retardo = 1, da como resultado un número comprendido entre -1 y 1. El signo de este número indica una correlación negativa o positiva, respectivamente. Un valor cercano a cero indica una correlación débil, mientras que un valor cercano a -1 o 1 indica una correlación fuerte.

# create an autocorrelation plot
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
autocorrelation_plot(series)
pyplot.show()
pyplot.close()

Remuestreo e interpolación

Es posible que tenga observaciones con una frecuencia incorrecta. Tal vez sean demasiado granulares o no lo suficientemente granulares. La librería Pandas en Python proporciona la capacidad de cambiar la frecuencia de sus datos de series temporales. En este tutorial, descubrirá cómo utilizar Pandas en Python tanto para aumentar como para disminuir la frecuencia de muestreo de los datos de series temporales.

Conjunto de datos

En esta lección, utilizaremos como ejemplo el conjunto de datos Ventas de champú. Este conjunto de datos describe el número mensual de ventas de champú durante un periodo de 3 años.

Remuestreo de datos

Las observaciones de las ventas de champú son mensuales. Imaginemos que quisiéramos obtener información sobre las ventas diarias. Tendríamos que remuestrear la frecuencia de mensual a diaria y utilizar un esquema de interpolación para obtener la nueva frecuencia diaria. La librería Pandas proporciona una función llamada resample() en los objetos Series y DataFrame. Esta función se puede utilizar para agrupar registros cuando se reduce el muestreo y para hacer espacio para nuevas observaciones cuando se aumenta el muestreo.

# upsample to daily intervals
from pandas import read_csv
from datetime import datetime

def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
upsampled = series.resample('D').mean()
print(upsampled.head(32))

Month
1901-01-01    266.0
1901-01-02      NaN
1901-01-03      NaN
1901-01-04      NaN
1901-01-05      NaN
1901-01-06      NaN
1901-01-07      NaN
1901-01-08      NaN
1901-01-09      NaN
1901-01-10      NaN
1901-01-11      NaN
1901-01-12      NaN
1901-01-13      NaN
1901-01-14      NaN
1901-01-15      NaN
1901-01-16      NaN
1901-01-17      NaN
1901-01-18      NaN
1901-01-19      NaN
1901-01-20      NaN
1901-01-21      NaN
1901-01-22      NaN
1901-01-23      NaN
1901-01-24      NaN
1901-01-25      NaN
1901-01-26      NaN
1901-01-27      NaN
1901-01-28      NaN
1901-01-29      NaN
1901-01-30      NaN
1901-01-31      NaN
1901-02-01    145.9
Freq: D, Name: Sales, dtype: float64

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\1241687185.py:7: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

Un buen punto de partida es utilizar una interpolación lineal. Esto dibuja una línea recta entre los datos disponibles, en este caso el primero del mes, y rellena los valores con la frecuencia elegida a partir de esta línea.

# upsample to daily intervals with linear interpolation
from pandas import read_csv
from datetime import datetime
from matplotlib import pyplot
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
upsampled = series.resample('D').mean()
interpolated = upsampled.interpolate(method='linear')
print(interpolated.head(32))
interpolated.plot()
pyplot.show()
pyplot.close()

Month
1901-01-01    266.000000
1901-01-02    262.125806
1901-01-03    258.251613
1901-01-04    254.377419
1901-01-05    250.503226
1901-01-06    246.629032
1901-01-07    242.754839
1901-01-08    238.880645
1901-01-09    235.006452
1901-01-10    231.132258
1901-01-11    227.258065
1901-01-12    223.383871
1901-01-13    219.509677
1901-01-14    215.635484
1901-01-15    211.761290
1901-01-16    207.887097
1901-01-17    204.012903
1901-01-18    200.138710
1901-01-19    196.264516
1901-01-20    192.390323
1901-01-21    188.516129
1901-01-22    184.641935
1901-01-23    180.767742
1901-01-24    176.893548
1901-01-25    173.019355
1901-01-26    169.145161
1901-01-27    165.270968
1901-01-28    161.396774
1901-01-29    157.522581
1901-01-30    153.648387
1901-01-31    149.774194
1901-02-01    145.900000
Freq: D, Name: Sales, dtype: float64

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\2755792297.py:7: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

Otro método común de interpolación es utilizar un polinomio o un spline para conectar los valores. Esto crea más curvas y puede parecer más natural en muchos conjuntos de datos. El uso de una interpolación spline requiere que especifique el orden (número de términos en el polinomio); en este caso, un orden de 2 está bien.

# upsample to daily intervals with spline interpolation
from pandas import read_csv
from datetime import datetime
from matplotlib import pyplot
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
upsampled = series.resample('D').mean()
interpolated = upsampled.interpolate(method='spline', order=2)
print(interpolated.head(32))
interpolated.plot()
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\3377822813.py:7: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

Month
1901-01-01    266.000000
1901-01-02    258.630160
1901-01-03    251.560886
1901-01-04    244.720748
1901-01-05    238.109746
1901-01-06    231.727880
1901-01-07    225.575149
1901-01-08    219.651553
1901-01-09    213.957094
1901-01-10    208.491770
1901-01-11    203.255582
1901-01-12    198.248529
1901-01-13    193.470612
1901-01-14    188.921831
1901-01-15    184.602185
1901-01-16    180.511676
1901-01-17    176.650301
1901-01-18    173.018063
1901-01-19    169.614960
1901-01-20    166.440993
1901-01-21    163.496161
1901-01-22    160.780465
1901-01-23    158.293905
1901-01-24    156.036481
1901-01-25    154.008192
1901-01-26    152.209039
1901-01-27    150.639021
1901-01-28    149.298139
1901-01-29    148.186393
1901-01-30    147.303783
1901-01-31    146.650308
1901-02-01    145.900000
Freq: D, Name: Sales, dtype: float64

Muestreo descendente de datos

Los datos de ventas son mensuales, pero quizá preferiríamos que los datos fueran trimestrales. El año puede dividirse en 4 trimestres comerciales, de 3 meses cada uno. En lugar de crear nuevas filas entre las observaciones existentes, la función resample() de Pandas agrupará todas las observaciones por la nueva frecuencia.

# downsample to quarterly intervals
from pandas import read_csv
from datetime import datetime
from matplotlib import pyplot
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
resample = series.resample('Q')
quarterly_mean_sales = resample.mean()
print(quarterly_mean_sales.head())
quarterly_mean_sales.plot()
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\3452663358.py:7: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\3452663358.py:8: FutureWarning:

'Q' is deprecated and will be removed in a future version, please use 'QE' instead.

Month
1901-03-31    198.333333
1901-06-30    156.033333
1901-09-30    216.366667
1901-12-31    215.100000
1902-03-31    184.633333
Freq: QE-DEC, Name: Sales, dtype: float64

Quizás queramos ir más allá y convertir los datos mensuales en datos anuales, y quizás utilizarlos más tarde para modelizar el año siguiente. Podemos reducir la muestra de los datos utilizando el alias A para la frecuencia de fin de año y esta vez utilizar la suma para calcular las ventas totales de cada año.

series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
resample = series.resample('Y')
yearly_mean_sales = resample.sum()
print(yearly_mean_sales.head())
yearly_mean_sales.plot()
pyplot.show()
pyplot.close()

Month
1901-12-31    2357.5
1902-12-31    3153.5
1903-12-31    5742.6
Freq: YE-DEC, Name: Sales, dtype: float64

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\867565110.py:1: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\867565110.py:2: FutureWarning:

'Y' is deprecated and will be removed in a future version, please use 'YE' instead.

Transformaciones de potencia

El objetivo de las transformaciones de datos es eliminar el ruido y mejorar la señal en la predicción de series temporales. Puede resultar muy difícil seleccionar una buena, o incluso la mejor, transformación para un determinado problema de predicción. Hay muchas transformaciones entre las que elegir y cada una tiene una intuición matemática diferente. En este tutorial, descubrirá cómo explorar diferentes transformadas basadas en potencia para la predicción de series temporales con Python.

Conjunto de datos de pasajeros de líneas aéreas

En esta seccion, utilizaremos como ejemplo el conjunto de datos Pasajeros de líneas aéreas. Este conjunto de datos describe el número total de pasajeros de líneas aéreas a lo largo del tiempo.

# load and plot a time series
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(series)
# histogram
pyplot.subplot(212)
pyplot.hist(series)
pyplot.show()
pyplot.close()

Transformación de raíz cuadrada

Una serie temporal con una tendencia de crecimiento cuadrática puede convertirse en lineal sacando la raíz cuadrada. Vamos a demostrarlo con un ejemplo rápido. Consideremos una serie de números del 1 al 99 al cuadrado. El gráfico de líneas de esta serie mostrará una tendencia de crecimiento cuadrática y un histograma de los valores mostrará una distribución exponencial con una cola larga. El siguiente fragmento de código crea y representa gráficamente esta serie.

# contrive a quadratic time series
from matplotlib import pyplot
series = [i**2 for i in range(1,100)]
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(series)
# histogram
pyplot.subplot(212)
pyplot.hist(series)
pyplot.show()
pyplot.close()

Si observa una estructura como ésta en su propia serie temporal, es posible que tenga una tendencia de crecimiento cuadrática. El siguiente ejemplo realiza una transformación sqrt() en la serie temporal y representa gráficamente el resultado.

# square root transform a contrived quadratic time series
from matplotlib import pyplot
from numpy import sqrt
series = [i**2 for i in range(1,100)]
# sqrt transform
transform = series = sqrt(series)
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(transform)
# histogram
pyplot.subplot(212)
pyplot.hist(transform)
pyplot.show()
pyplot.close()

Es posible que el conjunto de datos de pasajeros de líneas aéreas muestre un crecimiento cuadrático. Si este es el caso, podríamos esperar que una transformación de raíz cuadrada redujera la tendencia de crecimiento a lineal y cambiara la distribución de las observaciones para que fuera quizás casi gaussiana. El siguiente ejemplo realiza una raíz cuadrada del conjunto de datos y representa gráficamente los resultados.

# square root transform a time series
from pandas import read_csv
from pandas import DataFrame
from numpy import sqrt
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'] = sqrt(dataframe['passengers'])
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()

Transformada logarítmica

Una clase de tendencias más extremas son las exponenciales, a menudo representadas gráficamente como un palo de hockey. Las series temporales con una distribución exponencial pueden hacerse lineales tomando el logaritmo de los valores. Esto se denomina transformación logarítmica. Como en el caso anterior del cuadrado y la raíz cuadrada, podemos demostrarlo con un ejemplo rápido. El código siguiente crea una distribución exponencial elevando los números del 1 al 99 al valor e, que es la base de los logaritmos naturales o número de Euler.

# create and plot an exponential time series
from matplotlib import pyplot
from math import exp
series = [exp(i) for i in range(1,100)]
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(series)
# histogram
pyplot.subplot(212)
pyplot.hist(series)
pyplot.show()
pyplot.close()

De nuevo, podemos volver a transformar esta serie en lineal tomando el logaritmo natural de los valores. Esto haría que la serie fuera lineal y la distribución uniforme. El ejemplo siguiente lo demuestra para completar.

# log transform a contrived exponential time series
from matplotlib import pyplot
from math import exp
from numpy import log
series = [exp(i) for i in range(1,100)]
transform = log(series)
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(transform)
# histogram
pyplot.subplot(212)
pyplot.hist(transform)
pyplot.show()

Nuestro conjunto de datos Pasajeros de líneas aéreas tiene una distribución de este tipo, pero quizá no tan extrema. El siguiente ejemplo muestra una transformación logarítmica del conjunto de datos Pasajeros de líneas aéreas.

# log transform a time series
from pandas import read_csv
from pandas import DataFrame
from numpy import log
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'] = log(dataframe['passengers'])
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()

Transformada de Box-Cox

La transformada de raíz cuadrada y la transformada logarítmica pertenecen a una clase de transformadas denominadas transformadas de potencia. La transformación de Box-Cox2 es un método de transformación de datos configurable que admite la raíz cuadrada y la transformación logarítmica, así como un conjunto de transformaciones relacionadas.

# manually box-cox transform a time series
from pandas import read_csv
from pandas import DataFrame
from scipy.stats import boxcox
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'] = boxcox(dataframe['passengers'], lmbda=0.0)
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()

Podemos establecer el parámetro lambda en Ninguno (el valor predeterminado) y dejar que la función y un valor ajustado estadísticamente. El siguiente ejemplo demuestra este uso, devolviendo tanto el conjunto de datos transformado como el valor lambda elegido.

# automatically box-cox transform a time series
from pandas import read_csv
from pandas import DataFrame
from scipy.stats import boxcox
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'], lam = boxcox(dataframe['passengers'])
print('Lambda: %f' % lam)
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()

Lambda: 0.148023

Suavizado de medias móviles

El suavizado de medias móviles es una técnica ingenua y eficaz de previsión de series temporales. Se puede utilizar para la preparación de datos, ingeniería de características, e incluso directamente para hacer predicciones. En este tutorial, descubrirá cómo utilizar el suavizado de medias móviles para la previsión de series temporales con Python.

La media móvil puede utilizarse como técnica de preparación de datos para crear una versión suavizada del conjunto de datos original. La suavización es útil como técnica de preparación de datos, ya que puede reducir la variación aleatoria de las observaciones y exponer mejor la estructura de los procesos causales subyacentes. subyacentes.

# moving average smoothing as data preparation
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# tail-rolling average transform
rolling = series.rolling(window=3)
rolling_mean = rolling.mean()
print(rolling_mean.head(10))
# plot original and transformed dataset
series.plot()
rolling_mean.plot(color='red')
pyplot.show()
pyplot.close()
# zoomed plot original and transformed dataset
series[:100].plot()
rolling_mean[:100].plot(color='red')
pyplot.show()
pyplot.close()

Date
1959-01-01          NaN
1959-01-02          NaN
1959-01-03    32.333333
1959-01-04    31.000000
1959-01-05    35.000000
1959-01-06    34.666667
1959-01-07    39.333333
1959-01-08    39.000000
1959-01-09    42.000000
1959-01-10    36.000000
Name: Births, dtype: float64

Uso de la media móvil como ingeniería

La media móvil puede utilizarse como fuente de nueva información al modelizar una previsión de series temporales como un problema de aprendizaje supervisado. En este caso, la media móvil se calcula y se añade como una nueva característica de entrada utilizada para predecir el siguiente paso temporal. En primer lugar, hay que desplazar una copia de la serie un paso adelante. Esto representará la entrada a nuestro problema de predicción, o una versión lag=1 de la serie. Se trata de una visión estándar de aprendizaje supervisado del problema de las series temporales.

# moving average smoothing as feature engineering
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
df = DataFrame(series.values)
width = 3
lag1 = df.shift(1)
lag3 = df.shift(width - 1)
window = lag3.rolling(window=width)
means = window.mean()
dataframe = concat([means, lag1, df], axis=1)
dataframe.columns = ['mean', 't', 't+1']
print(dataframe.head(10))

        mean     t  t+1
0        NaN   NaN   35
1        NaN  35.0   32
2        NaN  32.0   30
3        NaN  30.0   31
4  32.333333  31.0   44
5  31.000000  44.0   29
6  35.000000  29.0   45
7  34.666667  45.0   43
8  39.333333  43.0   38
9  39.000000  38.0   27

La media móvil como predicción

El valor de la media móvil también puede utilizarse directamente para hacer predicciones. Es un modelo ingenuo y supone que los componentes de tendencia y estacionalidad de la serie temporal ya se han eliminado o ajustado. El modelo de media móvil para predicciones puede utilizarse fácilmente de forma progresiva. A medida que se disponga de nuevas observaciones (por ejemplo, diariamente), se puede actualizar el modelo y hacer una predicción para el día siguiente. Podemos implementar esto manualmente en Python.

# moving average smoothing as a forecast model
#!pip install scikit-learn
from math import sqrt
from pandas import read_csv
from numpy import mean
from sklearn.metrics import mean_squared_error
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# prepare situation
X = series.values
X
window = 3
history = [X[i] for i in range(window)]
history
test = [X[i] for i in range(window, len(X))]
test
predictions = list()
# walk forward over time steps in test
for t in range(len(test)):
    length = len(history)
    yhat = mean([history[i] for i in range(length-window,length)])
    obs = test[t]
    predictions.append(yhat)
    history.append(obs)
    #print('predicted=%f, expected=%f' % (yhat, obs))
rmse = sqrt(mean_squared_error(test, predictions))
print('Test RMSE: %.3f' % rmse)
# plot
pyplot.plot(test)
pyplot.plot(predictions, color='red')
pyplot.show()
# zoom plot
pyplot.plot(test[:100])
pyplot.plot(predictions[:100], color='red')
pyplot.show()

Test RMSE: 7.834

Introducción al ruido blanco

El ruido blanco es un concepto importante en la previsión de series temporales. Si una serie temporal es ruido blanco, se trata de una secuencia de números aleatorios y no puede predecirse. Si la serie de errores de previsión no es ruido blanco, sugiere que se podrían introducir mejoras en el modelo de predicción. En este tutorial, descubrirá las series temporales de ruido blanco con Python.

En esta sección, crearemos una serie de ruido blanco gaussiano en Python y realizaremos algunas comprobaciones. Es útil crear y revisar una serie temporal de ruido blanco en la práctica. Proporcionará el marco de referencia y gráficos de ejemplo y pruebas estadísticas para utilizar y comparar en sus propios proyectos de series temporales para comprobar si son ruido blanco. En primer lugar, podemos crear una lista de 1.000 variables aleatorias gaussianas utilizando la función gauss() del módulo aleatorio1. Extraeremos variables de una distribución gaussiana con una media (mu) de 0,0 y una desviación típica (sigma) de 1,0. Una vez creada, podemos envolver la lista en una Serie Pandas por conveniencia.

##  10.4 Example of White Noise Time Series
from random import gauss
from random import seed
from pandas import Series
from pandas.plotting import autocorrelation_plot
from matplotlib import pyplot
# seed random number generator
seed(1)
# create white noise series
series = [gauss(0.0, 1.0) for i in range(1000)]
series = Series(series)
# summary stats
print(series.describe())
# line plot
series.plot()
pyplot.show()
pyplot.close()
# histogram plot
series.hist()
pyplot.show()
pyplot.close()
# autocorrelation
autocorrelation_plot(series)
pyplot.show()
pyplot.close()

count    1000.000000
mean       -0.013222
std         1.003685
min        -2.961214
25%        -0.684192
50%        -0.010934
75%         0.703915
max         2.737260
dtype: float64

Introducción al paseo aleatorio

¿Cómo saber si un problema de series temporales es predecible? Se trata de una pregunta difícil en la previsión de series temporales. Existe una herramienta llamada paseo aleatorio que puede ayudarte a entender la predictibilidad de tu problema de previsión de series temporales. En este tutorial, descubrirá el paseo aleatorio y sus propiedades en Python.

Series aleatorias

La biblioteca estándar de Python contiene el módulo random1 que proporciona acceso a un conjunto de funciones para generar números aleatorios. La función randrange()2 puede utilizarse para generar un número entero aleatorio entre 0 y un límite superior. Podemos utilizar la función randrange() para generar una lista de 1.000 números enteros aleatorios entre 0 y 10.

#!pip install statsmodels
# 11.1 Random Series
# create and plot a random series
from random import seed
from random import randrange
from matplotlib import pyplot
seed(1)
series = [randrange(10) for i in range(1000)]
pyplot.plot(series)
pyplot.show()
pyplot.close()

Paseo aleatorio

Un paseo aleatorio es diferente de una lista de números aleatorios porque el siguiente valor de la secuencia es una modificación del valor anterior de la secuencia. El proceso utilizado para generar la serie fuerza la dependencia de un paso temporal al siguiente. Esta dependencia proporciona cierta coherencia de un paso a otro en lugar de los grandes saltos que proporciona una serie de números aleatorios independientes. Es esta dependencia la que da nombre al proceso como paseo aleatorio o paseo del borracho.

# create and plot a random walk
from random import seed
from random import random
from matplotlib import pyplot
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)

for i in range(1,1000):
    movement = -1 if random() < 0.5 else 1
    value = random_walk[i-1] + movement
    random_walk.append(value)

pyplot.plot(random_walk)
pyplot.show()
pyplot.close()

Recorrido aleatorio y autocorrelación

Podemos calcular la correlación entre cada observación y las observaciones de los pasos anteriores. Un gráfico de estas correlaciones se denomina gráfico de autocorrelación o correlograma3. Dada la forma en que se construye el paseo aleatorio, cabría esperar una fuerte autocorrelación con la observación anterior y una caída lineal a partir de ahí con los valores de retardo anteriores. Podemos utilizar la función autocorrelation plot() en Pandas para trazar el correlograma del paseo aleatorio.

# plot the autocorrelation of a random walk
from random import seed
from random import random
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1000):
    movement = -1 if random() < 0.5 else 1
    value = random_walk[i-1] + movement
    random_walk.append(value)
autocorrelation_plot(random_walk)
pyplot.show()
pyplot.close()

Recorrido aleatorio y estacionariedad

Una serie temporal estacionaria es aquella en la que los valores no son una función del tiempo (la estacionariedad se trata con más detalle en el capítulo 15). Dada la forma en que se construye el paseo aleatorio y los resultados de revisar la autocorrelación, sabemos que las observaciones en un paseo aleatorio dependen del tiempo. La observación actual es un paso aleatorio desde la observación anterior.

# calculate the stationarity of a random walk
from random import seed
from random import random
from statsmodels.tsa.stattools import adfuller
# generate random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1,1000):
    movement = -1 if random() < 0.5 else 1
    value = random_walk[i-1] + movement
    random_walk.append(value)
print()
# statistical test
result = adfuller(random_walk)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))

ADF Statistic: 0.341605
p-value: 0.979175
Critical Values:
    1%: -3.437
    5%: -2.864
    10%: -2.568

Podemos hacer que el paseo aleatorio sea estacionario tomando la primera diferencia. Es decir, sustituyendo cada observación por la diferencia entre ésta y el valor anterior. Dada la forma en que se construyó este paseo aleatorio, esperaríamos que el resultado fuera una serie temporal de valores -1 y 1. Esto es exactamente lo que vemos. Esto es exactamente lo que vemos.

# calculate and plot a differenced random walk
from random import seed
from random import random
from matplotlib import pyplot
# create random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1000):
  movement = -1 if random() < 0.5 else 1
  value = random_walk[i-1] + movement
  random_walk.append(value)
# take difference
diff = list()
for i in range(1, len(random_walk)):
  value = random_walk[i] - random_walk[i - 1]
  diff.append(value)
# line plot
pyplot.plot(diff)
pyplot.show()
pyplot.close()
autocorrelation_plot(diff)
pyplot.show()

Predecir un paseo aleatorio

Un paseo aleatorio es impredecible; no se puede predecir razonablemente. Dada la forma en que se construye el paseo aleatorio, podemos esperar que la mejor predicción que podríamos hacer sería utilizar la observación en el paso de tiempo anterior como lo que sucederá en el siguiente paso de tiempo. Simplemente porque sabemos que el siguiente paso temporal será una función del paso temporal anterior. Esto suele denominarse previsión ingenua o modelo de persistencia.

# persistence forecasts for a random walk
from random import seed
from random import random
from sklearn.metrics import mean_squared_error
from math import sqrt
# generate the random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1, 1000):
  movement = -1 if random() < 0.5 else 1
  value = random_walk[i-1] + movement
  random_walk.append(value)
# prepare dataset
train_size = int(len(random_walk) * 0.66)
train, test = random_walk[0:train_size], random_walk[train_size:]
# persistence
predictions = list()
history = train[-1]
for i in range(len(test)):
  yhat = history
  predictions.append(yhat)
  history = test[i]
rmse = sqrt(mean_squared_error(test, predictions))
print('Persistence RMSE: %.3f' % rmse)

Persistence RMSE: 1.000

Otro error que cometen los principiantes del paseo aleatorio es suponer que si se conoce el rango de error (varianza), entonces podemos hacer predicciones utilizando un proceso del tipo de generación del paseo aleatorio. Es decir, si sabemos que el error es -1 o 1, ¿por qué no hacer predicciones añadiendo un -1 o 1 elegido al azar al valor anterior.

# random predictions for a random walk
from random import seed
from random import random
from sklearn.metrics import mean_squared_error
from math import sqrt
# generate the random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1000):
  movement = -1 if random() < 0.5 else 1
  value = random_walk[i-1] + movement
  random_walk.append(value)
# prepare dataset
train_size = int(len(random_walk) * 0.66)
train, test = random_walk[0:train_size], random_walk[train_size:]
# random prediction
predictions = list()
history = train[-1]
for i in range(len(test)):
  yhat = history + (-1 if random() < 0.5 else 1)
  predictions.append(yhat)
  history = test[i]
rmse = sqrt(mean_squared_error(test, predictions))
print('Random RMSE: %.3f' % rmse)

Random RMSE: 9.827

Descomponer datos de series temporales

La descomposición de series temporales consiste en considerar una serie como una combinación de componentes de nivel, tendencia, estacionalidad y ruido. La descomposición proporciona un modelo abstracto útil para pensar en las series temporales en general y para comprender mejor los problemas durante el análisis y la previsión de series temporales. En este tutorial, descubrirá la descomposición de series temporales y cómo dividir automáticamente una serie temporal en sus componentes con Python.

Descomposición automática de series temporales

Existen métodos para descomponer automáticamente una serie temporal. La biblioteca Statsmodels proporciona una implementación del método de descomposición ingenuo, o clásico, en una función denominada decompose() estacional. Requiere que se especifique si el modelo es aditivo o multiplicativo.

Descomposición aditiva

Podemos crear una serie temporal compuesta por una tendencia lineal creciente de 1 a 99 y algo de ruido aleatorio y descomponerla como un modelo aditivo. Dado que la serie temporal es artificial y se proporciona como una matriz de números, debemos especificar la frecuencia de las observaciones (argumento freq=1).

# additive decompose a contrived additive time series
from random import randrange
from matplotlib import pyplot
from statsmodels.tsa.seasonal import seasonal_decompose
series = [i+randrange(10) for i in range(1,100)]
result = seasonal_decompose(series, model='additive', period=1)
result.plot()
pyplot.show()
pyplot.close()

Descomposición multiplicativa

Podemos concebir una serie temporal cuadrática como un cuadrado del paso temporal de 1 a 99, y luego descomponerla asumiendo un modelo multiplicativo.

# multiplicative decompose a contrived multiplicative time series
from matplotlib import pyplot
from statsmodels.tsa.seasonal import seasonal_decompose
series = [i**2.0 for i in range(1,100)]
result = seasonal_decompose(series, model='multiplicative', period=1)
result.plot()
pyplot.show()
pyplot.close()

Airline Passengers Dataset

En esta lección, utilizaremos como ejemplo el conjunto de datos Pasajeros de líneas aéreas. Este conjunto de datos describe el número total de pasajeros de líneas aéreas a lo largo del tiempo.

# multiplicative decompose time series
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.tsa.seasonal import seasonal_decompose
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
result = seasonal_decompose(series, model='multiplicative')
result.plot()
pyplot.show()
pyplot.close()

Utilizar y eliminar tendencias

Nuestro conjunto de datos de series temporales puede contener una tendencia. Una tendencia es un aumento o disminución continuado de la serie a lo largo del tiempo. Puede ser beneficioso identificar, modelar e incluso eliminar la información de tendencias de su conjunto de datos de series temporales. En este tutorial, descubrirá cómo modelar y eliminar la información de tendencias de los datos de series temporales en Python.

Desfase por diferenciación

Tal vez el método más sencillo para detraer la tendencia de una serie temporal sea la diferenciación. Específicamente, se construye una nueva serie en la que el valor en el paso temporal actual se calcula como la diferencia entre la observación original y la observación en el paso temporal anterior.

# detrend a time series using differencing
from pandas import read_csv
from datetime import datetime
from matplotlib import pyplot
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
X = series.values
diff = list()
for i in range(1, len(X)):
  value = X[i] - X[i - 1]
  diff.append(value)
pyplot.plot(diff)
pyplot.show()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\1556502088.py:7: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

Detrend por ajuste de modelos

Una tendencia suele visualizarse fácilmente como una línea que atraviesa las observaciones. Las tendencias lineales pueden resumirse mediante un modelo lineal, y las tendencias no lineales pueden resumirse mejor utilizando un polinomio u otro método de ajuste de curvas. Debido a la naturaleza subjetiva y específica de la identificación de tendencias, este enfoque puede ayudar a identificar si una tendencia está presente. Incluso puede ser útil ajustar un modelo lineal a una tendencia que es claramente superlineal o exponencial.

# use a linear model to detrend a time series
from pandas import read_csv
from datetime import datetime
from sklearn.linear_model import LinearRegression
from matplotlib import pyplot
import numpy
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
# fit linear model
X = [i for i in range(0, len(series))]
X = numpy.reshape(X, (len(X), 1))
y = series.values
model = LinearRegression()
model.fit(X, y)
# calculate trend
trend = model.predict(X)
# plot trend
pyplot.plot(y)
pyplot.plot(trend)
pyplot.show()
# detrend
detrended = [y[i]-trend[i] for i in range(0, len(series))]
# plot detrended
pyplot.plot(detrended)
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\926202335.py:9: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

Utilizar y eliminar la estacionalidad

Los conjuntos de datos de series temporales pueden contener un componente estacional. Se trata de un ciclo que se repite a lo largo del tiempo, por ejemplo mensual o anualmente. Este ciclo repetitivo puede oscurecer la señal que deseamos modelizar al hacer previsiones y, a su vez, puede proporcionar una señal fuerte a nuestros modelos predictivos. En este tutorial, descubrirá cómo identificar y corregir la estacionalidad en los datos de series temporales con Python.

Ajuste estacional con diferenciación

Una forma sencilla de corregir un componente estacional es utilizar la diferenciación. Si existe un componente estacional a nivel de una semana, podemos eliminarlo en una observación de hoy restando el valor de la semana pasada. En el caso del conjunto de datos de las temperaturas mínimas diarias, parece que tenemos un componente estacional cada año que muestra el paso del verano al invierno.

# deseasonalize a time series using differencing
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
X = series.values
diff = list()
days_in_year = 365
for i in range(days_in_year, len(X)):
  value = X[i] - X[i - days_in_year]
  diff.append(value)
pyplot.plot(diff)
pyplot.show()
pyplot.close()

Otra opción es considerar que la temperatura en cualquier periodo del año es probablemente estable. Quizá a lo largo de varias semanas. Podemos abreviar esta idea y considerar que todas las temperaturas dentro de un mes natural son estables. Un modelo mejorado puede consistir en restar la temperatura media del mismo mes natural del año anterior, en lugar de la del mismo día. Podemos empezar o remuestreando el conjunto de datos a una temperatura mínima media mensual.

# calculate and plot monthly average
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
resample = series.resample('M')
monthly_mean = resample.mean()
print(monthly_mean.head(13))
monthly_mean.plot()
pyplot.show()
pyplot.close()

Date
1981-01-31    17.712903
1981-02-28    17.678571
1981-03-31    13.500000
1981-04-30    12.356667
1981-05-31     9.490323
1981-06-30     7.306667
1981-07-31     7.577419
1981-08-31     7.238710
1981-09-30    10.143333
1981-10-31    10.087097
1981-11-30    11.890000
1981-12-31    13.680645
1982-01-31    16.567742
Freq: ME, Name: Temp, dtype: float64

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\3501869392.py:6: FutureWarning:

'M' is deprecated and will be removed in a future version, please use 'ME' instead.

Podemos probar el mismo método de diferenciación en los datos mensuales y confirmar que el conjunto de datos desestacionalizados elimina efectivamente los ciclos anuales.

# deseasonalize monthly data by differencing
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
resample = series.resample('M')
monthly_mean = resample.mean()
X = series.values
diff = list()
months_in_year = 12
for i in range(months_in_year, len(monthly_mean)):
  value = monthly_mean[i] - monthly_mean[i - months_in_year]
  diff.append(value)
pyplot.plot(diff)
pyplot.show()
pyplot.close()

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\966464281.py:6: FutureWarning:

'M' is deprecated and will be removed in a future version, please use 'ME' instead.

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\966464281.py:12: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

A continuación, podemos utilizar las temperaturas mínimas medias mensuales del mismo mes del año anterior para ajustar el conjunto de datos de temperaturas mínimas diarias. Una vez más, nos limitamos a omitir el primer año de datos, pero la corrección utilizando los datos mensuales en lugar de los diarios puede ser un enfoque más estable.

# deseasonalize a time series using month-based differencing
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
X = series.values
diff = list()
days_in_year = 365
for i in range(days_in_year, len(X)):
  month_str = str(series.index[i].year-1)+'-'+str(series.index[i].month)
  month_mean_last_year = series[month_str].mean()
  value = X[i] - month_mean_last_year
  diff.append(value)
pyplot.plot(diff)
pyplot.show()
pyplot.close()

Ajuste estacional mediante modelización

Podemos modelizar directamente el componente estacional y, a continuación, sustraerlo de las observaciones. El componente estacional de una serie temporal dada es probablemente una onda sinusoidal con un periodo y una amplitud generalmente fijos. Puede aproximarse fácilmente mediante un método de ajuste de curvas.

# model seasonality with a polynomial model
from pandas import read_csv
from matplotlib import pyplot
from numpy import polyfit
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# fit polynomial: x^2*b1 + x*b2 + ... + bn
X = [i%365 for i in range(0, len(series))]
y = series.values
degree = 4
coef = polyfit(X, y, degree)
print('Coefficients: %s' % coef)
# create curve
curve = list()
for i in range(len(X)):
  value = coef[-1]
  for d in range(degree):
    value += X[i]**(degree-d) * coef[d]
  curve.append(value)
# plot curve over original data
pyplot.plot(series.values)
pyplot.plot(curve, color='red', linewidth=3)
pyplot.show()
pyplot.close()

Coefficients: [-1.17308000e-08  9.30253946e-06 -2.15977594e-03  1.19147966e-01
  1.38980178e+01]

La curva parece ajustarse bien a la estructura estacional del conjunto de datos. Ahora podemos utilizar este modelo para crear una versión desestacionalizada del conjunto de datos. El ejemplo completo figura a continuación.

# deseasonalize by differencing with a polynomial model
from pandas import read_csv
from matplotlib import pyplot
from numpy import polyfit
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# fit polynomial: x^2*b1 + x*b2 + ... + bn
X = [i%365 for i in range(0, len(series))]
y = series.values
degree = 4
coef = polyfit(X, y, degree)
# create curve
curve = list()
for i in range(len(X)):
  value = coef[-1]
  for d in range(degree):
    value += X[i]**(degree-d) * coef[d]
  curve.append(value)
# create seasonally adjusted
values = series.values
diff = list()
for i in range(len(values)):
  value = values[i] - curve[i]
  diff.append(value)
pyplot.plot(diff)
pyplot.show()
pyplot.close()

Estacionariedad de los datos de series temporales

Las series temporales son diferentes de los problemas más tradicionales de modelización predictiva de clasificación y regresión. La estructura temporal añade un orden a las observaciones. Este orden impuesto significa que deben manejarse de forma específica supuestos importantes sobre la coherencia de dichas observaciones. de forma específica. Por ejemplo, al modelizar, hay que suponer que los estadísticos de resumen de las observaciones son coherentes. En terminología de series temporales, nos referimos a esta expectativa como que las series temporales son estacionarias. Estas suposiciones pueden ser fácilmente violadas en series de tiempo por la adición de una tendencia, estacionalidad, y otras estructuras dependientes del tiempo. En este tutorial, descubrirá cómo comprobar si su serie temporal es estacionaria con Python.

Series temporales estacionarias

Las observaciones de una serie temporal estacionaria no dependen del tiempo. Las series temporales son estacionarias si no tienen tendencia ni efectos estacionales. Los estadísticos de resumen calculados sobre las series temporales son consistentes a lo largo del tiempo, como la media o la varianza de las observaciones. Cuando una serie temporal es estacionaria, puede ser más fácil de modelizar. Los métodos de modelización estadística presuponen o exigen que las series temporales sean estacionarias para ser eficaces. A continuación se muestra un ejemplo del conjunto de datos Nacimientos Femeninos Diarios que es estacionario.

# load time series data
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.plot()
pyplot.show()
pyplot.close()

Series temporales no estacionarias

Las observaciones de una serie temporal no estacionaria muestran efectos estacionales, tendencias y otras estructuras que dependen del índice temporal. Los estadísticos de resumen, como la media y la varianza, cambian con el tiempo, lo que proporciona una deriva en los conceptos que un modelo puede intentar captar. Los métodos clásicos de análisis y previsión de series temporales tratan de hacer estacionarios los datos de series temporales no estacionarios identificando y eliminando las tendencias y suprimiendo los efectos estacionales. A continuación se muestra un ejemplo de un conjunto de datos de pasajeros de líneas aéreas que no es estacionario y muestra componentes de tendencia y estacionales.

# load time series data
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.plot()
pyplot.show()
pyplot.close()

Estadísticas resumen

Una comprobación rápida y sucia para ver si su serie temporal es no estacionaria es revisar los estadísticos de resumen. Puede dividir su serie temporal en dos (o más) particiones y comparar la media y la varianza de cada grupo. Si difieren y la diferencia es estadísticamente significativa, es probable que la serie temporal no sea estacionaria.

Conjunto de datos de nacimientos diarios

Como estamos viendo la media y la varianza, suponemos que los datos se ajustan a una distribución gaussiana (también llamada curva de campana o normal). También podemos comprobarlo rápidamente observando un histograma de nuestras observaciones.

# plot a histogram of a time series
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.hist()
pyplot.show()
pyplot.close()

# calculate statistics of partitioned time series data
from pandas import read_csv
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
split = int(len(X) / 2)
X1, X2 = X[0:split], X[split:]
mean1, mean2 = X1.mean(), X2.mean()
var1, var2 = X1.var(), X2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))
print('variance1=%f, variance2=%f' % (var1, var2))

mean1=39.763736, mean2=44.185792
variance1=49.213410, variance2=48.708651

Conjunto de datos de pasajeros de líneas aéreas

Para ir al grano, podemos dividir nuestro conjunto de datos y calcular la media y la varianza de cada grupo.

# calculate statistics of partitioned time series data
from pandas import read_csv
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
split = int(len(X) / 2)
X1, X2 = X[0:split], X[split:]
mean1, mean2 = X1.mean(), X2.mean()
var1, var2 = X1.var(), X2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))
print('variance1=%f, variance2=%f' % (var1, var2))

# plot a histogram of a time series
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.hist()
pyplot.show()

mean1=182.902778, mean2=377.694444
variance1=2244.087770, variance2=7367.962191

Repasando de nuevo el gráfico de la serie temporal, podemos ver que hay un componente de estacionalidad obvio, y parece que el componente de estacionalidad está creciendo. Esto puede sugerir un crecimiento exponencial de una temporada a otra. Se puede utilizar una transformación logarítmica para aplanar el cambio exponencial y volver a una relación lineal. A continuación se muestra el mismo histograma con una transformación logarítmica de la serie temporal.

# histogram and line plot of log transformed time series
from pandas import read_csv
from matplotlib import pyplot
from numpy import log
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
X = log(X)
pyplot.hist(X)
pyplot.show()
pyplot.plot(X)
pyplot.show()

# calculate statistics of partitioned log transformed time series data
from pandas import read_csv
from numpy import log
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
X = log(X)
split = int(len(X) / 2)
X1, X2 = X[0:split], X[split:]
mean1, mean2 = X1.mean(), X2.mean()
var1, var2 = X1.var(), X2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))
print('variance1=%f, variance2=%f' % (var1, var2))

mean1=5.175146, mean2=5.909206
variance1=0.068375, variance2=0.049264

Prueba Dickey-Fuller aumentada

Las pruebas estadísticas se basan en hipótesis sólidas sobre los datos. Sólo pueden utilizarse para informar del grado en que puede rechazarse (o no rechazarse) una hipótesis nula. El resultado debe interpretarse para que un problema determinado tenga sentido. No obstante, pueden proporcionar una comprobación rápida y una prueba confirmatoria de que su serie temporal es estacionaria o no estacionaria.

La prueba Dickey-Fuller aumentada es un tipo de prueba estadística denominada prueba de raíz unitaria1. La intuición que subyace a una prueba de raíz unitaria es que determina hasta qué punto una serie temporal está definida por una tendencia.

Existen varias pruebas de raíz unitaria y la de Dickey-Fuller aumentada puede ser una de las más utilizadas. Utiliza un modelo autorregresivo y optimiza un criterio de información a través de múltiples valores de retardo diferentes. La hipótesis nula de la prueba es que la serie temporal puede ser representada por una raíz unitaria, que no es estacionaria (tiene alguna estructura dependiente del tiempo). La hipótesis alternativa (que rechaza la hipótesis nula) es que la serie temporal es estacionaria.

# calculate stationarity test of time series data
from pandas import read_csv
from statsmodels.tsa.stattools import adfuller
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
result = adfuller(X)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
  print('\t%s: %.3f' % (key, value))

ADF Statistic: -4.808291
p-value: 0.000052
Critical Values:
    1%: -3.449
    5%: -2.870
    10%: -2.571

# calculate stationarity test of time series data
from pandas import read_csv
from statsmodels.tsa.stattools import adfuller
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
result = adfuller(X)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
  print('\t%s: %.3f' % (key, value))

ADF Statistic: 0.815369
p-value: 0.991880
Critical Values:
    1%: -3.482
    5%: -2.884
    10%: -2.579

# calculate stationarity test of log transformed time series data
from pandas import read_csv
from statsmodels.tsa.stattools import adfuller
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
X = log(X)
result = adfuller(X)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
  print('\t%s: %.3f' % (key, value))

ADF Statistic: -1.717017
p-value: 0.422367
Critical Values:
    1%: -3.482
    5%: -2.884
    10%: -2.579

Modelos de predicción de backtest

El objetivo de la previsión de series temporales es realizar predicciones precisas sobre el futuro. Los métodos rápidos y potentes en los que confiamos en el aprendizaje automático, como el uso de divisiones entrenamiento-prueba y la validación cruzada k-fold, no funcionan en el caso de los datos de series temporales. Esto se debe a que ignoran los componentes temporales inherentes al problema. En este tutorial, descubrirá cómo evaluar modelos de aprendizaje automático en datos de series temporales con Python. En el campo de la previsión de series temporales, esto se denomina backtesting o hindcasting.

Data de entrenamiento-prueba

Puede dividir su conjunto de datos en subconjuntos de entrenamiento y de prueba. Su modelo puede prepararse en el conjunto de datos de entrenamiento y las predicciones pueden hacerse y evaluarse para el conjunto de datos de prueba. Para ello, seleccione un punto de división arbitrario en la lista ordenada de observaciones y cree dos nuevos conjuntos de datos. Dependiendo de la cantidad de datos de que disponga y de la cantidad de datos necesaria, puede utilizar divisiones de 50-50, 70-30 y 90-10. Es muy sencillo dividir datos en Python.

Después de cargar el conjunto de datos como una Serie Pandas, podemos extraer el array NumPy de valores de datos. El punto de división se puede calcular como un índice específico en el array. Todos los registros hasta el punto de división se toman como el conjunto de datos de entrenamiento y todos los registros desde el punto de división hasta el final de la lista de observaciones se toman como el conjunto de prueba. A continuación se muestra un ejemplo en Python con una división de 66-34

# calculate a train-test split of a time series dataset
from pandas import read_csv
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
train_size = int(len(X) * 0.66)
train, test = X[0:train_size], X[train_size:len(X)]
print('Observations: %d' % (len(X)))
print('Training Observations: %d' % (len(train)))
print('Testing Observations: %d' % (len(test)))

Observations: 2820
Training Observations: 1861
Testing Observations: 959

# plot train-test split of time series data
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
train_size = int(len(X) * 0.66)
train, test = X[0:train_size], X[train_size:len(X)]
print('Observations: %d' % (len(X)))
print('Training Observations: %d' % (len(train)))
print('Testing Observations: %d' % (len(test)))
pyplot.plot(train)
pyplot.plot([None for i in train] + [x for x in test])
pyplot.show()

Observations: 2820
Training Observations: 1861
Testing Observations: 959

Múltiples divisiones entrenamiento-prueba

Podemos repetir varias veces el proceso de dividir las series temporales en conjuntos de entrenamiento y de prueba. Esto requerirá el entrenamiento y la evaluación de múltiples modelos, pero este gasto computacional adicional proporcionará una estimación más robusta del rendimiento esperado del método y la configuración elegidos en datos no vistos. Podríamos hacerlo manualmente repitiendo el proceso descrito en la sección anterior con diferentes puntos de división.

# calculate repeated train-test splits of time series data
from pandas import read_csv
from sklearn.model_selection import TimeSeriesSplit
from matplotlib import pyplot
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
splits = TimeSeriesSplit(n_splits=3)
pyplot.figure(1)
index = 1
for train_index, test_index in splits.split(X):
  train = X[train_index]
  test = X[test_index]
  #print('Observations: %d' % (len(train) + len(test)))
  #print('Training Observations: %d' % (len(train)))
  #print('Testing Observations: %d' % (len(test)))
  pyplot.subplot(310 + index)
  pyplot.plot(train)
  pyplot.plot([None for i in train] + [x for x in test])
  index += 1
pyplot.show()

Validación progresiva

En la práctica, es muy probable que volvamos a entrenar nuestro modelo a medida que dispongamos de nuevos datos. Esto daría al modelo la mejor oportunidad de hacer buenas previsiones en cada paso temporal. Podemos evaluar nuestros modelos de aprendizaje automático bajo este supuesto. Dado que esta metodología implica moverse a lo largo de la serie temporal paso a paso, a menudo se denomina prueba de avance o validación de avance.

# walk forward evaluation model for time series data
from pandas import read_csv
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
n_train = 500
n_records = len(X)
for i in range(n_train, n_records):
  train, test = X[0:i], X[i:i+1]
print('train=%d, test=%d' % (len(train), len(test)))

train=2819, test=1

Medidas de rendimiento de las predicciones

Las medidas de rendimiento de la predicción de series temporales proporcionan un resumen de la habilidad y la capacidad del modelo de predicción que realizó las predicciones. Hay muchas medidas de rendimiento entre las que elegir. Puede resultar confuso saber qué medida utilizar y cómo interpretar los resultados. En este tutorial, descubrirá las medidas de rendimiento para evaluar las predicciones de series temporales con Python. Las series temporales generalmente se centran en la predicción de valores reales, llamados problemas de regresión. Por lo tanto, las medidas de rendimiento de este tutorial se centrarán en métodos para evaluar predicciones de valores reales.

# Forecast Error (or Residual Forecast Error)
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
print('Forecast Errors: %s' % forecast_errors)

# Mean Forecast Error (or Forecast Bias)
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
bias = sum(forecast_errors) * 1.0/len(expected)
print('Bias: %f' % bias)

# Mean Absolute Error
from sklearn.metrics import mean_absolute_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mae = mean_absolute_error(expected, predictions)
print('MAE: %f' % mae)

# Mean Squared Error
from sklearn.metrics import mean_squared_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
print('MSE: %f' % mse)

# Root Mean Squared Error
from sklearn.metrics import mean_squared_error
from math import sqrt
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
rmse = sqrt(mse)
print('RMSE: %f' % rmse)

Forecast Errors: [-0.2, 0.09999999999999998, -0.1, -0.09999999999999998, -0.2]
Bias: -0.100000
MAE: 0.140000
MSE: 0.022000
RMSE: 0.148324

Modelo de persistencia

Establecer una línea de base es esencial en cualquier problema de previsión de series temporales. Una línea de base en el rendimiento le da una idea de lo bien que todos los demás modelos funcionarán realmente en su problema. En este tutorial, descubrirá cómo desarrollar un pronóstico de persistencia que puede utilizar para calcular un nivel de referencia de rendimiento en un conjunto de datos de series temporales con Python.

Un modelo de persistencia supone que el valor futuro de una serie temporal se calcula bajo el supuesto de que nada cambia entre el tiempo actual y el tiempo previsto.

# evaluate a persistence forecast model
from pandas import read_csv
from datetime import datetime
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
from sklearn.metrics import mean_squared_error
from math import sqrt
# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
print(dataframe.head(5))
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
def model_persistence(x):
  return x
# walk-forward validation
predictions = list()
for x in test_X:
  yhat = model_persistence(x)
  predictions.append(yhat)
rmse = sqrt(mean_squared_error(test_y, predictions))
print('Test RMSE: %.3f' % rmse)
# plot predictions and expected results
pyplot.plot(train_y)
pyplot.plot([None for i in train_y] + [x for x in test_y])
pyplot.plot([None for i in train_y] + [x for x in predictions])
pyplot.show()

       t    t+1
0    NaN  266.0
1  266.0  145.9
2  145.9  183.1
3  183.1  119.3
4  119.3  180.3
Test RMSE: 133.156

C:\Users\WINDOWS 11\AppData\Local\Temp\ipykernel_3536\1367142013.py:12: FutureWarning:

The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.

Visualizar los errores de las predicciones

Los errores de predicción en problemas de regresión de series temporales se denominan residuos o errores residuales. Una exploración cuidadosa de los errores residuales en su problema de predicción de series temporales puede decirle mucho sobre su modelo de predicción e incluso sugerirle mejoras. En este tutorial, descubrirá cómo visualizar los errores residuales de las predicciones de series temporales.

Modelo de prediccion de persistencia

La previsión más sencilla que podemos hacer es prever que lo que ocurrió en el paso temporal anterior será lo mismo que ocurrirá en el paso temporal siguiente. Esto se denomina previsión ingenua o modelo de previsión de persistencia. Podemos implementar el modelo de persistencia en Python.

# calculate residuals from a persistence forecast
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
print(residuals.head())

      0
0   9.0
1 -10.0
2   3.0
3  -6.0
4  30.0

Gráfico de líneas residuales

El primer gráfico consiste en observar los errores de previsión residuales a lo largo del tiempo en forma de línea. Esperamos que el gráfico sea aleatorio alrededor del valor 0 y que no muestre ninguna tendencia o estructura cíclica. La matriz de errores residuales se puede envolver en un Pandas DataFrame y trazar directamente. El siguiente código proporciona un ejemplo.

# line plot of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
# plot residuals
residuals.plot()
pyplot.show()

Estadísticas resumen de los residuos

Podemos calcular estadísticas de resumen sobre los errores residuales. En primer lugar, nos interesa el valor medio de los errores residuales. Un valor cercano a cero sugiere que no hay sesgo en las previsiones, mientras que los valores positivos y negativos sugieren un sesgo positivo o negativo en las previsiones realizadas. Es útil conocer un sesgo en las previsiones, ya que puede corregirse directamente en las previsiones antes de su utilización o evaluación.

A continuación se muestra un ejemplo de cálculo de estadísticas resumidas de la distribución de errores residuales. Se incluyen la media y la desviación típica de la distribución, así como los percentiles y los errores mínimo y máximo observados.

# summary statistics of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
# summary statistics
print(residuals.describe())

                0
count  125.000000
mean     0.064000
std      9.187776
min    -28.000000
25%     -6.000000
50%     -1.000000
75%      5.000000
max     30.000000

Histograma y gráficos de densidad de residuos

Los gráficos pueden utilizarse para comprender mejor la distribución de los errores más allá de las estadísticas resumidas. Es de esperar que los errores de previsión se distribuyan normalmente en torno a una media cero. Los gráficos pueden ayudar a descubrir los sesgos de esta distribución. Podemos utilizar tanto histogramas como diagramas de densidad para comprender mejor la distribución de los errores residuales.

# density plots of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
# histogram plot
residuals.hist()
pyplot.show()
# density plot
residuals.plot(kind='kde')
pyplot.show()

Gráfico Q-Q residual

Un gráfico Q-Q, o gráfico de cuantiles, compara dos distribuciones y puede utilizarse para ver lo similares o diferentes que son. Podemos crear un gráfico Q-Q utilizando la función qqplot()1 de la biblioteca Statsmodels. El gráfico Q-Q puede utilizarse para comprobar rápidamente la normalidad de la distribución de los errores residuales. Los valores se ordenan y se comparan con una distribución idealizada de Gauss. La comparación se muestra como un gráfico de dispersión (teórica en el eje x y observada en el eje y) en el que una coincidencia entre las dos distribuciones se muestra como una línea diagonal desde la parte inferior izquierda hasta la parte superior derecha del gráfico. El gráfico es útil para detectar desviaciones evidentes de esta expectativa. A continuación se muestra un ejemplo de gráfico Q-Q de los errores residuales. El eje x muestra los cuantiles teóricos y el eje y muestra los cuantiles de la muestra.

# qq plot of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
import numpy
from statsmodels.graphics.gofplots import qqplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = numpy.array(residuals)
qqplot(residuals, line='r')
pyplot.show()

Gráfico de autocorrelación de residuos

La autocorrelación calcula la fuerza de la relación entre una observación y las observaciones en pasos de tiempo anteriores. Podemos calcular la autocorrelación de la serie temporal del error residual y representar gráficamente los resultados. Esto se denomina gráfico de autocorrelación. No es de esperar que exista correlación entre los residuos. Los resultados de la autocorrelación estarían por debajo del umbral de significación (líneas horizontales discontinuas en el gráfico).

Una autocorrelación significativa en el gráfico de residuos sugiere que el modelo podría estar haciendo un mejor trabajo a la hora de incorporar la relación entre las observaciones y las observaciones retardadas, lo que se denomina autoregresión. Pandas proporciona una función integrada para calcular un gráfico de autocorrelación, llamada autocorrelation plot().

# autoregression plot of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
autocorrelation_plot(residuals)
pyplot.show()

Modelos de autorregresión

La autorregresión es un modelo de series temporales que utiliza observaciones de pasos temporales anteriores como entrada de una ecuación de regresión para predecir el valor en el siguiente paso temporal. Se trata de una idea muy simple que puede dar lugar a previsiones precisas en una serie de problemas de series temporales. En este tutorial, descubrirá cómo implementar un modelo autorregresivo para la predicción de series temporales con Python.

Comprobación de la autocorrelación

Podemos realizar una comprobación visual rápida para ver si existe autocorrelación en nuestro conjunto de datos de series temporales. Podemos trazar la observación en el paso de tiempo actual (t) con la observación en el paso de tiempo anterior (t-1) como un gráfico de dispersión. Esto podría hacerse manualmente creando primero una versión con retardo del conjunto de datos de series temporales y utilizando una función de gráfico de dispersión incorporada en la biblioteca Pandas. Pero hay una manera más fácil. Pandas proporciona un gráfico integrado para hacer exactamente esto, llamado la función lag plot()1. A continuación se muestra un ejemplo de creación de un gráfico de retardo del conjunto de datos de Temperaturas Mínimas Diarias.

# lag plot of time series
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import lag_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
lag_plot(series)
pyplot.show()

La correlación puede calcularse fácilmente utilizando la función corr()2 en el DataFrame del conjunto de datos retardados. El siguiente ejemplo crea una versión retardada del conjunto de datos Temperaturas mínimas diarias y calcula una matriz de correlación de cada columna con otras columnas, incluida ella misma.

# correlation of lag=1
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
result = dataframe.corr()
print(result)

           t      t+1
t    1.00000  0.77487
t+1  0.77487  1.00000

Gráficos de autocorrelación

Podemos trazar el coeficiente de correlación de cada variable de retardo. Esto puede dar rápidamente una idea de qué variables de retardo pueden ser buenas candidatas para su uso en un modelo predictivo y de cómo cambia con el tiempo la relación entre la observación y sus valores históricos. Podríamos calcular manualmente los valores de correlación para cada variable de retardo y representar gráficamente el resultado. Afortunadamente, Pandas proporciona un gráfico integrado llamado autocorrelation plot() function.

# autocorrelation plot of time series
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
autocorrelation_plot(series)
pyplot.show()

# autocorrelation plot of time series
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_acf
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
plot_acf(series, lags=31)
pyplot.show()

Modelo de persistencia

Supongamos que queremos desarrollar un modelo para predecir los últimos 7 días de temperaturas mínimas en el conjunto de datos dadas todas las observaciones anteriores. El modelo más sencillo que podríamos utilizar para hacer predicciones sería persistir la última observación. Podemos llamarlo modelo de persistencia y proporciona una línea de base de rendimiento para el problema que podemos utilizar para la comparación con un modelo de autoregresión.

Podemos desarrollar un arnés de prueba para el problema dividiendo las observaciones en conjuntos de entrenamiento y de prueba, con sólo las 7 últimas observaciones del conjunto de datos asignadas al conjunto de prueba como datos no vistos que deseamos predecir. Las predicciones se realizan utilizando un modelo de validación walk-forward, de forma que podamos persistir en las observaciones más recientes para el día siguiente. Esto significa que no estamos haciendo una predicción de 7 días, sino 7 predicciones de 1 día.

# evaluate a persistence model
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
from sklearn.metrics import mean_squared_error
from math import sqrt
# load dataset
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train, test = X[1:len(X)-7], X[len(X)-7:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
def model_persistence(x):
  return x
# walk-forward validation
predictions = list()
for x in test_X:
  yhat = model_persistence(x)
  predictions.append(yhat)
rmse = sqrt(mean_squared_error(test_y, predictions))
print('Test RMSE: %.3f' % rmse)
# plot predictions vs expected
pyplot.plot(test_y)
pyplot.plot(predictions, color='red')
pyplot.show()

Test RMSE: 1.850

Modelo de autorregresión

Un modelo de autoregresión es un modelo de regresión lineal que utiliza variables retardadas como variables de entrada. Podríamos calcular el modelo de regresión lineal manualmente utilizando la clase LinearRegession en scikit-learn y especificar manualmente las variables de entrada retardadas a utilizar. Alternativamente, la biblioteca Statsmodels proporciona un modelo de autoregresión que selecciona automáticamente un valor de retardo apropiado utilizando pruebas estadísticas y entrena un modelo de regresión lineal. Se proporciona en la clase AR. Podemos utilizar este modelo creando primero el modelo AR() y luego llamando a fit() para entrenarlo en nuestro conjunto de datos. Esto devuelve un objeto AR Result. Una vez hecho esto, podemos utilizar el modelo para hacer una predicción llamando a la función predict() para un número de observaciones en el futuro. Esto crea 1 predicción de 7 días, que es diferente del ejemplo de persistencia anterior. A continuación se muestra el ejemplo completo.

# evaluate a persistence model
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
from sklearn.metrics import mean_squared_error
from math import sqrt
from statsmodels.tsa.ar_model import AutoReg
# load dataset
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# split dataset
X = series.values
train, test = X[1:len(X)-7], X[len(X)-7:]
# train autoregression
model = AutoReg(train, lags=7)
model_fit = model.fit()

print('Coefficients: %s' % model_fit.params)
# make predictions
predictions = model_fit.predict(start=len(train), end=len(train)+len(test)-1, dynamic=False)
for i in range(len(predictions)):
  print('predicted=%f, expected=%f' % (predictions[i], test[i]))
rmse = sqrt(mean_squared_error(test, predictions))
print('Test RMSE: %.3f' % rmse)
# plot results
pyplot.plot(test)
pyplot.plot(predictions, color='red')
pyplot.show()

Coefficients: [ 1.1171522   0.62645825 -0.07504654  0.07384666  0.06208368  0.06561174
  0.04428808  0.10242029]
predicted=11.549744, expected=12.900000
predicted=12.495221, expected=14.600000
predicted=12.703069, expected=14.000000
predicted=12.449430, expected=13.600000
predicted=12.226329, expected=13.500000
predicted=12.180301, expected=15.700000
predicted=11.893612, expected=13.000000
Test RMSE: 1.871

Forecasting de series temporales

A partir de esta seccion se describe cómo utilizar modelos de regresión de Scikit-learn para realizar forecasting de series temporales. En concreto, se hace uso de Skforecast, una librería que contiene las clases y funciones necesarias para adaptar cualquier modelo de regresión de Scikit-learn a problemas de forecasting.

Skforecast: forecasting series temporales con Python, Machine Learning y Scikit-learn por Joaquín Amat Rodrigo y Javier Escobar Ortiz.

import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme

# Modelado y Forecasting
# ==============================================================================
import sklearn
from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn.preprocessing import StandardScaler
import skforecast
from skforecast.recursive import ForecasterRecursive
from skforecast.direct import ForecasterDirect
from skforecast.model_selection import TimeSeriesFold, grid_search_forecaster, backtesting_forecaster
from skforecast.preprocessing import RollingFeatures
from skforecast.utils import save_forecaster, load_forecaster
from skforecast.metrics import calculate_coverage
from skforecast.plot import plot_prediction_intervals
import shap

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

C:\Users\WINDOWS 11\Documents\.virtualenvs\r-tensorflow\lib\site-packages\tqdm\auto.py:21: TqdmWarning:

IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html

Datos

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. Los datos empleados en los ejemplos de este documento se han obtenido del magnífico libro Forecasting: Principles and Practice by Rob J Hyndman and George Athanasopoulos.

# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True)

# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()
datos.head()

print(f'Número de filas con missing values: {datos.isnull().any(axis=1).mean()}')

# Verificar que un índice temporal está completo
# ==============================================================================
fecha_inicio = datos.index.min()
fecha_fin = datos.index.max()
date_range_completo = pd.date_range(start=fecha_inicio, end=fecha_fin, freq=datos.index.freq)
print(f"Índice completo: {(datos.index == date_range_completo).all()}")


# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]
print(f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})")
print(f"Fechas test  : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})")

set_dark_theme()
fig, ax = plt.subplots(figsize=(8, 4))

datos_train['y'].plot(ax=ax, label='train', fontsize=12)
datos_test['y'].plot(ax=ax, label='test', fontsize=12)
ax.legend();
fig

h2o_exog
--------
Monthly expenditure ($AUD) on corticosteroid drugs that the Australian health
system had between 1991 and 2008. Two additional variables (exog_1, exog_2) are
simulated.
Hyndman R (2023). fpp3: Data for Forecasting: Principles and Practice (3rd
Edition). http://pkg.robjhyndman.com/fpp3package/,
https://github.com/robjhyndman/fpp3package, http://OTexts.com/fpp3.
Shape of the dataset: (195, 4)
Número de filas con missing values: 0.0
Índice completo: True
Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00  (n=159)
Fechas test  : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00  (n=36)

Forecasting autorregresivo recursivo

Se crea y entrena un modelo ForecasterRecursive a partir de un regresor RandomForestRegressor y una ventana temporal de 6 lags. Esto último significa que, el modelo, utiliza como predictores los 6 meses anteriores.

# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
                regressor = RandomForestRegressor(random_state=123),
                lags = 6
             )
forecaster.fit(y=datos_train['y'])
forecaster

# Predicciones
# ==============================================================================
steps = 36
predicciones = forecaster.predict(steps=steps)
predicciones.head(5)

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

Error de test (mse): 0.07326833976120374

Ajuste de hiperparámetros (tuning)

El ForecasterRecursive 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. La librería Skforecast proporciona varias estrategias de búsqueda para encontrar la mejor combinación de hiperparámetros y lags. En este caso, se utiliza la función grid_search_forecaster, que compara los resultados obtenidos con cada combinación de hiperparámetros y lags, e identifica la mejor.

# Búsqueda de hiperparámetros: grid search
# ==============================================================================
forecaster = ForecasterRecursive(
                regressor = RandomForestRegressor(random_state=123),
                lags      = 12 # Este valor será remplazado en el grid search
             )

# Particiones de entrenamiento y validación
cv = TimeSeriesFold(
      steps              = 36,
      initial_train_size = int(len(datos_train) * 0.5),
      refit              = False,
      fixed_train_size   = False,

    )

# Valores candidatos de lags
lags_grid = [10, 20]

# Valores candidatos de hiperparámetros del regresor
param_grid = {
      'n_estimators': [100, 250],
      'max_depth': [3, 5, 10]
}

resultados_grid = grid_search_forecaster(
                        forecaster  = forecaster,
                        y           = datos_train['y'],
                        cv          = cv,
                        param_grid  = param_grid,
                        lags_grid   = lags_grid,
                        metric      = 'mean_squared_error',  
                        return_best = True,
                        n_jobs      = 1,
                        verbose     = False
                  )

# Resultados de la búsqueda de hiperparámetros
# ==============================================================================
resultados_grid

# Crear y entrenar forecaster con mejores hiperparámetros
# ==============================================================================
regressor = RandomForestRegressor(n_estimators=250, max_depth=3, random_state=123)
forecaster = ForecasterRecursive(
                regressor = regressor,
                lags      = 20
             )
forecaster.fit(y=datos_train['y'])

# Predicciones
# ==============================================================================
predicciones = forecaster.predict(steps=steps)

# Gráfico de predicciones vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos_train['y'].plot(ax=ax, label='train')
datos_test['y'].plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
plt.show()
# Error de test
# ==============================================================================
error_mse = mean_squared_error(
                y_true = datos_test['y'],
                y_pred = predicciones
            )
print(f"Error de test (mse) {error_mse}")

lags grid:   0%|          | 0/2 [00:00<?, ?it/s]
params grid:   0%|          | 0/6 [00:00<?, ?it/s]
params grid:  17%|█▋        | 1/6 [00:00<00:01,  3.56it/s]
params grid:  33%|███▎      | 2/6 [00:00<00:02,  1.97it/s]
params grid:  50%|█████     | 3/6 [00:01<00:01,  2.47it/s]
params grid:  67%|██████▋   | 4/6 [00:01<00:01,  1.95it/s]
params grid:  83%|████████▎ | 5/6 [00:02<00:00,  2.31it/s]
params grid: 100%|██████████| 6/6 [00:02<00:00,  1.92it/s]
                                                          lags grid:  50%|█████     | 1/2 [00:02<00:02,  2.89s/it]lags grid: 100%|██████████| 2/2 [00:05<00:00,  2.92s/it]lags grid: 100%|██████████| 2/2 [00:05<00:00,  2.92s/it]

`Forecaster` refitted using the best-found lags and 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] 
  Parameters: {'max_depth': 3, 'n_estimators': 250}
  Backtesting metric: 0.02177319540541341

Error de test (mse) 0.004356831371529945

Backtesting

Para obtener una estimación robusta de la capacidad predictiva del modelo, se lleva a cabo un proceso de backtesting. El backtesting consiste en evaluar el comportamiento de un modelo predictivo al aplicarlo de forma retrospectiva sobre datos históricos. Por lo tanto, es una estrategia de validación que permite cuantificar la capacidad predictiva de un modelo.

# Backtesting partitions using TimeSeriesFold
# ==============================================================================
cv = TimeSeriesFold(
         steps              = 12 * 3, 
         initial_train_size = len(datos['y']) - 12 * 9,  # Last 9 years are separated for the backtest
         window_size        = 20,
         fixed_train_size   = False,
         refit              = True,
)

cv.split(X=datos['y'], as_pandas=True)

# Backtesting
# ==============================================================================
cv = TimeSeriesFold(
         steps              = 12 * 3, 
         initial_train_size = len(datos) - 12 * 9,
         fixed_train_size   = False,
         refit              = True,
     )

metrica, predicciones_backtest = backtesting_forecaster(
                                    forecaster = forecaster,
                                    y          = datos['y'],
                                    cv         = cv,
                                    metric     = 'mean_squared_error',
                                    verbose    = True,
                                    n_jobs      = 1
                                 )
metrica

# Gráfico de predicciones de backtest vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos.loc[predicciones_backtest.index, 'y'].plot(ax=ax, label='test')
predicciones_backtest.plot(ax=ax, label='predicciones')
ax.legend();
plt.show()

Information of folds
--------------------
Number of observations used for initial training: 87
Number of observations used for backtesting: 108
    Number of folds: 3
    Number skipped folds: 0 
    Number of steps per fold: 36
    Number of steps to exclude between last observed data (last window) and predictions (gap): 0

Fold: 0
    Training:   1992-04-01 00:00:00 -- 1999-06-01 00:00:00  (n=87)
    Validation: 1999-07-01 00:00:00 -- 2002-06-01 00:00:00  (n=36)
Fold: 1
    Training:   1992-04-01 00:00:00 -- 2002-06-01 00:00:00  (n=123)
    Validation: 2002-07-01 00:00:00 -- 2005-06-01 00:00:00  (n=36)
Fold: 2
    Training:   1992-04-01 00:00:00 -- 2005-06-01 00:00:00  (n=159)
    Validation: 2005-07-01 00:00:00 -- 2008-06-01 00:00:00  (n=36)

Information of folds
--------------------
Number of observations used for initial training: 87
Number of observations used for backtesting: 108
    Number of folds: 3
    Number skipped folds: 0 
    Number of steps per fold: 36
    Number of steps to exclude between last observed data (last window) and predictions (gap): 0

Fold: 0
    Training:   1992-04-01 00:00:00 -- 1999-06-01 00:00:00  (n=87)
    Validation: 1999-07-01 00:00:00 -- 2002-06-01 00:00:00  (n=36)
Fold: 1
    Training:   1992-04-01 00:00:00 -- 2002-06-01 00:00:00  (n=123)
    Validation: 2002-07-01 00:00:00 -- 2005-06-01 00:00:00  (n=36)
Fold: 2
    Training:   1992-04-01 00:00:00 -- 2005-06-01 00:00:00  (n=159)
    Validation: 2005-07-01 00:00:00 -- 2008-06-01 00:00:00  (n=36)

  0%|          | 0/3 [00:00<?, ?it/s] 33%|███▎      | 1/3 [00:00<00:00,  3.97it/s] 67%|██████▋   | 2/3 [00:00<00:00,  2.95it/s]100%|██████████| 3/3 [00:01<00:00,  2.65it/s]100%|██████████| 3/3 [00:01<00:00,  2.79it/s]

Explicabilidad del modelo

Debido a la naturaleza compleja de muchos de los actuales modelos de machine learning, a menudo funcionan como cajas negras, lo que dificulta entender por qué han hecho una predicción u otra. Las técnicas de explicabilidad pretenden desmitificar estos modelos, proporcionando información sobre su funcionamiento interno y ayudando a generar confianza, mejorar la transparencia y cumplir los requisitos normativos en diversos ámbitos. Mejorar la explicabilidad de los modelos no sólo ayuda a comprender su comportamiento, sino también a identificar sesgos, mejorar su rendimiento y permitir a las partes interesadas tomar decisiones más informadas basadas en los conocimientos del machine learning.

Skforecast es compatible con algunos de los métodos de explicabilidad más populares: model-specific feature importances, SHAP values, and partial dependence plots.

Importancia model-specific

# Importancia predictores
# ==============================================================================
importancia = forecaster.get_feature_importances()
importancia.head(10)

	feature	importance
11	lag_12	0.815564
1	lag_2	0.086286
13	lag_14	0.019047
9	lag_10	0.013819
2	lag_3	0.012943
14	lag_15	0.009637
0	lag_1	0.009141
10	lag_11	0.008130
7	lag_8	0.007377
8	lag_9	0.005268

Shap values

Los valores SHAP (SHapley Additive exPlanations) son un método muy utilizado para explicar los modelos de machine learning, ya que ayudan a comprender cómo influyen las variables y los valores en las predicciones de forma visual y cuantitativa.

Se puede obtener un análisis SHAP a partir de modelos skforecast con sólo dos elementos:

• El regresor interno del forecaster.

• Las matrices de entrenamiento creadas a partir de la serie temporal y variables exógenas, utilizadas para ajustar el pronosticador.

Aprovechando estos dos componentes, los usuarios pueden crear explicaciones interpretables para sus modelos de skforecast. Estas explicaciones pueden utilizarse para verificar la fiabilidad del modelo, identificar los factores más significativos que contribuyen a las predicciones y comprender mejor la relación subyacente entre las variables de entrada y la variable objetivo.

# Matrices de entrenamiento utilizadas por el forecaster para entrenar el regresor
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(y=datos_train['y'])

# Crear SHAP explainer (para modelos basados en árboles)
# ==============================================================================
explainer = shap.TreeExplainer(forecaster.regressor)
# Se selecciona una muestra del 50% de los datos para acelerar el cálculo
rng = np.random.default_rng(seed=785412)
sample = rng.choice(X_train.index, size=int(len(X_train) * 0.5), replace=False)
X_train_sample = X_train.loc[sample, :]
shap_values = explainer.shap_values(X_train_sample)

# Shap summary plot (top 10)
# ==============================================================================
shap.initjs()
shap.summary_plot(shap_values, X_train_sample, max_display=15, show=False)
fig, ax = plt.gcf(), plt.gca()
ax.set_title("SHAP Summary plot")
ax.tick_params(labelsize=8)
fig.set_size_inches(8, 4)

plt.show()

# Backtesting indicando que se devuelvan los predictores
# ==============================================================================
cv = TimeSeriesFold(
         steps              = 12 * 3, 
         initial_train_size = len(datos) - 12 * 9,
         fixed_train_size   = False,
         refit              = True,
     )

_, predicciones_backtest = backtesting_forecaster(
                                    forecaster        = forecaster,
                                    y                 = datos['y'],
                                    cv                = cv,
                                    return_predictors = True,
                                    metric            = 'mean_squared_error',
                                    n_jobs      = 1
                                 )

predicciones_backtest.head()

# Waterfall para una predicción concreta
# ==============================================================================
iloc_predicted_date = predicciones_backtest.index.get_loc('2005-04-01')
shap_values_single = explainer(predicciones_backtest.iloc[:, 2:])
shap.plots.waterfall(shap_values_single[iloc_predicted_date], show=False)
fig, ax = plt.gcf(), plt.gca()
fig.set_size_inches(7, 4)
plt.show()

  0%|          | 0/3 [00:00<?, ?it/s] 33%|███▎      | 1/3 [00:00<00:00,  2.03it/s] 67%|██████▋   | 2/3 [00:01<00:00,  1.72it/s]100%|██████████| 3/3 [00:01<00:00,  1.61it/s]100%|██████████| 3/3 [00:01<00:00,  1.66it/s]

Forecasting 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. Esto mísmo es aplicable a múltiples variables exógenas.

# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True, verbose=False)

# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()

fig, ax = plt.subplots(figsize=(8, 4))
datos['y'].plot(ax=ax, label='y')
datos['exog_1'].plot(ax=ax, label='variable exógena')
ax.legend();
plt.show()
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]
print(
    f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})"
)
print(
    f"Fechas test  : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})"
)

# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
                regressor = RandomForestRegressor(random_state=123),
                lags      = 8
             )
forecaster.fit(y=datos_train['y'], exog=datos_train['exog_1'])

# Predicciones
# ==============================================================================
predicciones = forecaster.predict(steps=steps, exog=datos_test['exog_1'])

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

Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00  (n=159)
Fechas test  : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00  (n=36)

Error de test (mse): 0.03989087922533575

Predictores custom y window features

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. El argumento window_features permite incorporar al modelo predictores adicionales, creados a partir de los valores pasados de la serie temporal.

La clase RollingFeatures disponible en skforecast permite crear algunos de los predictores más utilizados:

• ‘mean’: la media de los n valores anteriores.

• ‘std’: la desviación estándar de los n valores anteriores.

• ‘min’: el mínimo de los n valores anteriores.

• ‘max’: el máximo de los n valores anteriores.

• ‘sum’: la suma de los n valores anteriores.

• ‘median’: la mediana de los n valores anteriores.

• ‘ratio_min_max’: la relación entre el mínimo y el máximo de los n valores anteriores.

• ‘coef_variation’: el coeficiente de variación de los n valores anteriores.

El usuario puede especificar un tamaño de ventana diferente para cada uno de ellos o el mismo para todos.

# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True, verbose=False)

# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()

# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]
print(f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})")
print(f"Fechas test  : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})")

# Window features
# ==============================================================================
window_features = RollingFeatures(
    stats = ['mean', 'std', 'min', 'max'],
    window_sizes = 20
)

# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
                regressor       = RandomForestRegressor(random_state=123),
                lags            = 10,
                window_features = window_features,
             )
forecaster.fit(y=datos_train['y'])
forecaster

# Matrices de entrenamiento
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(y=datos_train['y'])
X_train.head(5)
y_train.head(5)

# Predicciones
# ==============================================================================
steps = 36
predicciones = forecaster.predict(steps=steps)

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

Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00  (n=159)
Fechas test  : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00  (n=36)

Error de test (mse): 0.04180143590431811

Direct multi-step forecasting

Para conseguir predicciones de varios steps a futuro, los modelos ForecasterRecursive 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 varios modelos. Sin embargo, en algunos escenarios, consigue mejores resultados. Este tipo de modelos pueden obtenerse con la clase ForecasterDirect y pueden incluir también window features y variables exógenas.

# Error test

# ==============================================================================

error_mse = mean_squared_error(y_true = datos_test[‘y’], y_pred = predicciones)

print(f”Error de test (mse) {error_mse}“)

# Crear forecaster
# ==============================================================================
forecaster = ForecasterDirect(
                regressor     = Ridge(random_state=123),
                transformer_y = StandardScaler(),
                steps         = 36,
                lags          = 8 # Este valor será remplazado en el grid search
             )
# Búsqueda de hiperparámetros
# ==============================================================================
from skforecast.exceptions import LongTrainingWarning
warnings.simplefilter('ignore', category=LongTrainingWarning)

cv = TimeSeriesFold(
         steps              = 36, 
         initial_train_size = int(len(datos_train) * 0.5),
         fixed_train_size   = False,
         refit              = False,
     )

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

lags_grid = [5, 12, 20]

resultados_grid = grid_search_forecaster(
                    forecaster  = forecaster,
                    y           = datos_train['y'],
                    cv          = cv,
                    param_grid  = param_grid,
                    lags_grid   = lags_grid,
                    metric      = 'mean_squared_error',
                    return_best = True,
                    n_jobs      = 1
                )
# Resultados de la búsqueda de hiperparámetros
# ==============================================================================
resultados_grid.head()
# Predicciones
# ==============================================================================
predicciones = forecaster.predict()

lags grid:   0%|          | 0/3 [00:00<?, ?it/s]
params grid:   0%|          | 0/10 [00:00<?, ?it/s]
params grid:  30%|███       | 3/10 [00:00<00:00, 27.98it/s]
params grid:  60%|██████    | 6/10 [00:00<00:00, 28.33it/s]
params grid:  90%|█████████ | 9/10 [00:00<00:00, 28.40it/s]
                                                           lags grid:  33%|███▎      | 1/3 [00:00<00:00,  2.82it/s]lags grid:  67%|██████▋   | 2/3 [00:00<00:00,  2.85it/s]lags grid: 100%|██████████| 3/3 [00:01<00:00,  2.83it/s]lags grid: 100%|██████████| 3/3 [00:01<00:00,  2.83it/s]

`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12] 
  Parameters: {'alpha': 0.2782559402207126}
  Backtesting metric: 0.027413948265204574

plt.show()

# Gráfico predicciones vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos_train['y'].plot(ax=ax, label='train')
datos_test['y'].plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
plt.show()

# Error test
# ==============================================================================
error_mse = mean_squared_error(y_true = datos_test['y'], y_pred = predicciones)
print(f"Error de test (mse) {error_mse}")

Error de test (mse) 0.011792965469623133

Forecasting probabilístico

El forecasting probabilístico es una familia de técnicas que permite predecir la distribución esperada del de la serie en lugar de un único valor. Este tipo de pronóstico proporciona mucha más información, ya que permite la creación de intervalos de predicción.

Un intervalo de predicción define el espacio dentro del cual es de esperar que se encuentre el verdadero valor de yy con una determinada probabilidad. Por ejemplo, es de esperar que el intervalo de predicción (1, 99) contenga el verdadero valor de la predicción con un 98% de probabilidad.

Skforecast implementa varios métodos para el pronóstico probabilístico:

• Forecasting probabilístico: Bootstrapping

• Forecasting probabilístico: Predicción conformal

• Forecasting probabilístico: Regresión cuantílica

# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True, verbose=False)

# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()

# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test  = datos[-steps:]
print(f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})")
print(f"Fechas test  : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})")
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
               regressor = Ridge(alpha=0.1, random_state=765),
               lags = 15
             )
forecaster.fit(y=datos_train['y'], store_in_sample_residuals=True)

# Intervalos de predicción
# ==============================================================================
predicciones = forecaster.predict_interval(
                  steps    = steps,
                  interval = [5, 95],
                  method   = 'bootstrapping',
                  n_boot   = 150
               )
predicciones.head(5)
# Error de predicción
# ==============================================================================
error_mse = mean_squared_error(
                y_true = datos_test['y'],
                y_pred = predicciones['pred']
            )
print(f"Error de test (mse): {error_mse}")

# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
plot_prediction_intervals(
    predictions     = predicciones,
    y_true          = datos_test,
    target_variable = "y",
    ax              = ax
)
plt.show()
# Cobertura del intervalo predicho
# ==============================================================================
cobertura = calculate_coverage(
                y_true      = datos.loc[predicciones.index, 'y'],
                lower_bound = predicciones['lower_bound'],
                upper_bound = predicciones['upper_bound'],
            )
print(f"Cobertura del intervalo predicho: {round(100 * cobertura, 2)} %")

Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00  (n=159)
Fechas test  : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00  (n=36)
Error de test (mse): 0.010465086161791183

Cobertura del intervalo predicho: 83.33 %

Métrica custom

En los procesos de backtesting (backtesting_forecaster) y optimización de hiperparámetros (grid_search_forecaster), además de las métricas mean_squared_error, mean_absolute_error y mean_absolute_percentage_error, el usuario puede utilizar cualquier función que desee siempre y cuando cumpla lo siguiente:

• Tiene como argumentos:

  • y_true: verdaderos valores de la serie.

  • y_pred: valores predichos.

• Devuelve un valor numérico (float o int).

• El modelo es mejor cuanto menor es la métrica. Esto únicamente es necesario si se quiere que la función grid_search_forecaster reentrene automáticamente el mejor modelo encontrado.

Gracias a esta flexibilidad, es posible evaluar la capacidad predictiva del modelo con métricas aplicables a escenarios muy diversos. Por ejemplo:

• Considerar únicamente determinados meses, días u horas.

• Considerar únicamente fechas que sean festivos.

• Considerar únicamente el último step del horizonte predicho.

Véase un ejemplo en el que se quiere predecir un horizonte de 12 meses, pero únicamente considerar los últimos 3 meses de cada año para calcular la métrica de interés.

# Métrica custom 
# ==============================================================================
def custom_metric(y_true, y_pred):
    '''
    Calcular el mean_absolute_error utilizando únicamente las predicciones de
    los últimos 3 meses del año.
    '''
    mask = y_true.index.month.isin([10, 11, 12])
    metric = mean_absolute_error(y_true[mask], y_pred[mask])
    
    return metric
# Backtesting 
# ==============================================================================
metrica, predicciones_backtest = backtesting_forecaster(
                                    forecaster = forecaster,
                                    y          = datos['y'],
                                    cv         = cv,
                                    metric     = custom_metric,
                                    verbose    = True
                                 )
metrica

Information of folds
--------------------
Number of observations used for initial training: 79
Number of observations used for backtesting: 116
    Number of folds: 4
    Number skipped folds: 0 
    Number of steps per fold: 36
    Number of steps to exclude between last observed data (last window) and predictions (gap): 0
    Last fold only includes 8 observations.

Fold: 0
    Training:   1992-04-01 00:00:00 -- 1998-10-01 00:00:00  (n=79)
    Validation: 1998-11-01 00:00:00 -- 2001-10-01 00:00:00  (n=36)
Fold: 1
    Training:   No training in this fold
    Validation: 2001-11-01 00:00:00 -- 2004-10-01 00:00:00  (n=36)
Fold: 2
    Training:   No training in this fold
    Validation: 2004-11-01 00:00:00 -- 2007-10-01 00:00:00  (n=36)
Fold: 3
    Training:   No training in this fold
    Validation: 2007-11-01 00:00:00 -- 2008-06-01 00:00:00  (n=8)

  0%|          | 0/4 [00:00<?, ?it/s]100%|██████████| 4/4 [00:00<00:00, 532.61it/s]

	custom_metric
0	0.171414

Guardar y cargar modelos

Los modelos generados con Skforecast se pueden cargar y guardar usando las librerías Pickle o Joblib. Para facilitar el proceso, dos funciones están disponibles: save_forecaster y load_forecaster. A continuación, se muestra un sencillo ejemplo. Para más información cosultar: skforecast save and load forecaster.

# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(RandomForestRegressor(random_state=123), lags=3)
forecaster.fit(y=datos['y'])
forecaster.predict(steps=3)
# Guardar modelo
# ==============================================================================
save_forecaster(forecaster, file_name='forecaster.joblib', verbose=False)
# Cargar modelo
# ==============================================================================
forecaster_cargado = load_forecaster('forecaster.joblib')
# Predicciones
# ==============================================================================
forecaster_cargado.predict(steps=3)

=================== 
ForecasterRecursive 
=================== 
Regressor: RandomForestRegressor 
Lags: [1 2 3] 
Window features: None 
Window size: 3 
Series name: y 
Exogenous included: False 
Exogenous names: None 
Transformer for y: None 
Transformer for exog: None 
Weight function included: False 
Differentiation order: None 
Training range: [Timestamp('1992-04-01 00:00:00'), Timestamp('2008-06-01 00:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: MS 
Regressor parameters: 
    {'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'squared_error', 'max_depth':
    None, 'max_features': 1.0, 'max_leaf_nodes': None, 'max_samples': None,
    'min_impurity_decrease': 0.0, 'min_samples_leaf': 1, 'min_samples_split': 2,
    'min_weight_fraction_leaf': 0.0, 'monotonic_cst': None, 'n_estimators': 100,
    'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0,
    'warm_start': False} 
fit_kwargs: {} 
Creation date: 2025-05-29 11:50:35 
Last fit date: 2025-05-29 11:50:35 
Skforecast version: 0.16.0 
Python version: 3.10.11 
Forecaster id: None 

2008-07-01    0.751967
2008-08-01    0.826505
2008-09-01    0.879444
Freq: MS, Name: pred, dtype: float64

Uso de modelos en producción

En los proyectos relacionados con forecasting es frecuente que, como resultado de la etapa de experimentación y desarrollo, se genere un modelo. Para que este modelo consiga un impacto real en el negocio, se tiene que poder poner en producción y generar predicciones cada cierto tiempo, con las que tomar decisiones. Esta necesidad ha guiado en gran medida el desarrollo de la librería Skforecast.

Supóngase un caso de uso en el que se han de generar predicciones de forma semanal, por ejemplo, cada lunes el modelo tiene que predecir el resto de la semana. Una forma de conseguir este comportamiento es reentrenando el modelo semanalmente justo antes de que se ejecute la primera predicción y llamar a continuación al método predict del objeto forecaster.

# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
                regressor = RandomForestRegressor(random_state=123),
                lags = 6
             )
forecaster.fit(y=datos_train['y'])
# Predecir con last_window
# ==============================================================================
last_window = datos_test['y'][-6:]
forecaster.predict(last_window=last_window, steps=4)

2008-07-01    0.757750
2008-08-01    0.836313
2008-09-01    0.877668
2008-10-01    0.911734
Freq: MS, Name: pred, dtype: float64