8 de mayo de 2018
El motivo de la introducción de los modelos ARIMA nace del hecho de que no se puede trabajar con una serie temporal no estacionaria. Se dice que una serie es estacionaria cuando su media, varianza y autocovarianza son invariantes en el tiempo.
La mayoría de series temporales económicas no son estacionarias pero diferenciándolas un número determinado de veces la serie original se transforma en estacionaria, con lo cual ya se podría aplicar la metodología de los modelos ARIMA.
En general, se dice que una serie temporal Yt admite una representación autoregresiva integrada y de medias móviles de órdenes p , d y q respectivamente, y se denota por ARIMA(p,d,q).
\[ Y_t = c + \phi_1y_d {_ {t-1}} + \phi_p y_d {_ {tp}} + ... + \theta_1 e_ {t-1} + \theta_q e_ {tq} + e_t \]
Donde p denota el número de términos autoregresivos, d el número de veces que la serie debe ser diferenciada para hacerla estacionaria y q el número de términos de la media móvil invertible.
La construcción de los modelos se lleva de manera iterativa mediante un proceso en el que se puede distinguir cuatro etapas:
Identificación: Utilizando los datos ordenados cronológicamente se intentara sugerir un modelo que merezca la pena ser investigada. El objetivo es determinar los valores que sean apropiados para reproducir la serie de tiempo. Comprende los pasos 1, 2 y 3 del ejemplo próximo.
Análisis y diferenciación de la serie temporal. Consiste en examinar la estacionariedad, los diagramas de autocorrelación, también conocidos como ACF y PACF, y elección del orden del modelo. Comprende los pasos 4 y 5.
Ajuste de un modelo ARIMA. Obtención de los coeficientes de determinación. Paso 6.
Predicción. Una vez seleccionado el mejor modelo candidato se pueden hacer pronósticos en términos probabilísticos de los valores futuros.
El objetivo del siguiente análisis es: Trazar, examinar y preparar series para modelar. Extrae el componente de estacionalidad de la serie temporal. Probar la estacionalidad y aplicar las transformaciones apropiadas. Elija el orden de un modelo ARIMA. Pronostica la serie.
La base de datos con la que hemos realizado el análisis está extraída del INE(Instituto Nacional de Estadística) y trata sobre viajeros que pernoctan en apartamentos turísticos de diferentes Comunidades Autónomas durante 17 años.
Para realizar el ejemplo, comenzamos instalando y ejecutando las librerías necesarias:
library('ggplot2')
## Warning: package 'ggplot2' was built under R version 3.4.4
library('forecast')
## Warning: package 'forecast' was built under R version 3.4.4
library('tseries')
## Warning: package 'tseries' was built under R version 3.4.4
Lectura y preparación de los datos.
datos = read.csv2('Datos.csv', header=TRUE,
stringsAsFactors=FALSE)
View(datos)
datos$tiempo = as.Date(datos$tiempo, '%d/%m/%Y')
ggplot(datos, aes(tiempo, apartamento)) + geom_line()
+scale_x_date('month')+ylab("Daily Apartamento Checkouts")
+ xlab("")
Para eliminar valores atípico, también puede ingresar valores perdidos en la serie, si los hay.
count_ts = ts(datos[, c('apartamento')])
datos$clean_apartamento = tsclean(count_ts)
Trazamos la serie limpia usando ggplot:
ggplot() + geom_line(data = datos, aes(x = tiempo,
y = clean_apartamento))+ylab('Cleaned Apartamento Count')
Podemos tomar la media móvil semanal o mensual, suavizando la serie en algo más estable y, por lo tanto, más predecible:
datos$apartamento_ma=ma(datos$clean_apartamento,order=7) datos$apartamento_maa=ma(datos$clean_apartamento,order=30)
ggplot() +
geom_line(data=datos,aes(x=tiempo,y=clean_apartamento,
colour="Counts")) + geom_line(data=datos,aes(x=tiempo,
y=apartamento_ma, colour="Weekly Moving Average"))+
geom_line(data=datos,aes(x=tiempo,y=apartamento_maa,
colour="Monthly Moving Average"))+ylab('Apartamento Count')
Primero, calculamos el componente estacional del uso de datos stl(). A continuación,hallamos el componente estacional de la serie mediante suavizado y ajusta la serie original restando la estacionalidad.
count_ma = ts(na.omit(datos$apartamento_ma), frequency=30) decomp = stl(count_ma, s.window="periodic") deseasonal_apartamento <- seasadj(decomp)
plot(decomp)
Deducimos que en la serie se puede observar estacionalidad pero no una clara tendencia puesto que los valores crecen y decrecen continuamente.
La instalación de un modelo ARIMA requiere que la serie sea estacionaria . Se dice que una serie es estacionaria cuando su media, varianza y autocovarianza son invariantes en el tiempo.
Esta suposición tiene un sentido intuitivo: dado que ARIMA usa retardos previos de series para modelar su comportamiento.
adf.test(count_ma, alternative = "stationary")
## ## Augmented Dickey-Fuller Test ## ## data: count_ma ## Dickey-Fuller = -2.9654, Lag order = 5, p-value = 0.1713 ## alternative hypothesis: stationary
Contraste de hipótesis:
H0: No estacionaria H1: Estacionaria
Tras realizar la prueba aumentada de Dickey-Fuller (ADF), obtenemos un p-valor = 0.1713. Como el p-valor > 0.05, no rechazamos H0.
Podemos concluir que nuestra serie temporal no es estacionaria.
Continuamos con las autocorrelaciones.
Las ACF proporcionan información sobre cómo una observación influye en las siguientes.
Acf(count_ma, main='')
Las PACF proporcionan la relación directa existente entre observaciones separadas por k retardos.
Pacf(count_ma, main='')
Para realizar un modelo ARIMA, la serie temporal debe ser estacionaria. Para conseguir esta estacionariedad, la diferenciaremos.
count_d1 = diff(deseasonal_apartamento, differences = 1) plot(count_d1)
Para comprobar que la serie es, efectivamente, estacionaria, hacemos de nuevo el test aumentado de Dickey-Fuller.
adf.test(count_d1, alternative = "stationary")
## Warning in adf.test(count_d1, alternative = "stationary"): p-value smaller ## than printed p-value
## ## Augmented Dickey-Fuller Test ## ## data: count_d1 ## Dickey-Fuller = -20.091, Lag order = 5, p-value = 0.01 ## alternative hypothesis: stationary
El p-valor obtenido es de 0.01 < 0.05. Rechazamos H0, la serie temporal es estacionaria. Al trazar la serie diferenciada, vemos un patron oscilante al rededor de 0, sin una tendencia fuerte visible. Esto sugiere que la diferenciacion de los terminos del orden 1 es suficiente y debe incluirse en el modelo.
A continuación, los picos en rezagos particulares de la serie diferenciada pueden ayudar a informar la elección de p o q para nuestro modelo:
Acf(count_d1, main='ACF for Differenced Series')
Pacf(count_d1, main='PACF for Differenced Series')
Ajustamos el modelo.
auto.arima(deseasonal_apartamento, seasonal=FALSE)
## Series: deseasonal_apartamento ## ARIMA(3,1,1) ## ## Coefficients: ## ar1 ar2 ar3 ma1 ## 1.5678 -0.7359 -0.1323 -0.7964 ## s.e. 0.0757 0.1292 0.0755 0.0334 ## ## sigma^2 estimated as 368699881: log likelihood=-2358.36 ## AIC=4726.73 AICc=4727.02 BIC=4743.44
Usando la notación ARIMA presentada anteriormente, el modelo ajustado se puede escribir como:
\[ \hat Y_ {d_t} = 0.15678 Y_ {t-1} - 0.7359 Y_ {t-2} - 0.1323 Y_ {t-3} - 0.07964 e_ {t-1} + E \] donde E es un error y la serie original se diferencia con la orden 1.
Si los parámetros de orden del modelo y la estructura se especifican correctamente, no esperaríamos que existieran autocorrelaciones significativas.
Hemos ajustado un modelo que puede producir un pronóstico, pero comenzamos examinando los gráficos ACF y PACF para los residuos del modelo. Si los parámetros de orden del modelo y la estructura se especifican correctamente, no esperaríamos que existieran autocorrelaciones significativas.
modeloarima<-auto.arima(deseasonal_apartamento, seasonal=FALSE) modeloarima
## Series: deseasonal_apartamento ## ARIMA(3,1,1) ## ## Coefficients: ## ar1 ar2 ar3 ma1 ## 1.5678 -0.7359 -0.1323 -0.7964 ## s.e. 0.0757 0.1292 0.0755 0.0334 ## ## sigma^2 estimated as 368699881: log likelihood=-2358.36 ## AIC=4726.73 AICc=4727.02 BIC=4743.44
tsdisplay(residuals(modeloarima), lag.max=10, main='(3,1,1) Model Residuals')
Hay un claro patrón presente en ACF que se repite en el retardo 4. Esto sugiere que nuestro modelo puede ser mejor con una especificación diferente como p=4 o q=4.
Existen varios criterios de este tipo para comparar la calidad del ajuste en múltiples modelos. Dos de los más utilizados son los criterios de información de Akaike (AIC) y los criterios de información de Baysian (BIC). Estos criterios están estrechamente relacionados y pueden interpretarse como una estimación de cuánta información se perdería si se elige un modelo dado.
En este caso, el AIC obtenido es 4726.73.
Podemos especificar el horizonte de pronóstico h periodos por delante para que se realicen las predicciones, y usar el modelo ajustado para generar dichas predicciones:
prediccion <- forecast(modeloarima, h=30) plot(prediccion)
