Cómo realizar una búsqueda en cuadrícula de hiperparámetros SARIMA para el pronóstico de series temporales

El modelo de media móvil integrada autorregresiva estacional, o SARIMA, es un enfoque para modelar datos de series temporales univariadas que pueden contener componentes estacionales y de tendencia.

Es un enfoque eficaz para el pronóstico de series de tiempo, aunque requiere un análisis cuidadoso y experiencia en el dominio para configurar los siete o más hiperparámetros del modelo.

Un enfoque alternativo para configurar el modelo que utiliza hardware moderno rápido y paralelo es buscar en red un conjunto de configuraciones de hiperparámetros para descubrir cuál funciona mejor. A menudo, este proceso puede revelar configuraciones del modelo no intuitivas que dan como resultado un error de pronóstico menor que aquellas configuraciones especificadas mediante un análisis cuidadoso.

En este tutorial, descubrirá cómo desarrollar un marco para la búsqueda en cuadrícula de todos los hiperparámetros del modelo SARIMA para el pronóstico de series temporales univariadas.

Después de completar este tutorial, sabrá:

• Cómo desarrollar un marco para la búsqueda de cuadrículas en modelos SARIMA desde cero mediante validación directa.
• Cómo realizar una búsqueda en cuadrícula de hiperparámetros del modelo SARIMA para datos de series temporales diarias de nacimientos.
• Cómo realizar una búsqueda en cuadrícula de hiperparámetros del modelo SARIMA para datos de series temporales mensuales de ventas de champú, ventas de automóviles y temperatura.

Pon en marcha tu proyecto con mi nuevo libro Deep Learning for Time Series Forecasting, que incluye tutoriales paso a paso y los archivos de código fuente de Python. para todos los ejemplos.

Empecemos.

• Actualizado en abril de 2019: se actualizó el enlace al conjunto de datos.

Descripción general del tutorial

Este tutorial se divide en seis partes; ellos son:

1. SARIMA para pronóstico de series temporales
2. Desarrollar un marco de búsqueda de cuadrícula
3. Estudio de caso 1: Sin tendencia ni estacionalidad
4. Estudio de caso 2: Tendencia
5. Estudio de caso 3: estacionalidad
6. Estudio de caso 4: Tendencia y estacionalidad

SARIMA para pronóstico de series temporales

La media móvil integrada autorregresiva estacional, SARIMA o ARIMA estacional, es una extensión de ARIMA que admite explícitamente datos de series temporales univariadas con un componente estacional.

Agrega tres nuevos hiperparámetros para especificar la autorregresión (AR), la diferenciación (I) y la media móvil (MA) para el componente estacional de la serie, así como un parámetro adicional para el período de la estacionalidad.

Un modelo ARIMA estacional se forma incluyendo términos estacionales adicionales en el ARIMA […] La parte estacional del modelo consta de términos que son muy similares a los componentes no estacionales del modelo, pero implican retrocesos del período estacional.

— Página 242, Previsión: principios y práctica, 2013.

Configurar un SARIMA requiere seleccionar hiperparámetros para los elementos de tendencia y estacionales de la serie.

Hay tres elementos de tendencia que requieren configuración.

Son iguales que el modelo ARIMA; específicamente:

• p: Orden de autorregresión de tendencia.
• d: Orden de diferencia de tendencia.
• q: Orden de media móvil de tendencia.

Hay cuatro elementos estacionales que no forman parte de ARIMA y que deben configurarse; ellos son:

• P: Orden autorregresivo estacional.
• D: Orden de diferencia estacional.
• Q: Orden de media móvil estacional.
• m: el número de pasos de tiempo para un único período estacional.

En conjunto, la notación para un modelo SARIMA se especifica como:

SARIMA(p,d,q)(P,D,Q)m

El modelo SARIMA puede incluir los modelos ARIMA, ARMA, AR y MA a través de los parámetros de configuración del modelo.

Los hiperparámetros de tendencia y estacionales del modelo se pueden configurar analizando gráficos de autocorrelación y autocorrelación parcial, y esto puede requerir algo de experiencia.

Un enfoque alternativo es buscar en una cuadrícula un conjunto de configuraciones de modelos y descubrir qué configuraciones funcionan mejor para una serie temporal univariada específica.

Los modelos ARIMA estacionales pueden tener potencialmente una gran cantidad de parámetros y combinaciones de términos. Por lo tanto, es apropiado probar una amplia gama de modelos al ajustar los datos y elegir el modelo que mejor se ajuste utilizando un criterio apropiado...

