library(readr)
library(gridExtra)
library(forecast)
library(fpp2)
## Loading required package: ggplot2
## Loading required package: fma
## Loading required package: expsmooth
library(lmtest)
## Loading required package: zoo
## 
## 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)
library(ggplot2)
library(rmarkdown)
library(zoo)
library(tseries)
library(vars)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following objects are masked from 'package:fma':
## 
##     cement, housing, petrol
## Loading required package: strucchange
## Loading required package: sandwich

variable original y declarada serie de tiempo

pasajeros <- read_csv("C:/Users/family/Desktop/pasajeros.csv")
## Parsed with column specification:
## cols(
##   Periodo = col_character(),
##   Ptrans = col_double()
## )
stp<-ts(pasajeros$Ptrans,frequency = 12, start = c(1986,1))

variable elegida y declarada serie de tiempo

trenes <- read_csv("C:/Users/family/Desktop/trenes.csv")
## Parsed with column specification:
## cols(
##   Periodo = col_character(),
##   trenes = col_integer()
## )
stt<-ts(trenes$trenes,frequency = 12, start = c(1986,1))

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

Serie:“Pasajeros transportados en metro de la CDMX”" Es el promedio diario de pasajeros transportado que 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: septiembre/2018

-Ultima actualizacion 2018/11/21

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

2. Graficar la serie de tiempo, analizar patrones y observaciones atipicas. Apoye su analisis con una descomposicion clasica.

Nuestra serie presenta patrones de estacioanalidad con una tendencia ciclo creciente

autoplot(stp, main="Pasajeros transportados", xlab = "Años", ylab = "Pasajeros")

En el siguiente gráfico se pudede visualizar cada componente de la serie, los datos originales, la tendencia, la estacionalidad, y los residuos, se utilizo la descomposicipon aditva

fit<- decompose(stp, type = "additive")
autoplot(fit)

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

Se utilizó logaritmo para estabilizar la varianza y se aplicó la transformación Box-Cox

logstp <-log(stp)
autoplot(logstp,main="Pasajeros transportados", xlab = "Años", ylab = "Pasajeros")

BoxCox.ar(logstp)

4. Modelar la serie de tiempo elegida como funcion de una tendencia, tendencia cuadratica y/o medias estacionales. Presente el diagnostico del modelo elegido.

Se modela la serie en funcion del tiempo

mes.<-season(stp)
modelo1<-lm(stp~mes.+time(stp)+ I(time(stp)^2))
summary(modelo1)
## 
## Call:
## lm(formula = stp ~ mes. + time(stp) + I(time(stp)^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1373.45  -124.28     7.69   106.74   647.21 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     6.628e+06  5.387e+05  12.303  < 2e-16 ***
## mes.February    1.724e+02  5.240e+01   3.290  0.00110 ** 
## mes.March       8.672e+01  5.240e+01   1.655  0.09875 .  
## mes.April      -1.471e+01  5.240e+01  -0.281  0.77915    
## mes.May         9.599e+01  5.240e+01   1.832  0.06775 .  
## mes.June        1.497e+02  5.240e+01   2.856  0.00452 ** 
## mes.July        4.859e+01  5.240e+01   0.927  0.35436    
## mes.August      2.146e+02  5.241e+01   4.095 5.16e-05 ***
## mes.September   1.398e+02  5.241e+01   2.667  0.00798 ** 
## mes.October     3.243e+02  5.281e+01   6.141 2.07e-09 ***
## mes.November    2.454e+02  5.282e+01   4.647 4.65e-06 ***
## mes.December    1.647e+01  5.282e+01   0.312  0.75540    
## time(stp)      -6.631e+03  5.381e+02 -12.323  < 2e-16 ***
## I(time(stp)^2)  1.659e+00  1.344e-01  12.350  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 212.9 on 379 degrees of freedom
## Multiple R-squared:  0.512,  Adjusted R-squared:  0.4953 
## F-statistic: 30.59 on 13 and 379 DF,  p-value: < 2.2e-16

