e.g ARIMA
Dante Bernedo
5th January 2022
Synopsis
This objective of this analysis and modelling is to review time series theory and experiment with R packages.
We will be following an ARIMA modeling procedure of the AirPassengers dataset as follows:
1. Perform exploratory data analysis
2. Decomposition of data
3. Test the stationarity
4. Fit a model used an automated algorithm
5. Calculate forecasts
Revision del DATAFRAME
La base de datos utilizada es: Airpassangers
AirPassengers## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432AP<-AirPassengers
class(AP)## [1] "ts"The AirPassenger dataset in R provides monthly totals of a US airline passengers, from 1949 to 1960. This dataset is already of a time series class therefore no further class or date manipulation is required.
Analisis descritivo de la serie:
sum(is.na(AP))## [1] 0frequency(AP)## [1] 12cycle(AP)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 1 2 3 4 5 6 7 8 9 10 11 12
## 1950 1 2 3 4 5 6 7 8 9 10 11 12
## 1951 1 2 3 4 5 6 7 8 9 10 11 12
## 1952 1 2 3 4 5 6 7 8 9 10 11 12
## 1953 1 2 3 4 5 6 7 8 9 10 11 12
## 1954 1 2 3 4 5 6 7 8 9 10 11 12
## 1955 1 2 3 4 5 6 7 8 9 10 11 12
## 1956 1 2 3 4 5 6 7 8 9 10 11 12
## 1957 1 2 3 4 5 6 7 8 9 10 11 12
## 1958 1 2 3 4 5 6 7 8 9 10 11 12
## 1959 1 2 3 4 5 6 7 8 9 10 11 12
## 1960 1 2 3 4 5 6 7 8 9 10 11 12summary(AP)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 104.0 180.0 265.5 280.3 360.5 622.0Grafico de la serie
Graficando los datos para ver los patrones de la serie, se puede observar que desde 1940 empieza a subir la cantidad de viajeros con el tiempo, tambien se hace mas variable esa cantidad de personas que viajan pero siempre incrementandose y se observa picos que se ubica cada cierto tiempo a lo largo de la serie, lo que nos indica seguramente la presencia de estacionalidad.
plot(AP,xlab="Date",type="o", ylab = "Passengers",main="Air Passenger numbers from 1949 to 1961")
hist(AP, probability = TRUE)
lines(density(AP), col = "red")
qqnorm(AP)
qqline(AP, col = "red")
boxplot(AP~cycle(AP),xlab="Date", ylab = "Passenger Numbers (1000's)" ,main ="Monthly Air Passengers Boxplot from 1949 to 1961")
Revisando la Autocorrelacion
El acf, le decimos que lo haga con un lag (retraso)de 12,24,48 el que sea adecuado dependiendo de la longitud de la serie. Si tenemos una serie de 24 elementos le ponemos un lag de 12 o de 6, pero si tenemos una serie mas amplia podemos usar lags de 48. Si son unidades de tiempo e.g horas y minutos nos interesa verlo cada 60. Colocamos el valor de 35 solo para ejemplificar.
En el acf se observa que todas las autocorrelaciones son positivas y significativas en los diferentes lags y que van decayendo conforme avanza la serie, lo cual nos indica que no es una serie estacionaria sino, que hay una tendencia que de acuerdo al grafico de la serie es una tendencia creciente.
Algunas series no presenta la tendencia tan claramente pero el grafico acf ayuda a indicar si existe la tendencia (la tendencia subyacente)
Ojo la DF ha sido acondicionada y especificado varable tiempo, ya que los lags en el eje x salen unidades de tiempo (año) pq el objeto de serie de tiempo se le indico frecuencia años con 12 meses, si no se hace eso sale como numero de lags
La pacf son los lags de interes pero controlado por los lags de los otros valores q se dan en la serie.
El patron del grafico muestra el “lag 1” significativo,esperable, debido a que la serie tiene tendencia y los otros lags no significativos. Aunque hay algunos picos significativos en el valor 1.2 que en este caso seria cada 12 meses, lo que se notaba en los datos como unos picos estacionales. Para poder seguir buscando el modelo adecuado como vimos que hay tendencia lo que podemos hacer es una diferenciacion de un lag en los datos. Para inspeccionar graficamos nuevamente los datos solamente que antes le aplicamo el comando “diff” e.g plot(diff(AP, xlab=“tiempo”,ylab=“valores”)) que es (Yt-Yt-1).
require(tseries)## Loading required package: tseries## Warning: package 'tseries' was built under R version 4.0.5## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoorequire(forecast)## Loading required package: forecast## Warning: package 'forecast' was built under R version 4.0.5acf(AP,lag=35, col="2", main="", xlab="lag", ylab="valores")
pacf(AP,lag=35, col="2", main="", xlab="lag", ylab="valores")
## Diferenciando la serie
Se observa que ya no se muestra el incremento en la serie, pq se ha eliminado al quitarlo con la diferenciacion que hace q el promedio sea constante a lo largo de la serie. Pero notamos que la varianza se va incrementando a lo largo del tiempo eso no se corrige con la diferenciacion. La acf de la serie diferenciada, se observa que ya no todas las autocorrealciones en los diferentes rezagos son postivas y significativas. Si no que ahora estan espaciodos en el tiempo los que tienen significancia positivos y negativos.
La pacf diferenciada presentaba un patron donde la varianza incrementaba con el tiempo es algo que se ve en el grafico y no necesariamente va a verse en la acf o pacf.
plot(diff(AP, xlab="tiempo",ylab="valores"))
acf(diff(AP),lag=35, col="2", main="", xlab="lag", ylab="valores")
pacf(diff(AP),lag=35, col="2", main="", xlab="lag", ylab="valores")
# Aplicando logaritmo a la serie Como la varianza se incrementa en el tiempo debemos sacar logaritmo natural a la serie. y volverla a graficar nuevamente pero antes aplicamos logaritmo natural a la serie original lo que va homegeinizar la varianza. “plot(diff(log(AP)), xlab=”tiempo“,ylab=”valores“)” Luego de aplicado el comando se nota ya que no se nota el patron de incremento de la varianza en el tiempo, que tambien podria ser un patron de creciemitno de la varianza en el tiempo o a media serie. Quitamos la heterocedisticidad de los datos. Las acp que se obtienen siguen presentando los picos que se deben a la estacionalidad que hay en los datos. La pacf tiene el mismo comportamiento tenemmos que analizarlo para ver el aspecto de la autocorrelacion acf y pacf para ver en cual de los modelos calza si es Regresivo, Media moviles o ARIMA.
plot(diff(log(AP)), xlab="tiempo",ylab="valores")
acf(diff(log(AP)),lag=35, col="2", main="", xlab="lag", ylab="valores")
pacf(diff(log(AP)),lag=35, col="2", main="", xlab="lag", ylab="valores")
# Test Dickey - Fuller El test sirve para ver no estacionariedad, devuelve la probabilidad de de ser no estacionario. Si es > 0.05 (La serie es no estacionaria, ie tiene tendencia) Si es < 0.05 (La serie es estacionaria,ie No hay tendencia)
- Chequeamos que nuestra serie con logaritmo y diferenciada sea estacionaria, ie que no tenga tendencia.
The null hypothesis H0 : that the time series is non stationary The alternative hypothesis H1 : that the time series is stationary
Como p=0.01 < 0.005 por lo tanto rechazamos la hipotesis nula de no estacionariedad
adf.test(diff(log(AP)),alternative = "stationary") ## Warning in adf.test(diff(log(AP)), alternative = "stationary"): p-value smaller
## than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff(log(AP))
## Dickey-Fuller = -6.4313, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryrenombrando la serie diferenciada y logaritmo
aerial<-diff(log(AP))Ajuste automatico de los coeficientes
Usamos el comando AUTOARIMA. 1. Hacemos el codigo con la serie con logaritmo El rsultado del modelo es: ARIMA(0,1,1)(0,1,1)[12], que nos indica que es lo que debe aplicarse a la serie. A este se le aplico logaritmo pq se sabia segun el grafico que no habia homogeneidad de varianzas. Parte regular(p,d,q),Parte estacional(P,D,Q)
ARIMA(0,1,1)(0,1,1)[12] ie. No hay componente autorregresivo, hay que aplicarle una difereciacion que ya dedejumimos del analisis grafico. Tiene una media movil, ademas tiene un componente estacional que seria el segundo parentesis a este componente estacional tambien se le hace una diferenciacion cada 12 meses y se le hace una media movil.
[12] son los valores de los coeficientes que estan en la formula del modelo con su error estandard. Lo importante aqui es cuando se quiere comprar varios modelos la verosmilitud (loglikehood)o el AIC q sirve tambien para comparaciones que se quieran hacer con otros modelos.
ARIMAmodelLOG<-auto.arima(log(AP))
ARIMAmodelLOG## Series: log(AP)
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.4018 -0.5569
## s.e. 0.0896 0.0731
##
## sigma^2 estimated as 0.001371: log likelihood=244.7
## AIC=-483.4 AICc=-483.21 BIC=-474.77Serie Diferenciada
Vamos a hacer el modelo tambien para demostrar que vamos a obtener el mismo tipo de resultado, solo que con la serie que estraba diferencia(log y diferenciada). El modelo resultante es: ARIMA(0,0,1)(0,1,1)[12] El modelo solo ocupa una media movil y no hay valor de diferenciacion pq ya se le habia aplicado a la serie anterior la estacionalidad cada 12 meses por una media movil y si se ve los valores de los coeficientes pues son los mismos con el mismo error estandar y el de verosimilitud y AIC igual la misma magnitud. Se ha realizado para que se note que se puede hacer con la serie sin diferenciar y obtener los mismos valores.
Pero no nos quita que antes hayamos hecho los graficos para entender como se esta comportando los datos en el cambioque tienen en el tiempo.
ARIMAmodelDIFF<-auto.arima(aerial)
ARIMAmodelDIFF## Series: aerial
## ARIMA(0,0,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.4018 -0.5569
## s.e. 0.0896 0.0731
##
## sigma^2 estimated as 0.001369: log likelihood=244.7
## AIC=-483.39 AICc=-483.2 BIC=-474.77Raices unitarias
Estos nos va ayudar a ver si los modelos son los apropiados, la funcion grafica para los componente AR y los componentes de MA. En este caso como el modelo no tiene componente AR no sale el grafico, pero para MA si sale el grafico mas los componente estacionales que tambien tenia media movil, todos tienen que quedar dentro del circulo, pero es mas pro ejemplificar pq siempre que se utilice el comando AUTOARIMA no va a devolvernos un modelo que no cumpla esa situacion. Si se hace en forma manual el modelo si se tiene que correr estos graficos para asegurar que esta bien. ARIMA(0,1,1)(0,1,1)[12]
require(forecast)
autoplot(ARIMAmodelDIFF, title = "Raices invertidas sobre AR y MA")
Diagnostico del modelo (metodo 1)
Como los modelos lo que busca es contar con uan serie que tenga ruido blanco, tenemos que hacer un diagnostico del modelo. Los graficos que genera explican como se comportaron los residuos a lo largo del tiempo, una autocorrelacion de los residuales donde vemos que solo el lag 0, ie donde no hay lag es significativo pero el resto donde ya tenemos rezago no hay ningun pico que sea significativo lo caul no indica que tenemos una serie con ruido blanco en los residuales y para asegurar se hace una prueba Ljung-Box donde todos sus valores tiene que estar fuera del intervalo donde se hace significativo por lo cual no hay presencia de autocorrelacion en esta serie de residuales cumpliendose lo que necesitamos de supuestos para que el modelo sea posible usar para proyeccion. Aqui corremos la prueba de Ljung-Box generalizada no por lag si no que tome en cuenta 10 lags en general para ver como sale la probabilidad. Como se oberva en los resultados las probabilidades son muy altas que indican que la serie no tiene ningun tipo de autocorrelacion y que es de ruido blanco.
# Model diagnostic
# Method 1
tsdiag(ARIMAmodelLOG)
Box.test(ARIMAmodelLOG$residuals,lag = 10,type = "Ljung-Box") ##
## Box-Ljung test
##
## data: ARIMAmodelLOG$residuals
## X-squared = 8.8097, df = 10, p-value = 0.5503Diagnostico del modelo (metodo 2)
Este metodo es chekresiduals, que basicamente va a hacer lo mismo, nos muestra la informacion del modelo que estamos probando del grafico de los residuales. Hace la autocorrelacion de los residuales. Se observa que hay un lag significativo pero esta ubicado bastante lejos en los rezagos y tiene un valor no tan alto de correlacion. Tambien genera un histograma de los residuales que vemos que se acerca a la normalidad bastante bien. Tambien genera otro grafico de Ljung-Box basada en otro metodo pero igual resulta no significativa por lo tanto la serie de los residuos muestra ruido blanco lo que nos indica que es apropiado para proyecciones.
# Method 2
checkresiduals(ARIMAmodelLOG)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 26.446, df = 22, p-value = 0.233
##
## Model df: 2. Total lags used: 24Observado vs esperados del modelo
Podemos visualizar como nuestro modelo se ve a la par de los datos observados que teniamos y esto lo hacemos con el comando plot para x que seria los datos de nuestro modelo y con color rojo para diferenciarlo. Acordarse que este es la serie que presenta logaritmo y ponemos una linea adicional para los datos que el modelo predijo tambien cuando se le aplicaba logaritmos se observa como caen los datos en azul esperados versus los datos observados que estaban en color rojo y los valores que obviamente se ve que no calzan son los que analizamos como excedentes o deficit re4specto al modelo en los residuales pero sabemos que cumplen ya un patron de reuido blanco. Entonces no nos afecta mucho en el uso del modelo simplemente lo vemos para ejemplificarlo un poco.
plot(ARIMAmodelDIFF$x,col="red")
lines(fitted(ARIMAmodelDIFF), col="blue")
predicciones
Especificamos la cantidad de fechas que queremos proyectar, en este caso lo hacemos para un año, entonces le decimos 12 fechas que serian los meses. el comando “forecast” nos devuelve los valores de la proyeccion que en este caso serian valores diferenciados con un logaritmo aplicado y nos muestra intervalos de confianza para los datos a 80% y 95%. En este ejemplo seria para 1961 por que nuestra serie termina en 1960. Para graficarlo usamos el comando autoplot. La salida de estos comandos se puede ver como se comporta el modelo y como se proyecta a un año siguiente, lo que va a ocurrir con esos viajeros.
Tambien hemos corrido para ejemplificar el codigo con el modelo sin diferenciar pero al que se le aplico logaritmo. EL resultado muestra un grafico del logaritmo de los pasajeros y como este varia el tiempo este es mas facil para despues convertir los datos en cantida de pasajeros por que al que simplemente se le aplicaria antilogaritmo natural o e^al valor, para que se convierta en numero de pasajeros.
# prediccion
# 1. modelo diferenciado
forecast(ARIMAmodelDIFF, h=12)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1961 0.041760005 -0.005657453 0.089177463 -0.03075874 0.11427875
## Feb 1961 -0.056410776 -0.107513080 -0.005308471 -0.13456501 0.02174346
## Mar 1961 0.117938868 0.066836564 0.169041173 0.03978464 0.19609310
## Apr 1961 0.027586656 -0.023515649 0.078688960 -0.05056758 0.10574089
## May 1961 0.033255644 -0.017846661 0.084357949 -0.04489859 0.11140987
## Jun 1961 0.136222424 0.085120119 0.187324729 0.05806819 0.21437665
## Jul 1961 0.138515503 0.087413198 0.189617808 0.06036127 0.21666973
## Aug 1961 -0.004387551 -0.055489856 0.046714754 -0.08254178 0.07376668
## Sep 1961 -0.178208533 -0.229310837 -0.127106228 -0.25636276 -0.10005430
## Oct 1961 -0.115689890 -0.166792195 -0.064587585 -0.19384412 -0.03753566
## Nov 1961 -0.145520903 -0.196623208 -0.094418598 -0.22367513 -0.06736667
## Dec 1961 0.104537303 0.053434998 0.155639608 0.02638307 0.18269153autoplot(forecast::forecast(ARIMAmodelDIFF,h=12))
# 2. modelo sin diferenciar
forecast(ARIMAmodelLOG, h=12)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1961 6.110186 6.062729 6.157642 6.037607 6.182764
## Feb 1961 6.053775 5.998476 6.109074 5.969203 6.138347
## Mar 1961 6.171715 6.109555 6.233874 6.076650 6.266779
## Apr 1961 6.199300 6.130966 6.267635 6.094792 6.303809
## May 1961 6.232556 6.158560 6.306552 6.119388 6.345724
## Jun 1961 6.368779 6.289524 6.448033 6.247569 6.489988
## Jul 1961 6.507294 6.423109 6.591479 6.378544 6.636044
## Aug 1961 6.502906 6.414064 6.591749 6.367034 6.638779
## Sep 1961 6.324698 6.231431 6.417965 6.182058 6.467338
## Oct 1961 6.209008 6.111516 6.306500 6.059908 6.358109
## Nov 1961 6.063487 5.961947 6.165028 5.908195 6.218780
## Dec 1961 6.168025 6.062591 6.273459 6.006778 6.329272autoplot(forecast::forecast(ARIMAmodelLOG,h=12))
### Convertiendo el modelo de proyeccion en dataframe Le colocamos el nombre al modelo de proyeccion como airpass nada mas para simplificar y si queremos solo el dato de los proyectado como cantida de individuos no como logaritmo, llamamos este punto de proyeccion en la tabla y aca seria el vector de la cantidad de viajeros en los siguientes 12 meses que se proyectaron y de ahi se podria hacer el grafico de forma en cantidad de personas se pueden transformar tambien los intervalos de confianza. Todo esto es el cuidado que hay que tener si la serie se le aplica o no el logaritmo, pq si no se hace la inspeccion se podria correr el modelo autoarima sin aplicar el logaritmo y tendriamos mas dificultades que el modelo ajuste puede ser que eso ocurra o no ocurra en ciertos casos. Si eso es el caso el modelo debe ser corrido con el logaritmo para homogenizar la varianza, ya despues los datos se pueden volver a la escala original.
airpass<-as.data.frame(forecast(ARIMAmodelLOG,h=12))
exp(airpass$`Point Forecast`)## [1] 450.4224 425.7172 479.0068 492.4045 509.0550 583.3449 670.0108 667.0776
## [9] 558.1894 497.2078 429.8720 477.2426