— Páginas 143-144, Series temporales introductorias con R, 2009.

Este enfoque puede ser más rápido en las computadoras modernas que un proceso de análisis y puede revelar hallazgos sorprendentes que podrían no ser obvios y dar como resultado un menor error de pronóstico.

Desarrollar un marco de búsqueda de cuadrícula

En esta sección, desarrollaremos un marco para la búsqueda de cuadrículas de hiperparámetros del modelo SARIMA para un problema de pronóstico de series de tiempo univariado determinado.

Usaremos la implementación de SARIMA proporcionada por la biblioteca statsmodels.

Este modelo cuenta con hiperparámetros que controlan la naturaleza del modelo realizado para la serie, tendencia y estacionalidad, específicamente:

• orden: Tupla de parámetros p, d y q para el modelado de la tendencia.
• sesonal_order: una tupla de parámetros P, D, Q y m para modelar la estacionalidad.
• tendencia: parámetro para controlar un modelo de tendencia determinista como uno de 'n','c','t','ct' sin tendencia, constante, lineal y constante con tendencia lineal. , respectivamente.

Si sabe lo suficiente sobre su problema como para especificar uno o más de estos parámetros, entonces debería especificarlos. De lo contrario, puede intentar buscar en la cuadrícula estos parámetros.

Podemos comenzar definiendo una función que se ajuste a un modelo con una configuración determinada y haga un pronóstico de un solo paso.

El sarima_forecast() a continuación implementa este comportamiento.

La función toma una matriz o lista de observaciones anteriores contiguas y una lista de parámetros de configuración utilizados para configurar el modelo, específicamente dos tuplas y una cadena para el orden de tendencia, la tendencia de orden estacional y el parámetro.

También intentamos hacer que el modelo sea robusto relajando restricciones, como que los datos deben ser estacionarios y que la transformada MA sea invertible.

# one-step sarima forecast
def sarima_forecast(history, config):
	order, sorder, trend = config
	# define model
	model = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)
	# fit model
	model_fit = model.fit(disp=False)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

A continuación, necesitamos crear algunas funciones para ajustar y evaluar un modelo repetidamente mediante validación directa, incluida la división de un conjunto de datos en conjuntos de tren y prueba y la evaluación de pronósticos de un solo paso.

Podemos dividir una lista o una matriz NumPy de datos usando un segmento dado un tamaño específico de división, p. el número de pasos de tiempo que se utilizarán a partir de los datos del conjunto de prueba.

La siguiente función train_test_split() implementa esto para un conjunto de datos proporcionado y un número específico de pasos de tiempo para usar en el conjunto de prueba.

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

Una vez realizados los pronósticos para cada paso del conjunto de datos de prueba, es necesario compararlos con el conjunto de prueba para calcular una puntuación de error.

Existen muchas puntuaciones de error populares para el pronóstico de series temporales. En este caso usaremos el error cuadrático medio (RMSE), pero puedes cambiarlo a tu medida preferida, p.e. MAPE, MAE, etc.

La función measure_rmse() a continuación calculará el RMSE dada una lista de valores reales (el conjunto de prueba) y pronosticados.

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

Ahora podemos implementar el esquema de validación walk-forward. Este es un enfoque estándar para evaluar un modelo de pronóstico de series temporales que respeta el orden temporal de las observaciones.

En primer lugar, un conjunto de datos de series temporales univariadas proporcionado se divide en conjuntos de entrenamiento y de prueba mediante la función train_test_split(). Luego se enumera el número de observaciones en el conjunto de prueba. Para cada uno ajustamos un modelo en toda la historia y hacemos un pronóstico de un paso. Luego, la observación verdadera para el paso de tiempo se agrega al historial y se repite el proceso. La función sarima_forecast() se llama para ajustar un modelo y hacer una predicción. Finalmente, se calcula una puntuación de error comparando todos los pronósticos de un paso con el conjunto de pruebas real llamando a la función measure_rmse().

La siguiente función walk_forward_validation() implementa esto, tomando una serie de tiempo univariada, una cantidad de pasos de tiempo para usar en el conjunto de prueba y una matriz de configuración del modelo.

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = sarima_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

Si está interesado en realizar predicciones de varios pasos, puede cambiar la llamada a predict() en la función sarima_forecast() y también cambiar el cálculo del error en el < funciónmeasure_rmse().

