Ejemplo práctico (parte I)

Img1

Comenzamos trabajando con la serie de tiempo almacenada en Airpassengers (número de viajes mensuales en aerolíneas en EU). Lo primero que haremos será visualizar los datos

autoplot(AirPassengers)

donde notamos un período de 12 y vemos que la varianza es creciente. De tal manera podemos realizar una transformación logarítmica para estabilizar la varianza

autoplot(log(AirPassengers))

Caso 1: Trabajar con la ST sin transformar

Podemos aplicar de forma directa la función auto.arima como sigue:

fit_auto <- auto.arima(AirPassengers)
fit_auto
## Series: AirPassengers 
## ARIMA(2,1,1)(0,1,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.5960  0.2143  -0.9819
## s.e.  0.0888  0.0880   0.0292
## 
## sigma^2 = 132.3:  log likelihood = -504.92
## AIC=1017.85   AICc=1018.17   BIC=1029.35

donde el modelo “propuesto” es \(ARIMA(2,1,1)(0,1,0)_{12}\).

Recordemos que uno de los supuestos hechos es que estamos trabajando con un ruido blanco \(\{Z_{t}\}\), lo cual implica que tenemos varianza constante, normalidad y no correlación en dicho proceso. Para analizar este proceso trabajaremos con los residuales del modelo (\(e_{i}=y_{i}-\hat{y_{i}}\)). Donde puede decirse que los residuales son una estimación del ruido blanco.

De tal manera, procedemos a analizar los residuales:

ggplot(fit_auto, aes(x=fit_auto$residuals)) + geom_histogram(bins=12)

y/o analizamos la qqplot correspondiente

qqnorm(fit_auto$residuals)
qqline(fit_auto$residuals, col="red")

donde parece ser que los residuales si siguen una distribución normal aunque vemos ciertos sectores donde los puntos no siguen la recta roja referente a la normalidad. Podemos hacer también pruebas de hipótesis:

Prueba de Anderson-Darling

ad.test(fit_auto$residuals)
## 
##  Anderson-Darling normality test
## 
## data:  fit_auto$residuals
## A = 1.2105, p-value = 0.003601

donde se rechaza \(H_{0}\), es decir, es plausible asumir que los residuales no siguen una distribución normal.

Veamos ahora pruebas sobre la varianza constante:

# Preparamos los valores
Y <- as.numeric(fit_auto$residuals)
X <- 1:length(fit_auto$residuals)

# Hacemos una prueba no paramétrica
bptest(Y ~ X)
## 
##  studentized Breusch-Pagan test
## 
## data:  Y ~ X
## BP = 12.75, df = 1, p-value = 0.0003559

donde \(H_{0}\) representa que la varianza es constante. Vemos que rechazamos \(H_{0}\), por lo que la varianza tampoco es constante. Podemos realizar una prueba alternativa

ncvTest(lm(Y~X))
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 22.81394, Df = 1, p = 1.7847e-06

donde también se rechaza \(H_{0}\).

Para verificar la independencia graficamos los residuales

ggtsdisplay(fit_auto$residuals)

vemos en los correlogramas que hay observaciones que caen fuera del intervalo de confianza, por lo que éstas son significativas, de este modo podemos pensar que no hay independencia. Para tener más evidencia respecto a tal afirmación, realizamos una prueba de hipótesis

Box.test(fit_auto$residuals)
## 
##  Box-Pierce test
## 
## data:  fit_auto$residuals
## X-squared = 0.00022445, df = 1, p-value = 0.988

donde no se rechaza \(H_{0}\), es decir, no se rechaza que los datos no se correlacionan, en otras palabras, no se encontró evidencia en contra de que los residuales son “independientes”. Podemos realizar la prueba en lags específicos para ver su correlación correspondiente

Box.test(fit_auto$residuals, lag = 12)
## 
##  Box-Pierce test
## 
## data:  fit_auto$residuals
## X-squared = 10.208, df = 12, p-value = 0.5978
Box.test(fit_auto$residuals, lag = 24)
## 
##  Box-Pierce test
## 
## data:  fit_auto$residuals
## X-squared = 32.803, df = 24, p-value = 0.1083

donde en ambos no se rechaza \(H_{0}\).

Notamos que nuestro modelo es deficiente pues no se están cumpliendo los supuestos sobre los residuos (i.e. que no tenemos un ruido blanco).

Finalmente podemos hacer una prueba de estacionariedad:

adf.test(fit_auto$residuals)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  fit_auto$residuals
## Dickey-Fuller = -4.5795, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

donde no se rechaza \(H_{0}\) (no se rechaza que la serie sea estacionaria). Por ende, debemos regresar a la etapa 1 e identificar de nuevo los parámetros del modelo, para lo cual aplicaremos transformación logarítmo.

Caso 2: Transformar la ST

Una vez hecha la transformación logramos estabilizar la varianza. Realizamos nuevamente

fit_auto_t <- auto.arima(log(AirPassengers))
fit_auto_t
## Series: log(AirPassengers) 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 = 0.001371:  log likelihood = 244.7
## AIC=-483.4   AICc=-483.21   BIC=-474.77

