Modelos Arima

1. Breve referencia de la serie de tiempo elegida, así como frecuencia, unidad de medida, fuente, etc.
2. Graficar los datos, analizar patrones y observaciones atìpicas. Apoye su anàlisis con una descomposiciòn clàsica.
3. Si es necesario, utilizar transformaciòn Box-Cox / logaritmos para estabilizar la varianza.
- Transformación Box-Cox
- Transformación Logaritmica
4. En caso de estacionalidad, aplicar diferencias estacionales.
- Diferencias Estacionales
5. Usar prueba Dickey-Fuller para evaluar el orden de integraciòn de la serie. Diferenciar hasta que la serie sea estacionaria.
6. Usar herramientas ACF, PACF, EACF y criterios de Akaike/ Bayes para construir propuestas de modelos.
7. Como parte del diagnòstico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.
- Analizar residuos
- Propuesta Ljung-Box
8. Usar la funciòn auto.arima() y compare sus resultados con el inciso anterior.
9. Presente la ecuaciòn final y describa. 10. Crear un pronòstico a dos años.

1. Breve referencia de la serie de tiempo elegida, así como frecuencia, unidad de medida, fuente, etc.

Se trabaja con una base de datos para estimar el comportamiento de movimiento de pasajeros por vía aérea en aeropuertos nacionales

Frecuencia:
- La frecuencia que presenta la serie es mensual
Unidad de medida:
- Número de pasajeros

Enlace en línea

Ruta temática:
- Comunicaciones y transportes > Movimiento de pasajeros por vía aérea en aeropuertos > Nacionales > Total
Periodicidad:
- Mensual
Unidad de medida:
- Número de pasajeros
Fuente:
- Aeropuertos y Servicios Auxiliares (ASA). Dirección General Adjuna de Finanzas y Operación. Gerencia Informática.
Fecha inicial:
- 2005/01
Fecha final:
- 2018/05
Última actualización:
- 2018/06/28
Enlace en línea

library(forecast)
library(TSA)

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

library(urca)
library(ggplot2)

pasajeros = read.csv ("C:/Users/ACER/Documents/Series de tiempo/pasajeros.csv")
View(pasajeros)

Declarar base de datos en serie de tiempo ST

pas<-ts(pasajeros$Numero.de.Pasajeros, frequency=12, start=c(2005,1))

autoplot(pas,frequency=12,xlab="Años",ylab="Volumen",main="Pasajeros")

2. Graficar los datos, analizar patrones y observaciones atìpicas. Apoye su anàlisis con una descomposiciòn clàsica.

library(gridExtra)
fitA<-decompose(pas,type="additive")
autoplot(fitA, xlab="Años",ylab="Volumen",main="Pasajeros")

Se muestra en la serie la presencia de patrones como son:estacionalidad no se puede ser visible con exactitud pero presenta influencia de factores estacionales, tendecia ciclo debido a que en el largo plazo se visuliza un incremento en los datos.

library(forecast)
seasonplot(pas, col=rainbow(12), year.labels=TRUE)

La serie de tiempo presenta estacionalidad en los meses como se observa en la gráfica, los meses de Julio son con mayor incremento de vuelos.

ggseasonplot(pas,polar=TRUE,main="Pasajeros" )

Observemos en estas dos graficas con más claridad lo mencionado anteriormente.

3. Si es necesario, utilizar transformaciòn Box-Cox / logaritmos para estabilizar la varianza.

Transformación Box-Cox

BoxCox.ar(pas)

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

Transformación Logaritmica

Pasajeroslog <- log(pas)
plot(Pasajeroslog,xlab="Años",ylab="Volumen",type="o",main="Pasajeros",col="darkred")

4. En caso de estacionalidad, aplicar diferencias estacionales.

Diferencias Estacionales

Existe un pequeño aumento en la varianza, por lo que se toman logaritmos para estabilizar.

plot(log(pas),xlab="Años",ylab=" ", main='Logatirmo de Pasajeros',col="darkblue")

baseP<-diff(Pasajeroslog)
autoplot(baseP, xlab="Años",ylab=" ",main="Primera diferencia en log de Pasajeros")

Podemos observar la presencia de la estacionalidad aun cuando se aplica la primera diferencia

5. Usar prueba Dickey-Fuller para evaluar el orden de integraciòn de la serie. Diferenciar hasta que la serie sea estacionaria.

library(urca)
summary(ur.df(baseP, type = "none", selectlags = "AIC"))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.234571 -0.056689 -0.007543  0.096784  0.249740 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.72983    0.11300 -15.309  < 2e-16 ***
## z.diff.lag  0.41161    0.07209   5.709 5.57e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1007 on 156 degrees of freedom
## Multiple R-squared:  0.6833, Adjusted R-squared:  0.6793 
## F-statistic: 168.3 on 2 and 156 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -15.3085 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Se rechaza la Ho, la serie ya es estacionaria, -15.3085 es mayor a los valores critico taul.

6. Usar herramientas ACF, PACF, EACF y criterios de Akaike/ Bayes para construir propuestas de modelos.

