Box y Jenkins popularizaron un enfoque que combina el promedio móvil y los enfoques autorregresivos en el libro “Análisis de series temporales: pronóstico y control” (Box, Jenkins y Reinsel, 1994).
Aunque ya se conocían enfoques autorregresivos y de promedio móvil (y fueron investigados originalmente por Yule), la contribución de Box y Jenkins fue desarrollar una metodología sistemática para identificar y estimar modelos que pudieran incorporar ambos enfoques. Esto hace que los modelos Box-Jenkins sean una clase poderosa de modelos.
El modelo ARMA de Box-Jenkins es una combinación de los modelos AR y MA
Por lo general, el ajuste efectivo de los modelos Box-Jenkins requiere al menos una serie moderadamente larga. Chatfield (1996) recomienda al menos 50 observaciones. Muchos otros recomendarían al menos 100 observaciones.
Los datos deben ser estacionarios; por estacionario, significa que las propiedades de la serie no dependen del momento en que se capturan. Una serie de ruido blanco y series con comportamiento cíclico también pueden considerarse como series estacionarias.
Los datos deben ser univariados: ARIMA funciona en una sola variable. La regresión automática tiene que ver con la regresión con los valores pasados.
Análisis exploratorio de datos
Descomposición de la serie
Identificación del Modelo.
Ajustar el modelo.
Diagnostico del modelo.
Pronósticos.
Demanda mensual de gasolina Ontario millones de galones 1960-1975.
sales_ts Jan Feb Mar Apr May Jun Jul Aug Sep Oct
1960 87695 86890 96442 98133 113615 123924 128924 134775 117357 114626
1961 92188 88591 98683 99207 125485 124677 132543 140735 124008 121194
1962 101007 94228 104255 106922 130621 125251 140318 146174 122318 128770
1963 108497 100482 106140 118581 132371 132042 151938 150997 130931 137018
1964 109894 106061 112539 125745 136251 140892 158390 148314 144148 140138
1965 109895 109044 122499 124264 142296 150693 163331 165837 151731 142491
1966 116963 118049 137869 127392 154166 160227 165869 173522 155828 153771
1967 124046 121260 138870 129782 162312 167211 172897 189689 166496 160754
1968 139625 137361 138963 155301 172026 165004 185861 190270 163903 174270
1969 146182 137728 148932 156751 177998 174559 198079 189073 175702 180097
1970 154277 144998 159644 168646 166273 190176 205541 193657 182617 189614
1971 158167 156261 176353 175720 193939 201269 218960 209861 198688 190474
1972 166286 170699 181468 174241 210802 212262 218099 229001 203200 212557
1973 188992 175347 196265 203526 227443 233038 234119 255133 216478 232868
1974 194784 189756 193522 212870 248565 221532 252642 255007 206826 233231
1975 199024 191813 195997 208684 244113 243108 255918 244642 237579 237579
Nov Dec
1960 107677 108087
1961 111634 111565
1962 117518 115492
1963 121271 123548
1964 124075 136485
1965 140229 140463
1966 143963 143898
1967 155582 145936
1968 160272 165614
1969 155202 174508
1970 174176 184416
1971 194502 190755
1972 197095 193693
1973 221616 209893
1974 212678 217173
1975 217775 227621
El primer paso en la identificación del modelo es determinar si la serie es estacionaria, es decir, si la serie de tiempo parece variar alrededor de un nivel fijo. Es útil observar una gráfica de la serie junto con la función de autocorrelación de la muestra. Examinemos el comportamiento de la serie a través del tiempo mediante un gráfico.
En la figura aparece el gráfico de la serie original, este parece indicar que la serie tiene una tendencia aditiva y estacional.
Estos dos componentes se definen de la siguiente manera:
La estacionalidad es corroborada por la prueba de Dickey-Fuller.
library(tseries)
adf.test(sales_ts)
Augmented Dickey-Fuller Test
data: sales_ts
Dickey-Fuller = -10.096, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary
La descomposición de series de tiempo implica pensar en una serie como una combinación de componentes de nivel, tendencia, estacionalidad y ruido. La descomposición proporciona un modelo abstracto útil para pensar en series temporales en general y para comprender mejor los problemas durante el análisis y la previsión de series temporales.
Según el diagrama de descomposición anterior, los datos tienen una tendencia creciente y tienen estacionalidad.
sales_test <- tail(sales_ts, 24)
sales_train <- head(sales_ts, length(sales_ts) - length(sales_test))Debido a que nuestros datos de series de tiempo tienen un componente de tendencia y estacional, el modelo se construirá utilizando SARIMA.
¿QUE ES EL SARIMA? La media móvil integrada autorregresiva (ARIMA) es una combinación de dos métodos, que son la media móvil autorregresiva (AR) y la media móvil (MA). Un modelo ARIMA se puede entender delineando cada uno de sus componentes de la siguiente manera:
Autorregresión (AR): se refiere a un modelo que muestra una variable cambiante que retrocede en sus propios valores retrasados o previos.
Integrado (I): representa la diferenciación de las observaciones sin procesar para permitir que la serie temporal se vuelva estacionaria (es decir, los valores de los datos se reemplazan por la diferencia entre los valores de los datos y los valores anteriores).
Media móvil (MA): incorpora la dependencia entre una observación y un error residual de un modelo de media móvil aplicado a observaciones retrasadas.
El modelo ARIMA tiene 3 parámetros:
p: el número de observaciones de retraso en el modelo; también conocido como el orden de retardo. d: el número de veces que se diferencian las observaciones sin procesar; también conocido como el grado de diferenciación. q: el tamaño de la ventana de media móvil; también conocido como el orden de la media móvil.
Si los datos de la serie temporal tienen un patrón estacional, entonces podemos agregar el término estacional a ARIMA y se convierte en Promedio Móvil Integrado Auto Regresivo Estacional (SARIMA). El requisito para que los datos se procesen con ARIMA es que los datos deben ser estacionarios. Los datos son estacionarios si su valor fluctúa alrededor de su media. Podemos diferenciar para que los datos sean estacionarios.
Los datos son estacionarios usando 1 diff y 1 diff estacional.
adf.test(diff(sales_train, lag = 12))
Augmented Dickey-Fuller Test
data: diff(sales_train, lag = 12)
Dickey-Fuller = -3.1454, Lag order = 5, p-value = 0.09954
alternative hypothesis: stationary
adf.test(diff(sales_train, lag = 12) %>%
diff())
Augmented Dickey-Fuller Test
data: diff(sales_train, lag = 12) %>% diff()
Dickey-Fuller = -8.0316, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary
A pesar de poder usar auto.arima, podemos construir un modelo manualmente, definiendo los parametros p,d,q y P,D,Q usando el gráfico ACF y PACF.
Recordamos la estructura de SARIMA.
\[\begin{equation} SARIMA_ (p,d,q)(P,D,Q)(M) \end{equation}\] Donde, (p,d,q)<- procesos arima, (P,D,Q)M <- procesos estacionales y M = 12; datos mensuales (M frequency data).
Para la data general (p,d,q) vemos el retrazo desde la linea 1,2,3,… Podemos observar que el gráfico del ACF se corta y el gráfico del PACF decrece. Por lo tanto, se ajusta usando MA model. Luego vemos la línea de cada parcela que sea significativa o cruce la línea azul.
Para el efecto estacional (P,D,Q) vemos el retraso desde la línea 12, 24, 32,… y podemos observar que el gráfico ACF y PACF decrece. Por lo tanto, se ajusta usando AR-MA model. Luego vemos la línea de cada parcela que sea significativa o cruce la línea azul.
Del resultado anterior, la combinación de modelos para construir:
ARIMA(0,1,1)(0,1,2)[12]
ARIMA(0,1,3)(0,1,2)[12]
ARIMA(0,1,1)(0,1,3)[12]
ARIMA(0,1,3)(0,1,3)[12]
model_sarima_1 <- arima(x = sales_train, order = c(0,1,1),
seasonal = list(order = c(0,1,2), period = 12))
model_sarima_2 <- arima(x = sales_train, order = c(0,1,3),
seasonal = list(order = c(0,1,2), period = 12))
model_sarima_3 <- arima(x = sales_train, order = c(0,1,1),
seasonal = list(order = c(0,1,3), period = 12))
model_sarima_4 <- arima(x = sales_train, order = c(0,1,3),
seasonal = list(order = c(0,1,3), period = 12))forecast_model1 <- forecast(model_sarima_1, h = 24)
forecast_model2 <- forecast(model_sarima_2, h = 24)
forecast_model3 <- forecast(model_sarima_3, h = 24)
forecast_model4 <- forecast(model_sarima_4, h = 24)MAPE(y_pred = forecast_model1$mean, y_true = sales_test)*100[1] 5.183373
MAPE(y_pred = forecast_model2$mean, y_true = sales_test)*100[1] 5.950886
MAPE(y_pred = forecast_model3$mean, y_true = sales_test)*100[1] 5.478983
MAPE(y_pred = forecast_model4$mean, y_true = sales_test)*100[1] 6.202979
Call:
arima(x = sales_train, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 2),
period = 12))
Coefficients:
ma1 sma1 sma2
-0.8183 -0.6322 -0.3678
s.e. 0.0387 0.2392 0.1128
sigma^2 estimated as 2.7e+07: log likelihood = -1559.79, aic = 3127.59
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 824.2185 4990.742 3780.052 0.3594269 2.341505 0.3443899 -0.2881937
plot(model_sarima_1$residuals)
Shapiro-Wilk normality test
data: model_sarima_1$residuals
W = 0.98792, p-value = 0.1588
# Ljung-Box test
Box.test(model_sarima_1$residuals)
Box-Pierce test
data: model_sarima_1$residuals
X-squared = 13.953, df = 1, p-value = 0.0001874