Podemos llamar a walk_forward_validation() repetidamente con diferentes listas de configuraciones de modelos.

Un posible problema es que es posible que algunas combinaciones de configuraciones del modelo no se soliciten para el modelo y generarán una excepción, p. especificando algunos pero no todos los aspectos de la estructura estacional de los datos.

Además, algunos modelos también pueden generar advertencias sobre algunos datos, p. de las bibliotecas de álgebra lineal llamadas por la biblioteca statsmodels.

Podemos detectar excepciones e ignorar advertencias durante la búsqueda en la cuadrícula envolviendo todas las llamadas a walk_forward_validation() con un try-except y un bloque para ignorar advertencias. También podemos agregar soporte de depuración para desactivar estas protecciones en caso de que queramos ver qué está pasando realmente. Finalmente, si ocurre un error, podemos devolver un resultado Ninguno; de lo contrario, podemos imprimir información sobre la habilidad de cada modelo evaluado. Esto resulta útil cuando se evalúa una gran cantidad de modelos.

La siguiente función score_model() implementa esto y devuelve una tupla de (clave y resultado), donde la clave es una versión de cadena de la configuración del modelo probado.

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

A continuación, necesitamos un bucle para probar una lista de diferentes configuraciones de modelos.

Esta es la función principal que impulsa el proceso de búsqueda de la cuadrícula y llamará a la función score_model() para cada configuración del modelo.

Podemos acelerar drásticamente el proceso de búsqueda de la red evaluando las configuraciones del modelo en paralelo. Una forma de hacerlo es utilizar la biblioteca Joblib.

Podemos definir un objeto paralelo con la cantidad de núcleos a usar y configurarlo según la cantidad de puntuaciones detectadas en su hardware.

executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')

Luego podemos crear una lista de tareas para ejecutar en paralelo, que será una llamada a la función score_model() para cada configuración de modelo que tengamos.

tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)

Finalmente, podemos usar el objeto Parallel para ejecutar la lista de tareas en paralelo.

scores = executor(tasks)

Eso es todo.

También podemos proporcionar una versión no paralela para evaluar todas las configuraciones del modelo en caso de que queramos depurar algo.

scores = [score_model(data, n_test, cfg) for cfg in cfg_list]

El resultado de evaluar una lista de configuraciones será una lista de tuplas, cada una con un nombre que resume una configuración de modelo específica y el error del modelo evaluado con esa configuración como RMSE o Ninguno si hubo un error.

Podemos filtrar todas las puntuaciones con Ninguno.

scores = [r for r in scores if r[1] != None]

Luego podemos ordenar todas las tuplas de la lista por puntuación en orden ascendente (las mejores son las primeras) y luego devolver esta lista de puntuaciones para su revisión.

La función grid_search() siguiente implementa este comportamiento dado un conjunto de datos de series temporales univariadas, una lista de configuraciones de modelo (lista de listas) y la cantidad de pasos de tiempo que se usarán en el conjunto de prueba. Un argumento paralelo opcional permite activar o desactivar la evaluación de modelos en todos los núcleos y está activado de forma predeterminada.

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

Ya casi hemos terminado.

Lo único que queda por hacer es definir una lista de configuraciones de modelo para probar en un conjunto de datos.

Podemos definir esto genéricamente. El único parámetro que quizás queramos especificar es la periodicidad del componente estacional de la serie, si existe. De forma predeterminada, no asumiremos ningún componente estacional.

La función sarima_configs() a continuación creará una lista de configuraciones del modelo para evaluar.

Las configuraciones suponen que cada uno de los componentes AR, MA e I para tendencia y estacionalidad son de orden inferior, p. apagado (0) o en [1,2]. Es posible que desee ampliar estos rangos si cree que el orden puede ser mayor. Se puede especificar una lista opcional de períodos estacionales e incluso podría cambiar la función para especificar otros elementos que pueda conocer sobre su serie temporal.

En teoría, hay 1296 configuraciones de modelo posibles para evaluar, pero en la práctica, muchas no serán válidas y darán como resultado un error que atraparemos e ignoraremos.

# create a set of sarima configs to try
def sarima_configs(seasonal=[0]):
	models = list()
	# define config lists
	p_params = [0, 1, 2]
	d_params = [0, 1]
	q_params = [0, 1, 2]
	t_params = ['n','c','t','ct']
	P_params = [0, 1, 2]
	D_params = [0, 1]
	Q_params = [0, 1, 2]
	m_params = seasonal
	# create config instances
	for p in p_params:
		for d in d_params:
			for q in q_params:
				for t in t_params:
					for P in P_params:
						for D in D_params:
							for Q in Q_params:
								for m in m_params:
									cfg = [(p,d,q), (P,D,Q,m), t]
									models.append(cfg)
	return models

