Previsiones con modelos ARIMA
El ARIMA es un tipo de modelo de previsión que se basa en los procesos AR y MA. Es considerado como el enfoque más avanzado y sólido en comparación con otros modelos tradicionales como la regresión lineal y el alisamiento exponencial. En este capítulo, se describirán los elementos fundamentales del modelo, incluyendo los procesos AR y MA, así como la componente de diferenciación. Además, se presentarán los métodos y técnicas para ajustar los parámetros del modelo a través del uso de la diferenciación, el autocorrelador y las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF).
El proceso estacionario
Una de las principales suposiciones de la familia de modelos ARIMA es que la serie de entrada sigue una estructura de proceso estacionario. Esta suposición se basa en el teorema de representación de Wold, el cual establece que cualquier proceso estacionario puede ser representado como una combinación lineal de ruido blanco. Por lo tanto, antes de profundizar en los componentes del modelo ARIMA, detengámonos un momento y hablemos sobre el proceso estacionario. El proceso estacionario, en el contexto de los datos de series temporales, describe un estado estocástico de la serie.
Por ejemplo, en el siguiente ejemplo, simularemos un proceso AR con un retardo (es decir,p = 1) y 500 observaciones con la función arima. sim. Antes de ejecutar la simulación estableceremos el valor de la semilla en 12345:
set.seed(12345)
stationary_ts <- arima.sim(model = list(order = c(1,0,0),
ar = 0.5 ),
n = 500)Ahora, vamos a trazar la serie temporal simulada con la función ts_plot:
library(TSstudio)
ts_plot(stationary_ts,
title = "Stationary Time Series",
Ytitle = "Value",
Xtitle = "Index")En este caso, puede ver que, en general, la media de la serie, a lo largo del tiempo, se mantiene en torno a la línea cero. Además, la varianza de la serie no cambia con el tiempo. Utilicemos la función arima.sim para crear un ejemplo de series no estacionarias:
set.seed(12345)
non_stationary_ts <- arima.sim(model = list(order = c(1,1,0),ar = 0.3 ),n = 500)
ts_plot(non_stationary_ts,
title = "Non-Stationary Time Series",
Ytitle = "Value",
Xtitle = "Index")La serie clásica de pasajeros aéreos (el número mensual de pasajeros de líneas aéreas entre 1949 y 1960) de la paquete de datos es un buen ejemplo de una serie que viola las dos condiciones del proceso estacionario. Dado que la serie tiene tanto una fuerte tendencia lineal como un componente estacional multiplicativo, la media y la varianza cambian con el tiempo:
data(AirPassengers)
ts_plot(AirPassengers,
title = "Monthly Airline Passenger Numbers 1949-1960",
Ytitle = "Thousands of Passengers",
Xtitle = "Year")Transformando series no estacionarias en series estacionarias
En la mayoría de los casos, a menos que tengas mucha suerte, tus datos en bruto probablemente vendrán con una tendencia u otra forma de oscilación que viola las suposiciones del proceso estacionario. Por lo tanto, para manejar esto, tendrás que aplicar algunas transformaciones para llevar la serie a un estado estacionario. Los métodos de transformación comunes son la diferenciación de la serie (o la eliminación de la tendencia) y la transformación logarítmica (o ambas). Revisemos las aplicaciones de estos métodos.
Diferenciación de series temporales
El enfoque más común para lograr que los datos de una serie temporal no estacionaria se transformen en un estado estacionario es a través de la diferenciación de la serie con respecto a sus retardos. La principal ventaja de diferenciar una serie es que se elimina la tendencia (o se “desdiferencia” la serie), lo que ayuda a estabilizar la media de la serie. El nivel o orden de diferenciación de la serie se mide por la cantidad de veces que se aplica la diferenciación con respecto a sus retardos.
Obtenemos el siguiente resultado:
ts_plot(diff(AirPassengers, lag = 1),
title = "AirPassengers Series - First Differencing",
Xtitle = "Year",
Ytitle = "Differencing of Thousands of Passengers")Se puede observar que la primera diferencia de la serie AirPassenger eliminó la tendencia de la serie y que la media de la serie se mantiene constante a lo largo del tiempo. Por otro lado, hay evidencia clara de que la variación de la serie está aumentando con el tiempo, lo que indica que la serie aún no es estacionaria. Además de la diferencia de primer orden, tomar la diferencia estacional de la serie podría resolver este problema. Añadamos la diferencia estacional a la diferencia de primer orden y tracemos la serie nuevamente:
ts_plot(diff(diff(AirPassengers, lag = 1), 12),
title = "AirPassengers Series - First and Seasonal Differencing",
Xtitle = "Year",
Ytitle = "Differencing of Thousands of Passengers")Transformación logarítmica
En el capítulo 5, Descomposición de datos de series temporales, vimos las aplicaciones de la transformación logarítmica para manejar datos de series temporales multiplicativas. De manera similar, podemos utilizar este enfoque para estabilizar una oscilación estacional multiplicativa, si existe. Este enfoque no reemplaza la diferenciación, sino que se añade a ella. Por ejemplo, en el ejemplo de AirPassenger en la sección anterior, vimos que la primera diferenciación está haciendo un gran trabajo en estabilizar la media de la serie, pero no es suficiente para estabilizar la varianza de la serie. Por lo tanto, podemos aplicar una transformación logarítmica para transformar la estructura estacional de multiplicativa a aditiva y luego aplicar la diferencia de primer orden para estacionarizar la serie:
ts_plot(diff(log(AirPassengers), lag = 1),
title = "AirPassengers Series - First Differencing with Log Transformation",
Xtitle = "Year",
Ytitle = "Differencing/Log of Thousands of Passengers")La combinación de la transformación logarítmica y la diferenciación de primer orden logra una mejor estacionarización de la serie en comparación con el método anterior de doble diferenciación (primera diferencia y luego diferencia estacional).
El proceso del camino aleatorio
El paseo aleatorio, en el contexto de las series temporales, describe un proceso estocástico de un objeto en el tiempo.
En el siguiente código se simulan 20 caminatas aleatorias distintas con 500 pasos cada una, todas partiendo desde el mismo punto de inicio. Se grafican las trayectorias de estas caminatas, así como su primera diferencia, utilizando el paquete plotly.
Obtenemos el siguiente resultado.
set.seed(12345)
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':
##
## layoutp1 <- plot_ly()
p2 <- plot_ly()
for(i in 1:20){
rm <- NULL
rw <- arima.sim(model = list(order = c(0, 1, 0)), n = 500)
p1 <- p1 %>% add_lines(x = time(rw), y = as.numeric(rw))
p2 <- p2 %>% add_lines(x = time(diff(rw)), y = as.numeric(diff(rw)))
}Aquí, p1 representa el gráfico de la simulación del paseo aleatorio:
p1 %>% layout(title = "Simulate Random Walk",
yaxis = list(title = "Value"),
xaxis = list(title = "Index")) %>%
hide_legend()Aquí, p2 representa el gráfico correspondiente de la diferenciación de primer orden del azar:
p2 %>% layout(title = "Simulate Random Walk with First-Order Differencing",
yaxis = list(title = "Value"),
xaxis = list(title = "Index")) %>%
hide_legend()El proceso AR
El proceso AR define el valor actual de la serie Y, como una combinación lineal de los p rezagos anteriores de la serie
set.seed(12345)
ar2 <- arima.sim(model = list(order = c(2,0,0),
ar = c(0.9, -0.3)),
n = 500) Revisemos nuevamente la simulación de series temporales.
ts_plot(ar2,
title = "Simulate AR(2) Series",
Ytitle = "Value",
Xtitle = "Index")La función ar del paquete stats nos permite ajustar un modelo AR a los datos de una serie temporal y pronosticar sus valores futuros. Esta función identifica automáticamente el orden AR en función del criterio de información de Akaike (AIC). El argumento “method” permite definir el método de estimación de los coeficientes, como el método de mínimos cuadrados ordinarios (OLS) (que vimos en el Capítulo 9, Pronóstico con Regresión Lineal), la estimación de máxima verosimilitud (MLE) y Yule-Walker (predeterminado). Apliquemos la función ar para identificar el orden AR y estimar sus coeficientes correspondientemente.
md_ar <- ar(ar2)
md_ar##
## Call:
## ar(x = ar2)
##
## Coefficients:
## 1 2
## 0.8851 -0.2900
##
## Order selected 2 sigma^2 estimated as 1.049El uso de la función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF), que se presentaron en el capítulo 7, Análisis de correlación, nos permite clasificar el tipo de proceso y identificar su orden. Si la salida de la ACF disminuye gradualmente y la salida de la PACF se corta en el rezago p, esto indica que la serie es un proceso AR(p). Calculemos y grafiquemos la ACF y la PACF para la serie AR(2) simulada que creamos anteriormente con las funciones ACF y PACF. Primero, usaremos la función par para graficar ambos gráficos uno al lado del otro, estableciendo el argumento mfrow en c (1,2) (una fila, dos columnas):
par(mfrow=c(1,2))
acf(ar2)
pacf(ar2)
El proceso de medias móviles
En algunas situaciones, el modelo de pronóstico no puede capturar todos los patrones de la serie, por lo que alguna información queda en los residuos del modelo (o error de pronóstico) e. El objetivo del proceso de media móvil es capturar patrones en los residuos, si existen, mediante la modelización de la relación entre Y, el término de error, e, y los g términos de error pasados del modelo (por ejemplo, © 1 ¥1–2+ ++ +1 £i—4). La estructura del proceso de media móvil es bastante similar a la de los procesos AR. La siguiente ecuación define un proceso MA con un orden g:
set.seed(12345)
ma2 <- arima.sim(model = list(order = c(0, 0, 2),
ma = c(0.5, -0.3)),
n = 500) De forma similar, ya que hemos simulado la serie AR(2) en la sección anterior, utilizaremos la función arima. sim para simular una serie con una estructura MA(2).
ts_plot(ma2,
title = "Simulate MA(2) Series",
Ytitle = "Value",
Xtitle = "Index")md_ma <- arima(ma2, order = c(0,0,2), method = "ML")
md_ma##
## Call:
## arima(x = ma2, order = c(0, 0, 2), method = "ML")
##
## Coefficients:
## ma1 ma2 intercept
## 0.530 -0.3454 0.0875
## s.e. 0.041 0.0406 0.0525
##
## sigma^2 estimated as 0.9829: log likelihood = -705.81, aic = 1419.62Identificar el proceso de MA y sus características
Podemos identificar un proceso de media móvil y su orden utilizando las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF), de manera similar al proceso AR. Si el ACF se corta en el retraso g y la función PACF tiene una cola, podemos concluir que el proceso es una MA(g). Realicemos el mismo proceso que aplicamos en la serie ar2 con la serie ma2:
par(mfrow=c(1,2))
acf(ma2)
pacf(ma2)
####En la serie ma2, si observamos la trama de ACF, se corta en el segundo rezago, y si observamos la trama de PACF, disminuye. Basándonos en esto, podemos concluir que la serie ma2 es un proceso MAA(2).
El modelo ARMA
Hasta ahora, hemos aprendido cómo aplicar los modelos AR y MA por separado. Pero en algunos casos, es útil combinar ambos para manejar datos de series de tiempo más complicados. El modelo ARMA es una combinación de los modelos AR(p) y MA(q).
set.seed(12345)
arma <- arima.sim(model = list(order(1,0,2), ar = c(0.7), ma = c(0.5,-0.3)), n = 500)ts_plot(arma,
title = "Simulate ARMA(1,2) Series",
Ytitle = "Value",
Xtitle = "Index")La función arima permite ajustar un modelo ARMA de manera sencilla. Para hacerlo, debemos definir los parámetros p y q en el argumento de orden, tal como lo hacemos con la simulación. Una vez ajustado el modelo, podemos observar que los valores de los coeficientes son similares a los que obtuvimos en la simulación.
arma_md <- arima(arma, order = c(1,0,2))
arma_md##
## Call:
## arima(x = arma, order = c(1, 0, 2))
##
## Coefficients:
## ar1 ma1 ma2 intercept
## 0.7439 0.4785 -0.3954 0.2853
## s.e. 0.0657 0.0878 0.0849 0.1891
##
## sigma^2 estimated as 1.01: log likelihood = -713.18, aic = 1436.36Identificación de una función ARMA
La identificación del proceso ARMA sigue el mismo planteamiento que utilizamos anteriormente con los procesos AR y MA. Existe un proceso ARMA en los datos de series temporales si tanto los gráficos ACF como PACE llegan a la cola, como podemos ver en el siguiente ejemplo:
par(mfrow=c(1,2))
acf(arma)
pacf(arma)
Ajuste manual del proceso ARMA
El ajuste manual del modelo ARMA se basa principalmente en la experimentación, la intuición, el sentido común y la experiencia.
arima(arma, order = c(1, 0, 1), include.mean = FALSE)##
## Call:
## arima(x = arma, order = c(1, 0, 1), include.mean = FALSE)
##
## Coefficients:
## ar1 ma1
## 0.4144 0.8864
## s.e. 0.0432 0.0248
##
## sigma^2 estimated as 1.051: log likelihood = -723, aic = 1452arima(arma, order = c(2, 0, 1), include.mean = FALSE)##
## Call:
## arima(x = arma, order = c(2, 0, 1), include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ma1
## 0.3136 0.1918 0.9227
## s.e. 0.0486 0.0484 0.0183
##
## sigma^2 estimated as 1.019: log likelihood = -715.26, aic = 1438.52arima(arma, order = c(1, 0, 2), include.mean = FALSE)##
## Call:
## arima(x = arma, order = c(1, 0, 2), include.mean = FALSE)
##
## Coefficients:
## ar1 ma1 ma2
## 0.7602 0.4654 -0.4079
## s.e. 0.0626 0.0858 0.0832
##
## sigma^2 estimated as 1.014: log likelihood = -714.27, aic = 1436.54arima(arma, order = c(2, 0, 2), include.mean = FALSE)##
## Call:
## arima(x = arma, order = c(2, 0, 2), include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.7239 0.0194 0.4997 -0.3783
## s.e. 0.2458 0.1257 0.2427 0.2134
##
## sigma^2 estimated as 1.014: log likelihood = -714.26, aic = 1438.51La salida de “checkresiduals” muestra que los residuos del modelo ARMAC(1,2) cumplen con las suposiciones del modelo que definimos previamente.
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoobest_arma <- arima(arma, order = c(1, 0, 2), include.mean = FALSE)
checkresiduals(best_arma)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2) with zero mean
## Q* = 5.3129, df = 7, p-value = 0.6218
##
## Model df: 3. Total lags used: 10Previsión de modelos AR, MA y ARMA
Hacer predicciones con cualquiera de los modelos que hemos visto hasta ahora fue fácil: utilizamos la función “forecast” del paquete “forecast” de forma similar a cómo lo hicimos en el capítulo anterior. Por ejemplo, el siguiente código muestra cómo predecir las próximas 100 observaciones del modelo AR que entrenamos anteriormente en la sección “El proceso AR”.
ar_fc <- forecast(md_ar, h = 100)plot_forecast(ar_fc,
title = "Forecast AR(2) Model",
Ytitle = "Value",
Xtitle = "Year")El proceso ARIMA
Las limitaciones de los modelos AR, MA y ARMA radican en que no pueden manejar datos de series de tiempo no estacionarias. Para resolver este problema, se requiere un paso de preprocesamiento para transformar la serie de datos de un estado no estacionario a un estado estacionario. El modelo ARIMA proporciona una solución a este problema al agregar un proceso integrado para el modelo ARMA. El proceso integrado (I) es simplemente la diferencia entre la serie y sus lags, donde el grado de diferencia está representado por el parámetro d. La diferencia, como vimos anteriormente, es una de las formas en que se pueden transformar los métodos de una serie de no estacionarios a estacionarios. Por ejemplo, Y, - Y,, representa la primera diferencia de la serie, mientras que (Y;- Y,,) - (Y,, - Y,,) representa la segunda diferencia.
Determinación del grado de diferenciación del modelo
De manera similar a como se establecen los parámetros p y q, el parámetro d (que representa el grado de diferenciación de la serie) puede ser fijado observando los gráficos ACF y PACF. En el siguiente ejemplo, se utilizarán los precios mensuales del café Robusta desde el año 2000. Estos precios forman parte del objeto de series de tiempo múltiples llamado “Coffee_Prices” que se encuentra en el paquete “TSstudio”. Comenzaremos cargando la serie de precios del café y realizando la diferencia entre los precios mensuales de Robusta desde enero del 2010 utilizando la función “window”.
data(Coffee_Prices)
robusta_price <- window(Coffee_Prices[,1], start = c(2000, 1))ts_plot(robusta_price,
title = "The Robusta Coffee Monthly Prices",
Ytitle = "Price in USD",
Xtitle = "Year") Como puede observarse, los precios del café Robusta a lo largo del tiempo muestran una tendencia al alza, lo que indica que no está en un estado estacionario. Además, dado que esta serie representa precios continuos, es probable que tenga una fuerte correlación con sus retrasos pasados, ya que los cambios en el precio suelen estar cerca del precio anterior. Utilizaremos de nuevo la función act para identificar el tipo de relación entre la serie y sus retrasos.
obtenemos el siguiente resultado:
acf(robusta_price)
####La correlación de la serie con sus rezagos se va reduciendo gradualmente a lo largo del tiempo en un patrón lineal, lo que se puede apreciar en la salida del gráfico ACF. Para eliminar la tendencia y la correlación entre la serie y sus rezagos, se procederá a diferenciar la serie. Se empezará con la primera diferenciación utilizando la función diff.
robusta_price_d1 <- diff(robusta_price)par(mfrow=c(1,2))
acf(robusta_price_d1)
pacf(robusta_price_d1)
####El gráfico ACF y PACE de la primera diferencia de la serie indica que es apropiado utilizar un proceso AR(1) en la serie diferenciada ya que el ACF se está reduciendo gradualmente y el PACF se corta en el primer retraso. Por lo tanto, aplicaremos un modelo ARIMA(1,1,0) en la serie “robusta_price” para incluir la primera diferencia.
robusta_md <- arima(robusta_price, order = c(1, 1, 0))summary(robusta_md)##
## Call:
## arima(x = robusta_price, order = c(1, 1, 0))
##
## Coefficients:
## ar1
## 0.2780
## s.e. 0.0647
##
## sigma^2 estimated as 0.007142: log likelihood = 231.38, aic = -458.76
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.002595604 0.08432096 0.06494772 0.08104715 4.254984 1.001542
## ACF1
## Training set 0.001526295Con base a los resultados de los anteriores codigos, obtenemos la siguiente grafica:
checkresiduals(robusta_md)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)
## Q* = 26.896, df = 23, p-value = 0.2604
##
## Model df: 1. Total lags used: 24En general, la trama de los residuos del modelo y la prueba de Ljung-Box indican que los residuos son ruido blanco. La trama ACF indica que hay algunos rezagos correlacionados, pero están en el borde de ser significativos, por lo que podemos ignorarlos.
El modelo ARIMA estacional
El modelo ARIMA estacional (SARIMA), como su nombre indica, es una versión designada del modelo ARIMA para series temporales con un componente estacional. Como vimos en el Capítulo 6, Análisis de la estacionalidad, y en el Capítulo 7, Análisis de la correlación, una serie temporal con un componente estacional tiene una fuerte relación con sus desfases estacionales. El modelo SARIMA utiliza los desfases estacionales de manera similar a cómo el modelo ARIMA utiliza los desfases no estacionales con los procesos AR y MA y el diferenciado.
Previsión del consumo mensual de gas natural en EE.UU. con el modelo SARIMA - un estudio de caso
En esta sección, aplicaremos lo aprendido en todo este capítulo y pronosticaremos el consumo mensual de gas natural en los Estados Unidos utilizando el modelo SARIMA.
data(USgas)ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")En este caso, como vimos en los capítulos anteriores, la serie Usgas tiene un patrón estacional fuerte y, por lo tanto, entre los modelos de la familia ARIMA, el modelo SARIMA es el más apropiado para utilizar. Además, como la serie muestra una tendencia al alza, ya podemos concluir que no es estacionaria y se requiere cierto grado de diferenciación de la serie.
USgas_split <- ts_split(USgas, sample.out = 12)
train <- USgas_split$train
test <- USgas_split$testpar(mfrow=c(1,2))
acf(train, lag.max = 60)
pacf(train, lag.max = 60)
El gráfico de ACF anterior indica que la serie tiene una fuerte correlación con tanto los desfases estacionales como no estacionales. Además, la disminución lineal de los desfases estacionales indica que la serie no es estacionaria y que se requiere una diferenciación estacional. Comenzaremos con una diferenciación estacional de la serie y graficaremos la salida para identificar si la serie está en un estado estacionario.
USgas_d12 <- diff(train, 12)
ts_plot(USgas_d12,
title = "US Monthly Natural Gas consumption - First Seasonal Difference",
Ytitle = "Billion Cubic Feet (First Difference)",
Xtitle = "Year")Aunque hemos eliminado la tendencia de la serie, la variación de la serie aún no es estable. Por lo tanto, intentaremos también tratar de tomar la primera diferencia de la serie:
USgas_d12_1 <- diff(diff(USgas_d12, 1))
ts_plot(USgas_d12_1,
title = "US Monthly Natural Gas consumption - First Seasonal and Non-Seasonal Differencing",
Ytitle = "Billion Cubic Feet (Difference)",
Xtitle = "Year")Después de aplicar una diferencia de orden 1 y una diferencia estacional de orden 1, la serie parece estabilizarse alrededor de la línea del eje x en cero (o cercano a una estabilidad razonable). Una vez que se ha transformado la serie en un estado estacionario, podemos volver a revisar las funciones ACF y PACF para identificar el proceso requerido.
Obtenemos la siguiente grafica:
par(mfrow=c(1,2))
acf(USgas_d12_1, lag.max = 60)
pacf(USgas_d12_1, lag.max = 60)
p <- q <- P <- Q <- 0:2arima_grid <- expand.grid(p,q,P,Q)
names(arima_grid) <- c("p", "q", "P", "Q")
arima_grid$d <- 1
arima_grid$D <- 1arima_grid$k <- rowSums(arima_grid)USgas_best_md <- arima(train, order = c(1,1,1), seasonal = list(order = c(2,1,1)))
USgas_best_md##
## Call:
## arima(x = train, order = c(1, 1, 1), seasonal = list(order = c(2, 1, 1)))
##
## Coefficients:
## ar1 ma1 sar1 sar2 sma1
## 0.4247 -0.9180 0.0132 -0.2639 -0.7449
## s.e. 0.0770 0.0376 0.0894 0.0834 0.0753
##
## sigma^2 estimated as 10160: log likelihood = -1292.96, aic = 2597.91USgas_test_fc <- forecast(USgas_best_md, h = 12)accuracy(USgas_test_fc, test)## ME RMSE MAE MPE MAPE MASE
## Training set 6.081099 97.85701 73.36854 0.1298714 3.517097 0.6371821
## Test set 42.211254 104.79281 83.09943 1.4913412 3.314280 0.7216918
## ACF1 Theil's U
## Training set 0.004565601 NA
## Test set -0.049999869 0.3469228test_forecast(USgas,
forecast.obj = USgas_test_fc,
test = test)En este fragmento se muestra que el modelo SARIMA logra capturar con éxito el patrón estacional y de tendencia de la serie. Sin embargo, el modelo encuentra dificultades para capturar los picos estacionales (mes de enero) en la partición de entrenamiento y tiene un error absoluto del 6,7% para el mes de enero (pico anual) en la partición de prueba. Esto se debe a la alta fluctuación durante el invierno. Podemos manejar esta incertidumbre del modelo durante los momentos de mayor demanda con los intervalos de confianza del modelo o las simulaciones de ruta que vimos en el capítulo 8, Estrategias de Pronóstico.
final_md <- arima(USgas, order = c(1,1,1), seasonal = list(order = c(2,1,1)))Antes de hacer una predicción para los próximos 12 meses, debemos verificar que los residuos del modelo cumplan con la condición del modelo.
checkresiduals(final_md)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(2,1,1)[12]
## Q* = 30.173, df = 19, p-value = 0.04964
##
## Model df: 5. Total lags used: 24El resultado de la prueba de Ljung-Box sugiere que los residuos del modelo son ruido blanco.
Al observar el gráfico de residuos anterior, se puede ver que los residuos son ruido blanco y están distribuidos normalmente. Además, la prueba de Ljung-Box confirma que no queda autocorrelación en los residuos: con un valor de p de 0,12, no podemos rechazar la hipótesis nula de que los residuos son ruido blanco. ¡Estamos listos para seguir adelante! Utilicemos la función de pronóstico para pronosticar los próximos 12 meses de la serie Usgas.
USgas_fc <- forecast(final_md, h = 12)Y al final obtenemos el siguiente resultado:
plot_forecast(USgas_fc,
title = "US Natural Gas Consumption - Forecast",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")La función auto.arima
Una de las principales dificultades en la realización de pronósticos con la familia de modelos ARIMA es el engorroso proceso de ajuste de los modelos. Como se vio en este capítulo, este proceso incluye numerosos pasos manuales para verificar la estructura de la serie (estacionaria o no estacionaria), transformaciones de datos, análisis descriptivo con los gráficos ACF y PACF para identificar el tipo de proceso y ajustar finalmente los parámetros del modelo. Si bien puede tomar unos pocos minutos entrenar un modelo ARIMA para una sola serie, esto puede no ser escalable si se tienen docenas de series para pronosticar. La función “auto.arima” del paquete “forecast” proporciona una solución a este problema. Este algoritmo automatiza el proceso de ajuste del modelo ARIMA con el uso de métodos estadísticos para identificar tanto la estructura de la serie (estacionaria o no estacionaria) como el tipo (estacional o no) y ajusta los parámetros del modelo en consecuencia.
USgas_auto_md1 <- auto.arima(train)
USgas_auto_md1## Series: train
## ARIMA(2,1,1)(2,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 sma1
## 0.4301 -0.0372 -0.9098 0.0118 -0.2673 -0.7431
## s.e. 0.0790 0.0752 0.0452 0.0893 0.0829 0.0751
##
## sigma^2 = 10446: log likelihood = -1292.83
## AIC=2599.67 AICc=2600.22 BIC=2623.2USgas_auto_md2 <- auto.arima(train,
max.order = 5,
D = 1,
d = 1,
ic = "aic",
stepwise = FALSE,
approximation = FALSE)
USgas_auto_md2## Series: train
## ARIMA(1,1,1)(2,1,1)[12]
##
## Coefficients:
## ar1 ma1 sar1 sar2 sma1
## 0.4247 -0.9180 0.0132 -0.2639 -0.7449
## s.e. 0.0770 0.0376 0.0894 0.0834 0.0753
##
## sigma^2 = 10405: log likelihood = -1292.96
## AIC=2597.91 AICc=2598.32 BIC=2618.08df <- ts_to_prophet(AirPassengers) %>% setNames(c("date", "y"))
df$lag12 <- dplyr::lag(df$y, n = 12)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, uniondf$month <- factor(month(df$date, label = TRUE), ordered = FALSE)
df$trend <- 1:nrow(df)par <- ts_split(ts.obj = AirPassengers, sample.out = 12)
train <- par$train
test <- par$testtrain_df <- df[1:(nrow(df) - 12), ]
test_df <- df[(nrow(df) - 12 + 1):nrow(df), ]md1 <- tslm(train ~ season + trend + lag12, data = train_df)checkresiduals(md1)
##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals from Linear regression model
## LM test = 83.147, df = 24, p-value = 1.903e-08En el gráfico de autocorrelación anterior, podemos observar que la serie de residuos presenta una fuerte correlación con sus retrasos anteriores, lo que indica que la serie no cumple con el criterio de ruido blanco. Esta conclusión sugiere que el modelo de regresión no ha logrado capturar todos los patrones de la serie, lo cual es común ya que la variación en la serie puede verse influenciada por otras variables como los precios del petróleo y de los boletos, la tasa de desempleo, entre otros factores externos.
Para resolver este problema, una solución simple es modelar los términos de error con un modelo ARIMA y agregarlos a la regresión.
md2 <- auto.arima(train,
xreg = cbind(model.matrix(~ month,train_df)[,-1],
train_df$trend,
train_df$lag12),
seasonal = TRUE,
stepwise = FALSE,
approximation = FALSE)summary(md2)## Series: train
## Regression with ARIMA(2,0,0)(2,0,0)[12] errors
##
## Coefficients:
## ar1 ar2 sar1 sar2 monthFeb monthMar monthApr monthMay
## 0.5849 0.3056 -0.4421 -0.2063 -2.7523 0.8231 0.2066 2.8279
## s.e. 0.0897 0.0903 0.1050 0.1060 1.8417 2.3248 2.4162 2.5560
## monthJun monthJul monthAug monthSep monthOct monthNov monthDec
## 6.6057 11.2337 12.1909 3.8269 0.6350 -2.2723 -0.9918
## s.e. 3.2916 4.2324 4.1198 2.9243 2.3405 2.4211 1.9172
##
## 0.2726 1.0244
## s.e. 0.1367 0.0426
##
## sigma^2 = 72.82: log likelihood = -426.93
## AIC=889.86 AICc=896.63 BIC=940.04
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4117148 8.353629 6.397538 0.2571543 2.549682 0.2100998
## ACF1
## Training set 0.005966715checkresiduals(md2)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(2,0,0)(2,0,0)[12] errors
## Q* = 25.816, df = 20, p-value = 0.172
##
## Model df: 4. Total lags used: 24Después de aplicar el modelo de regresión lineal con errores ARIMA, se observa un cambio en el gráfico ACF, ya que la mayoría de los retrasos no son estadísticamente significativos. Esto es diferente del modelo de regresión lineal anterior (md1), que no logró capturar todos los patrones de la serie y mostró una fuerte correlación en la serie de residuos con sus retrasos anteriores. La incorporación de los errores del modelo ARIMA en la regresión ha ayudado a capturar patrones que antes se perdían y ha mejorado el rendimiento del modelo.
fc1 <- forecast(md1, newdata = test_df)
fc2 <- forecast(md2, xreg = cbind(model.matrix(~ month,test_df)[,-1],
test_df$trend,
test_df$lag12))accuracy(fc1, test)## ME RMSE MAE MPE MAPE MASE
## Training set -9.478240e-16 13.99206 11.47136 -0.1200023 4.472154 0.3541758
## Test set 1.079611e+01 20.82782 18.57391 1.9530865 3.828115 0.5734654
## ACF1 Theil's U
## Training set 0.7502578 NA
## Test set 0.1653119 0.4052175accuracy(fc2, test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.4117148 8.353629 6.397538 0.2571543 2.549682 0.2100998
## Test set -11.2123706 17.928195 13.353910 -2.4898843 2.924174 0.4385521
## ACF1 Theil's U
## Training set 0.005966715 NA
## Test set -0.325044933 0.3923795