Estrategias de Forecasting
Lic. Miguel Angel
2024-03-09
Hablaremos sobre el tercer componente del análisis de series de tiempo: la predicción. Se presenta una introducción a los principales elementos del flujo de trabajo de la predicción, incluyendo enfoques de entrenamiento, evaluación del rendimiento y métodos de referencia. El capítulo cubre temas como el entrenamiento y prueba de modelos de predicción, la evaluación del rendimiento, la medición de errores, los métodos de referencia y la cuantificación de la incertidumbre de la predicción.
El proceso de pronóstico de series de tiempo es un conjunto de pasos que se siguen para hacer una predicción de una serie de tiempo en el futuro. Este proceso se divide en cuatro etapas principales: preparación de datos, entrenamiento del modelo, prueba del modelo y evaluación del modelo.
En la primera etapa, se preparan los datos para el entrenamiento y prueba del modelo. Esto implica dividir la serie en particiones de entrenamiento y prueba, crear nuevas características y aplicar transformaciones según sea necesario.
En la segunda etapa, se entrena el modelo utilizando la partición de entrenamiento. El objetivo principal de esta etapa es ajustar el modelo a los datos de entrenamiento y estimar los parámetros del modelo que minimizan los criterios de error seleccionados.
En la tercera etapa, se prueba el modelo utilizando la partición de prueba. Aquí, se utilizan las estimaciones del modelo para predecir los valores futuros de la serie de tiempo.
En la cuarta y última etapa, se evalúa el rendimiento del modelo. Se comparan las predicciones del modelo con los valores reales y se calculan los errores. Si el modelo cumple con ciertos criterios de desempeño, se utiliza para generar la predicción final. De lo contrario, se selecciona un nuevo modelo o se ajusta el parámetro de entrenamiento y se repite el proceso de entrenamiento.
El proceso de entrenamiento del modelo es fundamental en el proceso de pronóstico, ya que su calidad impacta directamente en la precisión del resultado. Durante este proceso, se busca formalizar la relación de la serie con otros factores, ajustar los parámetros del modelo y asegurar que sea escalable a nuevos datos y evite el sobreajuste. La serie se divide en particiones de entrenamiento y prueba en orden cronológico, para establecer una relación matemática entre la serie y sus valores pasados y errores.
Partición de entrenamiento y prueba
El enfoque de entrenamiento más común es el uso de particiones de entrenamiento y prueba simples. Este enfoque consiste en dividir la serie en particiones de entrenamiento y prueba, entrenar el modelo en la partición de entrenamiento y probar su rendimiento en el conjunto de prueba. El enfoque utiliza un parámetro único: la longitud de la muestra de prueba, que se deriva de reglas empíricas como que debe ser hasta el 30% de la longitud total de la serie. Además, la longitud del horizonte de pronóstico puede ser considerada como una guía adicional para determinar la longitud de la muestra de prueba. Un ejemplo práctico sería si tuviéramos una serie mensual de 72 observaciones y quisiéramos pronosticar el próximo año, utilizaríamos las primeras 60 observaciones para el entrenamiento y las últimas 12 para la prueba. La división de particiones puede ser realizada manualmente en R con la función window del paquete stats.
Cargamos las serie de USgas de los paquetes de TSstudio.
Esta base contiene datos sobre el consumo de gas natural en Estados Unidos. Representa el volumen de gas consumido en miles de millones de pies cúbicos. En esta serie se podría esperar ver patrones estacionales dados por los cambios en el uso del gas natural durante diferentes temporadas del año, como un mayor consumo en invierno debido a la calefacción. También podría haber una tendencia subyacente a largo plazo, que podría estar influenciada por factores como el crecimiento de la población, los cambios en la industria y la economía, o la introducción de tecnologías de eficiencia energética, etc.
library (TSstudio)
data (USgas)
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooSe puede observar las características principales de USgas con la función ts_info(USgas)
ts_info(USgas)## The USgas series is a ts object with 1 variable and 238 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2019 10Se usa la función window para dividir las series en particiones de entrenamiento y testeo
train <- window(USgas,
start = time (USgas)[1],
end = time(USgas)[length(USgas) - 12])
test <- window (USgas,
start = time(USgas)[length(USgas) - 12 + 1],
end = time(USgas)[length(USgas)])ts_info(train)## The train series is a ts object with 1 variable and 226 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2018 10ts_info(test)## The test series is a ts object with 1 variable and 12 observations
## Frequency: 12
## Start time: 2018 11
## End time: 2019 10Luego usamos la función ts_split de TSstudio para crear una forma de las particiones training y testing para datos de series de tiempo.
ts_split## function (ts.obj, sample.out = NULL)
## {
## if (!stats::is.ts(ts.obj) && !tsibble::is_tsibble(ts.obj)) {
## stop("The 'ts.obj' is not a valid 'ts' or 'tsibble' object")
## }
## l <- train <- test <- split <- NULL
## if (stats::is.ts(ts.obj)) {
## l <- base::length(ts.obj)
## }
## else if (tsibble::is_tsibble(ts.obj)) {
## l <- base::nrow(ts.obj)
## }
## if (base::is.null(sample.out)) {
## h <- base::round(l * 0.3)
## }
## else if (base::round(sample.out) != sample.out) {
## stop("The 'sample.out' parameter is not a valid number (must be an integer)")
## }
## else if (sample.out >= l) {
## warning("The length of the sample out period is", " longer than the length of the series, ",
## "using the default option (30% of the length of the series)")
## h <- base::round(l * 0.3)
## }
## else {
## h <- sample.out
## }
## if (stats::is.ts(ts.obj)) {
## split <- base::list(train <- stats::window(ts.obj, start = stats::time(ts.obj)[1],
## end = stats::time(ts.obj)[base::length(stats::time(ts.obj)) -
## h]), test <- stats::window(ts.obj, start = stats::time(ts.obj)[base::length(stats::time(ts.obj)) -
## h + 1], end = stats::time(ts.obj)[base::length(stats::time(ts.obj))]))
## }
## else if (tsibble::is_tsibble(ts.obj)) {
## split <- base::list(train = ts.obj[1:(base::nrow(ts.obj) -
## h), ], test = ts.obj[(base::nrow(ts.obj) - h + 1):base::nrow(ts.obj),
## ])
## }
## base::names(split) <- c("train", "test")
## return(split)
## }
## <bytecode: 0x000001eb1dc20fa0>
## <environment: namespace:TSstudio>Esta función se utiliza para dividir un objeto de series de tiempo, que puede ser de clase “ts” o “tsibble”, en particiones de entrenamiento y prueba. Se puede especificar el número de periodos que se quiere incluir en la partición de prueba mediante el argumento sample.out; si no se especifica, el valor predeterminado corresponderá al 30% de la longitud de la serie. Por ejemplo, si se desea establecer los últimos 12 meses de datos como partición de prueba y el resto como entrenamiento, puede hacerse indicando sample.out = 12
USgas_partitions <- ts_split(USgas, sample.out = 12)
train <- USgas_partitions$train
test <- USgas_partitions$testts_info(train)## The train series is a ts object with 1 variable and 226 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2018 10ts_info(test)## The test series is a ts object with 1 variable and 12 observations
## Frequency: 12
## Start time: 2018 11
## End time: 2019 10El método es simple, rápido y económico para entrenar y probar modelos, pero no es posible evaluar la estabilidad y escalabilidad del rendimiento del modelo con una única prueba. El enfoque de backtesting ayuda a mitigar este riesgo entrenando un modelo con múltiples particiones de entrenamiento y prueba.
Forecasting con backtesting
El enfoque de backtesting para entrenar un modelo de pronóstico es una versión avanzada del enfoque único fuera de muestra. Se basa en el uso de una ventana rodante para dividir la serie en múltiples pares de particiones de entrenamiento y prueba. El proceso básico de entrenamiento de backtesting incluye la preparación de datos, entrenamiento del modelo, prueba del modelo y evaluación del modelo. La metodología de puntuación nos permite evaluar la estabilidad del modelo examinando la tasa de error del modelo en los diferentes conjuntos de prueba. Este método es similar al enfoque de validación cruzada para entrenar modelos de aprendizaje automático, pero las particiones son ordenadas cronológicamente en lugar de ser basadas en un muestreo aleatorio. La elección de la ventana depende de las características de la serie de entrada y el equilibrio entre el costo computacional y la información obtenida de la evaluación del modelo.
Forecast evaluación
El objetivo principal de la etapa de evaluación es evaluar (valga la redundancia) la capacidad del modelo entrenado para predecir (o basado en otros criterios) las observaciones futuras de la serie con precisión. Este proceso incluye lo siguiente:
- Análisis residual
- Puntuaciones los pronosticos
Análisis Residual
El análisis residual evalúa cómo de bien el modelo ha capturado los patrones de la serie y proporciona información sobre las distribuciones de residuos necesarias para construir intervalos de confianza para la predicción. Se utiliza herramientas de visualización de datos y pruebas estadísticas para evaluar la bondad del ajuste del modelo frente a los valores reales, la autocorrelación de los residuos y la distribución de los mismos. Los residuos con oscilación aleatoria alrededor de cero y con variación constante indican que el modelo es capaz de capturar la mayoría de la variación de la serie, mientras que una oscilación residual que no tiene las características del ruido blanco indica que el modelo no logró capturar los patrones de la serie. Si los residuos no están distribuidos normalmente, no se pueden utilizar para crear intervalos de confianza.
Para demostrar esto, se puede entrenar un modelo ARIMA en la partición de entrenamiento que se creó antes para la serie de USgas.
Usamos la función auto.arima del paquete forecas:
library(forecast)
md <- auto.arima(train)
checkresiduals(md)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,1)(2,1,1)[12]
## Q* = 24.95, df = 18, p-value = 0.1263
##
## Model df: 6. Total lags used: 24Los resultados indican que el modelo no captura completamente todos los patrones de la serie y que puede haber valores atípicos durante algunos períodos. Se sugiere que se modifiquen los parámetros del modelo para mejorar su ajuste. Además, se señala que es difícil obtener residuos completamente independientes y se recomienda seleccionar un modelo que se acerque a ese objetivo.
Puntuaciones de forecast
Ésta parte se centra en la importancia de evaluar el rendimiento de los modelos de aprendizaje automático. Una vez que se completa la sintonización del modelo, es necesario evaluar su capacidad para hacer predicciones precisas en observaciones que nunca ha visto antes. La evaluación del rendimiento se lleva a cabo utilizando métricas de precisión o error, como el Error Cuadrático Medio (MSE), la Raíz del Error Cuadrático Medio (RMSE), el Error Absoluto Medio (MAE) y el Error Porcentual Absoluto Medio (MAPE).
El MSE es una medida de la distancia media al cuadrado entre los valores reales y pronosticados, lo que evita que los valores positivos y negativos se cancelen entre sí y penaliza la puntuación de error a medida que aumenta la tasa de error. El RMSE toma la raíz cuadrada de MSE y también es sensible a valores atípicos debido al efecto cuadrático del error.
El MAE mide la tasa de error absoluto del pronóstico y solo puede tener valores positivos, lo que la hace menos sensible a los valores atípicos. Por otro lado, el MAPE mide el error porcentual absoluto promedio, lo que facilita la comparación y la comunicación de resultados.
La elección de una métrica de error adecuada depende de los objetivos de precisión de la predicción y las limitaciones del modelo. La evaluación del modelo es un proceso iterativo y continuo que requiere ajustes y mejoras constantes para aumentar su precisión.
Una forma de ejemplificar esto es usando el modelo que se entrenó antes. Se usa la función forecast del paquete homologo.
fc <- forecast(md, h = 12)Luego se usa la función de accuracy para revisar el rendimiento del modelo.
accuracy(fc, test)## ME RMSE MAE MPE MAPE MASE
## Training set 5.844136 97.81626 73.42657 0.1170672 3.522348 0.6376860
## Test set 37.847885 103.22848 81.46603 1.3107987 3.261643 0.7075062
## ACF1 Theil's U
## Training set -0.004183172 NA
## Test set -0.046708926 0.3404092La función de accuracy será muy utilizada, ya que devuelve métricas de error para los valores ajustados y la previsión real en el conjunto de entrenamiento y prueba. Es importante destacar que el MAPE del modelo es del 3,52% en el conjunto de entrenamiento y del 3,26% en el conjunto de prueba. Sería normal que la tasa de error sea más alta en la partición de prueba ya que el modelo solo vio los datos de la partición de entrenamiento durante el proceso de entrenamiento. Si la tasa de error es baja en el conjunto de entrenamiento pero alta en el conjunto de prueba, esto es una señal de sobreajuste en el modelo.
Para evaluar el ajuste del modelo en el conjunto de entrenamiento y prueba, se puede usar la función test_forecast de TSstudio, que muestra la serie real, los valores ajustados en el conjunto de entrenamiento y los valores pronosticados en el conjunto de prueba. Al pasar el cursor sobre los valores ajustados o pronosticados, aparecen los resultados de RMSE y MAPE tanto en las particiones de entrenamiento como de prueba
test_forecast(actual = USgas,
forecast.obj = fc,
test = test)Training Set: MAPE (3.52%): Error Porcentual Absoluto Medio en el conjunto de entrenamiento, indicando que las predicciones se desvían en promedio un 3.52% de los valores reales.
RMSE (97.82): Raíz del Error Cuadrático Medio en el conjunto de entrenamiento, una medida de la magnitud de los errores del modelo.
Testing Set: MAPE (3.26%): Error Porcentual Absoluto Medio en el conjunto de prueba, sugiere que las predicciones son ligeramente más precisas que en el conjunto de entrenamiento.
RMSE (103.23): RMSE en el conjunto de prueba, ligeramente superior al del conjunto de entrenamiento, lo que podría indicar una disminución en el rendimiento del modelo cuando se aplica a nuevos datos.
Es más sencillo y rápido identificar conclusiones acerca de la calidad del ajuste tanto de los valores ajustados como pronosticados cuando se grafican dichos valores contra los valores reales de la serie. Al hacerlo, se puede notar de inmediato que el pico residual durante el 2006 se debe a valores atípicos (o consumo inferior al patrón normal de la serie).
Es más sencillo y rápido identificar conclusiones acerca de la calidad del ajuste tanto de los valores ajustados como pronosticados cuando se grafican dichos valores contra los valores reales de la serie. Al hacerlo, se puede notar de inmediato que el pico residual durante el 2006 se debe a valores atípicos (o consumo inferior al patrón normal de la serie).
Punto de referencia de pronóstico
De acuerdo con las métricas de error, el modelo entrenado logró una puntuación de MAPE del 3,52% o un RMSE de 97,82. Para determinar si estos resultados son satisfactorios, es común comparar el rendimiento del modelo con alguna previsión de referencia o algún método ya existente que deseamos mejorar. Una forma popular de establecer una referencia es utilizar un enfoque de previsión simplista como punto de partida. Por ejemplo, podríamos pronosticar la serie con un enfoque naive y compararlo con el pronóstico anterior creado con el modelo ARIMA.
El enfoque naive asume que el valor más reciente es el representante verdadero del futuro, por lo que se utiliza este valor como base para el pronóstico futuro. Podemos crear un pronóstico ingenuo utilizando la función naive del paquete forecast, y utilizar el conjunto de entrenamiento como entrada para el modelo. De esta manera, podemos determinar si el rendimiento del modelo es aceptable en comparación con la referencia establecida:
Primero llamamos la libreria y paquete de Forecast, luego usamamos la función test_forecast para testear y entrenar la función.
library(forecast)
naive_model <- naive(train, h = 12)
test_forecast(actual = USgas,
forecast.obj = naive_model,
test = test)También podemos usar accuracy para evaluar el rendimiento del modelo:
accuracy(naive_model, test)## ME RMSE MAE MPE MAPE MASE
## Training set -1.028444 285.6607 228.5084 -0.9218463 10.97123 1.984522
## Test set 301.891667 499.6914 379.1417 9.6798015 13.28187 3.292723
## ACF1 Theil's U
## Training set 0.3761105 NA
## Test set 0.7002486 1.499679En ésta parte se resalta el uso de modelos naive y modelos naive estacionales en la predicción de datos de series de tiempo, específicamente en el caso de USgas que tiene un patrón estacional fuerte. El modelo naive no requiere entrenamiento y establece los valores ajustados como los valores reales. En cambio, el modelo naive estacional toma en cuenta la variación estacional y utiliza el último punto estacional como pronóstico para todas las observaciones estacionales correspondientes en el futuro.
Por ejemplo, se usa una serie de meses, para evaluar el enero más resiente en la serie, con respecto a USgas.
snaive_model <- snaive(train, h = 12)
test_forecast(actual =USgas,
forecast.obj = snaive_model,
test = test)Volvemos a usar accuracy para revisar el rendimiento de naive estacional en el entrenamiento y testeo.
accuracy(snaive_model, test)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 33.99953 148.7049 115.1453 1.379869 5.494048 1.000000 0.4859501
## Test set 96.45000 164.6967 135.8833 3.612060 5.220458 1.180103 -0.2120929
## Theil's U
## Training set NA
## Test set 0.4289964Esto puede mostrar que el modelo estacional de naive, tiene un mejor ajuste al tipo de series de tiempo que se están pronosticando con la serie de USgas, gracias al fuerte patrón estacional. Se usa como un punto de referencia para el modelo ARIMA, comparando el MAPE y RMSE, de ambos modelos en el testeo. Dejando un claro resultado de que ARIMA da un ascenso con respecto al modelo del punto de referencia
Finalizando the forecast
El último paso consiste en recalibrar los pesos o coeficientes del modelo con la serie completa. Para hacerlo, existen dos enfoques, dependiendo de si el modelo fue ajustado manualmente o automáticamente por un algoritmo. Si se ajustó manualmente, se deben utilizar los mismos parámetros de ajuste que se usaron en el modelo entrenado. Si se ajustó automáticamente, se puede extraer la configuración de parámetros usada con la partición de entrenamiento o permitir que el algoritmo ajuste los parámetros del modelo de nuevo usando la serie completa, asumiendo que tiene la capacidad de hacerlo correctamente al entrenar con nuevos datos.
Para simplificar todo, se seguirá usando el modelo auto.arima para entrenar el modelo final:
md_final <- auto.arima (USgas)
fc_final <- forecast(md_final, h = 12)Se usa la función plot_forecast de TSstudio para trazar el pronostico final:
plot_forecast(fc_final,
title = "The US Natural Gas Consumption Forecast",
Xtitle = "year",
Ytitle = "Billion Cubic Feet")Manejo de la incertidumbre del pronóstico
El propocito principal de los procesos de pronostico, es minimizar el nivel de incertidumbre sobre valores de series futuras. A pesar de que no se puede eliminar totalmente esta incertidumbre, podemos cuantificar y dar un rango para estimar el pronostico. Se puede realizar usando intervalos de confianza, o usando un simulador.
Condifence interval
Son un metodo de aproximación estadistica, que es usada para expresar el rango de prosibles valores, que contienen el valor real en algún punto de confianza.
La función forecast genera un intervalo de predicción con un nivel de confianza del 80% al 95%, pero se puede modificar.
Por ejemplo, se usan las funciones del modelo entrenado md_final y forecast, para los siguientes 60 meses, usando un intervalo de confianza del 80% al 90%.
fc_final2 <- forecast(md_final,
h = 60,
level = c(80, 90))
plot_forecast(fc_final2,
title = "The US Natural Gas Counsumption Forecast",
Xtitle = "Year",
Ytitle = "Billion Cubic Feet")Simulacion
La función forecast_sim del paquete TSstudio da una función de construción, para que la simulación sea posible en campos de pronosticos. Este estimado se puede usar para calcular pronosticos de puntos estimados, o calcular probabilidades de tener diferentes valores.
Usaremos el mismo modelo para la función y haremos 100 repeticiones:
fc_final3 <- forecast_sim(model = md_final,
h = 60,
n = 100)library(plotly)
fc_final3$plot %>%
layout(title = "US Natural Gas Consumption - Forecasting Simulation",
yaxis = list(title = "Billion Cubic Feet"),
xaxis = list(title = "Year"))Enfoque de competencia (Horse race approach)
El acercamiento Horse race es basado en el entrenamiento, testeo y evaluación de multiples modelos de pronostico. En este ejemplo se aplica el acercamiento horse race para 7 modelos diferentes
Los modelos que se evaluarán son:
-auto.arima -bsts -ets -hybrid -nnetar -tbats -HoltWinters
Usamos la función set.seed, para que se puedan arrojar los resultados:
set.seed(1234)#Realizamos una competencia entre varios modelos
methods <- list(ets1 = list(method = "ets",
method_arg = list(opt.crit = "lik"),
notes = "ETS model with opt.crit = lik"),
ets2 = list(method = "ets",
method_arg = list(opt.crit = "amse"),
notes = "ETS model with opt.crit = amse"),
arima1 = list(method = "arima",
method_arg = list(order = c(2,1,0)),
notes = "ARIMA(2,1,0)"),
arima2 = list(method = "arima",
method_arg = list(order = c(2,1,2),
seasonal = list(order = c(1,1,1))),
notes = "SARIMA(2,1,2)(1,1,1)"),
hw = list(method = "HoltWinters",
method_arg = NULL,
notes = "HoltWinters Model"),
tslm = list(method = "tslm",
method_arg = list(formula = input ~ trend + season),
notes = "tslm model with trend and seasonal components"))
#Entrenando los modelos con validación retrospectiva (backtesting)
md <- train_model(input = USgas,
methods = methods,
train_method = list(partitions = 6,
sample.out = 12,
space = 3),
horizon = 12,
error = "MAPE")## # A tibble: 6 × 7
## model_id model notes avg_mape avg_rmse `avg_coverage_80%` `avg_coverage_95%`
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 hw HoltWi… Holt… 0.0482 144. 0.792 0.931
## 2 ets1 ets ETS … 0.0526 156. 0.833 0.972
## 3 arima2 arima SARI… 0.0546 163. 0.583 0.819
## 4 ets2 ets ETS … 0.0650 185. 0.5 0.792
## 5 tslm tslm tslm… 0.0854 242. 0.319 0.611
## 6 arima1 arima ARIM… 0.163 539. 0.861 0.958Esta tabla muestra los resultados de un análisis comparativo de diferentes modelos de pronóstico para series de tiempo. Se presenta la siguiente información para cada modelo:
model_id: Identificador numérico para cada modelo.
model: Nombre abreviado del modelo. notes: Notas adicionales o la especificación completa del modelo.
avg_mape: Error porcentual absoluto medio, una medida de la precisión del modelo, donde valores más bajos indican mayor precisión.
avg_rmse: Raíz del error cuadrático medio, que mide la diferencia entre los valores que el modelo de pronóstico predice y los valores reales.
avg_coverage_80%: El porcentaje promedio de veces que el intervalo de confianza al 80% incluyó el valor real.
avg_coverage_95%: El porcentaje promedio de veces que el intervalo de confianza al 95% incluyó el valor real.
Por ejemplo, el modelo HoltWinters (hw) tiene un MAPE promedio del 4.82%, lo que indica que sus predicciones, en promedio, se desvían un 4.82% de los valores reales, mientras que su cobertura promedio del intervalo de confianza al 95% es del 93.1%, lo que indica que el 93.1% de los valores reales cayeron dentro de su intervalo de confianza pronosticado al 95%.
El modelo earroja una lista ordenada en el criterio de error que se da. En este caso, HoltWinters tiene el rate de error más bajo, por lo que es el modelo recomendado para usar.
Podemos extraer una taza de error de la siguiente manera
library(forecast)
plot_error(md)