La R cuadrad nos demuestra que el modelo se ajusta en un 49% y el intercepto explica que en enero en promedio se transportan 6.775e+06 miles de pasajeros en metro

library(ggplot2) 
autoplot(stp) + geom_smooth(method="lm", se=FALSE)

autoplot(stp) + geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Se puedeapreciar a que a pesar de los cambios que tiene la serie origina, presenta una tendencia creciente

4.1 Presentar y analizar gráfico de los residuos: gráfica en el tiempo, grafico cuantil-cuantil, histograma y funcion de autocorrelacioon (correlograma).

a)Gráfico de los residuos

plot(y=rstandard(modelo1), x=time(stp), type='l')

###b)Gráfico cuantil cuantil No existe normalidad ya que no se encuentran la mayoría de los datos dentro de la línea de normalidad

qqnorm(rstandard(modelo1)); qqline(rstandard(modelo1))

c)Histograma

El siguente gráfico demuestra un sesgo

hist(rstandard(modelo1),type="o")
## Warning in plot.window(xlim, ylim, "", ...): graphical parameter "type" is
## obsolete
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## graphical parameter "type" is obsolete
## Warning in axis(1, ...): graphical parameter "type" is obsolete
## Warning in axis(2, ...): graphical parameter "type" is obsolete

d)Autocorrelacion

El siguiente gráfico nos demuestra que existe autocorrelacion entre los datos

acf(rstandard(modelo1))

e) Evaluar la normalidad de los residuos ocupando la prueba Shapiro-Wilk.

Ho=No hay normalidad Ha=Hay normalidad Se acepta la hipótesis nula de que no hay normalidad ya que el valor P es mayor a .05

shapiro.test(rstandard(modelo1))
## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(modelo1)
## W = 0.95958, p-value = 6.317e-09

5. En caso de estacionalidad de la serie de tiempo, aplicar diferencias estacionales. Presentar grafica.

Se aplica diferencia tipica a la diferencia estacional para eliminar la tendencia y la estacionalidad

ddlogstp<-diff(diff(logstp,12))
autoplot(ddlogstp)

#6. prueba Dickey-Fuller para evaluar el orden de integraci´on de la serie. Diferenciar hasta que la serie sea estacionaria.

summary(ur.df(logstp))
## 
## ############################################### 
## # 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.34440 -0.03128  0.00738  0.03520  0.20603 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## z.lag.1     9.027e-05  3.164e-04   0.285    0.776    
## z.diff.lag -4.404e-01  4.554e-02  -9.671   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05186 on 389 degrees of freedom
## Multiple R-squared:  0.1939, Adjusted R-squared:  0.1897 
## F-statistic: 46.78 on 2 and 389 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: 0.2853 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

H0=No es estacionaria H1=Es estacionaria

La serie es no estacionaria

7 y 8.Usar herramientas ACF, PACF, EACF y criterios de Akaike/ Bayes para construir propuestas de modelos. Como parte del diagnostico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

ggtsdisplay(diff(diff(logstp,12)))

Se crean las siguientes propuestas ##Propuesta 1: AR 1

propuesta1 <- Arima(logstp, order=c(1,1,0), seasonal=c(1,1,0)) # AR(1)
propuesta1
## Series: logstp 
## ARIMA(1,1,0)(1,1,0)[12] 
## 
## Coefficients:
##           ar1     sar1
##       -0.5054  -0.5896
## s.e.   0.0447   0.0470
## 
## sigma^2 estimated as 0.00242:  log likelihood=603.81
## AIC=-1201.63   AICc=-1201.56   BIC=-1189.81
checkresiduals(propuesta1)

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

Propuesta 2: MA1

propuesta2<-Arima(logstp, order=c(0,1,1), seasonal=c(0,1,1))#MA1
propuesta2
## Series: logstp 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.6573  -0.8928
## s.e.   0.0370   0.0306
## 
## sigma^2 estimated as 0.001616:  log likelihood=673.44
## AIC=-1340.89   AICc=-1340.82   BIC=-1329.07
checkresiduals(propuesta2)

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

