Hands-On-Time-Series-Analysis-with-R
Juan Narvaez
2024-04-24
library(htmltools)#Sección - 8
Estrategias de Pronóstico
Hasta ahora, nos hemos centrado principalmente en los dos primeros componentes del flujo de trabajo del análisis de series temporales: preprocesamiento de datos y análisis descriptivo. A partir de este capítulo, cambiaremos de enfoque y nos centraremos en el tercer y último componente del análisis: el pronóstico. Antes de adentrarnos en los diferentes modelos de pronóstico en los próximos capítulos, presentaremos los principales elementos del flujo de trabajo de pronóstico. Esto incluye enfoques para entrenar un modelo de pronóstico, evaluación del rendimiento y métodos de referencia. Esto le proporcionará un conjunto de herramientas para diseñar y construir un modelo de pronóstico según el objetivo del análisis. Este capítulo cubre los siguientes temas:
• Enfoques de entrenamiento y pruebas para un modelo de pronóstico • Métodos de evaluación del rendimiento y matrices de medición de errores • Métodos de referencia • Cuantificación de la incertidumbre del pronóstico con intervalos de confianza y simulación
El flujo de trabajo de pronóstico
El pronóstico tradicional de series temporales sigue el mismo flujo de trabajo que la mayoría de los campos de análisis predictivo, como regresión o clasificación, e incluye típicamente los siguientes pasos: 1. Preparación de datos: Aquí, preparamos los datos para el proceso de entrenamiento y pruebas del modelo. Este paso incluye dividir la serie en particiones de entrenamiento (dentro de la muestra) y pruebas (fuera de la muestra), crear nuevas características (cuando corresponda) y aplicar una transformación si es necesario (por ejemplo, transformación logarítmica, escalado, etc.). 2. Entrenar el modelo: Aquí, usamos 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 discutiremos métricas de error comunes en detalle). 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. Probar el modelo: Aquí, utilizamos el modelo entrenado para pronosticar las observaciones correspondientes de la partición de pruebas. 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, después de que el modelo haya sido entrenado y probado, es hora de evaluar el rendimiento general del modelo en ambas particiones de entrenamiento y pruebas. Basándonos en el proceso de evaluación del modelo, si el modelo cumple con cierto umbral o criterio, entonces o retenemos 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.
Por otro lado, este proceso tiene características únicas que lo distinguen de otros campos predictivos: ⚫ Las particiones de entrenamiento y pruebas deben estar ordenadas cronológicamente, en lugar de muestreo aleatorio. ⚫ Típicamente, una vez que hemos entrenado y probado el modelo utilizando las particiones de entrenamiento y pruebas, 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 podría ser impactante y aterrador para personas con experiencia en modelado de aprendizaje automático tradicional o aprendizaje supervisado, ya que típicamente conduce al sobreajuste y otros problemas. Discutiremos la razón detrás de esto y cómo evitar el sobreajuste más adelante.
Enfoques de entrenamiento
Uno de los elementos centrales del flujo de trabajo de pronóstico es el proceso de entrenamiento del modelo. La calidad del entrenamiento del modelo tendrá un impacto directo en la salida 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 retrasos pasados y variables externas de manera predictiva. • Ajustar los parámetros del modelo (cuando corresponda). • Que el modelo sea escalable en nuevos datos, o en otras palabras, evitar el sobreajuste. Como mencionamos anteriormente, antes del proceso de entrenamiento, la serie se divide en particiones de entrenamiento y pruebas, donde el modelo se entrena en la partición de entrenamiento y se prueba en la partición de pruebas. Estas particiones deben estar en orden cronológico, independientemente del enfoque de entrenamiento que se haya utilizado. La razón principal de esto es que la mayoría de los modelos de series temporales establecen una relación matemática entre la serie en términos de sus retrasos pasados y términos de error.
Entrenamiento con una sola partición de entrenamiento y pruebas
Uno de los enfoques de entrenamiento más comunes es el uso de particiones únicas de entrenamiento y pruebas (o una sola fuera de la muestra). Este enfoque simplista se basa en dividir la serie en particiones de entrenamiento y pruebas (o dentro y fuera de la muestra, respectivamente), entrenar el modelo en la partición de entrenamiento y probar su rendimiento en el conjunto de pruebas:
Como se puede ver en el diagrama anterior, las particiones de entrenamiento y pruebas están ordenadas y organizadas cronológicamente. Este enfoque tiene un solo parámetro: la longitud de la fuera de la muestra (o la longitud de la partición de pruebas). Típicamente, la longitud de la partición de pruebas se deriva de las siguientes reglas empíricas: • La longitud de la partición de pruebas 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 longitud del horizonte de pronóstico (siempre y cuando no viole el término anterior). Esto está principalmente relacionado 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 proporcionar una estimación más cercana del error esperado del pronóstico. Por ejemplo, si tenemos una serie mensual con 72 observaciones (o 6 años) y el objetivo es pronosticar el próximo año (o 12 meses), tendría sentido usar las primeras 60 observaciones para el entrenamiento y probar el rendimiento usando las últimas 12 observaciones. La creación de particiones en R se puede hacer manualmente con la función window del paquete stats. Por ejemplo, dividamos la serie USgas en particiones, dejando las últimas 12 observaciones de la serie como la partición de pruebas y el resto como entrenamiento:
#Sección - 8 Código
——– Código 1 ——–
- Primero, carguemos la serie USgas desde el paquete TSstudio en el entorno:
- Podemos observar las principales características de la serie USgas con la función ts_info:
- Usemos la función window para dividir la serie en particiones de entrenamiento y pruebas:
- El resumen de la partición de entrenamiento y pruebas se puede ver en la siguiente salida:
- Ahora, echemos un vistazo a la partición de pruebas:
- Alternativamente, la función ts_split del paquete TSstudio proporciona una forma personalizada de crear particiones de entrenamiento y pruebas para datos de series temporales:
- Se puede observar a partir de la siguiente salida que recibimos los mismos resultados de ejecución que recibimos anteriormente:
- Ahora, echemos un vistazo a la partición de pruebas:
library(TSstudio)
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# Using the series time index to set the start and end point of each partiton
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 10La simplicidad de este método es su principal ventaja, ya que es bastante rápido entrenar y probar un modelo mientras se utiliza (relativamente) poca potencia informática. Por otro lado, no es posible llegar a una conclusión sobre la estabilidad y escalabilidad del rendimiento del modelo basado en una sola unidad de prueba. Es posible que un modelo, solo por casualidad, tenga un rendimiento relativamente bueno en el conjunto de pruebas pero se desempeñe mal en el pronóstico real ya que no es estable con el tiempo. Una forma de mitigar ese riesgo es con el enfoque de backtesting, que se basa en entrenar un modelo con múltiples particiones de entrenamiento y pruebas.
Pronóstico con backtesting
El enfoque de backtesting para entrenar un modelo de pronóstico es una versión avanzada del enfoque de única fuera de la muestra que vimos anteriormente. Se basa en el uso de una ventana deslizante para dividir la serie en múltiples pares de particiones de entrenamiento y pruebas. Un proceso básico de entrenamiento de backtesting incluye los siguientes pasos: 1. Preparación de datos: Crear múltiples pares de particiones de entrenamiento y pruebas. 2. Entrenar un modelo: Esto se hace en cada una de las particiones de entrenamiento. 3. Probar el modelo: Evaluar su rendimiento en las particiones de pruebas correspondientes. 4. Evaluar el modelo: Evaluar la precisión, escalabilidad y estabilidad del modelo en función de la puntuación de las pruebas. Basándose en la evaluación, se realizaría una de las siguientes acciones: • Generar el pronóstico final para verificar si la puntuación del modelo cumple con un umbral o criterio específico • Aplicar ajustes y optimización adicionales para el modelo y repetir los pasos de entrenamiento y evaluaciones.
El uso de metodología de puntuación nos permite evaluar la estabilidad del modelo examinando la tasa de error del modelo en los diferentes conjuntos de pruebas. Consideraríamos un modelo como estable siempre que la distribución de errores del modelo en los conjuntos de pruebas sea bastante estrecha. En este caso, la tasa de error del pronóstico real debería estar dentro del mismo rango que los conjuntos de pruebas (suponiendo que no haya eventos anormales que impacten la tasa de error del pronóstico).
Este método es, conceptualmente, similar al enfoque de validación cruzada para entrenar modelos de aprendizaje automático. La principal distinción entre los dos enfoques está relacionada con cómo se establecen sus particiones. Mientras que las particiones de backtesting están ordenadas cronológicamente, las del enfoque de validación cruzada se basan en muestreo aleatorio.
Los principales parámetros de configuración de un modelo de entrenamiento de backtesting son los siguientes: ⚫ La longitud de las particiones de entrenamiento: Se establece mediante la configuración de ventana. Hay dos tipos comunes de ventanas deslizantes: • Ventana en expansión: Usando 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 aumenta a medida que la ventana se desplaza sobre la serie. ⚫ Ventana deslizante: Configure las primeras N observaciones como la partición de entrenamiento inicial y desplace la ventana en n observaciones para crear la siguiente partición de entrenamiento. La longitud de las particiones de entrenamiento permanece igual mientras la ventana se desplaza sobre la serie: • La longitud de las particiones de pruebas: un valor constante, típicamente alineado con la longitud del horizonte de pronóstico (bajo la limitación del número mínimo de observaciones requeridas para entrenar el modelo) ⚫ El espacio entre cada partición de entrenamiento: define el ritmo de la ventana deslizante • El número de particiones de entrenamiento y pruebas
El siguiente diagrama muestra la estructura del backtesting con una ventana en expansión:
Descripción de la imagen
Como se puede ver en el diagrama anterior, todas las particiones de entrenamiento del método de ventana en expansión comienzan en el mismo punto de índice, To, donde cada partición de entrenamiento es una combinación de la partición de entrenamiento anterior con las siguientes n observaciones (donde n es una constante que representa la tasa de expansión de la ventana). Tendría sentido usar este método cuando la serie tenga un patrón estacional fuerte y una tendencia estable. En este caso, las primeras observaciones de la serie podrían, potencialmente, tener información predictiva que pueda ser utilizada por el modelo. La desventaja de este método es entrenar el modelo con particiones de longitud diferente, ya que típicamente un modelo tiende a aprender más y, por lo tanto, a rendir mejor cuando hay más datos disponibles. Por lo tanto, podemos observar que el rendimiento de los modelos entrenados en las particiones más recientes es mejor que el de los que se entrenan con las primeras particiones.
Tendría más sentido usar la ventana deslizante siempre que la serie de entrada tenga cambios estructurales o alta volatilidad, o cuando la mayor parte del poder predictivo esté vinculado a la historia más reciente (o alta correlación con los rezagos más recientes de la serie). Por ejemplo, como vimos en el capítulo anterior, los precios mensuales del petróleo crudo tienen una fuerte relación con los rezagos más recientes de la serie, y la historia lejana no contiene información predictiva poderosa sobre la dirección futura de la serie. Si bien el enfoque de backtesting nos proporciona información intuitiva sobre la estabilidad y escalabilidad del modelo, tiene un mayor costo computacional en comparación con el enfoque de partición única de entrenamiento y prueba. Por lo tanto, este compromiso entre los dos enfoques debe tenerse en cuenta al seleccionar el enfoque de entrenamiento.
Evaluación del pronóstico
El objetivo principal del paso de evaluación es evaluar la capacidad del modelo entrenado para pronosticar (o basado en otros criterios) con precisión las observaciones futuras de la serie. Este proceso incluye lo siguiente: ⚫ Análisis residual: Esto se centra en la calidad del modelo, con valores ajustados en la partición de entrenamiento. ⚫ Puntuación del pronóstico: Esto se basa en la capacidad del modelo para pronosticar los valores reales del conjunto de pruebas.
Análisis residual
El análisis residual prueba qué tan bien capturó e identificó el modelo los patrones de la serie. Además, proporciona información sobre las distribuciones de los residuos, que se requieren para construir intervalos de confianza para el pronóstico. La definición matemática de un residuo 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: Et = Yt - Ŷt Aquí, Et, Y y Ŷt representan el residuo, el valor real y los valores ajustados, respectivamente, en el tiempo t. Este proceso incluye el uso de herramientas de visualización de datos y pruebas estadísticas para evaluar lo siguiente: • Probar la bondad del ajuste frente a los valores reales: Haces esto trazando los valores de los residuos 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 de cero y con variación constante (ruido blanco) indican que el modelo puede capturar la mayoría 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 pudo capturar los patrones de la serie. Aquí hay algunas interpretaciones potenciales de la posible salida: • 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.
——– Codigo 2 ——–
# The sample.out argument set the size of the testing partition
# (and therefore the training partition)
USgas_partitions <- ts_split(USgas, sample.out = 12)
train <- USgas_partitions$train
test <- USgas_partitions$test
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 10——– Código 3 ——–
• Picos aleatorios: Esto es una indicación de posibles valores atípicos en la partición de entrenamiento. • Autocorrelación residual: Esto indica qué tan bien el modelo pudo 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 capturó. • Distribución de los residuos: Esto es necesario para concluir sobre la confiabilidad del intervalo de confianza del pronóstico. Si los residuos no siguen una distribución normal, no podemos usarlos para crear intervalos de confianza ya que se basa en la suposición de que los residuos siguen una distribución normal. • 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 te preocupes si no estás familiarizado con el modelo ARIMA: lo discutiremos en detalle en el Capítulo 11, Pronóstico con el Modelo ARIMA. Entrenaremos el modelo con la función auto.arima del paquete forecast:
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoomd <- auto.arima(train)——– Código 4 ——–
Usaremos la función checkresiduals del paquete forecast para examinar los residuos del modelo entrenado en la partición de entrenamiento. Esta función devuelve los siguientes cuatro resultados: ⚫ Gráfico de la serie temporal de los residuos ⚫ Gráfico de ACF de los residuos ⚫ Gráfico de distribución de los residuos ⚫ El resultado de la prueba de 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 residuos) es diferente de cero y utiliza las siguientes hipótesis: ⚫ Ho: El nivel de correlación entre la serie y su rezago es igual a cero, y por lo tanto las observaciones de la serie son independientes. ⚫ H1: 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.
——– Codigo 4 ——–
Vamos a usar la función checkresiduals para evaluar el rendimiento del modelo entrenado en la partición de entrenamiento:
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: 24Obtenemos los siguientes residuos del modelo ARIMA:
Comenzando con la salida de la prueba de Ljung-Box, notarás que, basándonos en los resultados del valor p, podemos rechazar la hipótesis nula con un nivel de significancia del 0.01. Por lo tanto, hay una indicación de que la correlación entre la serie residual y sus rezagos es diferente de cero. El gráfico de ACF proporciona un apoyo adicional para eso también. Esto indica que el modelo no capturó completamente todos los patrones de la serie, y es posible que desees modificar los parámetros de ajuste del modelo. El gráfico de serie temporal de los residuos oscila alrededor del eje x, con la excepción de algunos residuos, que cruzan el valor de ± 250. Esto podría indicar que ocurren algunos valores atípicos durante estos períodos, y deberías revisar esos puntos de datos en la serie (más adelante en este capítulo, veremos un método más intuitivo para comparar los valores ajustados con los valores reales utilizando la función test_forecast del paquete TSstudio). Por último, pero no menos importante, está el gráfico de distribución de los residuos, que parece ser una representación bastante buena de una distribución normal.
En un mundo ideal, deberías terminar este proceso con ruido blanco y residuos independientes. Sin embargo, en la realidad, en algunos casos, será más difícil lograr este objetivo debido a la estructura de la serie (valores atípicos, cambios estructurales, etc.), y es posible que tengas que seleccionar un modelo que te acerque más a este objetivo.
Evaluación del pronóstico
Una vez que finalices la sintonización del modelo, es hora de probar la capacidad del modelo para predecir observaciones que el modelo no vio antes (a diferencia de los valores ajustados que el modelo vio durante todo el proceso de entrenamiento). El método más común para evaluar el éxito del pronóstico para predecir los valores reales es utilizar métricas de precisión o error. El método más común para evaluar el éxito del pronóstico es predecir los valores reales con el uso de una métrica de error para cuantificar la precisión general del pronóstico. La selección de una métrica de error específica depende de los objetivos de precisión del pronóstico. Un ejemplo de métricas de error comunes son las siguientes:
• Error Cuadrático Medio (MSE): Esto cuantifica la distancia cuadrada promedio entre los valores reales y los pronosticados:
Descripción de la imagen
El efecto cuadrático del error evita que los valores positivos y negativos se cancelen entre sí y penaliza la puntuación del error a medida que aumenta la tasa de error. • Error Cuadrático Medio (MSE): Esto cuantifica la distancia cuadrada promedio entre los valores reales y los pronosticados:
• Raíz del Error Cuadrático Medio (RMSE): Esta es la raíz de la distancia cuadrada promedio entre los valores reales y los pronosticados:
Descripción de la imagen
Descripción de la imagen
Descripción de la imagen
Es más fácil comparar, establecer un punto de referencia o comunicarse con personas no técnicas debido a la representación en porcentaje.
——– Codigo 5 ——–
Por ejemplo, usemos el modelo que entrenamos anteriormente para pronosticar las 12 observaciones que dejamos para la prueba y evaluar su rendimiento. Utilizaremos la función forecast del paquete forecast para pronosticar los siguientes 12 meses (con respecto al punto final de la partición de entrenamiento):
fc <- forecast(md, h = 12)Ahora que hemos asignado el pronóstico al objeto fc, usaremos la función accuracy del paquete forecast para evaluar el rendimiento del modelo con respecto a los valores reales en la partición de prueba: # ——– Code 6 ——–
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.3404092——– Codigo 7 ——–
La función accuracy, que utilizaremos intensivamente en los próximos capítulos, devuelve varias métricas de error tanto para los valores ajustados (fila de Conjunto de Entrenamiento) como para el pronóstico real (fila de Conjunto de Prueba). Notarás que los resultados de MAPE del modelo son del 3.52% y del 7.84% en las particiones de entrenamiento y prueba, respectivamente. No debería sorprender una tasa de error más alta en la partición de prueba en comparación con la de entrenamiento, ya que típicamente el modelo vio los datos de la partición de entrenamiento durante todo el proceso de entrenamiento. Una tasa de error bastante baja en el conjunto de entrenamiento, junto con una alta tasa de error en el conjunto de prueba, es una clara indicación de sobreajuste en el modelo.
Un enfoque alternativo para evaluar el ajuste del modelo tanto en el conjunto de entrenamiento como en el de prueba es con la función test_forecast del paquete TSstudio. Esta función visualiza la serie real, los valores ajustados en la partición de entrenamiento y los valores pronosticados en el conjunto de prueba. Al pasar el cursor sobre los valores ajustados o pronosticados, aparece un cuadro de texto emergente con los resultados de RMSE y MAPE tanto en las particiones de entrenamiento como de prueba (que, lamentablemente, no se pueden transferir al libro en sí):
test_forecast(actual = USgas,
forecast.obj = fc,
test = test) Es más fácil y rápido identificar información sobre la bondad del ajuste tanto de los valores ajustados como de los pronosticados al trazar esos valores contra los valores reales de la serie. Por ejemplo, inmediatamente notarás que el pico residual durante 2006 es causado por valores atípicos (o un consumo más bajo que el patrón normal de la serie). Además, el pronóstico real no capturó el pico anual de 2018. Estas observaciones no se pueden obtener con métricas de error. Benchmark del pronóstico
Según las métricas de error, el modelo entrenado obtuvo un MAPE del 7.84% o un RMSE de 208.01. ¿Cómo podemos evaluar si estos resultados son demasiado altos o bajos? El método más común es comparar el rendimiento del modelo con algún pronóstico de referencia o algún método heredado que deseamos reemplazar. Un enfoque de referencia popular sería utilizar un enfoque de pronóstico simplista como línea base. Por ejemplo, pronostiquemos la serie con un enfoque ingenuo y lo usemos como referencia para el pronóstico anterior que creamos con el modelo ARIMA.
——– Code 8 ——–
Un enfoque ingenuo simple típicamente asume que el valor observado más recientemente es el verdadero representante del futuro. Por lo tanto, continuará con el último valor hasta el infinito (o hasta el horizonte del pronóstico). Podemos crear un pronóstico ingenuo con la función naive del paquete forecast y usar 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:
library(forecast)naive_model <- naive(train, h = 12)
test_forecast(actual = USgas,
forecast.obj = naive_model,
test = test)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.499679——– Codigo 9 ——–
En el caso del modelo ingenuo, no hay proceso de entrenamiento, y los valores ajustados se establecen como los valores reales (como puedes ver en el gráfico anterior). Dado que USgas tiene un patrón estacional fuerte, tendría sentido utilizar un modelo ingenuo estacional que tenga en cuenta la variación estacional. El modelo snaive del paquete forecast utiliza el último punto estacional como pronóstico de todas las observaciones estacionales correspondientes. Por ejemplo, si estamos utilizando una serie mensual, el valor del último enero en la serie se utilizará como el pronóstico puntual para todos los futuros meses de enero:
snaive_model <- snaive(train, h = 12)
test_forecast(actual = USgas,
forecast.obj = snaive_model,
test = test)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.4289964Parece que el modelo ingenuo estacional tiene un mejor ajuste para el tipo de serie que estamos pronosticando, es decir, USgas, debido a su fuerte patrón estacional (en comparación con el modelo ingenuo). Por lo 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.
——– Código 10 ——–
Finalizando el pronóstico Ahora que el modelo ha sido entrenado, probado, ajustado (si es necesario) y evaluado con éxito, podemos pasar 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 enfoques para utilizar la configuración de parámetros del modelo: ⚫ Si el modelo fue ajustado manualmente, debes utilizar los parámetros de ajuste exactos que se utilizaron en el modelo entrenado. ⚫ Si el modelo fue ajustado automáticamente por un algoritmo (como la función auto.arima que utilizamos anteriormente), puedes hacer cualquiera de las siguientes opciones: • Extraer la configuración de parámetros que se utilizó con la partición de entrenamiento. • Permitir que el algoritmo vuelva a ajustar los parámetros del modelo utilizando la serie completa, bajo la suposición de que el algoritmo tiene la capacidad de ajustar correctamente los parámetros del modelo al entrenar el modelo con nuevos datos.
Se recomienda el uso de algoritmos para automatizar el proceso de ajuste del modelo cuando se prueba la capacidad del modelo para ajustar el modelo con backtesting. Esto te permite revisar si el algoritmo tiene la capacidad de ajustar correctamente los parámetros del modelo, según los resultados del backtesting. Por razones de simplicidad, seguiremos utilizando el modelo auto.arima para entrenar el modelo final.
md_final <- auto.arima(USgas)
fc_final <- forecast(md_final, h = 12)——– Codigo 11 ——–
plot_forecast(fc_final,
title = "The US Natural Gas Consumption Forecast",
Xtitle = "Year",
Ytitle = "Billion Cubic Feet")——– Code 12 ——–
Manejo de la incertidumbre en el pronóstico
El objetivo principal del proceso de pronóstico, como vimos anteriormente, es minimizar el nivel de incertidumbre en torno a los valores futuros de la serie. Aunque no podemos eliminar completamente esta incertidumbre, podemos cuantificarla y proporcionar un rango alrededor de la estimación puntual del pronóstico (que no es más que el valor esperado del modelo de cada punto en el futuro). Esto se puede hacer utilizando el intervalo de confianza (o un intervalo creíble, cuando se utiliza el modelo bayesiano) o mediante simulación.
Intervalo de confianza
El intervalo de confianza es un método de aproximación estadística que se utiliza para expresar el rango de valores posibles que contienen el valor verdadero con cierto grado de confianza (o probabilidad). Hay dos parámetros que determinan el rango del intervalo de confianza:
- El nivel de confianza o la probabilidad de que el valor verdadero esté en ese rango. Cuanto mayor sea el nivel de confianza, mayor será el rango del intervalo.
- La desviación estándar estimada del pronóstico en el tiempo T+i, donde T representa la longitud de la serie e i representa el valor pronosticado i. Cuanto menor sea la tasa de error, más corto será el rango del intervalo de predicción.
Por defecto, la función forecast genera un intervalo de predicción con un nivel de confianza del 80% y 95%, pero puedes modificarlo utilizando el argumento level. Por ejemplo, usemos la función md_final del modelo entrenado y la función forecast para los próximos 60 meses utilizando el intervalo de predicción con niveles de confianza del 80% y 90%:
fc_final2 <- forecast(md_final,
h = 60,
level = c(80, 90))
plot_forecast(fc_final2,
title = "The US Natural Gas Consumption Forecast",
Xtitle = "Year",
Ytitle = "Billion Cubic Feet")——– Code 13 ——–
Simulación
Un enfoque alternativo es utilizar la distribución del modelo para simular posibles trayectorias para el pronóstico. Este método solo se puede usar cuando la distribución del modelo está disponible. La función forecast_sim del paquete TSstudio proporciona una función incorporada para simular posibles trayectorias de pronóstico. Esta estimación se puede utilizar para calcular la estimación puntual del pronóstico (por ejemplo, utilizando la media o la mediana de todas las trayectorias), o para calcular probabilidades de obtener diferentes valores. Alimentaremos el mismo modelo a la función y ejecutaremos 100 iteraciones:
fc_final3 <- forecast_sim(model = md_final,
h = 60,
n = 500) ——– Codigo 14 ——–
La salida de la función anterior contiene todas las simulaciones calculadas y las trayectorias simuladas. Extraigamos el gráfico de simulación (y usemos el paquete plotly para agregar títulos al gráfico):
library(plotly)## 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 = "US Natural Gas Consumption - Forecasting Simulation",
yaxis = list(title = "Billion Cubic Feet"),
xaxis = list(title = "Year"))#——– Código 15 ——– Enfoque de carrera de caballos Por último, pero no menos importante, terminaremos este capítulo con un enfoque de pronóstico robusto que combina lo que hemos aprendido hasta ahora en este capítulo. El enfoque de carrera de caballos se basa en entrenar, probar y evaluar múltiples modelos de pronóstico y seleccionar el modelo que tenga el mejor desempeño en las particiones de prueba. En el siguiente ejemplo, aplicaremos una carrera de caballos entre siete modelos diferentes (revisaremos los modelos en los próximos capítulos; por ahora, no te preocupes si no estás familiarizado con ellos) usando seis períodos de backtesting. La función ts_backtesting del paquete TSstudio realiza el proceso completo de entrenamiento, prueba, evaluación y luego pronóstico, utilizando el modelo que tuvo el mejor desempeño en las particiones de prueba del backtesting. De forma predeterminada, el modelo probará los siguientes modelos:
auto.arima: Modelo ARIMA automatizado bsts: Modelo bayesiano de series temporales estructurales ets: Modelo de espacio de estados de suavización exponencial híbrido: Un conjunto de múltiples modelos nnetar: Modelo de series temporales de red neuronal tbats: Modelo de espacio de estados de suavización exponencial, junto con transformación Box-Cox, tendencia, errores ARMA y componentes estacionales Holt Winters: Filtrado Holt-Winters Antes de ejecutar la función, establezcamos el valor de la semilla con la función set.seed para que podamos reproducir los resultados.
set.seed(1234)——– Code 16 ——–
{r}# USgas_forecast <- ts_backtesting(ts.obj = USgas, periods = 6, models = "abehntw", error = "MAPE", window_size = 12, h = 60, plot = FALSE)
——– Code 17 ——–
El modelo proporciona una tabla de clasificación (como podemos ver en la salida anterior) que está ordenada según el criterio de error establecido. En este caso, el modelo bsts tuvo la tasa de error más baja, por lo que la función recomendó que lo usemos (aunque todos los modelos y sus pronósticos están disponibles para extraer de la salida de la función). Podemos trazar la tasa de error y el pronóstico sugerido utilizando el modelo almacenado en el objeto de salida:
USgas_forecast$summary_plotEstos gráficos de error proporcionan una visión general informativa del rendimiento de cada modelo. Consideramos que un modelo es bueno si la distribución del error es estrecha y baja con respecto al resto de los modelos probados (como el modelo bsts).
#Sección - 9 # Código del Capítulo 9
#Predicción con regresión lineal El modelo de regresión lineal es uno de los métodos más comunes para identificar y cuantificar la relación entre una variable dependiente y una variable independiente única (regresión lineal univariada) o múltiple (regresión lineal multivariada). Este modelo tiene una amplia gama de aplicaciones, desde la inferencia causal hasta el análisis predictivo y, en particular, la predicción de series temporales. El objetivo de este capítulo son los métodos y enfoques para pronosticar datos de series de tiempo con regresión lineal. Esto incluye métodos para descomponer y pronosticar los componentes de la serie (por ejemplo, la tendencia y los patrones estacionales), manejar eventos especiales (como valores atípicos y días festivos) 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. • Serie de pronósticos con multiestacionalidad. Requerimiento técnico Los siguientes paquetes se utilizarán en este capítulo: •TSstudio: Versión 0.1.4 y superior •plotly: Versión 4.8 y superior •dplyr: Versión 0.8.1 y superior •lubridate: Versión 1.7.4 y superior •previsión: Versión 8.5 y superior
La regresión lineal 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: Y; Bo+B1X1,i+€¡ • Para n variables independientes, la ecuación tiene el siguiente aspecto: Y;= Bo+B1X1,i+B2 X2,i+…+3n Xn‚í + €i Las variables del modelo para estas ecuaciones son las siguientes: ⚫ i representa el índice de observaciones, i = 1,…, N Y representa la i observación de la variable dependiente X representa el valor i de la variable independiente j, donde j = 1,…, n Be representa el valor del término constante (o intercepto) ⚫B; representa los parámetros correspondientes (o coeficientes) de las j variables independientes & define el término de error, que no es más que toda la información que no fue capturada por las variables independientes para la observación i
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, Bo, B1, …, Bn), que puede formalizarse mediante las siguientes ecuaciones: • Para el modelo de regresión lineal univariado, la ecuación es la siguiente: Yi = Bo+B1X1, ⚫ Para el modelo de regresión lineal multivariante, la ecuación es la siguiente: Yi = ẞo + 1X1,¿ + B2X2,i+…+nXni Las variables para estas ecuaciones son las siguientes: i representa el índice de observación, i = 1,…, N ⚫ representa la estimación de la variable dependiente i observación X, representa el valor i de la variable independiente j, donde j = 1,…, n Bo representa la estimación del término constante (o intercepto) $1,…, son la estimación de los parámetros (o coeficientes) correspondientes de las n variables independientes La estimación de los coeficientes del modelo se basa en los dos pasos siguientes: • Definir una función de costo (también conocida como función de pérdida), estableciendo alguna métrica de error para minimizar • Aplicar optimización matemática para minimizar la función de costos. El método de estimación más común es aplicar los mínimos cuadrados ordinarios (OLS). norte método como técnica de optimización para minimizar la suma de cuadrados residuales(). 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
#Estimación de coeficientes con el método MCO 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:
Descripción de la imagen
#Predicción con regresión lineal El modelo de regresión lineal, a diferencia de los modelos tradicionales de series de tiempo como ARIMA o Holt-Winters, no fue diseñado explícitamente para manejar y pronosticar datos de series de tiempo. Más bien, es un modelo genérico con una amplia gama de aplicaciones, desde la inferencia causal hasta el análisis predictivo. Por tanto, la previsión con un modelo de regresión lineal se basa principalmente en los dos pasos siguientes: 1. Identificar la estructura de la serie, las características clave, los patrones, los valores atípicos y otras características. 2. Transformar esas características en variables de entrada y hacerlas retroceder con la serie para crear un modelo de pronóstico. Las características principales de un modelo de pronóstico de regresión lineal son los componentes de tendencia y estacional. La siguiente sección se centra en identificar la tendencia de la serie y los componentes estacionales y luego transformarlos en variables de entrada del modelo de regresión.
#Previsión de la tendencia y los componentes estacionales. 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.
Recordemos que la descomposición de una serie se puede definir mediante una de las siguientes expresiones: Y=T+S+C+4, cuando la serie tiene estructura aditiva Y = TX St Ct It, cuando la serie tiene estructura multiplicativa La explicación es la siguiente: ⚫ Tendencia: Representa el crecimiento de la serie en el tiempo después de ajustar y eliminar los efectos estacionales. ⚫ Estacional: un patrón cíclico recurrente que se deriva directamente de las unidades de frecuencia de la serie (por ejemplo, el mes del año para una serie con una frecuencia mensual) ⚫ Ciclo: Un patrón cíclico que no está relacionado con la unidad de frecuencia en serie. • Irregular: cualquier otro patrón que no sea capturado por los componentes de tendencia, estacional y cíclico. En aras de la simplicidad, podemos eliminar el componente de ciclo, ya que normalmente se fusiona con el componente de tendencia (o ignoramos el componente de ciclo). Por lo tanto, podemos actualizar y reemplazar las ecuaciones anteriores con las siguientes: . Y1 = T + St+It, cuando la serie tiene estructura aditiva YT St It, cuando la serie tiene estructura multiplicativa Ahora podemos transformar esas ecuaciones para un modelo de regresión lineal, modificando la notación de las ecuaciones: Y = Bo+B1T+ B2 St + €t Dónde: ⚫ Y representa una serie de tiempo con n observaciones ⚫ T, una variable independiente con n observaciones, representa el componente de tendencia de la serie ⚫ S, una variable independiente con n observaciones, representa el componente estacional de la serie. • ⚫ Et, el término de error de regresión, representa el componente irregular o cualquier patrón que no sea capturado por la tendencia de la serie y el componente estacional. • Bo, B1 y B2 representan el origen del modelo y los coeficientes de los componentes tendencial y estacional, respectivamente.
La transformación de una serie con estructura multiplicativa en una formación de regresión lineal requirió una transformación de la serie en una estructura aditiva. Esto se puede hacer aplicando la transformación logarítmica para ambos lados de las ecuaciones: log(Y) = log(T;) + log(St) + log(lt) Una vez que la serie se transforma en una estructura aditiva, la transformación a una formación de regresión lineal es sencilla y sigue el mismo proceso descrito anteriormente: log(Yt) = Bo+B1 log(Ti) + B2 log(St) +&t
——– Código 1 ——–
#Características de la ingeniería de los componentes de la serie. Antes de crear los datos de entrada de la 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")——– Codigo 2 ——–
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——– Codigo 3 ——–
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)——– Codigo 4 ——–
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)——– Codigo 5 ——–
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.1——– Codigo 6 ——–
Después de transformar la serie en un objeto data.frame, podemos comenzar a crear las funciones 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: USgas_df$tendencia <- 1:nrow (USgas_df) 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.
USgas_df$trend <- 1:nrow(USgas_df)——– Codigo 7 ——–
Después de transformar la serie en un objeto data.frame, podemos comenzar a crear las funciones 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.
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionUSgas_df$seasonal <- month(USgas_df$ds, label = T)——– Codigo 8 ——–
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 funciones:
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 jun——– Codigo 9 ——–
Modelando la serie tendencia y estacional. componentes Primero modelaremos la tendencia de la serie haciendo una regresión de la serie con la variable de tendencia, en la partición de entrenamiento:
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), ]——– Codigo 10 ——–
md_trend <- lm(y ~ trend, data = train)
summary(md_trend)##
## 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-09——– Codigo 11 ——–
Como puede verse en el resultado de la regresión anterior, el coeficiente de la variable de tendencia es estadísticamente significativo hasta un nivel de 0,001. Sin embargo, el R cuadrado ajustado de la regresión es bastante bajo, lo que generalmente tiene sentido, ya que la mayor parte de la variación de la serie está relacionada con el patrón estacional, como vimos en los gráficos anteriores. i Como puede observar en el resultado de la regresión anterior, la cuarta columna representa el valor p de cada uno de los coeficientes del modelo. El valor p proporciona la probabilidad de que rechacemos la hipótesis nula dado que en realidad es cierta, o el error tipo I. Por lo tanto, para el valor p menor que α, el valor umbral, rechazaremos la hipótesis nula con un nivel de significancia de α, donde los valores típicos de a son 0,1, 0,05, 0,01, etc. Como siempre, se recomienda poner algo de contexto a los números con visualización de datos. Por lo tanto, usaremos el modelo que creamos para predecir los valores ajustados en la partición de entrenamiento y los valores pronosticados en la partición de prueba. La función de predicción del paquete de estadísticas, como su nombre lo indica, predice los valores de unos datos de entrada en función de un modelo determinado. Lo usaremos para predecir los valores ajustados y pronosticados del modelo de tendencia que entrenamos antes:
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)——– Codigo 12 ——–
Crearemos una función de utilidad que traza la serie y el resultado del modelo, utilizando el paquete plotly:
library(plotly)
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)
}——– Codigo 13 ——–
Los argumentos de la función son los siguientes:
datos: Los datos de entrada, un dato. objeto frame que sigue la misma estructura que la de USgas_df (incluida la variable yhat) tren: el conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo t 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")——– Codigo 14 ——–
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.1299951——– Codigo 15 ——–
El 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:
md_seasonal <- lm(y ~ seasonal, data = train)
summary(md_seasonal)##
## 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) 2030.88 14.43 140.747 < 2e-16 ***
## seasonal.L -480.00 50.24 -9.554 < 2e-16 ***
## seasonal.Q 1103.83 50.17 22.000 < 2e-16 ***
## seasonal.C 72.42 50.05 1.447 0.149389
## seasonal^4 174.60 50.07 3.487 0.000592 ***
## seasonal^5 288.01 50.13 5.745 3.13e-08 ***
## seasonal^6 -44.67 50.09 -0.892 0.373460
## seasonal^7 -187.91 49.96 -3.762 0.000218 ***
## seasonal^8 84.95 49.84 1.705 0.089706 .
## seasonal^9 46.16 49.78 0.927 0.354828
## seasonal^10 77.55 49.76 1.559 0.120587
## seasonal^11 -11.09 49.75 -0.223 0.823856
## ---
## 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-16——– Codigo 16 ——–
Dado 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")——– Codigo 17 ——–
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.21502100——– Codigo 18 ——–
La 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:
md1 <- lm(y ~ seasonal + trend, data = train)
summary(md1)##
## 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) 1733.7153 17.0794 101.509 < 2e-16 ***
## seasonal.L -498.1709 29.6354 -16.810 < 2e-16 ***
## seasonal.Q 1115.2951 29.5872 37.695 < 2e-16 ***
## seasonal.C 78.9835 29.5103 2.676 0.00802 **
## seasonal^4 175.6505 29.5196 5.950 1.09e-08 ***
## seasonal^5 285.0192 29.5568 9.643 < 2e-16 ***
## seasonal^6 -49.3611 29.5319 -1.671 0.09610 .
## seasonal^7 -192.3050 29.4540 -6.529 4.77e-10 ***
## seasonal^8 81.8787 29.3835 2.787 0.00581 **
## seasonal^9 44.4849 29.3480 1.516 0.13106
## seasonal^10 76.8636 29.3372 2.620 0.00943 **
## seasonal^11 -11.2755 29.3353 -0.384 0.70109
## 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-16——– Codigo 19 ——–
La 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")——– Codigo 20 ——–
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.09143290——– Codigo 21 ——–
Hacer una regresión de la serie con los componentes de tendencia y estacional proporciona un aumento significativo tanto en la calidad del ajuste del modelo como en la precisión del modelo. Sin embargo, al observar el gráfico del ajuste y pronóstico del modelo, se puede observar que la tendencia del modelo es demasiado lineal y le falta la ruptura estructural de la tendencia de la serie. Este es el punto en el que agregar un componente polinomial al modelo podría proporcionar una mejora adicional en la precisión del modelo.
Una técnica sencilla para capturar una tendencia no lineal es agregar un componente polinómico a la tendencia de la serie para capturar la curvatura de la tendencia a lo largo del tiempo. Usaremos el argumento I, que nos permite aplicar operaciones matemáticas en cualquiera de los objetos de entrada. Por lo tanto, usaremos este argumento para agregar un segundo grado del polinomio para la entrada de tendencia:
md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)
summary(md2)##
## 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) 1.882e+03 2.199e+01 85.568 < 2e-16 ***
## seasonal.L -4.917e+02 2.530e+01 -19.438 < 2e-16 ***
## seasonal.Q 1.121e+03 2.525e+01 44.381 < 2e-16 ***
## seasonal.C 8.247e+01 2.518e+01 3.275 0.00123 **
## seasonal^4 1.763e+02 2.519e+01 6.999 3.35e-11 ***
## seasonal^5 2.835e+02 2.522e+01 11.243 < 2e-16 ***
## seasonal^6 -5.175e+01 2.520e+01 -2.054 0.04123 *
## seasonal^7 -1.946e+02 2.513e+01 -7.741 3.97e-13 ***
## seasonal^8 8.030e+01 2.507e+01 3.203 0.00157 **
## seasonal^9 4.362e+01 2.504e+01 1.742 0.08293 .
## seasonal^10 7.651e+01 2.503e+01 3.057 0.00253 **
## seasonal^11 -1.137e+01 2.503e+01 -0.454 0.65005
## 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-16——– Codigo 22 ——–
Agregar 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:
Además, 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%:
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")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.04212618——– Code Chank 23 ——–
#La función tslm
Hasta ahora, hemos visto el proceso manual de transformar un objeto ts a un formato de modelo de pronóstico de regresión lineal. La función ts1m del paquete de pronóstico proporciona una función incorporada para transformar un objeto ts en un modelo de pronóstico de regresión lineal. Con la función tslm, puede configurar el componente de regresión junto con otras funciones. Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la serie USgas con la función tslm utilizando la tendencia, el cuadrado de la tendencia y los componentes estacionales. Primero, dividamos la serie en particiones de entrenamiento y prueba usando la función ts_split:
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$test——– Codigo 24 ——–
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-16——– Codigo 25 ——–
Como 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 a un dato. ingeniería de características y objetos de marco • 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)
#Modelado de eventos únicos y eventos no estacionales 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.
Al no expresarse en el modelo de regresión, este tipo de eventos sesgará los coeficientes estimados, ya que el modelo ponderará ese tipo de eventos junto con los eventos regulares de la serie. El uso de variables de codificación activa, binarias o de bandera podría ayudar al modelo a ignorar este tipo de eventos o ajustar los coeficientes del modelo en consecuencia. Por ejemplo, se puede observar en el gráfico descompuesto de la serie de gas estadounidense mostrada anteriormente que la tendencia de la serie tuvo una ruptura estructural alrededor del año 2010. Si bien el crecimiento antes del año 2010 fue relativamente plano, la pendiente de la tendencia cambió después, con un crecimiento positivo. . En este caso, podemos utilizar una variable binaria que sea igual a cero para las observaciones anteriores al año 2010 y un año después. Hacer una regresión en el modelo slm con variables externas requiere un objeto data.frame separado con las variables correspondientes. El siguiente ejemplo demuestra el proceso de creación de una variable binaria externa que es igual a 0 antes del año 2010 y a 1 después, utilizando 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] <- 1md3 <- 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-16Como se puede observar en el resumen del modelo anterior, la variable ruptura estructural es estadísticamente significativa, con un nivel de 0,01. Del mismo modo, en el caso de valores atípicos o días festivos, se puede aplicar codificación activa estableciendo una variable binaria que sea igual a 1 siempre que ocurra un evento atípico o recurrente no estacional, y 0 en caso contrario. i Tenga en cuenta que, una vez que haya entrenado un modelo de pronóstico con tslm función con el uso de variables externas, tendrá que producir el valores futuros de esas variables, ya que se utilizarán como entrada del pronóstico.
#Previsión de una serie con componentes multiestacionales: un estudio de caso Una 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
La serie UKgrid representa la demanda de electricidad de la red nacional en el Reino Unido y está disponible en el paquete UKgrid. Esta serie representa datos de series temporales de alta frecuencia con frecuencia de media hora. Utilizaremos la función de extracción de cuadrícula del paquete UKgrid para definir la serie y las características principales (por ejemplo, formato de datos, variables, frecuencia, etc.). Esta función de transformación nos permite agregar la frecuencia de la serie desde cada media hora a una frecuencia más baja, como por hora, por día o por mes. Como nuestro objetivo aquí es pronosticar la demanda diaria en los próximos 365 días, configuraremos la serie en frecuencia diaria utilizando la estructura 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 puedes ver, esta serie tiene dos variables: • TIMESTAMP: un objeto de fecha utilizado como marca de tiempo o índice de la serie. ND: La demanda neta de electricidad Usaremos la función ts_plot para trazar y revisar la estructura de la serie:
ts_plot(UKdaily,
title = "The UK National Demand for Electricity",
Ytitle = "MW",
Xtitle = "Year")——– Codigo 26 ——–
Como 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 Se puede ver evidencia de esos patrones en el siguiente mapa de calor de la serie desde 2016 usando la función ts_heatmap del paquete TSstudio:
ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
title = "UK the Daily National Grid Demand Heatmap")——– Codigo 27 ——–
Como puede ver en el mapa de calor de la serie, la demanda general aumenta a lo largo de las semanas de invierno (por ejemplo, las semanas 1 a 12 del calendario y las semanas 44 a 52). Además, se puede observar el cambio de la serie durante los días de la semana, ya que la demanda aumenta durante los días laborables de la semana, y disminuye durante el fin de semana.
#Preprocesamiento e ingeniería de funciones de la serie UKdaily 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
Además, como es razonable suponer (confirmaremos esta suposición con la función ACF una vez que hayamos transformado la serie en un objeto ts) que la serie tiene una fuerte correlación con los rezagos estacionales, crearemos una variable de rezago con un retraso de 365 observaciones. Utilizaremos el paquete dplyr para crear esas funciones:
library(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)Recuerde que el costo de utilizar una variable de rezago con una longitud de N es la pérdida de las primeras N observaciones (ya que los rezagos de esas observaciones no pueden generarse a partir de la serie). Por lo tanto, utilizamos las funciones de filtro para eliminando las filas de la tabla en las que falta la variable lag365 (es decir, las primeras 365 observaciones).
Repasemos la estructura de la tabla UKdaily después de agregar esas nuevas características:
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 ...——– Code Chank 28 ——–
start_date <- min(UKdaily$TIMESTAMP)
start <- c(year(start_date),yday(start_date))UK_ts <- ts(UKdaily$ND,
start = c(year(start_date), yday(start_date)),
frequency = 365)```{r}# ts_acf(UK_ts, lag.max = 365 * 4)
# -------- Codigo 29 --------
El gráfico ACF anterior confirma nuestra suposición y la serie tiene una fuerte relación con los rezagos estacionales, en particular el rezago 365, el primer rezago.
Como nota al margen, puede estar seguro de que la serie también tiene una fuerte correlación con los desfases semanales (es decir, desfases 7, 14, 21, etc.). Sin embargo, en general, no se recomienda utilizar retrasos que sean menores que el horizonte de pronóstico (por ejemplo, retraso 7, cuando el horizonte de pronóstico es 365), ya que también tendrá que pronosticar esos retrasos para poder usarlos. como entrada en el modelo. Esto implica un esfuerzo adicional, ya que habrá que pronosticar esos retrasos. Además, puede aumentar el sesgo del pronóstico ya que utilizamos valores pronosticados como datos de entrada.
Ahora, después de haber creado las nuevas características para el modelo y configurar el objeto ts, estamos listos para dividir la serie de entrada y el objeto de características externas correspondiente que creamos (UKdaily) en una partición de entrenamiento y prueba. Como nuestro objetivo es pronosticar las próximas 365 observaciones y la longitud de la serie es lo suficientemente grande (2540 observaciones), podemos permitirnos establecer la partición de prueba a la longitud del horizonte de pronóstico de 365 observaciones. Estableceremos h, una variable indicadora, en 365 y la usaremos para definir las particiones y, más adelante, el horizonte de previsión:
Como antes, dividiremos la serie en particiones de entrenamiento y prueba con la función ts_split:
De manera similar, tenemos que dividir las características que creamos para el modelo de regresión (las características estacionales y de retraso) en una partición de entrenamiento y prueba siguiendo exactamente el mismo orden que usamos para el objeto ts correspondiente. Usaremos los datos. Funcionalidad de índice de cuadros para configurar la tabla UKdaily para entrenar y probar particiones:
```r
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), ]——– Codigo 30 ——–
Entrenamiento y prueba del modelo de pronóstico. 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 referencia: Regresión de la serie con el componente estacional y de tendencia utilizando las funciones integradas de la función ts1m. 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 usando la puntuación MAPE • Visualización de los valores ajustados y pronosticados versus los valores reales de la serie usando la función de pronóstico de prueba
Comenzar con un modelo simplista (modelo de referencia) nos permitirá observar si agregar nuevas características contribuye al rendimiento del modelo o si debemos evitarlo, ya que agregar más características o complejidad al modelo no siempre produce mejores resultados. Comenzaremos con el modelo base, haciendo una regresión de la serie con sus componentes estacionales y de tendencia:
A continuación, utilizaremos el modelo entrenado, md_ts1m1, para pronosticar los próximos 365 días de la serie, correspondientes a las observaciones de la partición de prueba, usando la función de pronóstico:
Comparemos el rendimiento del modelo en los conjuntos de entrenamiento y prueba 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)Podemos observar en el gráfico de rendimiento anterior que el modelo de referencia está haciendo un gran trabajo al capturar tanto la tendencia de la serie como la estacionalidad del día del año. Por otro lado, no logra capturar la oscilación relacionada con el día de la semana. La puntuación MAPE del modelo, como podemos ver en el resultado de la siguiente función de precisión, es 6,09 % y 7,77 % en las particiones de entrenamiento y prueba, respectivamente:
accuracy(fc_tslm1, 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——– Codigo 31 ——–
md_tslm2 <- tslm(train_ts ~ season + trend + wday, data = train_df)Como ahora utilizamos funciones de una fuente de datos externa, tenemos que especificar los datos de entrada con el argumento de datos. 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)Al observar el gráfico de rendimiento anterior del segundo modelo, podemos observar la contribución de las características estacionales en el pronóstico, ya que el segundo modelo fue capaz de capturar tanto la tendencia como los patrones multiestacionales de la serie. Esto también se puede observar en la puntuación del modelo MAPE, que cayó al 2,87 % y al 5,23 % en las particiones de entrenamiento y prueba, respectivamente:
accuracy(fc_tslm2, 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.68445——– Codigo 32 ——–
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)