Ahora tenemos un marco para la búsqueda en cuadrícula de hiperparámetros del modelo SARIMA mediante una validación directa de un solo paso.

Es genérico y funcionará para cualquier serie temporal univariada en memoria proporcionada como una lista o matriz NumPy.

Podemos asegurarnos de que todas las piezas funcionen juntas probándolo en un conjunto de datos artificial de 10 pasos.

El ejemplo completo se enumera a continuación.

# grid search sarima hyperparameters
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error

# one-step sarima forecast
def sarima_forecast(history, config):
	order, sorder, trend = config
	# define model
	model = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)
	# fit model
	model_fit = model.fit(disp=False)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = sarima_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of sarima configs to try
def sarima_configs(seasonal=[0]):
	models = list()
	# define config lists
	p_params = [0, 1, 2]
	d_params = [0, 1]
	q_params = [0, 1, 2]
	t_params = ['n','c','t','ct']
	P_params = [0, 1, 2]
	D_params = [0, 1]
	Q_params = [0, 1, 2]
	m_params = seasonal
	# create config instances
	for p in p_params:
		for d in d_params:
			for q in q_params:
				for t in t_params:
					for P in P_params:
						for D in D_params:
							for Q in Q_params:
								for m in m_params:
									cfg = [(p,d,q), (P,D,Q,m), t]
									models.append(cfg)
	return models

if __name__ == '__main__':
	# define dataset
	data = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0]
	print(data)
	# data split
	n_test = 4
	# model configs
	cfg_list = sarima_configs()
	# grid search
	scores = grid_search(data, cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

Al ejecutar el ejemplo, primero se imprime el conjunto de datos de serie temporal artificial.

A continuación, se informan las configuraciones del modelo y sus errores a medida que se evalúan, truncadas a continuación para mayor brevedad.

Finalmente, se informan las configuraciones y el error de las tres configuraciones principales. Podemos ver que muchos modelos logran un rendimiento perfecto en este simple problema artificial de series de tiempo que aumenta linealmente.

[10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0]

...
 > Model[[(2, 0, 0), (2, 0, 0, 0), 'ct']] 0.001
 > Model[[(2, 0, 0), (2, 0, 1, 0), 'ct']] 0.000
 > Model[[(2, 0, 1), (0, 0, 0, 0), 'n']] 0.000
 > Model[[(2, 0, 1), (0, 0, 1, 0), 'n']] 0.000
done

[(2, 1, 0), (1, 0, 0, 0), 'n'] 0.0
[(2, 1, 0), (2, 0, 0, 0), 'n'] 0.0
[(2, 1, 1), (1, 0, 1, 0), 'n'] 0.0

Ahora que tenemos un marco sólido para la búsqueda de hiperparámetros del modelo SARIMA en cuadrículas, probémoslo en un conjunto de conjuntos de datos de series temporales univariadas estándar.

Los conjuntos de datos se eligieron con fines de demostración; No estoy sugiriendo que un modelo SARIMA sea el mejor enfoque para cada conjunto de datos; quizás en algunos casos sería más apropiado un ETS u otra cosa.

Estudio de caso 1: Sin tendencia ni estacionalidad

El conjunto de datos sobre “nacimientos femeninos diarios” resume el total diario de nacimientos femeninos en California, EE. UU., en 1959.

El conjunto de datos no tiene ninguna tendencia obvia ni componente estacional.

Descargue el conjunto de datos directamente desde aquí:

• total-diario-nacimientos-femeninos.csv

Guarde el archivo con el nombre 'daily-total-female-births.csv' en su directorio de trabajo actual.

Podemos cargar este conjunto de datos como una serie de Pandas usando la función read_csv().

series = read_csv('daily-total-female-births.csv', header=0, index_col=0)

El conjunto de datos tiene un año o 365 observaciones. Usaremos los primeros 200 para entrenamiento y los 165 restantes como conjunto de prueba.

A continuación se muestra el ejemplo completo de cuadrícula que busca el problema de pronóstico de series temporales univariadas femeninas diarias.

# grid search sarima hyperparameters for daily female dataset
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error
from pandas import read_csv

# one-step sarima forecast
def sarima_forecast(history, config):
	order, sorder, trend = config
	# define model
	model = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)
	# fit model
	model_fit = model.fit(disp=False)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = sarima_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of sarima configs to try