Propuesta 3: MA2

propuesta3<-Arima(logstp, order=c(0,1,2), seasonal=c(0,1,2)) #MA2
propuesta3
## Series: logstp 
## ARIMA(0,1,2)(0,1,2)[12] 
## 
## Coefficients:
##           ma1     ma2     sma1    sma2
##       -0.7119  0.0894  -0.9762  0.1004
## s.e.   0.0524  0.0562   0.0553  0.0540
## 
## sigma^2 estimated as 0.0016:  log likelihood=676.55
## AIC=-1343.1   AICc=-1342.94   BIC=-1323.4
checkresiduals(propuesta3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,1,2)[12]
## Q* = 23.191, df = 20, p-value = 0.2795
## 
## 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

9.Usar la funcion auto.arima() y compare sus resultados con el inciso anterior

propuesta4<-auto.arima(logstp)
propuesta4
## Series: logstp 
## ARIMA(0,1,1)(0,0,2)[12] 
## 
## Coefficients:
##           ma1    sma1    sma2
##       -0.6758  0.1396  0.3393
## s.e.   0.0364  0.0501  0.0561
## 
## sigma^2 estimated as 0.001954:  log likelihood=666.1
## AIC=-1324.2   AICc=-1324.1   BIC=-1308.31
checkresiduals(propuesta4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,0,2)[12]
## Q* = 30.239, df = 21, p-value = 0.08728
## 
## 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

10.Presente la ecuacion final y describa.

\[ Y_t=-0.7119y_{t-1} + 0.894y_{t-2} -0.9762y_{t-12} +0.1004y_{t-13} \]

11. Crear un pronostico a dos años.

Se realiza el pronostico a dos años, donde se puede apreciar que mantiene el patron de estacionalidad constante

autoplot(forecast(logstp, h=24), xlab="años", ylab="cantidad de pasajeros transportados", main= "Pronostico a dos años")

12. Utilizando metodos de pronosticos simples y/o suavizamiento exponencial generar un pronostico para dos anos. Elegir el m´etodo que presente la ra´iz del error medio al cuadrado m´as pequeno.

stp2 <- window(stp,start=1986,end=c(2017,9))

Se crea la grafica

autoplot(stp) +
forecast::autolayer(meanf(stp2, h=11), PI=FALSE, series="Mean") +
forecast::autolayer(naive(stp2, h=11), PI=FALSE, series="Naïve") +
forecast::autolayer(snaive(stp2, h=11), PI=FALSE, series="Seasonal naïve") +
ggtitle("Pasajeros trasnsportados en metro de laCDMX") +
xlab("Años") + ylab("Pasaeros") +
guides(colour=guide_legend(title="Forecast"))

PrOnosticos simples

stpfit1 <- meanf(stp2,h=10)
stpfit2 <- rwf(stp2,h=10)
stpfit3 <- snaive(stp2,h=10)
autoplot(window(stp, start=1986)) + forecast::autolayer(stpfit1, series="Mean", PI=FALSE) + forecast::autolayer(stpfit2, series="Naïve", PI=FALSE) + forecast::autolayer(stpfit3, series="Seasonal naïve", PI=FALSE) + xlab("Año") + ylab("Pasajeros") + guides(colour=guide_legend(title="Forecast"))

stp3 <- window(stp, start=1986)

accuracy (stpfit1,stp3)#mean
##                        ME     RMSE      MAE        MPE     MAPE     MASE
## Training set 1.204340e-13 295.1918 226.1572 -0.5391072 5.663983 1.194219
## Test set     3.650237e+02 391.0060 365.0237  8.3067246 8.306725 1.927500
##                     ACF1 Theil's U
## Training set  0.71909061        NA
## Test set     -0.03574802  1.938152
accuracy (stpfit2,stp3)#naive
##                       ME      RMSE       MAE        MPE      MAPE     MASE
## Training set   -1.207932  215.7945  163.0011 -0.1878968  4.137075 0.860724
## Test set     1246.943352 1254.7953 1246.9434 28.6288975 28.628897 6.584460
##                     ACF1 Theil's U
## Training set -0.37824403        NA
## Test set     -0.03574802  6.418325
accuracy (stpfit3,stp3)#seasonal naive
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  18.68603 249.2875 189.3767  0.2581491 4.752692 1.0000000
## Test set     -47.52562 158.5519 112.8833 -1.1366098 2.611300 0.5960781
##                    ACF1 Theil's U
## Training set  0.4739069        NA
## Test set     -0.4476061 0.8506243

