Metodología de Box-Jenkins
Brian Gómez
26 de enero de 2020
Para las series analizadas, estimar modelos ARIMA apropiados basados en la metodología de Box-Jenkins
El acrónimo ARIMA significa media móvil integrada autorregresiva. Los retrasos de las series estacionarias en la ecuación de pronóstico se denominan términos “autorregresivos”, los retrasos de los errores de pronóstico se denominan términos de “promedio móvil” y se dice que una serie de tiempo que debe diferenciarse para hacerse estacionaria es “integrada” versión de una serie estacionaria. Los modelos de caminata aleatoria y de tendencia aleatoria, los modelos autorregresivos y los modelos de suavizado exponencial son casos especiales de los modelos ARIMA.
Un modelo ARIMA no estacional se clasifica como un modelo “ARIMA (p, d, q)”, donde: p es el número de términos autorregresivos, d es el número de diferencias no estacionales necesarias para la estacionariedad, y q es el número de errores de pronóstico rezagados en la ecuación de predicción.
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.
Un modelo autorregresivo es simplemente una regresión lineal del valor actual de la serie contra uno o más valores anteriores de la serie. El valor de p se llama el orden del modelo AR.
Los modelos AR se pueden analizar con uno de varios métodos, incluidas las técnicas estándar de mínimos cuadrados lineales. También tienen una interpretación directa.
Un modelo de promedio móvil (MA) es conceptualmente una regresión lineal del valor actual de la serie contra el ruido blanco o choques aleatorios de uno o más valores anteriores de la serie. Se supone que los choques aleatorios en cada punto provienen de la misma distribución, típicamente una distribución normal, con ubicación en cero y escala constante. La distinción en este modelo es que estos choques aleatorios se propagan a valores futuros de las series de tiempo. Ajustar las estimaciones de MA es más complicado que con los modelos AR porque los términos de error no son observables.
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.
Supuestos del modelo ARIMA
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.
Pasos a seguir para el modelado ARIMA:
Identificación del Modelo.
Ajustar el modelo.
Diagnostico del modelo.
Pronósticos.
Previamente analizamos las series IBM y Y1. Ahora estimaremos modelos ARIMA apropiados para las series IBM y Y1 basados en la metodología de Box-Jenkins.
IBM
1. Identificación del Modelo
1.1. Gráfico de la Serie
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.
Y=read.csv("TsData.csv")
ibm<-ts(Y$ibm,start=1990,frequency=12, end=c(1999,12))
plot.ts(ibm)
En la figura aparece el gráfico de la serie original, el grafico nos indica cambios en el nivel de la serie así como tendencia, por tanto, es no estacionaria. Esto es corroborado por el correlograma muestral así como las pruebas de Dickey-Fuller.
1.2. Análisis de la Autocorrelación
La función de autocorrelación (acf ()) proporciona la autocorrelación en todos los retrasos posibles. La autocorrelación en el retraso 0 se incluye por defecto, que siempre toma el valor 1, ya que representa la correlación entre los datos y ellos mismos.
acf(ibm,main="",sub="Figura 1: Función de Autocorrelación Simple")
Como podemos deducir del gráfico anterior, la autocorrelación continúa disminuyendo a medida que aumenta el retraso, lo que confirma que no existe una asociación lineal entre las observaciones separadas por retrasos más grandes.
pacf () en el rezago \(k\) es una función de autocorrelación que describe la correlación entre todos los puntos de datos que son exactamente \(k\) pasos separados, después de tener en cuenta su correlación con los datos entre esos \(k\) pasos. Ayuda a identificar el número de coeficientes de autorregresión (AR) (valor p) en un modelo ARIMA.
pacf(ibm, main="", sub="Figura 2: Función de Autocorrelación Parcial")
Confirmaremos lo dicho anteriormente con la prueba de Dickey - Fuller.
1.3 Prueba de Dickey - Fuller
Hipótesis:
\(H_0:\) La serie es no estacionaria: Tiene raíz unitaria
\(H_1:\) La serie es estacionaria: No tiene raíz unitaria
Estadístico de prueba:
adfTest(ibm)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: 1.8653
## P VALUE:
## 0.9837
##
## Description:
## Fri Feb 14 12:25:39 2020 by user: mamaEl valor del estadístico es 0,987, con un \(p\)-valor mayor a 0,05 por lo que la hipótesis nula no se rechaza, la serie original no es estacionaria.
Tomando la primera diferencia
El siguiente paso es tomar la primera diferencia a la serie.
plot.ts(diff(ibm), main="", sub="Figura 3: Serie en primera diferencia", xlab="Tiempo",ylab="Primera diferencia")
Al tomar las primeras diferencias observamos que la serie se estabiliza en torno a un valor medio (Figura 3).
acf(diff(ibm),main="",sub="Figura 4: AFC en Primera Diferencia")
pacf(diff(ibm),main="",sub="Figura 5: PAFC en Primera Diferencia")
Así mismo el correlograma confirma un declinamiento rápido en el tiempo de la AFC y PAFC. (Figuras 4 y 5).
Prueba de Dickey-Fuller a la serie en primera diferencia
Hipótesis:
\(H_0:\) La serie es no estacionaria: Tiene raíz unitaria
\(H_1:\) La serie es estacionaria: No Tiene raíz unitaria
Estadístico de prueba:
adfTest(diff(ibm))## Warning in adfTest(diff(ibm)): p-value smaller than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -7.6162
## P VALUE:
## 0.01
##
## Description:
## Fri Feb 14 12:25:40 2020 by user: mamaEl valor del estadístico es -7,6162, con un \(p\) valor menor a 0,05 por lo que la hipótesis nula se rechaza. La serie en primera diferencia es estacionaria.
Una vez que se ha obtenido una serie estacionaria, debemos identificar la forma del modelo que se utilizará. La identificación de la forma del modelo se lleva a cabo comparando las autocorrelaciones y las autocorrelaciones parciales calculadas con los datos de la autocorrelaciones y las autocorrelaciones parciales teóricas de los diferentes modelos ARIMA.
En los correlogramas podemos observar lo siguiente:
- En la ACF, la primera autocorrelación es significativa y además parece ser que los valores de la PACF disminuyen en forma oscilatoria. Esto nos lleva a pensar en un primer modelo ARIMA(2,1,1) y otro modelo ARIMA(1,1,2).
2. Ajuste del modelo
Una vez que se han seleccionado los modelos tentativos, se deben estimar los parámetros para estos modelos.
Modelo 1:
La salida de arima() incluye los coeficientes ajustados y el error estándar (s.e.) para cada coeficiente. Observando los coeficientes podemos excluir los insignificantes. Podemos usar una función confint() para este propósito.
arima1<-arima(diff(ibm),c(2,1,1), method="ML")
arima1##
## Call:
## arima(x = diff(ibm), order = c(2, 1, 1), method = "ML")
##
## Coefficients:
## ar1 ar2 ma1
## 0.2635 -0.1800 -1.000
## s.e. 0.0914 0.0915 0.039
##
## sigma^2 estimated as 25.2: log likelihood = -360.17, aic = 726.34coeftest(arima1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.263524 0.091408 2.8829 0.00394 **
## ar2 -0.180016 0.091452 -1.9684 0.04902 *
## ma1 -0.999999 0.039042 -25.6137 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(arima1)## 2.5 % 97.5 %
## ar1 0.08436662 0.4426805111
## ar2 -0.35925735 -0.0007740059
## ma1 -1.07651953 -0.9234794188Modelo 2:
arima2<-arima(diff(ibm),c(1,1,2), method="ML")
arima2##
## Call:
## arima(x = diff(ibm), order = c(1, 1, 2), method = "ML")
##
## Coefficients:
## ar1 ma1 ma2
## -0.0765 -0.6561 -0.3439
## s.e. 0.2355 0.2134 0.2113
##
## sigma^2 estimated as 25.65: log likelihood = -361.07, aic = 728.13coeftest(arima2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.076546 0.235518 -0.3250 0.745173
## ma1 -0.656052 0.213421 -3.0740 0.002112 **
## ma2 -0.343946 0.211300 -1.6278 0.103576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(arima2)## 2.5 % 97.5 %
## ar1 -0.5381529 0.38506114
## ma1 -1.0743502 -0.23775448
## ma2 -0.7580867 0.07019455Para escoger el mejor modelo se debe encontrar el Criterio de información (AIC) de Akaike para un conjunto de modelos y comparar los modelos con los valores de AIC más bajos. Los resultados son similares en cuanto al valor de AIC y \(\sigma^{2}\) por lo que decidimos por el modelo ARIMA(2,1,1).
3. Validación del Modelo
3.1. Estacionariedad de los residuos
La siguiente figura muestra una posible estacionariedad entre los residuales del modelo. Esto lo confirmaremos con un correlograma.
plot.ts(arima1$residuals, main="", sub="Figura 6: Residuos del modelo ARIMA(2,1,1)", xlab="Tiempo",ylab="Residuos")
3.2. Autocorrelación de los residuos:
acf(arima1$residuals,main="", sub="Figura 7: Autocorrelaciones de los residuos del modelo ARIMA(2,1,1)", xlab="Tiempo",ylab="Autocorrelación")
pacf(arima1$residuals,main="", sub="Figura 8: Autocorrelaciones Parciales de los residuos del modelo ARIMA(2,1,1)", xlab="Tiempo",ylab="Autocorrelación")
Los correlogramas muestran que no existe autocorrelación significativa en los residuos.
3.3. Normalidad de los residuos
Los gráficos Q-Q son una herramienta eficaz para evaluar la normalidad. Aquí las aplicamos a los residuos.
qqnorm(arima1$residuals,main="", sub="Figura 9: Gráfico Q para evaluar normalidad"); qqline(arima1$residuals)
El diagrama Q-Q de los residuos del modelo ARIMA(2,1,1) estimado para la serie se muestra en la Figura 9. Los puntos parecen seguir la línea recta bastante de cerca. Este gráfico no nos llevará a rechazar la normalidad de los términos de error en este modelo.
Además, la prueba de normalidad de Shapiro-Wilk aplicada a los residuos produce un estadístico de prueba de W = 0,97148, lo que corresponde a un valor de p de 0,1445, y no rechazará la normalidad basado en esta prueba.
shapiro.test(arima1$residuals)##
## Shapiro-Wilk normality test
##
## data: arima1$residuals
## W = 0.98839, p-value = 0.4078Como todos los gráficos respaldan la suposición de que no hay un patrón en los residuos, podemos seguir adelante y calcular pronósticos.
Pronósticos
Los parámetros del modelo ARIMA se pueden usar como modelo predictivo para hacer pronósticos de valores futuros de la serie una vez que se selecciona el modelo más adecuado para los datos de la serie temporal.
El valor d afecta los intervalos de predicción; los intervalos de predicción aumentan de tamaño con valores más altos de “d”. Todos los intervalos de predicción serán esencialmente los mismos cuando d = 0 porque la desviación estándar del pronóstico a largo plazo irá a la desviación estándar de los datos históricos.
predict(arima1,n.ahead = 5)## $pred
## Jan Feb Mar Apr May
## 2000 3.0388604 0.2078499 0.5349105 1.1307249 1.2288600
##
## $se
## Jan Feb Mar Apr May
## 2000 5.040701 5.223500 5.247507 5.258124 5.258112Necesitamos asegurarnos de que los errores de pronóstico no están correlacionados, normalmente distribuidos con media cero y varianza constante. Podemos usar la medida de diagnóstico para encontrar el modelo apropiado con los mejores valores de pronóstico posibles.
plot(forecast(auto.arima(ibm)), main="", sub = "Figura 10: Valores pronosticados")
Y1
Para estimar el modelo ARIMA apropiado para la serie Y1 aplicaremos la misma metodología que la realizada a la serie IBM previamente.
1. Identificación del Modelo
1.1. Gráfico de la Serie
El primer paso en la identificación del modelo es determinar si la serie es estacionaria.
Y=read.csv("4Series.csv")
y1<-ts(Y$Y1,start=1,frequency=12, end=c(22,12))
plot.ts(y1)
En la figura aparece el gráfico de la serie original, la cual muestra un comportamiento primero creciente y luego decreciente, no podemos afirmar que la serie sea o no estacionaria.
1.2. Análisis de la Autocorrelación
acf(y1,main="",sub="Figura 11: Función de Autocorrelación Simple")
Como podemos deducir del gráfico anterior, la autocorrelación continúa disminuyendo a medida que aumenta el retraso, lo que hace parecer que no existe una asociación lineal entre las observaciones separadas por retrasos más grandes.
pacf () Ayuda a identificar el número de coeficientes de autorregresión (AR) (valor p) en un modelo ARIMA.
pacf(y1,main="", sub="Figura 12: Función de Autocorrelación Parcial")
1.3 Prueba de Dickey - Fuller
Hipótesis:
\(H_0:\) La serie es no estacionaria: Tiene raíz unitaria.
\(H_1:\) La serie es estacionaria: no tiene raíz unitaria.
Estadístico de prueba:
adfTest(y1)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: 0.0966
## P VALUE:
## 0.6471
##
## Description:
## Fri Feb 14 12:25:46 2020 by user: mamaEl valor del estadístico es 0,0966, con un \(p\)-valor mayor a 0,05 por lo que la hipótesis nula no se rechaza, la serie original no es estacionaria.
Tomando Primera Diferencia
El siguiente paso es tomar la primera diferencia a la serie.
plot.ts(diff(y1), main="", sub="Figura 13: Serie en primera diferencia", xlab="Tiempo",ylab="Primera diferencia")
Al tomar las primeras diferencias observamos que la serie se estabiliza en torno a un valor medio. (Figura 3).
acf(diff(y1),main="",sub="Figura 14: AFC en Primera Diferencia")
pacf(diff(y1),main="",sub="Figura 15: PAFC en Primera Diferencia")
Así mismo el correlograma confirma un declinamiento rápido en el tiempo de la AFC y PAFC. (Figuras 4 y 5).
Prueba de Dickey-Fuller a la serie en Primera Diferencia.
Hipótesis:
\(H_0:\) La serie es no estacionaria: Tiene raíz unitaria.
\(H_1:\) La serie es estacionaria: no tiene raíz unitaria.
Estadístico de prueba:
adfTest(diff(y1))## Warning in adfTest(diff(y1)): p-value smaller than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -12.4194
## P VALUE:
## 0.01
##
## Description:
## Fri Feb 14 12:25:47 2020 by user: mamaEl valor del estadístico es -12,4194, con un \(p\) valor menor a 0,05 por lo que la hipótesis nula se rechaza. La serie en primera diferencia es estacionaria.
Una vez que se ha obtenido una serie estacionaria, debemos identificar la forma del modelo que se utilizará. En los correlogramas podemos observar lo siguiente:
- En la ACF, ninguna autocorrelación es significativa y además parece ser que los valores de la PACF disminuyen en forma oscilatoria. Esto nos lleva a pensar en los modelos ARIMA(1,1,0), ARIMA(0,1,1) y ARIMA(0,1,0).
2. Ajuste del modelo
Una vez que se han seleccionado los modelos tentativos, se deben estimar los parámetros para estos modelos.
Modelo 1:
La salida de arima() incluye los coeficientes ajustados y el error estándar (s.e.) para cada coeficiente.
arima3<-arima(diff(y1),c(1,1,0), method="ML")
arima3##
## Call:
## arima(x = diff(y1), order = c(1, 1, 0), method = "ML")
##
## Coefficients:
## ar1
## -0.4126
## s.e. 0.0562
##
## sigma^2 estimated as 0.3364: log likelihood = -229.15, aic = 460.31coeftest(arima3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.412604 0.056169 -7.3458 2.046e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(arima3)## 2.5 % 97.5 %
## ar1 -0.5226927 -0.3025148Modelo 2:
arima4<-arima(diff(y1),c(0,1,1), method="ML")
arima4##
## Call:
## arima(x = diff(y1), order = c(0, 1, 1), method = "ML")
##
## Coefficients:
## ma1
## -0.995
## s.e. 0.020
##
## sigma^2 estimated as 0.2154: log likelihood = -172.88, aic = 347.76coeftest(arima4)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.994976 0.019964 -49.838 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(arima4)## 2.5 % 97.5 %
## ma1 -1.034105 -0.955847Modelo 3:
arima5<-arima(diff(y1),c(0,1,0), method="ML")
arima5##
## Call:
## arima(x = diff(y1), order = c(0, 1, 0), method = "ML")
##
##
## sigma^2 estimated as 0.406: log likelihood = -253.66, aic = 507.33Para escoger el mejor modelo se debe encontrar el Criterio de información (AIC) de Akaike para un conjunto de modelos y comparar los modelos con los valores de AIC más bajos. Observamos que el modelo 2 ARIMA(0,1,1) tiene el menor AIC con 347,76 y al mismo tiempo la menor varianza por lo tanto nos quedamos con este modelo.
Validación del Modelo
3.1. Estacionariedad de los residuos
La siguiente figura muestra una posible estacionariedad entre los residuales del modelo. Esto lo confirmaremos con un correlograma.
plot.ts(arima4$residuals, main="", sub="Figura 16: Residuos del modelo ARIMA(0,1,1)", xlab="Tiempo",ylab="Residuos")
Autocorrelación de los residuos:
acf(arima4$residuals, main="", sub="Figura 17: Autocorrelaciones de los residuos del modelo ARIMA(0,1,1)", xlab="Tiempo",ylab="Autocorrelación")
pacf(arima4$residuals, main="", sub="Figura 18: Autocorrelaciones Parciales de los residuos del modelo ARIMA(0,1,1)", xlab="Tiempo",ylab="Autocorrelación")
Los correlogramas muestran que no existe autocorrelación significativa en los residuos.
3.2. Normalidad de los residuos
Los gráficos Q-Q son una herramienta eficaz para evaluar la normalidad. Aquí las aplicamos a los residuos.
qqnorm(arima4$residuals, main="", sub="Figura 19: Gráfico Q para evaluar normalidad");qqline(arima4$residuals)
El diagrama Q-Q de los residuos del modelo ARIMA(0,1,1) estimado para la serie se muestra en la Figura 9. Los puntos parecen seguir la línea recta bastante de cerca. Este gráfico no nos llevará a rechazar la normalidad de los términos de error en este modelo.
Además, la prueba de normalidad de Shapiro-Wilk aplicada a los residuos produce un estadístico de prueba de W = 0,97148, lo que corresponde a un valor de \(p\) de 0,1445, y no rechazará la normalidad basado en esta prueba.
shapiro.test(arima4$residuals)##
## Shapiro-Wilk normality test
##
## data: arima4$residuals
## W = 0.97637, p-value = 0.0002344Como todos los gráficos respaldan la suposición de que no hay un patrón en los residuos, podemos seguir adelante y calcular pronósticos.
Pronósticos
Los parámetros del modelo ARIMA se pueden usar como modelo predictivo para hacer pronósticos de valores futuros de la serie una vez que se selecciona el modelo más adecuado para los datos de la serie temporal.
El valor d afecta los intervalos de predicción; los intervalos de predicción aumentan de tamaÃño con valores más altos de “d”. Todos los intervalos de predicción serán esencialmente los mismos cuando d = 0 porque la desviación estándar del pronóstico a largo plazo irá a la desviación estándar de los datos históricos.
predict(arima4,n.ahead = 5)## $pred
## Jan Feb Mar Apr May
## 23 0.002648955 0.002648955 0.002648955 0.002648955 0.002648955
##
## $se
## Jan Feb Mar Apr May
## 23 0.4642431 0.4642489 0.4642548 0.4642606 0.4642665Necesitamos asegurarnos de que los errores de pronóstico no están correlacionados, normalmente distribuidos con media cero y varianza constante. Podemos usar la medida de diagnóstico para encontrar el modelo apropiado con los mejores valores de pronóstico posibles.
plot(forecast(auto.arima(y1)), main="", sub = "Figura 20: Valores pronosticados")