def sarima_configs(seasonal=[0]):
	models = list()
	# define config lists
	p_params = [0, 1, 2]
	d_params = [0, 1]
	q_params = [0, 1, 2]
	t_params = ['n','c','t','ct']
	P_params = [0, 1, 2]
	D_params = [0, 1]
	Q_params = [0, 1, 2]
	m_params = seasonal
	# create config instances
	for p in p_params:
		for d in d_params:
			for q in q_params:
				for t in t_params:
					for P in P_params:
						for D in D_params:
							for Q in Q_params:
								for m in m_params:
									cfg = [(p,d,q), (P,D,Q,m), t]
									models.append(cfg)
	return models

if __name__ == '__main__':
	# load dataset
	series = read_csv('daily-total-female-births.csv', header=0, index_col=0)
	data = series.values
	print(data.shape)
	# data split
	n_test = 165
	# model configs
	cfg_list = sarima_configs()
	# grid search
	scores = grid_search(data, cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

La ejecución del ejemplo puede tardar unos minutos en hardware moderno.

Las configuraciones del modelo y el RMSE se imprimen a medida que se evalúan los modelos. Las tres configuraciones principales del modelo y su error se informan al final de la ejecución.

Nota: Sus resultados pueden variar dada la naturaleza estocástica del algoritmo o procedimiento de evaluación, o diferencias en la precisión numérica. Considere ejecutar el ejemplo varias veces y comparar el resultado promedio.

Podemos ver que el mejor resultado fue un RMSE de unos 6,77 nacimientos con la siguiente configuración:

• Orden: (1, 0, 2)
• Orden estacional: (1, 0, 1, 0)
• Parámetro de tendencia: 't' para tendencia lineal

Es sorprendente que una configuración con algunos elementos estacionales diera como resultado el error más bajo. No habría adivinado esta configuración y probablemente me habría quedado con un modelo ARIMA.

...
> Model[[(2, 1, 2), (1, 0, 1, 0), 'ct']] 6.905
> Model[[(2, 1, 2), (2, 0, 0, 0), 'ct']] 7.031
> Model[[(2, 1, 2), (2, 0, 1, 0), 'ct']] 6.985
> Model[[(2, 1, 2), (1, 0, 2, 0), 'ct']] 6.941
> Model[[(2, 1, 2), (2, 0, 2, 0), 'ct']] 7.056
done

[(1, 0, 2), (1, 0, 1, 0), 't'] 6.770349800255089
[(0, 1, 2), (1, 0, 2, 0), 'ct'] 6.773217122759515
[(2, 1, 1), (2, 0, 2, 0), 'ct'] 6.886633191752254

Estudio de caso 2: Tendencia

El conjunto de datos sobre "champú" resume las ventas mensuales de champú durante un período de tres años.

El conjunto de datos contiene una tendencia obvia pero ningún componente estacional obvio.

Descargue el conjunto de datos directamente desde aquí:

• champú.csv

Guarde el archivo con el nombre 'shampoo.csv' en su directorio de trabajo actual.

Podemos cargar este conjunto de datos como una serie de Pandas usando la función read_csv().

# parse dates
def custom_parser(x):
	return datetime.strptime('195'+x, '%Y-%m')

# load dataset
series = read_csv('shampoo.csv', header=0, index_col=0, date_parser=custom_parser)

El conjunto de datos tiene tres años, o 36 observaciones. Usaremos los primeros 24 para entrenamiento y los 12 restantes como conjunto de prueba.

A continuación se muestra la cuadrícula de ejemplo completa que busca el problema de pronóstico de series temporales univariadas de ventas de champú.

# grid search sarima hyperparameters for monthly shampoo sales dataset
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error
from pandas import read_csv
from pandas import datetime

# one-step sarima forecast
def sarima_forecast(history, config):
	order, sorder, trend = config
	# define model
	model = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)
	# fit model
	model_fit = model.fit(disp=False)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = sarima_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of sarima configs to try
def sarima_configs(seasonal=[0]):
	models = list()
	# define config lists
	p_params = [0, 1, 2]
	d_params = [0, 1]
	q_params = [0, 1, 2]
	t_params = ['n','c','t','ct']
	P_params = [0, 1, 2]
	D_params = [0, 1]
	Q_params = [0, 1, 2]
	m_params = seasonal
	# create config instances
	for p in p_params:
		for d in d_params:
			for q in q_params:
				for t in t_params:
					for P in P_params:
						for D in D_params:
							for Q in Q_params:
								for m in m_params:
									cfg = [(p,d,q), (P,D,Q,m), t]
									models.append(cfg)
	return models