La raíz del promedio de los errores al cuadrado (RMSE) de menor valor= 158 fue la de Seasonal Naive, por lo que es el pronóstico que más se adapta a nuestro modelo

14. Analizar un VAR, necesita incluir los siguientes puntos:

a)Buscar otra variable(incluso más variable)y redactar la teoría que se encuentra detrás de la relación de variables que propone.

Cifra relativa al número de trenes en operación en día laborable y horario de máxima demanda

Ruta temática: Comunicaciones y transportes > Principales características del sistema de transporte colectivo metro de la Ciudad de México. > Trenes en servicio Periodicidad: Mensual Unidad de medida: Número de unidades Fuente: Gobierno de la Ciudad de México. Sistema de Transporte Colectivo Metro. Cifras preliminares a partir de: 2017/01 Nota:
Datos no acumulables.

Fecha inicial:
1986/01 Fecha final:
2018/09 Última actualización: 2018/11/21

La variable que se analiza es la cantidad de trenes del Sistema colectivo metro en la CDMX, para así ver que relación tienen con la cantidad de pasajeros transportados en el Metro de la CDMX y que tanto pudeden depender una de la otra

b)Evaluar la estacionalidad de la nueva variable

ts.stp <- diff (diff(logstp,12)) # diferencia de pasajeros
ts.stt <- diff(diff(stt,12)) # diferencia de trenes


base2 <- cbind.zoo(P = ts.stp, T = ts.stt)

View(base2)

autoplot(ts.stp)

autoplot(ts.stt)

Podemos observar que tanto la serie de trenes y de pasajeros transportados en el metro de la CDMX presentan estacionalidad , para esto se aplico el método de diferencias de las diferencia estacional porque además tienen tendencia ,
“No hubo necesidad de crear otra base, ya que, las variables tienen la misma longitud”

c)Proponga el orden de rezagos óptimo para el VAR y presentar la estimación.

VARselect(base2, lag.max = 12, type = "const")[["selection"]]
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     12     12      2     12
VARselect(base2, lag.max = 12, type = "const")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     12     12      2     12 
## 
## $criteria
##                 1          2          3          4         5          6
## AIC(n) -1.5263317 -1.6674556 -1.6731877 -1.6625388 -1.653958 -1.6453446
## HQ(n)  -1.5010169 -1.6252642 -1.6141198 -1.5865943 -1.561137 -1.5356470
## SC(n)  -1.4626129 -1.5612577 -1.5245107 -1.4713825 -1.420323 -1.3692301
## FPE(n)  0.2173316  0.1887273  0.1876497  0.1896606  0.191298  0.1929574
##                7          8          9         10         11         12
## AIC(n) -1.633488 -1.6467962 -1.6315378 -1.6321438 -1.7833246 -2.0626450
## HQ(n)  -1.506914 -1.5033455 -1.4712106 -1.4549400 -1.5892443 -1.8516881
## SC(n)  -1.314894 -1.2857233 -1.2279858 -1.1861126 -1.2948143 -1.5316555
## FPE(n)  0.195265  0.1926915  0.1956645  0.1955585  0.1681333  0.1271706

El número máximo de rezagos son 12.

Se crean las siguientes propuestas:

VARselect(base2[,1:2], lag.max=12, type="const")[["selection"]]
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     12     12      2     12

VAR1

