library(readr)
## Warning: package 'readr' was built under R version 3.3.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 3.3.3
library(forecast)
## Warning: package 'forecast' was built under R version 3.3.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 3.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.3.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(TSA)
##
## Attaching package: 'TSA'
## The following object is masked from 'package:readr':
##
## spec
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
library(urca)
## Warning: package 'urca' was built under R version 3.3.3
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:forecast':
##
## autolayer
library(rmarkdown)
library(zoo)
metro<- read_csv("C:/Users/Dkevog/Desktop/PASAJEROS.csv")
## Parsed with column specification:
## cols(
## Periodo = col_character(),
## Ptrans = col_double()
## )
View(metro)
st<-ts(metro$Ptrans,frequency = 12, start = c(1986,1))
Serie:“Pasajeros transportados en metro de la CDMX”" Es el promedio diario de pasajeros transportados. Incluye a los pasajeros con boleto pagado y usuarios de acceso gratutito (personas de la tercera edad, Discapacitados, menores de 5 a?os y derechohabientes del Sistema de Transporte Colectivo, principalmente).
-Unidad de medida: Miles de pasajeros -Periodicidad: Mensual -Fecha inicial: Enero/1986 -Fecha final: Junio/ 2018 -Ultima actualizacion 2018/08/20 Fuente:(http://www.inegi.org.mx/sistemas/bie/) -Ruta:Comunicaciones y transportes> Principales caracteristicas del sistema de transporte colectivo metro de la Ciudad de Mexico.> Pasajeros transportados
La serie de los pasajeros transportados en el metro de la CDMX nos presenta patrones de estacionalidad y tendencia ciclo creciente
autoplot(st, main="Pasajeros transportados", xlab = "A?os", ylab = "Pasajeros")
En el siguiente grafico se pudede visualizar cada componente de la serie, los datos originales, la tendencia, la estacionalidad, y los residuos
fit<- decompose(st, type = "additive")
autoplot(fit)
En esta grafica se aprecia mas el patron de estacionalidad, se visualizan los meses en los que se transportan mas pasajeros
ggseasonplot(st)
Se utilizo logaritmo para estabilizar la varianza y se aplica la transformacion Box-Cox
logst <-log(st)
autoplot(logst,main="Pasajeros transportados", xlab = "Años", ylab = "Pasajeros")
BoxCox.ar(logst)
Se aplica diferencia tipica a la diferencia estacional para eliminar la tendencia y la estacionalidad
ddlogst<-diff(diff(logst,12))
autoplot(ddlogst)
summary(ur.df(logst))
##
## ###############################################
## # 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.34395 -0.03125 0.00687 0.03525 0.20649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 7.981e-05 3.214e-04 0.248 0.804
## z.diff.lag -4.391e-01 4.622e-02 -9.501 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05206 on 380 degrees of freedom
## Multiple R-squared: 0.192, Adjusted R-squared: 0.1877
## F-statistic: 45.14 on 2 and 380 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: 0.2483
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
H0=No es estacionaria H1=Es estacionaria
En base al valor P se rechaza la hiptesis nula, por lo que nos arroja que la serie ya es estacionaria
ggtsdisplay(diff(diff(logst,12)))
Se crean las siguientes propuestas
propuesta1 <- Arima(logst, order=c(1,1,0), seasonal=c(1,1,0)) # AR(1)
propuesta1
## Series: logst
## ARIMA(1,1,0)(1,1,0)[12]
##
## Coefficients:
## ar1 sar1
## -0.5047 -0.5356
## s.e. 0.0448 0.0527
##
## sigma^2 estimated as 0.002421: log likelihood=589.99
## AIC=-1173.97 AICc=-1173.91 BIC=-1162.22
checkresiduals(propuesta1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(1,1,0)[12]
## Q* = 116.12, df = 22, p-value = 8.771e-15
##
## Model df: 2. Total lags used: 24
propuesta2<-Arima(logst, order=c(0,1,1), seasonal=c(0,1,1))#MA1
propuesta2
## Series: logst
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.6603 -0.8914
## s.e. 0.0369 0.0327
##
## sigma^2 estimated as 0.001634: log likelihood=655.31
## AIC=-1304.62 AICc=-1304.55 BIC=-1292.87
checkresiduals(propuesta2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 24.192, df = 22, p-value = 0.3372
##
## Model df: 2. Total lags used: 24
propuesta3<-Arima(logst, order=c(0,1,2), seasonal=c(0,1,2)) #MA2
propuesta3
## Series: logst
## ARIMA(0,1,2)(0,1,2)[12]
##
## Coefficients:
## ma1 ma2 sma1 sma2
## -0.7178 0.0886 -0.9752 0.0939
## s.e. 0.0536 0.0565 0.0660 0.0632
##
## sigma^2 estimated as 0.001621: log likelihood=657.79
## AIC=-1305.58 AICc=-1305.42 BIC=-1286
checkresiduals(propuesta3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(0,1,2)[12]
## Q* = 20.985, df = 20, p-value = 0.398
##
## Model df: 4. Total lags used: 24
Con el analisis de los residuos de las propuestas, la que mas se ajusta a nuestro modelo es lapropuesta 3 con un MA2
propuesta4<-auto.arima(logst)
propuesta4
## Series: logst
## ARIMA(0,1,1)(0,0,2)[12]
##
## Coefficients:
## ma1 sma1 sma2
## -0.6763 0.1702 0.3283
## s.e. 0.0366 0.0556 0.0618
##
## sigma^2 estimated as 0.00198: log likelihood=648.35
## AIC=-1288.7 AICc=-1288.59 BIC=-1272.91
checkresiduals(propuesta4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,0,2)[12]
## Q* = 30.619, df = 21, p-value = 0.08023
##
## Model df: 3. Total lags used: 24
Aunque el valor de AICc de la propuesta 4 es bajo, la propuesta 3 es la que mas se ajusta a nuestro modelo
\[ Y_t=-0.7178y_{t-1} + 0.886y_{t-2} -0.9752y_{t-12} +0.939y_{t-13} \]
Se realiza el pronostico a dos años, donde se puede apreciar que mantiene el patron de estacionalidad constante
autoplot(forecast(logst, h=24), xlab="años", ylab="cantidad de pasajeros transportados", main= "Pronostico a dos años")