Capitulo11seriedetiempo
Jesus Castro Solano y Margerith Escobar
2023-05-03
INTRODUCTION
El modelo de media móvil integrada autorregresiva (ARIMA) es el nombre genérico de un familia de modelos de pronóstico que se basan en el autorregresivo (AR) y en movimiento Procesos promedio (MA). Entre los modelos de pronóstico tradicionales (por ejemplo, lineal regresión, suavizado exponencial, etc.), el modelo ARfMA se considera como el más enfoque avanzado y robusto. En este capítulo presentaremos el modelo componentes, los procesos AR y MA y el componente de diferenciación. Además, nosotros se centrará en métodos y enfoques para ajustar los parámetros del modelo con el uso de diferenciación, la función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF).
El proceso estacionario
El proceso estacionario, en el contexto de datos de series de tiempo, describe un estado estocástico de la serie. Los datos de la serie temporal son estacionario si se dan las siguientes condiciones: * La media y la varianza de la serie no cambian con el tiempo * La estructura de correlación de la serie, junto con sus rezagos, permanece igual a lo largo tiempo
Por ejemplo, en el siguiente ejemplo, simularemos un proceso AR con un retraso (es decir, p = 1) y 500 observaciones con la función arima.sim. Antes de ejecutar la simulación, estableceremos el valor de s e e d en 12345
set.seed(12345)
stationary_ts <- arima.sim(model = list(order = c(1,0,0),
ar = 0.5 ),
n = 500)Ahora, tracemos la serie temporal simulada con la función ts_pIot We get the following output:
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 alrededor de 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 para series no estacionarias y obetendremos la siguuiente salida:
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")Por otro lado, el segundo ejemplo viola la condición estacionaria ya que tiene tendencia ,con el tiempo, lo que significa que está cambiando con el tiempo. Consideraríamos los datos de una serie de tiempo como no estacionarios siempre que esas condiciones no se sostienen.
La serie clásica Ai rPassenger (el número mensual de pasajeros de aerolíneas entre 1949 y 1960) del paquete de conjuntos de datos es un buen ejemplo de una serie que viola los dos condiciones del proceso estacionario. Dado que la serie tiene tanto una fuerte tendencia lineal como una componente estacional multiplicativo, la media y la varianza cambian con el tiempo
Obtenemos la siguiente salida
data(AirPassengers)
ts_plot(AirPassengers,
title = "Monthly Airline Passenger Numbers 1949-1960",
Ytitle = "Thousands of Passengers",
Xtitle = "Year")Transformar una serie no estacionaria en una serie estacionaria
para manejar esto, deberá aplicar algunos pasos de transformación para llevar la serie a un estado estacionario. Los métodos de transformación comunes son diferentes la serie (o eliminación de tendencia) y la transformación logarítmica (o ambas).
Diferenciar series de tiempo
El enfoque más común para transformar datos de series de tiempo no estacionarios en estado estacionario es diferenciando la serie con sus rezagos. El principal efecto de diferenciar un La serie es la eliminación de la tendencia de la serie (o eliminación de la serie), que ayuda a estabilizar la media de la serie.
Aquí, f representa la frecuencia de la serie e (Yt-f) representa el retraso estacional de la serie. Es común utilizar la diferenciación estacional cuando una serie tiene un componente estacional. La función diff del paquete base diferencia la serie de entrada con un retraso específico por estableciendo el argumento “lag” de la función en el retraso relevante.
Volvamos a la serie AirPassenger y veamos cómo el primer orden y la diferenciación estacional afectan el estructura de la serie. Comenzaremos con la diferencia de primer orden
ts_plot(diff(AirPassengers, lag = 1),
title = "AirPassengers Series - First Differencing",
Xtitle = "Year",
Ytitle = "Differencing of Thousands of Passengers")Puede ver que la primera diferencia de la serie AirPassenger eliminó la tendencia de la serie y que la media de la serie es, en general, constante a lo largo del tiempo. Por otra parte, existe una clara evidencia de que la variación de la serie va aumentando en el tiempo, y por tanto la 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. Agreguemos la diferencia estacional a la diferencia de primer orden y grafiquémosla de nuevo
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 de registro
La Descomposición de Datos de Series Temporales, vimos las aplicaciones del log transformación para manejar datos de series de tiempo multiplicativas. Podemos utilizar esto enfoque para estabilizar una oscilación estacional multiplicativa, Por ejemplo, en el ejemplo de la AirPassenger en la sección anterior, vimos que la primera diferenciación está haciendo un gran trabajo en la estabilización de la media de la serie, pero no es suficiente para estabilizar la varianza de la serie
por lo tanto, podemos aplicar una transformación log para transformar la estacionalidad estructura de multiplicable a aditivo y luego aplicar la diferencia de primer orden a 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")El proceso de paseo aleatorio
El proceso aleatorio, en el contexto de series de tiempo, describe un proceso estocástico de un objeto con el tiempo
Una caminata aleatoria se usa comúnmente para simular posibles caminos futuros de un serie. Por ejemplo, podemos simular una caminata aleatoria con el arima. función de simulación por estableciendo el parámetro d del argumento or de r en 1. Esto es equivalente a un no estacionario serie con una estructura de primera diferencia. El siguiente código demuestra la simulación de 2 O caminos aleatorios diferentes de 5 0 0 pasos, todos comenzando en el punto 0 en el tiempo 0. Crearemos dos gráficos: uno para los caminos aleatorios y otro para su diferencia de primer orden. Nosotros usará el paquete plotly para hacer esto:
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í, pl representa la trama de la simulación de paseo aleatorio:
p1 %>% layout(title = "Simulate Random Walk",
yaxis = list(title = "Value"),
xaxis = list(title = "Index")) %>%
hide_legend()Aquí, p2 representa la gráfica correspondiente de la diferenciación de primer orden de la distribución aleatoria simulación de caminata:
p2 %>% layout(title = "Simulate Random Walk with First-Order Differencing",
yaxis = list(title = "Value"),
xaxis = list(title = "Index")) %>%
hide_legend()The AR process
En el siguiente ejemplo, utilizaremos el arima .sim función de nuevo para simular un AR(2) procesa datos de series de tiempo de estructura con 5 0 0 observaciones, y luego los usa para ajustar un AR modelo. Usaremos el argumento del modelo para establecer el orden AR en 2 y establecer los retrasos coeficientes fl - 0,9 y f2 - -0’3:
set.seed(12345)
ar2 <- arima.sim(model = list(order = c(2,0,0),
ar = c(0.9, -0.3)),
n = 500) Revisemos la serie temporal simulada:
Obtenemos la siguiente salida:
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 en datos de series de tiempo y que pronosticar sus valores futuros. Esta función identifica la orden AR automáticamente en función de el criterio de información de Akaike (AIC). El argumento met hod le permite definir el método de estimación de coeficientes, como el COLS de mínimos cuadrados ordinarios) (que vimos en Capítulo 9, Pronósticos con regresión lineal), estimación de máxima verosimilitud (MLE), y Yule-Walker (predeterminado). Apliquemos la función a r para identificar el orden AR y estimar su coeficientes en consecuencia:
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.049Identificar el proceso AR y sus características
En el ejemplo anterior, simulamos una serie (2), y estaba claro que necesitamos aplicar un modelo AR sobre los datos. Sin embargo, cuando trabaja con datos de series en tiempo real, tendrá que identificar la estructura de la serie antes de modelarla.
Análisis de correlación, nos permite clasificar el proceso escriba e identifique su orden. Si la salida de ACF se desvanece y la salida de PACF se corta en el retraso p, esto indica que la serie es un proceso AR(p). Calculemos y tracemos el ACF y el PACF para la serie AR(2) simulada que creamos previamente con las funciones ACF y PACF. Primero, usaremos la función pa r para trazar los dos gráficos uno al lado del otro configurando la fila mf argumento a C (I , 2 ) (una fila, dos columnas
Y generamos los gráficos con las funciones acf y pacf:
par(mfrow=c(1,2))
acf(ar2)
pacf(ar2)
El proceso de la media móvil
En algunos casos, el modelo de pronóstico no puede capturar todos los patrones de la serie, y por lo tanto, queda algo de información en los residuos del modelo (o error de pronóstico) et . La meta del proceso de promedio móvil es capturar patrones en los residuos, si existen, mediante modelando la relación entre Yt, el término de error, et, y los términos de error pasados q del modelos
De manera similar, ya que simulamos la serie AR(2) en la sección anterior, vamos a utiliza el a r ima . s im función para simular una serie con una estructura (2). En este caso, nosotros establecerá el parámetro q en el argumento order r en 2 y establecerá los coeficientes MA en 01 - 0.5 y f1= 02 -f2= -0.3
set.seed(12345)
ma2 <- arima.sim(model = list(order = c(0, 0, 2),
ma = c(0.5, -0.3)),
n = 500) Usaremos la función ts_pIot para trazar la serie simulada
Obtenemos la siguiente salida:
ts_plot(ma2,
title = "Simulate MA(2) Series",
Ytitle = "Value",
Xtitle = "Index")Identificar el proceso de MA y sus características
Similar al proceso AR, podemos identificar un proceso MA y su orden con el ACF y Funciones PACF. Si el ACF se corta en el retraso q y la función PACF se desvanece, podemos concluir que el proceso es un (q). Repitamos el proceso que aplicamos en la serie ar2 con la serie ma2:
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.62par(mfrow=c(1,2))
acf(ma2)
pacf(ma2)
The ARMA model
Hasta ahora hemos visto como las solicitudes de AR y MA se tramitan por separado. Sin embargo, en algunos casos, la combinación de los dos nos permite manejar series temporales más complejas. datos. El modelo ARMA es una combinación de los procesos AR(p) y MA(q) y puede ser escrito de la siguiente manera:
Simulemos una serie de datos de tiempo con una estructura ARMA (1,2) con arima. función de simulación y revisar las características del modelo ARMA. Estableceremos los parámetros p y q de el argumento de orden a 1 y 2, respectivamente, y establecer el coeficiente AR como f1- 0,7, y el Coeficientes MA como 01 - O’5 y f1= 02 - f2 =-0,3:
set.seed(12345)
arma <- arima.sim(model = list(order(1,0,2), ar = c(0.7), ma = c(0.5,-0.3)), n = 500)Grafiquemos y revisemos la estructura de la serie con la función ts_plot:
ts_plot(arma,
title = "Simulate ARMA(1,2) Series",
Ytitle = "Value",
Xtitle = "Index")Ajustar un modelo A es sencillo con la función arima. En este caso, tenemos para establecer los parámetros p y q en el argumento de pedido
Puede observar a partir de la siguiente salida del modelo ajustado que los valores del modelo los coeficientes son bastante cercanos al que simulamos
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.36#Identifying an ARMA process
Identificar el proceso ARMA sigue el mismo enfoque que usamos previamente con el Procesos AR y MA. Existe un proceso ARMA en datos de series de tiempo si tanto el ACF como el Las parcelas PACF disminuyen, como podemos ver en el siguiente ejemplo
par(mfrow=c(1,2))
acf(arma)
pacf(arma)
Ajuste manual del modelo ARMA
El ajuste manual del modelo ARMA se basa principalmente en la experimentación, la intuición, el común sentido y experiencia.
usemos el Serie ARMA simulada que creamos anteriormente para aplicar este enfoque de ajuste.
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 = 1452Ahora comenzaremos a aumentar el valor de p y q mientras monitoreamos el cambio de la puntuación AIC. El siguiente modelo es el modelo ARMA(2,1):
arima(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.52El puntaje AIC de ARMA (1,2) es incluso más bajo que el de ARMA (2,1), y por lo tanto es ahora el modelo superior. Ahora, probaremos la última combinación de ARMA (2,2) y veremos si puede lograr mejoras adicionales:
arima(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.51Por ejemplo, podemos usar la función checkresidual s del paquete de pronóstico para para validar esos supuestos en el modelo seleccionado
Puede observar en el siguiente gráfico de salida de checkresidual que los residuos de la El modelo ARMA(1,2) satisface los supuestos 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
Pronosticar cualquiera de los modelos que vimos hasta ahora fue sencillo: usamos el función de pronóstico del paquete de pronóstico de manera similar a como lo usamos en el capítulo previo. Por ejemplo, el siguiente código demuestra el pronóstico de los próximos 100 observaciones del modelo AR que entrenamos previamente en la sección de proceso AR con el ar función:
ar_fc <- forecast(md_ar, h = 100)Podemos usar pIot_forecast para trazar el resultado del pronóstico
Obtenemos la siguiente salida:
plot_forecast(ar_fc,
title = "Forecast AR(2) Model",
Ytitle = "Value",
Xtitle = "Year")#El modelo ARIMA
Una de las limitaciones de los modelos AR, MA y ARMA es que no pueden manejar datos de series de tiempo no estacionarios. Por lo tanto, si la serie de entrada no es estacionaria, un preprocesamiento se requiere un paso para transformar la serie de un estado no estacionario a un estado estacionario. El modelo proporciona una solución para este problema al agregar el proceso integrado para El modelo ARMA. El proceso Integrado (I) es simplemente diferenciar la serie con sus rezagos, donde el grado de diferenciación está representado por el parámetro d. la diferenciación proceso, como vimos anteriormente, es una de las formas en que puede transformar los métodos de una serie de no estacionario a estacionario.
Identificación del grado de diferenciación del modelo
Similar a los parámetros p y q, establecer el parámetro d (el grado de diferenciación de la serie) se puede hacer con las gráficas ACF y PACF. En el siguiente ejemplo, utilizaremos el precios mensuales del café Robusta desde 2000. Esta serie es parte de Coffee_Prices varios objetos de series temporales del paquete TSstudio. Comenzaremos cargando el Serie Coffee_Prices y restando los precios mensuales de Robusta desde enero de 2010 usando la función de ventana
data(Coffee_Prices)
robusta_price <- window(Coffee_Prices[,1], start = c(2000, 1))Tracemos la serie robusta_price y revisemos su estructura con la función ts_plot y obetenmos la siguiente entrada:
ts_plot(robusta_price,
title = "The Robusta Coffee Monthly Prices",
Ytitle = "Price in USD",
Xtitle = "Year")Como puede ver, los precios del café Robusta a lo largo del tiempo tienen una tendencia alcista y, por lo tanto, no es en estado estacionario. Además, dado que esta serie representa precios continuos, es probable que la serie tiene una fuerte relación de correlación con sus rezagos pasados (ya que los cambios en el precio son normalmente cerca del precio anterior). Usaremos la función acf nuevamente para identificar el tipo de relación entre la serie y sus rezagos
acf(robusta_price)
Como puede ver en el resultado anterior de la gráfica ACF, la correlación de la serie con su lags decae lentamente con el tiempo de manera lineal. Eliminando tanto la tendencia de la serie tbe como la correlación entre la serie y sus rezagos se puede hacer diferenciando la serie. Lo haremos Comience con la primera diferenciación usando la función diff:
robusta_price_d1 <- diff(robusta_price)Repasemos la primera diferencia de la serie con las funciones acf y pacf, y obtenemos la siguiente entrada
par(mfrow=c(1,2))
acf(robusta_price_d1)
pacf(robusta_price_d1)
Los gráficos ACF y PACF de la primera diferencia de la serie indican que un proceso AR(1) es apropiado para usar en la serie diferenciada ya que el ACF está disminuyendo y el PACF corta en el primer retraso. Por lo tanto, ganamos aplicar un ARIMA (1,1,0)
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.001526295checkresiduals(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 este chunk podemos observar que los residuos del modelo y la prueba, arrojan que estos residuos son ruido blanco. No obstante, hay otros datos, pero se encuentran casi que insignificativos. Por lo tanto, se omiten.
data(USgas)ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Esta serie presenta un patron estacioanl de ARIMA. Tmanbien, la grafica muestra una tendencia creciente o al alza.
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)
En este grafico se puede ver una correlacion de los desfases
estacionales y no estacionales. Al igual que, una disminucion en los
desfases(lineal), la cual se puede asociar con que la serie no es
estacionaria.
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")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")La anterior grafica muestra una serie que oscila en un rango amplio. Por lo tanto,la serie no se encuentra estable.
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.211253 104.79281 83.09943 1.4913412 3.314280 0.7216918
## ACF1 Theil's U
## Training set 0.004565602 NA
## Test set -0.049999868 0.3469228test_forecast(USgas,
forecast.obj = USgas_test_fc,
test = test)El fragmento indica que el modelo SARIMA es capaz de identificar adecuadamente los patrones estacionales y de tendencia de la serie, pero tiene dificultades para capturar los picos estacionales (específicamente, en enero) durante la fase de entrenamiento. Además, en la fase de prueba, el modelo presenta un error absoluto del 6,7% para el mes de enero (que es el pico anual), debido a la alta fluctuación durante el invierno. Sin embargo, para hacer frente a esta incertidumbre en los momentos de mayor demanda, es posible utilizar los intervalos de confianza del modelo o las simulaciones de ruta, tal como se describe en el capítulo 8 de estrategias de pronóstico.
final_md <- arima(USgas, order = c(1,1,1), seasonal = list(order = c(2,1,1)))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: 24al analizar el gráfico de residuos, se puede observar que estos están distribuidos normalmente y se comportan como ruido blanco. Además, la prueba de Ljung-Box confirma que no hay autocorrelación en los residuos, ya que el valor de p es de 0,12 y no se puede rechazar la hipótesis nula.
USgas_fc <- forecast(final_md, h = 12)plot_forecast(USgas_fc,
title = "US Natural Gas Consumption - Forecast",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")la función “auto.arima” del paquete “forecast” automatiza el proceso de ajuste de los modelos ARIMA, eliminando la necesidad de realizar numerosos pasos manuales para verificar la estacionaridad de la serie y ajustar los parámetros del modelo. Esto facilita el proceso de pronóstico, especialmente cuando se tienen múltiples series para pronosticar.
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.0117 -0.2673 -0.7431
## s.e. 0.0794 0.0741 0.0452 0.0887 0.0830 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-08Claro, puedo parafrasear ese párrafo. En el gráfico de autocorrelación, podemos ver que los residuos de la serie tienen una alta correlación con los retrasos anteriores, lo que indica que la serie no sigue el patrón de ruido blanco. Esto implica que el modelo de regresión no ha capturado todos los patrones en la serie, lo cual es normal ya que los cambios en la serie pueden estar influenciados por diversas variables externas como los precios del petróleo y los boletos, la tasa de desempleo y otros factores. Es posible que se necesiten ajustes adicionales en el modelo para capturar estos patrones faltantes y que los residuos finalmente se asemejen al ruido blanco.
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 monthabr 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 monthago monthsept monthoct monthnov monthdic
## 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.005966712checkresiduals(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: 24fc1 <- forecast(md1, newdata = test_df)
fc2 <- forecast(md2, xreg = cbind(model.matrix(~ month,test_df)[,-1],
test_df$trend,
test_df$lag12))accuracy(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.005966712 NA
## Test set -0.325044929 0.3923795