# parse dates
def custom_parser(x):
	return datetime.strptime('195'+x, '%Y-%m')

if __name__ == '__main__':
	# load dataset
	series = read_csv('shampoo.csv', header=0, index_col=0, date_parser=custom_parser)
	data = series.values
	print(data.shape)
	# data split
	n_test = 12
	# model configs
	cfg_list = sarima_configs()
	# grid search
	scores = grid_search(data, cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

La ejecución del ejemplo puede tardar unos minutos en hardware moderno.

Las configuraciones del modelo y el RMSE se imprimen a medida que se evalúan los modelos. Las tres configuraciones principales del modelo y su error se informan al final de la ejecución.

Nota: Sus resultados pueden variar dada la naturaleza estocástica del algoritmo o procedimiento de evaluación, o diferencias en la precisión numérica. Considere ejecutar el ejemplo varias veces y comparar el resultado promedio.

Podemos ver que el mejor resultado fue un RMSE de unas 54,76 ventas con la siguiente configuración:

• Orden de tendencia: (0, 1, 2)
• Orden estacional: (2, 0, 2, 0)
• Parámetro de tendencia: 't' (tendencia lineal)

...
> Model[[(2, 1, 2), (1, 0, 1, 0), 'ct']] 68.891
> Model[[(2, 1, 2), (2, 0, 0, 0), 'ct']] 75.406
> Model[[(2, 1, 2), (1, 0, 2, 0), 'ct']] 80.908
> Model[[(2, 1, 2), (2, 0, 1, 0), 'ct']] 78.734
> Model[[(2, 1, 2), (2, 0, 2, 0), 'ct']] 82.958
done
[(0, 1, 2), (2, 0, 2, 0), 't'] 54.767582003072874
[(0, 1, 1), (2, 0, 2, 0), 'ct'] 58.69987083057107
[(1, 1, 2), (0, 0, 1, 0), 't'] 58.709089340600094

Estudio de caso 3: estacionalidad

El conjunto de datos de "temperaturas medias mensuales" resume las temperaturas medias mensuales del aire en el castillo de Nottingham, Inglaterra, de 1920 a 1939 en grados Fahrenheit.

El conjunto de datos tiene un componente estacional obvio y ninguna tendencia obvia.

Descargue el conjunto de datos directamente desde aquí:

• temperatura-media-mensual.csv

Guarde el archivo con el nombre de archivo 'monthly-mean-temp.csv' en su directorio de trabajo actual.

Podemos cargar este conjunto de datos como una serie de Pandas usando la función read_csv().

series = read_csv('monthly-mean-temp.csv', header=0, index_col=0)

El conjunto de datos tiene 20 años o 240 observaciones. Recortaremos el conjunto de datos a los últimos cinco años (60 observaciones) para acelerar el proceso de evaluación del modelo y usaremos el último año, o 12 observaciones, para el conjunto de prueba.

# trim dataset to 5 years
data = data[-(5*12):]

El período del componente estacional es de aproximadamente un año, o 12 observaciones. Usaremos esto como el período estacional en la llamada a la función sarima_configs() al preparar las configuraciones del modelo.

# model configs
cfg_list = sarima_configs(seasonal=[0, 12])

A continuación se enumera el ejemplo completo de cuadrícula que busca el problema de pronóstico de series temporales de temperatura media mensual.

# grid search sarima hyperparameters for monthly mean temp dataset
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error
from pandas import read_csv

# one-step sarima forecast
def sarima_forecast(history, config):
	order, sorder, trend = config
	# define model
	model = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)
	# fit model
	model_fit = model.fit(disp=False)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = sarima_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of sarima configs to try
def sarima_configs(seasonal=[0]):
	models = list()
	# define config lists
	p_params = [0, 1, 2]
	d_params = [0, 1]
	q_params = [0, 1, 2]
	t_params = ['n','c','t','ct']
	P_params = [0, 1, 2]
	D_params = [0, 1]
	Q_params = [0, 1, 2]
	m_params = seasonal
	# create config instances
	for p in p_params:
		for d in d_params:
			for q in q_params:
				for t in t_params:
					for P in P_params:
						for D in D_params:
							for Q in Q_params:
								for m in m_params:
									cfg = [(p,d,q), (P,D,Q,m), t]
									models.append(cfg)
	return models