donde el nuevo modelo ahora es \(ARIMA(0,1,1)(0,1,1)_{12}\) y repetimos el análisis de los supuestos:

  1. Normalidad:
# Q-Q plot
qqnorm(fit_auto_t$residuals)
qqline(fit_auto_t$residuals, col="red")

y prueba de hipótesis

ad.test(fit_auto_t$residuals)
## 
##  Anderson-Darling normality test
## 
## data:  fit_auto_t$residuals
## A = 0.68703, p-value = 0.07132

donde no se rechaza \(H_{0}\) (no rechazamos la normalidad).

  1. Varianza constante:
# Preparamos los valores
Y <- as.numeric(fit_auto_t$residuals)
X <- 1:length(fit_auto_t$residuals)

# Hacemos las pruebas de hipótesis
bptest(Y ~ X)
## 
##  studentized Breusch-Pagan test
## 
## data:  Y ~ X
## BP = 0.31915, df = 1, p-value = 0.5721
ncvTest(lm(Y ~ X))
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.4571643, Df = 1, p = 0.49895

en ambas prueba no rechazamos \(H_{0}\) (no rechazamos que la varianza sea constante).

  1. Correlación entre los datos:
ggtsdisplay(fit_auto_t$residuals)

notamos algunos valores en los correlogramas que salen del intervalo de confianza. Nos auxiliamos de la prueba de hipótesis siguiente:

# Utilizamos distintos lags

Box.test(fit_auto$residuals)
## 
##  Box-Pierce test
## 
## data:  fit_auto$residuals
## X-squared = 0.00022445, df = 1, p-value = 0.988
Box.test(fit_auto$residuals, lag = 1)
## 
##  Box-Pierce test
## 
## data:  fit_auto$residuals
## X-squared = 0.00022445, df = 1, p-value = 0.988
Box.test(fit_auto$residuals, lag = 24)
## 
##  Box-Pierce test
## 
## data:  fit_auto$residuals
## X-squared = 32.803, df = 24, p-value = 0.1083

donde no se rechaza \(H_{0}\), es decir que los datos no están correlacionados.

  1. Estacionariedad:
adf.test(fit_auto_t$residuals)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  fit_auto_t$residuals
## Dickey-Fuller = -4.7907, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

vemos que obtenemos el mismo \(p-value\) por lo que es lo mismo aplicar la prueba a los adots transformados y no transformados. Tenemos entonces que no se rechaza que la serie no sea estacionaria.

Tenemos entonces que el modelo \(ARIMA(0,1,1)(0,1,1))_{12}\) cumple los supuestos. Tenemos ahora que ver que los parámetros sean estadísticamente significativos. Para ello tenemos que verificar que el o los parámetros sean distintos de cero (o que el cero no esté en el intervalo de confianza):

# Calculamos los intervalos de confianza para las estimaciones de los parámetros
confint(fit_auto_t)
##           2.5 %     97.5 %
## ma1  -0.5775267 -0.2261293
## sma1 -0.7002176 -0.4136721

en ambos casos el cero no pertenece a los intervalos de confianza por lo que ambos parámetros son estadísticamente significativos.

Criterios de comparación

Consideremos la información del ajuste del modelo hecho a los datos transformados

fit_auto_t
## Series: log(AirPassengers) 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 = 0.001371:  log likelihood = 244.7
## AIC=-483.4   AICc=-483.21   BIC=-474.77

donde

  • \(AIC=-2 log(verosimilitud)+2p\), con \(p\) el número de parámetros estimados. Buscamos un valor menor de \(AIC\).

Dichos valores, como el AIC y BIC, sirven para elegir el mejor modelo (los cuales deben cumplir los supuestos).

También podemos ver

summary(fit_auto_t)
## Series: log(AirPassengers) 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 = 0.001371:  log likelihood = 244.7
## AIC=-483.4   AICc=-483.21   BIC=-474.77
## 
## Training set error measures:
##                        ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.0005730622 0.03504883 0.02626034 0.01098898 0.4752815 0.2169522
##                    ACF1
## Training set 0.01443892

donde \(MAE\) (mean absolute error). Estos valores también nos servirán para elegir el mejor modelo.

  • Un buen modelo: buscamos donde el error sea menor.

Comparación de modelos

El modelo que hallamos para los datos transformados no es el único, puede hallarse algunos otros modelos para dicha serie:

# Modelo hallado
fit1 <- arima(log(AirPassengers), c(0,1,1), seasonal = list(order = c(0,1,1), period = 12))

# Alternativas
fit2 <- arima(log(AirPassengers), c(0,1,1), seasonal = list(order = c(2,1,2), period = 12))
fit3 <- arima(log(AirPassengers), c(2,1,1), seasonal = list(order = c(0,1,0), period = 12))

Para los ajustes 2 y 3 debemos repetir el proceso que llevamos a cabo para el ajuste 1, es decir, debemos repetir el análisis de los supuesto, para posteriormente comparar y ver cual de los modelos es el mejor utilizando los criterios de comparación.