Box Jenkins en R
Análisis de la serie
En este ejercicio haremos un análisis de series de tiempo acerca de los precios del petróleo observados desde enero de 2013.
library(tseries)
library(forecast)
library(readxl)
library(knitr)
library(dplyr)
library(kableExtra)Arimar <- read_excel("C:/Users/IvanaTorresUgalde/Documents/Modelos Lineales/Arimar.xlsx")
Arimar.ts <- ts(Arimar, start=c(2013,1), frequency = 12)Arimar.ts## Jan Feb Mar Apr May Jun Jul
## 2013 97.98278 103.35294 100.77090 97.75538 97.90686 96.99315 99.67255
## 2014 88.31081 90.82224 91.39522 93.88863 95.08631 97.01105 93.32017
## 2015 40.45127 45.98606 46.01733 49.32555 52.03910 52.82573 45.16595
## 2016 23.10238 23.72933 28.32656 31.33582 36.72636 39.80715 37.97153
## 2017 44.50008 44.17497 41.93243 43.22118 43.85074 41.15105 43.87771
## 2018 57.45399 56.15676 57.22458 58.15973 62.89929 64.63686 66.42162
## 2019 54.05666 57.37981 59.46449 62.07703 60.33674 56.90537 57.87580
## 2020 49.79299 44.63505 28.92302
## Aug Sep Oct Nov Dec
## 2013 99.24414 98.14572 93.87329 88.48590 89.40966
## 2014 89.50745 84.12824 72.50936 63.47980 50.18138
## 2015 38.41927 36.29036 34.86077 32.50883 26.19610
## 2016 37.74090 36.82666 40.66836 38.36162 42.32596
## 2017 45.55484 48.16237 48.89114 53.34505 54.05576
## 2018 64.26567 68.36348 71.15330 59.82583 51.86867
## 2019 49.58496 55.05522 51.00463 50.68969 54.54901
## 2020Iniciaremos observando el comportamiento de los datos a través del tiempo y descomponiendo la serie en 4 factores.
plot(Arimar.ts, main="Precios de petroleo", ylab="Precio", col="red")
plot(decompose(Arimar.ts))
De la primera gráfica observamos que el comportamiento no es estacional, podemos comprobar esto a través del test Dickey Fuller.
df_ts<-adf.test(Arimar.ts)
df_ts##
## Augmented Dickey-Fuller Test
##
## data: Arimar.ts
## Dickey-Fuller = -1.6185, Lag order = 4, p-value = 0.7328
## alternative hypothesis: stationaryAl obtener un p-value de 0.7328 no rechazamos la hipótesis nula de que la serie no sea estacionaria.
Diferenciación
Una función útil para darnos una idea de cuántas veces será necesario hacer diferenciación es
ndiffs(Arimar.ts, alpha = 0.05)## [1] 1Sin embargo, para asegurar estacionalidad volveremos a aplicar el test de Diceky Fuller.
seriedif<-diff(Arimar.ts)
plot(seriedif)
df_d1<-adf.test(seriedif)
df_d1##
## Augmented Dickey-Fuller Test
##
## data: seriedif
## Dickey-Fuller = -3.3327, Lag order = 4, p-value = 0.07162
## alternative hypothesis: stationarySi establecemos un \(\alpha=0.05\), y lo comparamos con el p-value obtenido de 0.0716 nos encontramos con que, a pesar de que la serie se observa más o menos estacional, será necesario aplicar diferenciación una vez más.
seriedif2<-diff(Arimar.ts, differences =2)
plot(seriedif2)
adf.test(seriedif2)## Warning in adf.test(seriedif2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: seriedif2
## Dickey-Fuller = -5.1673, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryEn esta ocasión rechazamos la hipótesis nula.
Análisis de las gráficas ACF y PACF
# AR(p) I(d) MA(q)
acf(ts(seriedif2, frequency=1)) #parámetro q
pacf(ts(seriedif2, frequency=1)) #parámetro p
Con lo observado en la gráfica proponemos un modelo ARIMA(1,2,1), así como un modelo alternativo ARIMA(1,2,2).
Verificaremos adecuación del modelo mediante la prueba de Ljung-Box, buscaremos no rechazar la hipótesis nula.
Modelo 1: ARIMA(1,2,1)
modelo1<-arima(Arimar.ts,order=c(1,2,1))
summary(modelo1)##
## Call:
## arima(x = Arimar.ts, order = c(1, 2, 1))
##
## Coefficients:
## ar1 ma1
## 0.4173 -1.000
## s.e. 0.1100 0.047
##
## sigma^2 estimated as 17.51: log likelihood = -244.07, aic = 494.14
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1842755 4.136183 3.135096 -0.2061955 6.279221 0.9152235
## ACF1
## Training set -0.02472726Box.test(residuals(modelo1),type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo1)
## X-squared = 0.055051, df = 1, p-value = 0.8145Modelo 2: ARIMA(1,2,2)
modelo2<-arima(Arimar.ts,order=c(1,2,2))
summary(modelo2)##
## Call:
## arima(x = Arimar.ts, order = c(1, 2, 2))
##
## Coefficients:
## ar1 ma1 ma2
## 0.5405 -1.1421 0.1421
## s.e. 0.1992 0.2171 0.2128
##
## sigma^2 estimated as 17.45: log likelihood = -243.86, aic = 495.72
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.2039827 4.129552 3.12132 -0.2394872 6.262272 0.9112018
## ACF1
## Training set 0.001116518Box.test(residuals(modelo2),type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo2)
## X-squared = 0.00011224, df = 1, p-value = 0.9915Comparación de modelos
Por resultados de la prueba Ljung-Box podemos afirmar que ambos modelos se ajustan bien a los datos, usaremos el criterio AIC para elegir alguno de los dos
| Modelo 1 | Modelo 2 | |
|---|---|---|
| AIC | 494.1389 | 495.719 |
Elegiremos el Modelo 1 ya que tiene un menor AIC. Por lo que ahora podemos calcular los pronósticos.
Pronósticos
pronostico<-forecast::forecast(modelo1,h=10)
pronostico## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Apr 2020 21.856654 16.4633719 27.24994 13.608341 30.10497
## May 2020 18.398192 8.9934192 27.80296 4.014833 32.78155
## Jun 2020 16.445341 3.6136910 29.27699 -3.178974 36.06966
## Jul 2020 15.120795 -0.6776994 30.91929 -9.040917 39.28251
## Aug 2020 14.058445 -4.3652157 32.48211 -14.118113 42.23500
## Sep 2020 13.105513 -7.6904599 33.90149 -18.699183 44.91021
## Oct 2020 12.198241 -10.7797404 35.17622 -22.943549 47.34003
## Nov 2020 11.310024 -13.7033148 36.32336 -26.944576 49.56462
## Dec 2020 10.429758 -16.5031190 37.36264 -30.760523 51.62004
## Jan 2021 9.552811 -19.2060240 38.31165 -34.430032 53.53565plot(pronostico)