if __name__ == '__main__':
	# load dataset
	series = read_csv('monthly-mean-temp.csv', header=0, index_col=0)
	data = series.values
	# trim dataset to 5 years
	data = data[-(5*12):]
	# data split
	n_test = 12
	# model configs
	cfg_list = sarima_configs(seasonal=[0, 12])
	# grid search
	scores = grid_search(data, cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

La ejecución del ejemplo puede tardar unos minutos en hardware moderno.

Las configuraciones del modelo y el RMSE se imprimen a medida que se evalúan los modelos. Las tres configuraciones principales del modelo y su error se informan al final de la ejecución.

Nota: Sus resultados pueden variar dada la naturaleza estocástica del algoritmo o procedimiento de evaluación, o diferencias en la precisión numérica. Considere ejecutar el ejemplo varias veces y comparar el resultado promedio.

Podemos ver que el mejor resultado fue un RMSE de aproximadamente 1,5 grados con la siguiente configuración:

• Orden de tendencia: (0, 0, 0)
• Orden estacional: (1, 0, 1, 12)
• Parámetro de tendencia: 'n' (sin tendencia)

Como era de esperar, el modelo no tiene un componente de tendencia y sí un componente ARMA estacional de 12 meses.

...
> Model[[(2, 1, 2), (2, 1, 0, 12), 't']] 4.599
> Model[[(2, 1, 2), (1, 1, 0, 12), 'ct']] 2.477
> Model[[(2, 1, 2), (2, 0, 0, 12), 'ct']] 2.548
> Model[[(2, 1, 2), (2, 0, 1, 12), 'ct']] 2.893
> Model[[(2, 1, 2), (2, 1, 0, 12), 'ct']] 5.404
done

[(0, 0, 0), (1, 0, 1, 12), 'n'] 1.5577613610905712
[(0, 0, 0), (1, 1, 0, 12), 'n'] 1.6469530713847962
[(0, 0, 0), (2, 0, 0, 12), 'n'] 1.7314448163607488

Estudio de caso 4: Tendencia y estacionalidad

El conjunto de datos de "ventas mensuales de automóviles" resume las ventas mensuales de automóviles en Quebec, Canadá, entre 1960 y 1968.

El conjunto de datos tiene una tendencia obvia y un componente estacional.

Descargue el conjunto de datos directamente desde aquí:

• ventas-mensuales-de-autos.csv

Guarde el archivo con el nombre 'monthly-car-sales.csv' en su directorio de trabajo actual.

Podemos cargar este conjunto de datos como una serie de Pandas usando la función read_csv().

series = read_csv('monthly-car-sales.csv', header=0, index_col=0)

El conjunto de datos tiene 9 años o 108 observaciones. Usaremos el último año, o 12 observaciones, como conjunto de prueba.

El período del componente estacional podría ser de seis meses o de 12 meses. Probaremos ambos como período estacional en la llamada a la función sarima_configs() al preparar las configuraciones del modelo.

# model configs
cfg_list = sarima_configs(seasonal=[0,6,12])

A continuación se muestra el ejemplo completo de cuadrícula que busca el problema de pronóstico de series temporales de ventas mensuales de automóviles.

# grid search sarima hyperparameters for monthly car sales dataset
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error
from pandas import read_csv

# one-step sarima forecast
def sarima_forecast(history, config):
	order, sorder, trend = config
	# define model
	model = SARIMAX(history, order=order, seasonal_order=sorder, trend=trend, enforce_stationarity=False, enforce_invertibility=False)
	# fit model
	model_fit = model.fit(disp=False)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = sarima_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of sarima configs to try
def sarima_configs(seasonal=[0]):
	models = list()
	# define config lists
	p_params = [0, 1, 2]
	d_params = [0, 1]
	q_params = [0, 1, 2]
	t_params = ['n','c','t','ct']
	P_params = [0, 1, 2]
	D_params = [0, 1]
	Q_params = [0, 1, 2]
	m_params = seasonal
	# create config instances
	for p in p_params:
		for d in d_params:
			for q in q_params:
				for t in t_params:
					for P in P_params:
						for D in D_params:
							for Q in Q_params:
								for m in m_params:
									cfg = [(p,d,q), (P,D,Q,m), t]
									models.append(cfg)
	return models

if __name__ == '__main__':
	# load dataset
	series = read_csv('monthly-car-sales.csv', header=0, index_col=0)
	data = series.values
	print(data.shape)
	# data split
	n_test = 12
	# model configs
	cfg_list = sarima_configs(seasonal=[0,6,12])
	# grid search
	scores = grid_search(data, cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

La ejecución del ejemplo puede tardar unos minutos en hardware moderno.

Las configuraciones del modelo y el RMSE se imprimen a medida que se evalúan los modelos. Las tres configuraciones principales del modelo y su error se informan al final de la ejecución.

Nota: Sus resultados pueden variar dada la naturaleza estocástica del algoritmo o procedimiento de evaluación, o diferencias en la precisión numérica. Considere ejecutar el ejemplo varias veces y comparar el resultado promedio.

Podemos ver que el mejor resultado fue un RMSE de unas 1.551 ventas con la siguiente configuración:

• Orden de tendencia: (0, 0, 0)
• Orden estacional: (1, 1, 0, 12)
• Parámetro de tendencia: 't' (tendencia lineal)

> Model[[(2, 1, 2), (2, 1, 1, 6), 'ct']] 2246.248
> Model[[(2, 1, 2), (2, 0, 2, 12), 'ct']] 10710.462
> Model[[(2, 1, 2), (2, 1, 2, 6), 'ct']] 2183.568
> Model[[(2, 1, 2), (2, 1, 0, 12), 'ct']] 2105.800
> Model[[(2, 1, 2), (2, 1, 1, 12), 'ct']] 2330.361
> Model[[(2, 1, 2), (2, 1, 2, 12), 'ct']] 31580326686.803
done
[(0, 0, 0), (1, 1, 0, 12), 't'] 1551.8423920342414
[(0, 0, 0), (2, 1, 1, 12), 'c'] 1557.334614575545
[(0, 0, 0), (1, 1, 0, 12), 'c'] 1559.3276311282675

Extensiones

Esta sección enumera algunas ideas para ampliar el tutorial que quizás desee explorar.

• Transformaciones de datos. Actualice el marco para admitir transformaciones de datos configurables, como la normalización y la estandarización.
• Previsión de la trama. Actualice el marco para reajustar un modelo con la mejor configuración y pronosticar todo el conjunto de datos de prueba, luego trace el pronóstico en comparación con las observaciones reales en el conjunto de prueba.
• Sintonizar cantidad de historial. Actualice el marco para ajustar la cantidad de datos históricos utilizados para ajustarse al modelo (por ejemplo, en el caso de los 10 años de datos de temperatura máxima).

Si explora alguna de estas extensiones, me encantaría saberlo.

Lectura adicional

Esta sección proporciona más recursos sobre el tema si desea profundizar más.

Publicaciones

• Cómo crear un modelo ARIMA para pronóstico de series temporales con Python
• Cómo realizar una búsqueda en cuadrícula de hiperparámetros del modelo ARIMA con Python
• Una suave introducción a la autocorrelación y la autocorrelación parcial

Libros

• Capítulo 8 Modelos ARIMA, Previsión: principios y práctica, 2013.
• Capítulo 7, Modelos no estacionarios, Series temporales introductorias con R, 2009.

API

• Una suave introducción a SARIMA para el pronóstico de series temporales en Python
• Análisis de series temporales de modelos estadísticos mediante métodos de espacio de estados
• statsmodels.tsa.statespace.sarimax.SARIMAX API
• statsmodels.tsa.statespace.sarimax.SARIMAXResultados API
• Cuaderno Statsmodels SARIMAX
• Joblib: ejecutar funciones de Python como trabajos de canalización

Artículos

• Media móvil integrada autorregresiva en Wikipedia

Resumen

En este tutorial, descubrió cómo desarrollar un marco para la búsqueda en cuadrícula de todos los hiperparámetros del modelo SARIMA para el pronóstico de series temporales univariadas.

Específicamente, aprendiste:

• Cómo desarrollar un marco para la búsqueda de cuadrículas en modelos SARIMA desde cero mediante validación directa.
• Cómo realizar una búsqueda en cuadrícula de hiperparámetros del modelo SARIMA para datos de series temporales diarias de nacimientos.
• Cómo realizar una búsqueda en cuadrícula de hiperparámetros del modelo SARIMA para datos de series temporales mensuales para ventas de champú, ventas de automóviles y temperatura.

¿Tiene alguna pregunta?
Haga sus preguntas en los comentarios a continuación y haré todo lo posible para responder.