var1 <- VAR(base2[,1:2], p=1, type="const")
serial.test(var1, lags.pt=10, type="PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 88.063, df = 36, p-value = 2.973e-06

VAR2

var2 <- VAR(base2[,1:2], p=2, type="const")
serial.test(var2, lags.pt=10, type="PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var2
## Chi-squared = 46.316, df = 32, p-value = 0.04878

VAR3

var3 <- VAR(base2[,1:2], p=3, type="const")
serial.test(var3, lags.pt=10, type="PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var3
## Chi-squared = 33.086, df = 28, p-value = 0.2326

VAR12

var12 <- VAR(base2[,1:2], p=12, type="const")
serial.test(var12, lags.pt=10, type="PT.asymptotic")
## Warning in pchisq(STATISTIC, df = PARAMETER): NaNs produced

## Warning in pchisq(STATISTIC, df = PARAMETER): NaNs produced
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var12
## Chi-squared = 14.441, df = -8, p-value = NA

Se tomará el modelo del VAR3 Se analizan los rezagos

d)Probar la correlación serial y normalidad en los resisduos

correlacion serial

var3.serial <- serial.test(var3, lags.pt=12, type="PT.asymptotic")
plot(var3.serial, names = "Pasajeros")
## Warning in plot.varcheck(var3.serial, names = "Pasajeros"): 
## Invalid residual name(s) supplied, using residuals of first variable.

HO=No hay correlacion serial

No se puede rechazar la hipotesis nula de que no hay correlacion serial

normalidad

var.norm <- normality.test(var3, multivariate.only = TRUE)
var.norm
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object var3
## Chi-squared = 3131, df = 4, p-value < 2.2e-16
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object var3
## Chi-squared = 2.2728, df = 2, p-value = 0.321
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object var3
## Chi-squared = 3128.7, df = 2, p-value < 2.2e-16

En base al valor P nos que existe normalidad

e)Resultados de la prueba de causalidad de Granger

causality(var3, cause='P')$Granger
## 
##  Granger causality H0: P do not Granger-cause T
## 
## data:  VAR object var3
## F-Test = 1.1739, df1 = 3, df2 = 740, p-value = 0.3187

Claramente no es posible rechazar la hipótesis nula de no Granger causalidad.

causality(var3, cause='T')$Granger
## 
##  Granger causality H0: T do not Granger-cause P
## 
## data:  VAR object var3
## F-Test = 4.3935, df1 = 3, df2 = 740, p-value = 0.004481

Se rechaza la hipotesis nula de no Granger causalidad

f)Crear un pronóstico para dos años usando el VAR especificado.

predictions <- predict(var3, n.ahead = 10, ci = 0.95)
class(predictions)
## [1] "varprd"
plot(predictions, names = "Pasajeros ")
## Warning in plot.varprd(predictions, names = "Pasajeros "): 
## Invalid variable name(s) supplied, using first variable.

plot(predictions,names="P")

plot(predictions,names="T")

g)Interpretar las funciones impulso respuesta.

plot(irf(var3, impulse = "P", response = "T", ortho=T))

plot(irf(var3, impulse = "T", response = "P", ortho=T))

En el caso de los pasajeros, se muestra como la cantidad de trenes responden ante un impulso de la cantidad de pasajeros , por lo tanto, los pasajeros infiere en la cantidad de trenes.

En el caso de los trenes, nos muestra un resultado similar, entonces, un impulso en los pasajeros corresponde a una respuesta de la cantidad de los trenes

15. Evaluar la potencial cointegración usando la metodología de Engle-Granger y Johansen

par(mfrow=c(1,2))
plot(ts.stp,main="Pasajeros")
plot(ts.stt,main="Trenes")

PRUEBA

po.test(base2[,1:2])
## Warning in po.test(base2[, 1:2]): p-value smaller than printed p-value
## 
##  Phillips-Ouliaris Cointegration Test
## 
## data:  base2[, 1:2]
## Phillips-Ouliaris demeaned = -499.96, Truncation lag parameter =
## 3, p-value = 0.01

Con H0 : las variables no están cointegradas.

Se rechaza la hipotesis nula es decir las variables están cointegradas