En los capitulos previos, el enfoque fue fundamentalmente en los dos primeros componentes de flujo de trabajo para análisis de series temporales: preprocesamiento de datos y análisis descriptivo. Partiendo de este capitulo, se tomará una dirección distinta y se abordarán el tercer y último componente de análisis: el pronóstico y, asimismo, los elementos principales del flujo de trabajo de pronóstico. Esto incluye enfoques para entrenar un modelo de pronóstico, evaluación de desempeño y métodos de referencia. Esto es relevante para proporcionar un conjunto de herramientas para el diseño y construcción de modelos de pronóstico acorde al objetivo de análisis.
Este capitulo tratará las siguientes temáticas:
• Enfoques de capacitación.
• Métodos de evaluación de desempeño y matrices de medición de errores.
• Métodos de referencia.
• Cuantificación de incertidumbre del pronóstico con intervalos de confianza y simulación.
Los siguientes paquetes serán utilizados a lo largo de este capitulo.
• Forecast: Versión 8.5 y superior.
• TStudio: Versión 0.1.4 y superior.
• Plotly: Versión 4.8 y superior.
El pronóstico de series de tiempo tradicional sigue el mismo flujo de trabajo que la mayoría de los campos del análisis predictivo, como la regresión o la clasificación, y normalmente incluye los siguientes pasos:
1. Preparación de datos: Aquí preparamos los datos para el proceso de entrenamiento y prueba del modelo. Este paso incluye dividir la serie en particiones de entrenamiento (dentro de la muestra) y de prueba (fuera de la muestra), creando nuevas características (cuando corresponda) y aplicando una transformación si es necesario (por ejemplo, transformación de registros, escalado, etc.).
2. Entrene el modelo: aquí utilizamos la partición de entrenamiento para entrenar un modelo estadístico. El objetivo principal de este paso es utilizar el conjunto de entrenamiento para entrenar, ajustar y estimar los coeficientes del modelo que minimizan los criterios de error seleccionados (más adelante en este capítulo, analizaremos en detalle las métricas de error comunes). Los valores ajustados y la estimación del modelo de las observaciones de la partición de entrenamiento se utilizarán más adelante para evaluar el rendimiento general del modelo.
3. Pruebe el modelo: aquí utilizamos el modelo entrenado para pronosticar las observaciones correspondientes de la partición de prueba. La idea principal aquí es evaluar el rendimiento del modelo con un nuevo conjunto de datos (que el modelo no vio durante el proceso de entrenamiento).
4. Evaluación del modelo: por último, pero no menos importante, una vez entrenado y probado el modelo, es hora de evaluar el rendimiento general del modelo tanto en la partición de entrenamiento como en la de prueba.
Según el proceso de evaluación del modelo, si el modelo cumple con un cierto umbral o criterio, conservamos el modelo usando la serie completa para generar el pronóstico final o seleccionamos un nuevo parámetro de entrenamiento/modelo diferente y repetimos el proceso de entrenamiento.
⚠️ CONSEJO ⚠️
Uno de los principales errores al entrenar un modelo estadístico es sobreajustar el modelo al conjunto de entrenamiento. El sobreajuste ocurre cuando el modelo sobreaprende los datos con los que se entrenó y no logra generalizar el rendimiento del modelo en el conjunto de entrenamiento en otros conjuntos de datos. Uno de los principales signos de sobreajuste es una relación alta entre la tasa de error en el conjunto de prueba y el conjunto de entrenamiento (por ejemplo, un error alto en la partición de entrenamiento versus un error bajo en la partición de entrenamiento). Por lo tanto, la evaluación del desempeño del modelo tanto en la partición de entrenamiento como en la de prueba es esencial para identificar el sobreajuste. Generalmente, el sobreajuste se debe a un ajuste incorrecto de los parámetros del modelo o al uso de un modelo de alta complejidad.
Por otro lado, este proceso tiene características propias y singulares, que lo distinguen de otros campos predictivos:
• Las particiones de capacitación y pruebas deben ordenarse en orden cronológico, a diferencia del muestreo aleatorio.
• Normalmente, una vez que hayamos entrenado y probado el modelo usando las particiones de entrenamiento y prueba, volveremos a entrenar el modelo con todos los datos (o al menos la observación más reciente en orden cronológico). A primera vista, esto puede resultar impactante y aterrador para las personas con experiencia en aprendizaje automático tradicional o modelado de aprendizaje supervisado, ya que normalmente conduce a un sobreajuste y otros problemas. Discutiremos el motivo detrás de esto y cómo evitar el sobreajuste más adelante.
El siguiente diagrama demuestra el flujo de trabajo de previsión:
Enfoques de formación Uno de los elementos centrales del flujo de trabajo de previsión es el proceso de formación del modelo. La calidad del entrenamiento del modelo tendrá un impacto directo en el resultado del pronóstico. Los principales objetivos de este proceso son los siguientes:
• Formalizar la relación de la serie con otros factores, como patrones estacionales y de tendencia, correlación con rezagos pasados y variables externas de manera predictiva.
• Ajustar los parámetros del modelo (cuando corresponda).
• El modelo es escalable con datos nuevos o, en otras palabras, evita el sobreajuste.
Como se mencionó previamente, antes del proceso de entrenamiento, la serie se divide en conjuntos de entrenamiento y prueba. En el conjunto de entrenamiento, el modelo se ajusta, y luego se evalúa en el conjunto de prueba. Es fundamental que estas divisiones sigan un orden cronológico, sin importar la técnica de entrenamiento utilizada. Esto se debe a que la mayoría de los modelos de series temporales establecen relaciones matemáticas basadas en observaciones pasadas y errores anteriores de la serie.
Uno de los enfoques de capacitación más comunes es utilizar particiones únicas de capacitación y prueba (o particiones únicas fuera de muestra). Este enfoque simplista se basa en dividir la serie en particiones de entrenamiento y prueba (o particiones dentro y fuera de la muestra, respectivamente), entrenar el modelo en la partición de entrenamiento y probar su rendimiento en el conjunto de prueba:
En el gráfico previo se observa que las divisiones entre los conjuntos de entrenamiento y prueba siguen un patrón cronológico. Este método se caracteriza por un único factor: la extensión del segmento fuera del conjunto de prueba, o dicho de otra manera, la longitud de la partición de prueba. Por lo general, la longitud de esta partición se determina a partir de los parámetros generales que se describen a continuación:
• La longitud de la partición de prueba debe ser hasta el 30% de la longitud total de la serie para tener suficientes datos de observación para el proceso de entrenamiento.
• La duración del horizonte de previsión (siempre que no viole el plazo anterior). Esto se relaciona principalmente con el hecho de que el nivel de incertidumbre aumenta a medida que aumenta el horizonte de pronóstico. Por lo tanto, alinear el conjunto de pruebas con el horizonte de pronóstico podría, potencialmente, proporcionar una estimación más aproximada del error esperado del pronóstico.
Como ejemplo, si tenemos una serie mensual con 72 observaciones (equivalente a 6 años) y deseamos predecir el próximo año (12 meses), tiene sentido utilizar las primeras 60 observaciones para la capacitación y evaluar el rendimiento utilizando las 12 observaciones más recientes como conjunto de prueba. La creación de estas divisiones en R se puede llevar a cabo manualmente utilizando la función “ventana” del paquete de estadísticas. Por ejemplo, podemos dividir la serie USgas en particiones, utilizando las últimas 12 observaciones como conjunto de prueba y el resto como conjunto de entrenamiento.
#Es posible observar las principales características de la serie USgas con la función ts_info y obtendremos el siguiente resultado:
data(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 10# Usemos la función de ventana para dividir la serie en particiones de entrenamiento y prueba
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)])## The train series is a ts object with 1 variable and 226 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2018 10## The test series is a ts object with 1 variable and 12 observations
## Frequency: 12
## Start time: 2018 11
## End time: 2019 10#Alternativamente, la función ts_split del paquete TSstudio proporciona una forma personalizada de entrenamiento y prueba para datos de series temporales
USgas_partitions <- ts_split(USgas, sample.out=12)
train <- USgas_partitions$train
test <- USgas_partitions$test#Se puede observar en el siguiente resultado que recibimos los resultados de ejecución que recibimos anteriormente
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 10## The test series is a ts object with 1 variable and 12 observations
## Frequency: 12
## Start time: 2018 11
## End time: 2019 10La principal ventaja de este método radica en su sencillez, ya que permite entrenar y evaluar rápidamente un modelo con un costo computacional comparativamente bajo. Sin embargo, no es posible determinar la estabilidad y la capacidad de adaptación del modelo basándose únicamente en una única prueba. Existe la posibilidad de que un modelo tenga un desempeño aparentemente bueno en el conjunto de pruebas por pura casualidad, pero luego presente un rendimiento deficiente en la predicción real debido a su falta de estabilidad en el tiempo. Una manera de mitigar este riesgo es mediante la técnica de backtesting, que implica entrenar el modelo con múltiples divisiones de conjuntos de entrenamiento y prueba.
El enfoque de backtesting para entrenar un modelo de pronóstico es una versión avanzada del enfoque único fuera de muestra que vimos anteriormente. Se basa en el uso de una ventana móvil para dividir la serie en múltiples pares de particiones de entrenamiento y prueba. Un proceso básico de capacitación en backtesting incluye los siguientes pasos:
1. Preparación de datos: cree varios pares de particiones de entrenamiento y prueba.
2. Entrenamiento del modelo: Esto se realiza en cada una de las particiones de entrenamiento.
3. Prueba del modelo: califique su rendimiento en las particiones de prueba correspondientes.
4. Evaluación del modelo: Evalúe la precisión, escalabilidad y estabilidad del modelo en función de la puntuación de la prueba. Según la evaluación, haría una de los siguientes aspectos: 1. Generar el pronóstico final para verificar si la puntuación del modelo cumple con un umbral o criterio específico. 2. Aplique ajustes y optimizaciones adicionales para el modelo y repita los pasos de capacitación y evaluaciones.
Usar la técnica de calificación permite evaluar la estabilidad del modelo al observar la tasa de error en diferentes conjuntos de pruebas. Un modelo se considera estable si los errores en las pruebas tienen una distribución estrecha y si la tasa de error en la predicción real se encuentra dentro del rango de las pruebas, a menos que ocurran eventos inusuales.
Este enfoque es similar en concepto a la validación cruzada utilizada en el entrenamiento de modelos de aprendizaje automático, pero se diferencia en que en el enfoque de validación cruzada, las divisiones se realizan de manera aleatoria, mientras que en el backtesting siguen un orden cronológico.
Los principales parámetros de configuración de un modelo de entrenamiento de backtesting son los siguientes:
• La longitud de las particiones de entrenamiento: La establece la configuración de la ventana. Hay dos tipos comunes de ventanas enrollables:
Ventana en expansión: utilizando las primeras N observaciones como partición de entrenamiento inicial, agregue las siguientes n observaciones para crear la siguiente partición de entrenamiento. La longitud de la partición de entrenamiento crece a medida que la ventana pasa por la serie.
Ventana deslizante: establezca las primeras N observaciones como la partición de entrenamiento inicial y cambie la ventana en n observaciones para crear la siguiente partición de entrenamiento. La longitud de las particiones de entrenamiento sigue siendo la misma a medida que la ventana pasa por la serie:
La longitud de las particiones de prueba: un valor constante, generalmente alineado con la longitud del horizonte de pronóstico (bajo la limitación del número mínimo de observaciones necesarias para entrenar el modelo).
El espacio entre cada partición de entrenamiento define el ritmo de la ventana móvil
El número de particiones de entrenamiento y prueba.
El siguiente diagrama demuestra la estructura de backtesting con una ventana de expansión:
En el gráfico anterior, todas las divisiones de entrenamiento del método de ventana en expansión comienzan en el mismo punto, llamado “To”. Cada división de entrenamiento incluye la anterior y las siguientes “n” observaciones (donde “n” es una constante que representa la tasa de expansión de la ventana). Este enfoque es adecuado cuando la serie muestra un patrón estacional marcado y una tendencia constante. En tales casos, las primeras observaciones de la serie pueden contener información útil para el modelo. Sin embargo, un inconveniente es que las particiones tienen longitudes diferentes, y los modelos entrenados con particiones posteriores suelen funcionar mejor que los entrenados con las iniciales.
Este sesgo no existe en el segundo enfoque de entrenamiento: La ventana deslizante, como todas las particiones de entrenamiento, tiene la misma longitud, como se ve en el siguiente diagrama:
Es más apropiado utilizar la ventana deslizante cuando la serie de entrada experimenta cambios estructurales, alta volatilidad o cuando el poder predictivo se relaciona principalmente con los datos más recientes o con rezagos recientes de la serie. Por ejemplo, en el caso de los precios mensuales del petróleo crudo, la información más relevante para predecir su dirección futura se encuentra en los datos recientes, mientras que los datos lejanos tienen menos influencia. Aunque el enfoque de backtesting ofrece información sobre la estabilidad y escalabilidad del modelo, requiere más recursos computacionales en comparación con el enfoque de una sola partición de entrenamiento y prueba. Por lo tanto, se debe considerar este equilibrio al elegir el enfoque de capacitación.
El objetivo principal del paso de evaluación es evaluar la capacidad del modelo entrenado para pronosticar (o basándose en otros criterios) las observaciones futuras de la serie con precisión. Este proceso incluye hacer lo siguiente:
• Análisis residual: se centra en la calidad del modelo, con valores ajustados en la partición de entrenamiento.
• Puntuación del pronóstico: se basa en la capacidad del modelo para pronosticar los valores reales del conjunto de pruebas.
El análisis residual prueba qué tan bien el modelo capturó e identificó los patrones de la serie. Además, proporciona información sobre las distribuciones de residuos, que se requieren para construir intervalos de confianza para el pronóstico. La definición matemática de un residual es la diferencia entre la observación real y el valor ajustado correspondiente del modelo, que se observó en el proceso de entrenamiento, o como la siguiente ecuación:
Este proceso incluye el uso de herramientas de visualización de datos y pruebas estadísticas para evaluar lo siguiente:
■ Pruebe la bondad del ajuste frente a los valores reales: esto se hace trazando los valores residuales a lo largo del tiempo en orden cronológico. El gráfico indica qué tan bien el modelo pudo capturar la oscilación de la serie en la partición de entrenamiento.Los residuos con oscilación aleatoria alrededor del cero y con variación constante (ruido blanco) indican que el modelo es capaz de capturar la mayor parte de la variación de la serie. Por otro lado, si la oscilación residual no tiene las características de ruido blanco, es una indicación de que el modelo no logró capturar los patrones de la serie. A continuación se muestran algunas posibles interpretaciones del posible resultado:
Todos o la mayoría de los residuos están por encima de las líneas cero: esto es una indicación de que el modelo tiende a subestimar los valores reales.
Todos o la mayoría de los residuos están por debajo de las líneas cero: esto es una indicación de que el modelo tiende a sobreestimar los valores reales.
Picos aleatorios: Esto es una indicación de posibles valores atípicos en la partición del entrenamiento.
La autocorrelación residual: indica qué tan bien el modelo fue capaz de capturar los patrones de la serie. Los rezagos no correlacionados indican que no hubo patrones que el modelo no capturó. De manera similar, la existencia de rezagos correlacionados indica patrones que el modelo no captó.
La distribución de residuos: Se requiere para concluir sobre la confiabilidad del intervalo de confianza del pronóstico. Si los residuos no se distribuyen normalmente, no podemos usarlo para crear intervalos deconfianza, ya que se basa en el supuesto de que los residuos se distribuyen normalmente.
■ Para demostrar el proceso de análisis residual, entrenaremos un modelo ARIMA en la partición de entrenamiento que creamos anteriormente para la serie USgas. No se preocupe si no está familiarizado con el modelo ARIMA; lo analizaremos en detalle en el Capítulo 11, Pronósticos con el modelo ARIMA. Entrenaremos el modelo con la función auto.arima del paquete de pronóstico:
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooPara examinar los residuos, usaremos la función de verificación de residuos del paquete de pronóstico, que devuelve los siguientes cuatro resultados:
• Gráfico de series temporales de los residuos.
• Gráfico ACF de los residuos.
• Gráfico de distribución de los residuos.
• El resultado de la prueba Ljung-Box.
La prueba de Ljung-Box es un método estadístico para probar si la autocorrelación de una serie (que, en este caso, son los residuales) es diferente de cero y utiliza la siguiente hipótesis:
• Ho: El nivel de correlación entre la serie y su rezago es igual a cero, por lo que las observaciones de la serie son independientes.
• H₁: El nivel de correlación entre la serie y su rezago es diferente de cero, y por lo tanto las observaciones de la serie no son independientes
#Usemos la función para evaluar el rendimiento del modelo entrenado en la partición.
#Obtendremos los siguientes residuos del modelo ARIMA
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: 24Según los resultados de la prueba de Ljung-Box, se puede rechazar la hipótesis nula con un nivel de significancia del 0,01, lo que sugiere que existe correlación entre la serie residual y sus rezagos. Esto indica que el modelo no capturó todos los patrones de la serie, y es necesario ajustar los parámetros del modelo. La gráfica de la serie residual muestra oscilaciones alrededor del eje x, con algunos valores atípicos que cruzan 250. Esto sugiere la presencia de valores inusuales que deben ser investigados. En secciones posteriores, se explorará un método más intuitivo para comparar los valores ajustados con los reales. Por último, el gráfico de distribución de residuos se asemeja a una distribución normal.
En un escenario ideal, se buscaría obtener ruido blanco y residuos independientes, pero en la realidad, debido a la estructura de la serie (valores atípicos, rupturas estructurales, etc.), puede ser más desafiante. En tales casos, es importante seleccionar un modelo que se acerque a este objetivo.
Luego de ajustar el modelo, es momento de evaluar su capacidad para predecir observaciones que no ha visto previamente, en contraste con los valores ajustados que formaron parte del proceso de entrenamiento. Para medir la eficacia de estas predicciones, se emplean métricas de precisión o error, siendo la elección de la métrica específica dependiente de los objetivos de precisión en las predicciones. Ejemplos de métricas comunes incluyen:
• Error cuadrático medio (MSE): Cuantifica la distancia cuadrática promedio entre los valores reales y pronosticados:
El efecto cuadrado del error evita que los valores positivos y negativos se cancelen entre sí y panelicen la puntuación del error a medida que aumenta la tasa de error.
• Error cuadrático medio (RMSE): Esta es la raíz de la distancia cuadrática promedio de los valores reales y pronosticados:
Al igual que el MSE, el RMSE tiene una gran tasa de error debido al efecto al cuadrado y, por tanto, es sensible a los valores atípicos.
• Error absoluto medio (MAE): Mide la tasa de error absoluto del pronóstico:
De manera similar a MSE y RMSE, este método solo puede tener valores positivos. Esto es para evitar la cancelación de valores positivos y negativos. Por otro lado, no hay penalización por errores y, por tanto, este método no es sensible a valores atípicos.
• Error porcentual absoluto medio (MAPE): Mide el error absoluto porcentual medio:
Por ejemplo, usemos el modelo que entrenamos anteriormente para pronosticar las 12 observaciones que dejamos para probar y calificar su desempeño. Usaremos la función de pronóstico del paquete de pronóstico para pronosticar los siguientes 12 montajes (con respecto al punto final de la partición de entrenamiento):
Ahora que hemos asignado el pronóstico al objeto fc, usaremos la función de precisión del paquete de pronóstico para calificar el desempeño del modelo con respecto a los valores reales en la partición de prueba:
## 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 precisión proporciona varias métricas de error para los valores ajustados (entrenamiento) y las predicciones reales (prueba), siendo el MAPE del modelo 3,52 % para entrenamiento y 7,84 % para prueba. Es común que la tasa de error en prueba sea mayor, ya que el modelo conoce los datos de entrenamiento. No obstante, un bajo error en entrenamiento junto a un alto error en prueba señala un sobreajuste del modelo.
Una alternativa para evaluar el modelo en entrenamiento y prueba es la función “test_forecast” de TSstudio. Esta herramienta visualiza la serie real, los valores ajustados en entrenamiento y las predicciones en prueba, mostrando el RMSE y MAPE al pasar el cursor sobre estos valores (no transferibles al libro).
Es más fácil y rápido identificar ideas sobre la bondad del ajuste de los valores ajustados y pronosticados al graficar esos valores con los valores reales de la serie. Por ejemplo, se notará inmediatamente que el pico residual durante 2006 se debe a valores atípicos (o un consumo menor que el patrón normal de la serie). Además, el pronóstico real no alcanzó el pico anual de 2018. Esos conocimientos no se pueden observar con métricas de error.
Las métricas de error del modelo son MAPE 7,84% y RMSE 208,01. Para evaluar, comparamos con pronósticos de referencia. Un enfoque común es utilizar un método simple como punto de referencia. El enfoque ingenuo asume que el valor más reciente es el representante del futuro y lo sigue indefinidamente. Creamos un pronóstico ingenuo con la función ingenua del paquete de pronóstico, utilizando el conjunto de entrenamiento como entrada del modelo.
#Podemos revisar el rendimiento del modelo en las particiones de entrenamiento y prueba usando la siguiente función test_forecast:
naive_model <- naive (train, h = 12)
test_forecast( actual = USgas,
forecast.obj = naive_model,
test= test)#Utilizamos la función accuracy para evaluar el desempeño del modelo en tanto las particiones de entrenamiento como las de prueba:
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 el caso del modelo naive, no se entrena; en su lugar, los valores ajustados coinciden con los valores reales. Para series con un patrón estacional, es lógico usar un modelo ingenuo estacional que considera la variación estacional. Con “snaive_model” del paquete de pronóstico, el último valor estacional se usa como pronóstico para todas las observaciones estacionales futuras (por ejemplo, el valor más reciente de enero se proyecta para todos los futuros eneros).
snaive_model <- snaive(train, h=12)
test_forecast(actual = USgas,
forecast.obj = snaive_model,
test=test )Parece que el modelo ingenuo estacional se ajusta mejor al tipo de serie que estamos pronosticando, es decir, gas estadounidense, debido a su fuerte patrón estacional (en comparación con el modelo ingenuo). Por tanto, lo utilizaremos como punto de referencia para el modelo ARIMA. Al comparar tanto el MAPE como el RMSE de los dos modelos en la partición de prueba, queda claro que el modelo ARIMA proporciona una mejora (en términos de precisión) con respecto al modelo de referencia:
Ahora que el modelo ha sido entrenado, probado, ajustado (si es necesario) y evaluado exitosamente, podemos avanzar al último paso y finalizar el pronóstico. Este paso se basa en recalibrar los pesos o coeficientes del modelo con la serie completa. Hay dos métodos para emplear la configuración de parámetros del modelo:
■ Si el modelo se ajustó manualmente, debe utilizar los parámetros de ajuste exactos que se utilizaron en el modelo entrenado.
■ Si el modelo se ajustó automáticamente mediante un algoritmo (como la función auto.arima que usamos anteriormente), puede realizar cualquiera de las siguientes acciones:
Extraiga la configuración de parámetros que se utilizó con la partición deentrenamiento.
Deje que el algoritmo vuelva a ajustar los parámetros del modelo utilizando la serie completa, bajo el supuesto de que el algoritmo tiene la capacidad de ajustar los parámetros del modelo correctamente al entrenar el modelo con nuevos datos.
Se recomienda el uso de algoritmos para automatizar el proceso de ajuste del modelo cuando la capacidad del modelo para ajustarlo se prueba con pruebas retrospectivas. Esto le permite revisar si el algoritmo tiene la capacidad de ajustar los parámetros del modelo correctamente, según los resultados del backtesting. Por motivos de simplicidad, seguiremos usando el auto. Modelo Arima para entrenar el modelo final:
md_final <- auto.arima(USgas)
fc_final <- forecast(md_final, h=12)
#Utilización de función del paquete TSstudio
plot_forecast(fc_final,
title= "Pronóstico de consumo de gas natural en EE.UU",
Xtitle="year",
Ytitle="Miles de millones de pies cúbicos")El proceso de pronóstico tiene como objetivo principal reducir la incertidumbre en torno a los valores futuros de la serie. Aunque no podemos eliminar por completo esta incertidumbre, podemos medirla y establecer un rango alrededor de la estimación puntual del pronóstico, que representa el valor esperado del modelo para cada punto en el futuro. Esto se logra mediante el uso de un intervalo de confianza o, en el caso de un modelo bayesiano, un intervalo creíble, o a través de la simulación.
El intervalo de confianza, una técnica estadística, refleja un rango de valores posibles que contienen el valor verdadero con cierto grado de confianza. Dos factores clave lo determinan:
■ Nivel de confianza: Indica la probabilidad de que el valor verdadero esté en el rango. Un nivel de confianza mayor produce un rango más amplio.
■ Desviación estándar estimada: En el momento T+i, donde T es la longitud de la serie y i es el valor pronosticado. Menor error, rango más corto.
La función de pronóstico por defecto genera un intervalo con niveles de confianza del 80% y 95%, pero puedes personalizarlo con el argumento de nivel. Ejemplo: Utilicemos las funciones md_final y de pronóstico del modelo entrenado para los próximos 60 meses con intervalos de confianza del 80% y 90%.
fc_final2 <- forecast (md_final,
h=60,
level= c(80,90))
plot_forecast(fc_final2,
title="Pronóstico de consumo de gas natural en EE.UU")Un enfoque alternativo es simular trayectorias de pronóstico usando la distribución del modelo, lo cual requiere conocer dicha distribución. Con TSstudio, la función Forecast_sim permite realizar estas simulaciones y obtener estimaciones puntuales o evaluar probabilidades. Ejecutaremos 100 iteraciones en nuestro ejemplo.
## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutfc_final3$plot %>%
layout(title= "Consumo de gas natural - Simulación de pronóstico",
yaxis=list(title="Miles de millones de pies cúbicos"),
xaxis=list(title="año"))Finalizamos este capítulo con un enfoque sólido: el método de la Horse Race. Consiste en entrenar, probar y evaluar múltiples modelos de pronóstico, eligiendo el más eficaz en pruebas. En nuestro ejemplo, aplicamos este método a siete modelos diferentes en seis periodos de backtesting utilizando la función ts_backtesting del paquete TSstudio, que automatiza este proceso. Por defecto, se evalúan los siguientes modelos.
• auto.arima: modelo ARIMA automatizado
• bsts: modelo de series temporales estructurales bayesianas
• ets: modelo de espacio de estados de suavizado exponencial
• híbrido: un conjunto de múltiples modelos
• nnetar: modelo de series temporales de redes neuronales
• tbats: modelo de espacio de estados de suavizado exponencial, junto con transformación de Box-Cox, tendencia, errores ARMA y componentes estacionales
• HoltWinters: filtrado de Holt-Winters
El proceso de entrenamiento de un modelo de pronóstico es el paso final del análisis de series de tiempo. El objetivo de este capítulo fue presentar el principio del flujo de trabajo de previsión. Como vimos, existen varios métodos que podemos utilizar para entrenar un modelo de pronóstico, y el proceso de selección del método debe alinearse con los objetivos de pronóstico y los recursos disponibles. En los siguientes capítulos verá estas aplicaciones en la práctica. En el próximo capítulo, utilizaremos las aplicaciones del modelo de regresión lineal para pronosticar datos de series de tiempo.
El modelo de regresión lineal se usa comúnmente para cuantificar la relación entre una variable dependiente y una o varias variables independientes. Tiene aplicaciones que van desde la inferencia causal hasta la predicción de series temporales.
Este capítulo se centra en métodos de pronóstico de series de tiempo con regresión lineal, que incluyen descomponer y prever componentes como la tendencia y los patrones estacionales, manejar eventos especiales y utilizar variables externas como regresores.
Este capitulo cubre los siguientes topicos:
• Enfoques de pronóstico con modelos de regresión lineal.
• Extracción y estimación de los componentes de la serie.
• Manejo de rupturas estructurales, valores atípicos y eventos especiales.
• Series de pronóstico con multiestacionalidad.
Los siguientes paquetes se utilizarán en este capítulo.
• TSstudio: Versión 0.1.4 y superior
• trama: Versión 4.8 y superior
• dplyr: Versión 0.8.1 y superior
• lubricar: versión 1.7.4 y superior
• pronóstico: versión 8.5 y superior
El uso principal del modelo de regresión lineal es cuantificar la relación entre la variable dependiente Y (también conocida como variable de respuesta) y las variables independientes X (también conocidas como variables predictoras, impulsoras o regresadoras) de manera lineal. En otras palabras, el modelo expresa la variable dependiente como una combinación lineal de las variables independientes. Una relación lineal entre las variables dependientes e independientes se puede generalizar mediante las siguientes ecuaciones:
• En el caso de una sola variable independiente, la ecuación es la siguiente:
• Para n variables independientes, la ecuación queda de la siguiente:
El término lineal, en el contexto de la regresión, se refiere a los coeficientes del modelo, los cuales deben seguir una estructura lineal (ya que esto nos permite construir una combinación lineal a partir de las variables independientes). Por otro lado, las variables independientes pueden seguir una formación tanto lineal como no lineal.
Suponiendo que las ecuaciones anteriores representan la verdadera naturaleza de la relación lineal entre las variables dependientes e independientes, entonces el modelo de regresión lineal proporciona una estimación de esos coeficientes (es decir, Po, P1,…, Pn), que puede formalizarse mediante las siguientes ecuaciones:
• Para el modelo de regresión lineal univariado, la ecuación es la siguiente:
• Para el modelo de regresión lineal multivariado, la ecuación es la siguiente:
El enfoque de estimación más común es aplicar los mínimos cuadrados ordinarios (OLS) como técnica de optimización para minimizar la suma de cuadrados de los residuos (-1). Cuadrar los residuos tiene dos efectos:
■ Evita que los valores positivos y negativos se cancelen entre sí al sumarlos.
■ El efecto cuadrado proporciona una penalización exponencial para los residuos con mayor distancia, a medida que su costo aumenta.
Existen múltiples técnicas de estimación además del MCO, como la máxima verosimilitud, el método de momentos y el bayesiano. Aunque esos métodos no están dentro del alcance del libro, se recomienda encarecidamente leer y conocer enfoques alternativos.
El OLS es un método de optimización simple que se basa en álgebra y cálculo lineal básico, o cálculo matricial. (Esta sección es para conocimientos generales; si no está familiarizado con el cálculo matricial, puede omitir esta sección). El objetivo del MCO es identificar los coeficientes que minimizan la suma de cuadrados de los residuos. Supongamos que el residual de la i observación se define como lo siguiente:
Antes de aplicar el método MCO para minimizar la suma de cuadrados de los residuos, por razones de simplicidad, transformaremos el representante de la función de costos en una formación matricial:
Los principales supuestos del modelo OLS son los siguientes:
• Los coeficientes del modelo deben seguir una estructura lineal (por ejemplo, Y = B+Bex es un modelo lineal pero Y=+X no lo es).
• No existe una colinealidad perfecta entre las variables independientes X1, X2,…, Xn. En otras palabras, ninguna de las variables independientes es una combinación lineal de ninguna de las otras variables independientes.
• Todas las variables independientes deben tener una varianza distinta de cero (o no constante).
• El término de error €, condicionado a la matriz de variables independientes X, es una variable independiente e idénticamente distribuida (i.i.d) con media 0 y varianza constante o².
• Tanto la variable dependiente como la independiente se extraen de la población en una muestra aleatoria. Esta suposición no se cumple cuando se hacen regresiones de datos de series temporales, ya que normalmente las observaciones tienen cierto grado de correlación. Por lo tanto, este supuesto se relaja cuando se hace una regresión de datos de series de tiempo.
A diferencia de modelos como ARIMA o Holt-Winters específicamente diseñados para series de tiempo, el modelo de regresión lineal es más genérico y versátil, aplicándose en la inferencia causal y análisis predictivo. El proceso de pronóstico con regresión lineal implica:
Identificar la estructura, patrones y características clave de la serie, incluyendo valores atípicos.
Transformar estas características en variables de entrada para crear un modelo de pronóstico mediante regresión.
En este contexto, los componentes principales son la tendencia y las estacionales, que se exploran y convierten en variables del modelo de regresión en la siguiente sección.
En el Capítulo 5, Descomposición de datos de series temporales, introdujimos los componentes estructurales de la serie: los componentes de tendencia, ciclo, estacional e irregular, y los métodos para descomponerlos con la función de descomposición. Estos componentes desempeñan un papel crucial en la descripción y el análisis de las series de tiempo, ya que ayudan a comprender cómo los datos evolucionan y varían a lo largo del tiempo.
En ese capítulo, presentamos la función de descomposición, que permite separar la serie temporal en sus componentes estructurales individuales. La descomposición facilita la visualización de la tendencia a largo plazo, los patrones estacionales que se repiten y las fluctuaciones irregulares impredecibles presentes en los datos.A medida que avanzamos en el estudio del análisis de series temporales, profundizaremos en cómo modelar y pronosticar cada uno de estos componentes, lo que nos permitirá comprender mejor los procesos subyacentes que generan los datos y hacer proyecciones más precisas para el futuro.
Antes de crear los datos de regresión que representan la tendencia de la serie y los componentes estacionales, primero debemos comprender su estructura. En los siguientes ejemplos, demostraremos el proceso de creación de nuevas funciones a partir de series existentes utilizando la serie USgas. Carguemos la serie desde el paquete TSstudio nuevamente y grafiquemos con la función ts_plot:
library(TSstudio)
data(USgas)
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")## The USgas series is a ts object with 1 variable and 238 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2019 10#Como se puede ver en el gráfico de la serie, y como vimos en los capítulos anteriores, USgas es una serie mensual con un fuerte componente estacional mensual y una línea de tendencia bastante estable. Podemos explorar más a fondo la estructura de los componentes de la serie con la función ts_decompose:
ts_decompose(USgas)Puede ver en el gráfico anterior que la tendencia de la serie es bastante plana entre 2000 y 2010, y tiene un crecimiento bastante lineal en el futuro. Por tanto, la tendencia general entre 2000 y 2018 no es estrictamente lineal. Esta es una idea importante que nos ayudará a definir la entrada de tendencia para el modelo de regresión.
Antes de usar la función 1m, la función de regresión lineal R incorporada del paquete de estadísticas, tendremos que transformar la serie de un objeto ts a un dato. objeto de marco. Por lo tanto, utilizaremos la función ts_to_prophet del paquete TSstudio:
USgas_df <- ts_to_prophet(USgas)
#La función transforma el objeto ts en dos columnas de datos. marco, donde las dos columnas representan los componentes numéricos y de tiempo de la serie, respectivamente:
head(USgas_df)## ds y
## 1 2000-01-01 2510.5
## 2 2000-02-01 2330.7
## 3 2000-03-01 2050.6
## 4 2000-04-01 1783.3
## 5 2000-05-01 1632.9
## 6 2000-06-01 1513.1Después transformamos la serie en datos. objeto de marco, podemos comenzar a crear las características de entrada de regresión. La primera característica que crearemos es la tendencia de la serie. Un enfoque básico para construir la variable de tendencia es indexar las series de observaciones en orden cronológico:
La regresión de la serie con el índice de la serie proporciona una estimación del crecimiento marginal de mes a mes, ya que el índice está en orden cronológico con incrementos constantes.
La segunda característica que queremos crear es el componente estacional. Como queremos medir la contribución de cada unidad de frecuencia a la oscilación de la serie, usaremos una variable categórica para cada unidad de frecuencia. En el caso de la serie USgas, las unidades de frecuencia representan los meses del año, por lo que crearemos una variable categórica con 12 categorías, correspondiendo cada categoría a un mes específico del año. Usaremos la función mes del paquete lubridate para extraer el mes del año de la variable de fecha ds:
##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionUSgas_df$seasonal <- factor(month(USgas_df$ds, label = T), ordered=FALSE)
#Usamos la función de factor para convertir el resultado de la función de mes en una variable categórica no ordenada. Revisemos ahora el marco de datos después de agregar las nuevas características
head(USgas_df)## ds y trend seasonal
## 1 2000-01-01 2510.5 1 ene
## 2 2000-02-01 2330.7 2 feb
## 3 2000-03-01 2050.6 3 mar
## 4 2000-04-01 1783.3 4 abr
## 5 2000-05-01 1632.9 5 may
## 6 2000-06-01 1513.1 6 junPor último, pero no menos importante, antes de comenzar a hacer una regresión de la serie con esas características, dividiremos la serie en una partición de entrenamiento y prueba. Estableceremos los últimos 12 meses de la serie como partición de prueba:
h <- 12 # setting a testing partition length
train <- USgas_df[1:(nrow(USgas_df) - h), ]
test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ]Ahora, después de crear los marcos de datos de entrenamiento y prueba, revisemos cómo el modelo de regresión captura cada uno de los componentes por separado y todos juntos.
Primero modelaremos la tendencia de la serie haciendo una regresión de la serie con la variable de tendencia, en el partición de entrenamiento:
##
## Call:
## lm(formula = y ~ trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -547.2 -307.4 -153.2 333.1 1052.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1751.0074 52.6435 33.26 < 2e-16 ***
## trend 2.4489 0.4021 6.09 4.86e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 394.4 on 224 degrees of freedom
## Multiple R-squared: 0.1421, Adjusted R-squared: 0.1382
## F-statistic: 37.09 on 1 and 224 DF, p-value: 4.861e-09En la regresión anterior, el coeficiente de la variable de tendencia es estadísticamente significativo (nivel de significancia 0,001), pero el R cuadrado ajustado es bajo, lo cual es coherente ya que la variación principal está relacionada con el patrón estacional.
En el análisis de regresión, el valor p en la cuarta columna mide la probabilidad de cometer un error tipo I, determinando si rechazamos la hipótesis nula. Generalmente, se utiliza un nivel de significancia (α) como umbral, como 0,1, 0,05, 0,01, etc.
Para dar contexto a los números, se emplea la visualización de datos. Usaremos el modelo para predecir los valores ajustados en el conjunto de entrenamiento y los valores pronosticados en el conjunto de prueba. La función de predicción del paquete de estadísticas calcula estos valores en base al modelo. En nuestro caso, aplicaremos esto al modelo de tendencia que hemos entrenado anteriormente.
Crearemos una función de utilidad que traza la serie y el resultado del modelo, utilizando el paquete plotly:
library(plotly)
library(ggplot2)
plot_lm <- function(data, train, test, title = NULL){
p <- plot_ly(data = data,
x = ~ ds,
y = ~ y,
type = "scatter",
mode = "line",
name = "Actual") %>%
add_lines(x = ~ train$ds,
y = ~ train$yhat,
line = list(color = "red"),
name = "Fitted") %>%
add_lines(x = ~ test$ds,
y = ~ test$yhat,
line = list(color = "green", dash = "dot", width = 3),
name = "Forecasted") %>%
layout(title = title,
xaxis = list(title = ""),
yaxis = list(title = "Billion Cubic Feet"),
legend = list(x = 0.05, y = 0.95))
return(p)
}Los argumentos de la función son los siguientes:
• Datos: os datos de entrada, un dato. objeto frame que sigue la misma estructura que la de USgas_df.
• Entrenar: El conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo
• Prueba: Asimismo, el conjunto de pruebas correspondiente que se utilizó para evaluar el modelo de pronóstico
• Título: el título de la trama, de forma predeterminada, establecido en NULL
Configuremos las entradas de la función plot_1m con la salida del modelo:
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend Component of the Series")En general, el modelo pudo capturar el movimiento general de la tendencia, pero una tendencia lineal puede no capturar la ruptura estructural de la tendencia que ocurrió alrededor de 2010. Más adelante, en este capítulo, veremos un método avanzado para capturar un tendencia no lineal.
Por último, pero no menos importante, para el análisis comparativo, queremos medir la tasa de error del modelo tanto en el conjunto de entrenamiento como en el de prueba:
mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend## [1] 0.1644088 0.1299951El proceso de modelación y pronóstico del componente estacional sigue el mismo proceso que aplicamos con la tendencia, haciendo una regresión de la serie con la variable estacional que creamos antes:
##
## Call:
## lm(formula = y ~ seasonal, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -608.98 -162.34 -50.77 148.40 566.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2774.28 49.75 55.759 < 2e-16 ***
## seasonalfeb -297.92 70.36 -4.234 3.41e-05 ***
## seasonalmar -479.10 70.36 -6.809 9.77e-11 ***
## seasonalabr -905.28 70.36 -12.866 < 2e-16 ***
## seasonalmay -1088.42 70.36 -15.468 < 2e-16 ***
## seasonaljun -1105.49 70.36 -15.711 < 2e-16 ***
## seasonaljul -939.35 70.36 -13.350 < 2e-16 ***
## seasonalago -914.12 70.36 -12.991 < 2e-16 ***
## seasonalsept -1114.74 70.36 -15.843 < 2e-16 ***
## seasonaloct -1022.21 70.36 -14.527 < 2e-16 ***
## seasonalnov -797.53 71.33 -11.180 < 2e-16 ***
## seasonaldic -256.67 71.33 -3.598 0.000398 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 216.9 on 214 degrees of freedom
## Multiple R-squared: 0.7521, Adjusted R-squared: 0.7394
## F-statistic: 59.04 on 11 and 214 DF, p-value: < 2.2e-16Dado que hacemos una regresión de la variable dependiente con una variable categórica, el modelo de regresión crea coeficientes para 11 de las 12 categorías, que son aquellas integradas con los valores de pendiente. Como puede ver en el resumen de regresión del modelo estacional, todos los coeficientes del modelo son estadísticamente significativos. Además, se puede observar que el R cuadrado ajustado del modelo estacional es algo mayor con respecto al modelo de tendencia (0,78 frente a 0,1).
Antes de trazar el modelo ajustado y los valores pronosticados con la función plot_1m, actualizaremos los valores de yhat con la función de predicción:
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")Como puede verse en el gráfico anterior, el modelo está haciendo un trabajo bastante bueno al capturar la estructura del patrón estacional de la serie. Sin embargo, se puede observar que falta la tendencia de la serie. Antes de agregar los componentes de tendencia y estacional, calificaremos el desempeño del modelo:
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal## [1] 0.08628973 0.21502100La alta tasa de error en el conjunto de pruebas está relacionada con el componente de tendencia que no se incluyó en el modelo. El siguiente paso es unir los dos componentes en un modelo y pronosticar los valores característicos de la serie y repasemos el resumen del modelo después de hacer una regresión de la serie con los componentes tendencial y estacional.
##
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -514.73 -77.17 -17.70 85.80 336.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2488.8994 32.6011 76.344 < 2e-16 ***
## seasonalfeb -300.5392 41.4864 -7.244 7.84e-12 ***
## seasonalmar -484.3363 41.4870 -11.674 < 2e-16 ***
## seasonalabr -913.1334 41.4880 -22.010 < 2e-16 ***
## seasonalmay -1098.8884 41.4895 -26.486 < 2e-16 ***
## seasonaljun -1118.5855 41.4913 -26.960 < 2e-16 ***
## seasonaljul -955.0563 41.4936 -23.017 < 2e-16 ***
## seasonalago -932.4482 41.4962 -22.471 < 2e-16 ***
## seasonalsept -1135.6874 41.4993 -27.366 < 2e-16 ***
## seasonaloct -1045.7687 41.5028 -25.198 < 2e-16 ***
## seasonalnov -808.0016 42.0617 -19.210 < 2e-16 ***
## seasonaldic -269.7642 42.0635 -6.413 9.05e-10 ***
## trend 2.6182 0.1305 20.065 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 127.9 on 213 degrees of freedom
## Multiple R-squared: 0.9142, Adjusted R-squared: 0.9094
## F-statistic: 189.2 on 12 and 213 DF, p-value: < 2.2e-16La regresión de la serie con los componentes de tendencia y estacional juntos proporciona un impulso adicional al R cuadrado ajustado del modelo de 0,78 a 0,91. Esto se puede ver en el gráfico del resultado del modelo:
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")Midamos la puntuación MAPE del modelo tanto en la partición de entrenamiento como en la de prueba:
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1## [1] 0.04774945 0.09143290Incorporar una regresión con componentes de tendencia y estacionales mejora significativamente la calidad y precisión del modelo. No obstante, al analizar el ajuste y pronóstico del modelo, se evidencia que la tendencia es demasiado lineal y carece de la ruptura estructural de la serie. Para abordar esto, añadir un componente polinómico al modelo puede ofrecer mejoras adicionales en la precisión.
Una técnica simple para capturar una tendencia no lineal consiste en agregar un componente polinómico para capturar la curvatura a lo largo del tiempo. Utilizaremos el argumento “I” que permite realizar operaciones matemáticas en las entradas. De esta manera, añadiremos un polinomio de segundo grado a la tendencia de entrada.
##
## Call:
## lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -468.47 -54.66 -2.21 63.11 294.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
## seasonalfeb -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
## seasonalmar -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
## seasonalabr -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
## seasonalmay -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
## seasonaljun -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
## seasonaljul -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
## seasonalago -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
## seasonalsept -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
## seasonaloct -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
## seasonalnov -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
## seasonaldic -2.618e+02 3.590e+01 -7.293 5.95e-12 ***
## trend -1.270e+00 4.472e-01 -2.840 0.00494 **
## I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
## F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16Agregar el polinomio de segundo grado al modelo de regresión no condujo a una mejora significativa de la bondad de ajuste del modelo. En el otro modelo, como se puede ver en el siguiente gráfico de resultados del modelo, este simple cambio en la estructura del modelo nos permite capturar la ruptura estructural de la tendencia a lo largo del tiempo:
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")#Como podemos ver en el modelo que sigue la puntuación MAPE, la precisión del modelo mejoró significativamente al agregar la tendencia polinómica al modelo de regresión, ya que el error en el conjunto de pruebas cayó del 9,2% al 4,5%:
mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md2## [1] 0.03688770 0.04212618Hasta ahora, hemos examinado cómo transformar manualmente un objeto ts en un modelo de regresión lineal para el pronóstico. Sin embargo, la función “tslm” del paquete de pronóstico automatiza este proceso, permitiéndote configurar el componente de regresión y otros parámetros.
Repetiremos el ejemplo anterior pronosticando las últimas 12 observaciones de la serie USgas utilizando “tslm”. Para ello, dividiremos la serie en conjuntos de entrenamiento y prueba utilizando “ts_split”.
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$test
#A continuación, aplicaremos la misma fórmula que usamos para crear el modelo de pronóstico md2 anterior usando la función tslm:
library(forecast)
md3 <- tslm(train.ts ~ season + trend + I(trend^2))
summary(md3)##
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -468.47 -54.66 -2.21 63.11 294.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
## season2 -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
## season3 -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
## season4 -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
## season5 -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
## season6 -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
## season7 -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
## season8 -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
## season9 -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
## season10 -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
## season11 -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
## season12 -2.618e+02 3.590e+01 -7.293 5.95e-12 ***
## trend -1.270e+00 4.472e-01 -2.840 0.00494 **
## I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
## F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16Como puede observar en el resultado anterior, ambos modelos (md2 y md3) son idénticos.
Existen varias ventajas al utilizar la función ts1m, en lugar de configurar manualmente un modelo de regresión para la serie con la función 1m:
• Eficiencia: no requiere transformar la serie en un objeto de marco de datos ni ingeniería de características.
• El objeto de salida admite todas las funciones del pronóstico (como las funciones de precisión y verificación de residuos) y los paquetes de TSstudio (como las funciones test_forecast y plot_forecast).
En algunos casos, los datos de series temporales pueden contener patrones inusuales que se repiten o no con el tiempo. Los siguientes son ejemplos de tales eventos:
■ Valores atípicos: uno o varios eventos únicos que se salen de los patrones normales de la serie.
■ Ruptura estructural: Un hecho significativo que cambia los patrones históricos de la serie. Un ejemplo común es un cambio en el crecimiento de la serie.
■ Eventos recurrentes no estacionales: Un evento que se repite de un ciclo a otro, pero el momento en que ocurren cambia de un ciclo a otro. Un ejemplo común de este tipo de eventos son las vacaciones de Semana Santa, que ocurren cada año alrededor de marzo/abril.
Eventos inusuales pueden sesgar los coeficientes del modelo si no se consideran. Para corregir esto, se usan variables binarias que permiten al modelo ignorar o ajustar estos eventos. Por ejemplo, si se detecta una ruptura estructural en la tendencia de una serie temporal, como un cambio en el crecimiento, se puede agregar una variable binaria que distinga antes y después de ese cambio.
Para aplicar esta técnica en un modelo “tslm” con variables externas, se necesita un marco de datos con las variables adecuadas. Aquí hay un ejemplo de cómo crear una variable binaria para marcar el cambio antes y después de un año específico en la tabla “USgas_df”.
r <- which(USgas_df$ds == as.Date("2014-01-01"))
USgas_df$s_break <- ifelse(year(USgas_df$ds) >= 2010, 1, 0)
USgas_df$s_break[r] <- 1
md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)
summary(md3)##
## Call:
## tslm(formula = USgas ~ season + trend + I(trend^2) + s_break,
## data = USgas_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -469.25 -50.68 -2.66 63.63 275.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.661e+03 3.200e+01 83.164 < 2e-16 ***
## season2 -3.054e+02 3.448e+01 -8.858 2.61e-16 ***
## season3 -4.849e+02 3.448e+01 -14.062 < 2e-16 ***
## season4 -9.272e+02 3.449e+01 -26.885 < 2e-16 ***
## season5 -1.108e+03 3.449e+01 -32.114 < 2e-16 ***
## season6 -1.127e+03 3.450e+01 -32.660 < 2e-16 ***
## season7 -9.568e+02 3.450e+01 -27.730 < 2e-16 ***
## season8 -9.340e+02 3.451e+01 -27.061 < 2e-16 ***
## season9 -1.138e+03 3.452e+01 -32.972 < 2e-16 ***
## season10 -1.040e+03 3.453e+01 -30.122 < 2e-16 ***
## season11 -7.896e+02 3.497e+01 -22.577 < 2e-16 ***
## season12 -2.649e+02 3.498e+01 -7.571 9.72e-13 ***
## trend -1.928e+00 4.479e-01 -4.304 2.51e-05 ***
## I(trend^2) 1.862e-02 1.676e-03 11.113 < 2e-16 ***
## s_break 6.060e+01 2.836e+01 2.137 0.0337 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109 on 223 degrees of freedom
## Multiple R-squared: 0.9423, Adjusted R-squared: 0.9387
## F-statistic: 260.3 on 14 and 223 DF, p-value: < 2.2e-16Una de las principales ventajas del modelo de regresión, frente a los modelos tradicionales de series temporales como ARIMA o Holt-Winters, es que proporciona una amplia gama de opciones de personalización y nos permite modelar y pronosticar datos de series temporales complejas como series con multiestacionalidad. En los siguientes ejemplos, utilizaremos la serie UKgrid para demostrar el enfoque de pronóstico de una serie multiestacional con un modelo de regresión lineal.
La serie UKgrid, que registra la demanda eléctrica en el Reino Unido, se encuentra en el paquete UKgrid. Esta serie es de alta frecuencia, registrando datos cada media hora. Utilizaremos la función “extracción de cuadrícula” de UKgrid para definir la serie, especificando sus características, como formato, variables y frecuencia. Esta función nos permite cambiar la frecuencia de los datos de media hora a una frecuencia más baja, como por hora, por día o por mes. Nuestro objetivo es pronosticar la demanda diaria para los próximos 365 días, así que configuraremos la serie en una frecuencia diaria utilizando la estructura de datos “data.frame”.
library(UKgrid)
UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")
head(UKdaily)## TIMESTAMP ND
## 1 2005-04-01 1920069
## 2 2005-04-02 1674699
## 3 2005-04-03 1631352
## 4 2005-04-04 1916693
## 5 2005-04-05 1952082
## 6 2005-04-06 1964584Como se puede ver en el gráfico anterior, la serie tiene una clara tendencia bajista y un patrón estacional de cadena. Como vimos en el Capítulo 6, Análisis de estacionalidad, esta serie tiene múltiples patrones de estacionalidad:
■ Diario: Un ciclo de 365 días al año
■ Día de la semana: un ciclo de 7 días
■ Mensual: Efectuado por el clima
ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
title = "UK the Daily National Grid Demand Heatmap")El mapa de calor de la serie muestra que la demanda aumenta durante las semanas de invierno (semanas 1 a 12 y 44 a 52) y varía según los días de la semana, con picos en días laborables y disminuciones los fines de semana.
Para capturar los componentes estacionales de la serie, estableceremos la serie como frecuencia diaria y crearemos las dos funciones siguientes:
• Indicador del día de la semana
• Indicador del mes del año
Suponemos una fuerte correlación con los rezagos estacionales y creamos una variable de retraso de 365 observaciones utilizando el paquete dplyr.
##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionUKdaily <- UKdaily %>%
mutate(wday = wday(TIMESTAMP, label = TRUE),
month = month(TIMESTAMP, label = TRUE),
lag365 = dplyr::lag(ND, 365)) %>%
filter(!is.na(lag365)) %>%
arrange(TIMESTAMP)
str(UKdaily)## 'data.frame': 4939 obs. of 5 variables:
## $ TIMESTAMP: Date, format: "2006-04-01" "2006-04-02" ...
## $ ND : int 1718405 1691341 1960325 2023886 2026204 2008422 1981175 1770440 1749715 2012865 ...
## $ wday : Ord.factor w/ 7 levels "dom\\."<"lun\\."<..: 7 1 2 3 4 5 6 7 1 2 ...
## $ month : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 4 4 4 4 4 4 4 4 4 4 ...
## $ lag365 : int 1920069 1674699 1631352 1916693 1952082 1964584 1990895 2003982 1811436 1684720 ...Como la entrada de la función ts1m debe estar en formato ts (al menos para la serie), convertiremos la serie en un objeto ts. Usaremos la primera marca de tiempo de la serie y las funciones año e ydía (el día del año) del paquete lubridate para establecer el punto de partida del objeto:
start_date <- min(UKdaily$TIMESTAMP)
UK_ts <- ts(UKdaily$ND,
start = c(year(start_date), yday(start_date)),
frequency = 365)
#ts_acf(UK_ts, lag.max = 365 * 4)Ahora, tras crear nuevas características y configurar el objeto ts, procedemos a dividir la serie de entrada y el objeto de características externas (UKdaily) en conjuntos de entrenamiento y prueba. Dado que queremos pronosticar las próximas 365 observaciones y la serie es lo suficientemente extensa (2540 observaciones), establecemos la partición de prueba en un horizonte de pronóstico de 365 observaciones. Usaremos una variable llamada “h” con valor 365 para definir las particiones y el horizonte de pronóstico.
h <- 365
UKpartitions <- ts_split(UK_ts, sample.out = h)
train_ts <- UKpartitions$train
test_ts <- UKpartitions$test
train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]Una vez que hayamos creado las diferentes funciones para el modelo, estamos listos para iniciar el proceso de capacitación del modelo de pronóstico. Usaremos la partición de entrenamiento y comenzaremos a entrenar los siguientes tres modelos:
• Modelo de línea base: Regresión de la serie con el componente estacional y de tendencia utilizando las características integradas de la función tslm. Como establecemos la frecuencia de la serie en 365, la característica estacional de la serie se refiere a la estacionalidad diaria.
• Modelo multiestacional: Agregando los indicadores de día de la semana y mes del año para capturar la multiestacionalidad de la serie.
• Un modelo multiestacional con rezago estacional: Utilizando, además de los indicadores estacionales, la variable de rezago estacional.
La comparación de estos tres modelos se basará en los siguientes criterios:
• Rendimiento del modelo en el conjunto de entrenamiento y prueba utilizando la puntuación MAPE.
• Visualizar los valores ajustados y pronosticados versus los valores reales de la serie usando la función test_forecast.
md_tslm1 <- tslm(train_ts ~ season + trend)
fc_tslm1 <- forecast(md_tslm1, h = h)
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -4.781286e-12 121551.4 100439.64 -0.5721337 6.294111 0.8418677
## Test set -1.740215e+04 123156.6 96785.27 -1.8735336 7.160573 0.8112374
## ACF1 Theil's U
## Training set 0.5277884 NA
## Test set 0.5106681 1.071899#Mejoramiento de la precisión del modelo agregando la semana y el mes del año al promedio
md_tslm2 <- tslm(train_ts ~ season + trend + wday, data = train_df)
#Ahora repetiremos el mismo proceso que antes, utilizando un pronóstico con el modelo entrenado y evaluando el rendimiento del modelo
fc_tslm2 <- forecast(md_tslm2, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm2,
test = test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set 8.172823e-12 70245.98 52146.79 -0.1738708 3.167605 0.4370853
## Test set -1.764563e+04 80711.71 65373.21 -1.3715505 4.682071 0.5479470
## ACF1 Theil's U
## Training set 0.7513664 NA
## Test set 0.6075598 0.68445Por último, pero no menos importante, agreguemos la variable de retraso al modelo y repitamos el mismo proceso que antes.Obtenemos la trama como se ve en el siguiente gráfico:
md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365, data = train_df)
fc_tslm3 <- forecast(md_tslm3, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm3,
test = test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -9.836939e-13 69904.92 52006.75 -0.1717563 3.163385 0.4359116
## Test set -1.754061e+04 81783.55 65957.66 -1.3613252 4.722083 0.5528457
## ACF1 Theil's U
## Training set 0.7500146 NA
## Test set 0.6094666 0.6925767Ahora es el momento de seleccionar un modelo. En este punto, está claro que el segundo y tercer modelo funcionan mejor que el primer modelo. Como tanto el segundo como el tercer modelo lograron una puntuación MAPE bastante similar, con una pequeña ventaja para el tercer modelo, deberíamospreguntarnos si una mejora de menos del 0,2% en la tasa de error del conjunto de pruebas vale el costo de usar el variable de rezago (es decir, la pérdida de 365 observaciones y el costo adicional de un grado de libertad para el modelo).
No hay una respuesta sencilla a esta pregunta. Además, depende del objetivo del pronóstico. Se recomienda que consideres la siguiente prueba:
• La primera pregunta que debe hacerse en este caso: ¿Es la variable de rezago estadísticamente significativa? Si la variable no es estadísticamente significativa, no tiene sentido continuar la discusión y es mejor descartarla. En el caso del tercer modelo, podemos utilizar la función de resumen para observar el nivel de significancia de la variable lag365.
## Estimate Std. Error t value Pr(>|t|)
## lag365 -0.06038328 0.01545495 -3.907051 9.490321e-05#Aplicación de una última prueba ANOVA para comprobar si la variable adicional es significativa
anova(md_tslm3)## Analysis of Variance Table
##
## Response: train_ts
## Df Sum Sq Mean Sq F value Pr(>F)
## season 364 1.4279e+14 3.9227e+11 73.5348 < 2.2e-16 ***
## trend 1 7.2634e+13 7.2634e+13 13615.7078 < 2.2e-16 ***
## wday 6 4.5009e+13 7.5016e+12 1406.2214 < 2.2e-16 ***
## month 11 1.3721e+11 1.2473e+10 2.3382 0.007266 **
## lag365 1 8.1432e+10 8.1432e+10 15.2650 9.49e-05 ***
## Residuals 4190 2.2352e+13 5.3345e+09
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La prueba ANOVA también indicó que la variable lag365 es estadísticamente significativa.
• Backtesting: El backtesting ayuda a verificar si la variable de retraso contribuye a la precisión del modelo. Realizar pruebas retrospectivas para ambos modelos es un ejercicio útil para validar su contribución en varios períodos de prueba.
En aras de la simplicidad, seguiremos los criterios de precisión y seleccionaremos el tercer modelo para pronosticar la serie. El siguiente paso es volver a entrenar el modelo en todas las series y pronosticar los próximos 365 días:
Justo antes de finalizar el pronóstico, hagamos un análisis de los residuos del modelo seleccionado utilizando la función checkresiduals:
##
## Breusch-Godfrey test for serial correlation of order up to 730
##
## data: Residuals from Linear regression model
## LM test = 3301.1, df = 730, p-value < 2.2e-16El gráfico de resumen de residuos muestra que los residuos no son ruido blanco; muestran cierta autocorrelación, lo que sugiere que el modelo no capturó todos los patrones de la serie. Abordar esto implica encontrar variables adicionales que puedan explicar la variación en los residuos, pero identificar tales variables es un desafío, especialmente si deben preverse con antelación, como los patrones climáticos que afectan la demanda de electricidad. Un enfoque alternativo es tratar los residuos como una serie de tiempo separada y modelarlos con un modelo ARIMA de regresión de errores, como se detallará en el Capítulo 11 sobre Pronósticos con modelos ARIMA.
Finalicemos el proceso y utilicemos el modelo entrenado seleccionado para pronosticar las 365 observaciones futuras. Como usamos variables externas con la función tslm, tendremos que generar sus valores futuros. Esto es relativamente simple, ya que utilizamos variables deterministas. Por lo tanto, crearemos un dato. objeto de marco con los valores de wday, mes y lag365 para las próximas 365 observaciones futuras. Un enfoque simplista es crear primero las fechas correspondientes de las observaciones pronosticadas y luego extraer de este objeto el día de la semana y el mes del año:
UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) + days(1),
by = "day",
length.out = h))
UK_fc_df$wday <- factor(lubridate::wday(UK_fc_df$date, label = TRUE), ordered = FALSE)
UK_fc_df$month <- factor(month(UK_fc_df$date, label = TRUE), ordered = FALSE)
UK_fc_df$lag365 <- tail(UKdaily$ND, h)
UKgrid_fc <- forecast(final_md, h = h, newdata = UK_fc_df)
plot_forecast(UKgrid_fc,
title = "The UK National Demand for Electricity Forecast",
Ytitle = "MW",
Xtitle = "Year")En este capítulo, presentamos las aplicaciones de pronóstico del modelo de regresión lineal. Aunque el modelo de regresión lineal no fue diseñado para manejar datos de series de tiempo, con una ingeniería de características simple podemos transformar un problema de pronóstico en un problema de regresión lineal. La principal ventaja del modelo de regresión lineal con respecto a otros modelos tradicionales de series temporales es la capacidad del modelo para incorporar variables y factores externos. Sin embargo, este modelo puede manejar series temporales con patrones multiestacionales, como vimos con el pronóstico de la demanda de electricidad del Reino Unido. Por último, pero no menos importante, los enfoques de pronóstico que demostramos en este capítulo serán la base para el modelado avanzado con modelos de aprendizaje automático que discutiremos en el Capítulo 12, Pronósticos con modelos de aprendizaje automático.
En el próximo capítulo, presentaremos los métodos de suavizado exponencial, una familia de modelos de pronóstico basados en un enfoque de promedio móvil ponderado.