ACF

ggAcf(baseP)

Esta grafica muestra que sería mejor una MA(1)

PACF

ggPacf(baseP)

Esta grafica la posbilidad de utilizar un ar(2)

EACF

eacf(baseP)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o x x x o x o x o  x  o  x 
## 1 x x o o o x x x o x o  x  x  x 
## 2 x x o o o x x o o o o  x  x  x 
## 3 o x o x o x o x o o o  x  o  x 
## 4 o o o x o x o x o o o  x  x  o 
## 5 x x x x o o o x o o o  x  x  o 
## 6 o x o x o o o x o o o  x  x  x 
## 7 x x o x x o x x o o o  x  x  x

Esta grafica no nos muestra un mejor modelo

ggtsdisplay(baseP)

Criterios de Akaike/Bayes

res <- armasubsets(y=diff(log(pas)),nar=7,nma=7,y.name='test',ar.method='ols')

plot(res)

La anterior gràfica nos muestra al parece ser un ar(1)

PROPUESTAS DE MODELO MA(1) y AR(1)

PROPUESTA 1 ARIMA(1)

propuesta1<-arima(baseP,order = c(1,1,0), seasonal=c(1,1,0)) 
propuesta1

## 
## Call:
## arima(x = baseP, order = c(1, 1, 0), seasonal = c(1, 1, 0))
## 
## Coefficients:
##           ar1     sar1
##       -0.6484  -0.3842
## s.e.   0.0625   0.0740
## 
## sigma^2 estimated as 0.004589:  log likelihood = 185.91,  aic = -367.82

Tenemos un AIC=-367.82

checkresiduals(propuesta1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)(1,1,0)[12]
## Q* = 113.81, df = 22, p-value = 2.287e-14
## 
## Model df: 2.   Total lags used: 24

PROPUESTA 2 MA(1) probamos con MA(1) estacional y no estacional

propuesta2<-arima(baseP,order = c(0,1,1), seasonal=c(0,1,1)) 
propuesta2

## 
## Call:
## arima(x = baseP, order = c(0, 1, 1), seasonal = c(0, 1, 1))
## 
## Coefficients:
##          ma1     sma1
##       -1.000  -0.8565
## s.e.   0.048   0.0843
## 
## sigma^2 estimated as 0.001996:  log likelihood = 236.73,  aic = -469.46

Tenemos un AIC=-469.46 aun maás pequeño que el anterior

checkresiduals(propuesta2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 49.775, df = 22, p-value = 0.000629
## 
## Model df: 2.   Total lags used: 24

PROPUESTA 3 ARIMA(1) probamos con MA(1) estacional y no estacional

propuesta3<-arima(baseP,order = c(1,1,1), seasonal=c(1,1,1)) 
propuesta3

## 
## Call:
## arima(x = baseP, order = c(1, 1, 1), seasonal = c(1, 1, 1))
## 
## Coefficients:
##           ar1      ma1    sar1     sma1
##       -0.3034  -0.9735  0.0290  -0.7963
## s.e.   0.0824   0.0412  0.1134   0.0974
## 
## sigma^2 estimated as 0.001914:  log likelihood = 243.06,  aic = -478.13

Tenemos un AIC=-478.13 aun más pequeño que el anterior, pero no pasa la pruebas

checkresiduals(propuesta3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(1,1,1)[12]
## Q* = 33.616, df = 20, p-value = 0.02884
## 
## Model df: 4.   Total lags used: 24

No pasa las pruebas, se mantiene la propuesta 2

7. Como parte del diagnòstico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

Analizar residuos

ggAcf(residuals(arima((baseP), order = c(0,1,1), seasonal=c(0,1,1))))

Propuesta Ljung-Box

checkresiduals(propuesta2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 49.775, df = 22, p-value = 0.000629
## 
## Model df: 2.   Total lags used: 24

Se acepta Ho, existe correlacion serial en los residuos

8. Usar la funciòn auto.arima() y compare sus resultados con el inciso anterior.

propuesta.auto<-auto.arima(baseP)
propuesta.auto

## Series: baseP 
## ARIMA(2,0,1)(1,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2      ma1    sar1     sma1  drift
##       0.2195  0.0657  -0.5541  0.0415  -0.8255  0e+00
## s.e.  0.3378  0.1459   0.3242  0.1135   0.1034  1e-04
## 
## sigma^2 estimated as 0.001903:  log likelihood=250.13
## AIC=-486.26   AICc=-485.46   BIC=-465.28

checkresiduals(propuesta.auto)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1)(1,1,1)[12] with drift
## Q* = 29.412, df = 18, p-value = 0.04357
## 
## Model df: 6.   Total lags used: 24

Como lo habiamos dicho anteriormente en las propuestas, la mejor opción estaba entre un MA(1) con un AR(2); donde la función auto.arima nos sugiere un ARIMA(2,0,1). Se decide mantener la propuesta 2 aunque se tiene autocorrelación se preeve que incluso auto arima nos marca la autocorrelación