Series de Tiempo:Modelos de regresión/ARIMA/VAR/Cointegración
Pasajeros del STCM
Estefanía Rodríguez Torres & César Espindola Palestino
Introducción
Texto
MODELO ARIMA
1. Breve referencia a los datos, frecuencia, unidad de medida, fuente, etc.
Las series de tiempo que se emplearan corresponden al número de pasajeros del Sistema de Transporte Colectivo Metro (STCM) de la Ciudad de México y la distancia medidad en kilometros.
La serie comienza en enero de 1986 y culmina en febrero de 2018, con una frecuencia de datos mensual por lo que se cuentan con 386 observaciones. La serie se mide en miles de personas, a partir del cálculo del promedio de pasajeros transportados diariamente por dicho sistema de transporte, contempla en su cálculo a los pasajeros con boleto pagado y usuarios de acceso gratuito (personas de la tercera edad, discapacitados, menores de 5 años y derechohabientes del STCM).
Los datos fueron obtenidos del Banco de Información Estadística del Instituto Nacional de Estadística, Geografía e Informática, sin embargo, la fuente original de la serie es del Gobierno de la Ciudad de México, Sistema de Transporte Colectivo Metro.
Para comprender mejor el comportamiento de los datos que se tiene registro es necesario conocer un poco de la historia del sistema de transporte metropolitano más grande del país, pues su creación se explica en parte de la gran expansión de la mancha urbana en la Zona Metropolitana del Valle de México (ZMVM) a partir de los años 70´s y a su vez el nivel de demanda que ha alcanzado en la actualidad. El STCM comenzó actividades en 1969 con tan solo 16 estaciones y una línea, para el año en el que se comienza a tener registro (1986), ya se contaba con 7 líneas y 91 estaciones, para fines del siglo pasado se construyeron 4 líneas más y 58 estaciones, y para los años que van del siglo XXI el incremento ha sido de tan solo 28 estaciones y una línea, sumando un total de 12 líneas y 195 estaciones para satisfacer la demanda de 20,886,703 personas que habitan la ZMVM actualmente. (GOBIERNO CDMX, INEGI,2017).
2. Graficar los datos, analizar patrones y observaciones atípicas. Apoye su análisis con una descomposición clásica.
I Gráfica de la serie original
autoplot(ts.p)+ylab("Número de Pasajeros Promedio")+xlab("Años")+ggtitle("Pasajeros del STCM 1986-2018")
II Descomposición Clásica
fit.m<-decompose(ts.p,type = 'multiplicative')autoplot(fit.m)
En el gráfico I y II se puede analizar el comportamiento de la serie de pasajeros del STCM a durante 32 años, en dichos grafico logra apreciarse que la varianza de la serie original no muestra cambios abruptos, así mismo esto se aprecia en el comportamiento de los rezagos, con excepción en la observación 381 correspondiente al promedio de pasajeros de septiembre del 2017, cabe mencionar que este es el punto más bajo en la serie, debido a que el Gobierno de la Ciudad de México tomó como medida el uso gratuito de este transporte del 19 al 28 de septiembre ,después del terremoto ocurrido en el mismo mes.
La serie también da indicios de que existe estacionalidad, lo cual tiene sentido ya que el servicio de transporte disminuye su uso en periodo de vacaciones, así mismo hay presencia de tendencia, se muestra que a partir del año de 2010 la demanda de este medio de transporte aumenta, siendo 2013 el año con más demanda del transporte y siendo 1998 el más bajo en cuanto a demanda.
3. Si es necesario, utilizar transformación Box-Cox / logaritmos para estabilizar la varianza.
BoxCox.ar(ts.p)
Aplicando la prueba Box-Cox el valor máximo de lambda,??? es muy cercano a 1, lo cual indica que no es necesario transformar la serie y se puede trabajar con la serie original, ya que una transformación de tipo logarítmica no tendría significancia en la estabilidad de varianza, puesto que la serie original ya se distribuye normalmente.
4. Modelar la serie de tiempo elegida como función de una tendencia, tendencua cuadrática y/o medias estacionales. Presente el diagnostico del modelo elegido.
Modelo lineal
modelo1<-lm(pasajeros~time(pasajeros))
summary(modelo1)##
## Call:
## lm(formula = pasajeros ~ time(pasajeros))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1105.89 -171.31 -5.12 159.75 842.55
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3758.8551 27.2845 137.765 <2e-16 ***
## time(pasajeros) 1.1662 0.1219 9.568 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.9 on 385 degrees of freedom
## Multiple R-squared: 0.1921, Adjusted R-squared: 0.19
## F-statistic: 91.56 on 1 and 385 DF, p-value: < 2.2e-16Modelo Cuadratico
modelo2 <- lm(log(pasajeros) ~ time(pasajeros) + I(time(pasajeros)^2))
summary(modelo2)##
## Call:
## lm(formula = log(pasajeros) ~ time(pasajeros) + I(time(pasajeros)^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36697 -0.03624 0.00390 0.03445 0.16789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.306e+00 8.992e-03 923.723 < 2e-16 ***
## time(pasajeros) -8.451e-04 1.070e-04 -7.896 3.05e-14 ***
## I(time(pasajeros)^2) 2.901e-06 2.671e-07 10.859 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05866 on 384 degrees of freedom
## Multiple R-squared: 0.3729, Adjusted R-squared: 0.3696
## F-statistic: 114.2 on 2 and 384 DF, p-value: < 2.2e-165. En caso de estacionalidad, aplicar diferencias estacionales.
Gráficas para detectar estacionalidad
ggseasonplot(ts.p,polar = TRUE,main='Análisis de Estacionalidad')
ggseasonplot(ts.p,year.labels = TRUE,year.labels.left = TRUE)+ylab("Número de Pasajeros")+ggtitle("Pasajeros del STCM 1986-2018")+xlab("Meses")
En los gráficos anteriores se aprecia que en enero, diciembre, marzo y julio el número de pasajeros promedio disminuye en todos los años, estos meses corresponden a los periodos vacacionales. Por otro lado, en agosto y octubre el uso del servicio tiende a aumentar en mayor medida, lo mismo pasa en menor medida a los meses que les suceden a los periodos vacacionales, lo cual indica que el STCM puede tener como principal objetivo trasladar a trabajadores y estudiantes para la realización de sus actividades diarias.
Diferencia Típica de la Serie Original
ggtsdisplay(diff(ts.p,1),main="Diferencia Típica")
Al aplicar diferencia típica a la serie podemos comprobar que la serie tiene estacionalidad debido a que los rezagos salen de las medias cada 12 periodos en ACF y en PACF cada 11.
Diferencia estacional de la Serie Original
ggtsdisplay(diff(ts.p,12),main="Diferencia Estacional")
Al existir un alto componente estacional aplicamos diferencia estacional, pero observamos que la media no es estable y que al bajar de forma lenta los rezagos hay tendencia, para mejorar el comportamiento de los rezagos se aplicará la diferencia de la diferencia estacional.
Diferencia típica y diferencia estacional de la Serie Original
ggtsdisplay(diff(diff(ts.p,12)),main="Diferencia de la diferencia estacional")
Aplicando la diferencia de la diferencia estacional se estabilizan los rezagos respecto a la media, sin embargo los rezagos el primer rezago en ACF y PACF salen de las bandas de ajuste, por lo que se proponen modelos AR(1) y MA(1) estacional y no estacional.
6. 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(ts.p,type = "trend",selectlags = "AIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1238.54 -126.77 28.45 137.61 680.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 941.34841 156.62790 6.010 4.34e-09 ***
## z.lag.1 -0.24915 0.04131 -6.031 3.85e-09 ***
## tt 0.27966 0.10194 2.743 0.00637 **
## z.diff.lag -0.31340 0.04867 -6.439 3.63e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 196.4 on 381 degrees of freedom
## Multiple R-squared: 0.262, Adjusted R-squared: 0.2561
## F-statistic: 45.08 on 3 and 381 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -6.0313 12.1562 18.1925
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36summary(ur.df(ts.p,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
## -1273.78 -118.79 29.92 143.13 779.01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.0002016 0.0026167 -0.077 0.939
## z.diff.lag -0.4378278 0.0460360 -9.511 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 205 on 383 degrees of freedom
## Multiple R-squared: 0.1914, Adjusted R-squared: 0.1872
## F-statistic: 45.32 on 2 and 383 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -0.077
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62Utilizando la prueba Dickey Fuller se obtuvo un valor p de ???2.2e.???^(-16) se rechaza la presencia de raíz unitaria, es decir que la serie sea estacionaria, por lo cual es conveniente aplicar diferencias a la serie original. De acuerdo con esto lo correcto es aplicar diferencia estacional para corregir problemas de estacionalidad y a su vez diferenciar esto para convertir la serie en estacionaria y modelar adecuadamente.
7. Usar herramientas ACF, PACF, EACF y/o criterios de Akaike/ Bayes para construir propuestas de modelos.
Propuesta de modelo AR(1) estacional y no estacional
fit1<-Arima(ts.p,order = c(1,1,0),seasonal = c(1,1,0))fit1## Series: ts.p
## ARIMA(1,1,0)(1,1,0)[12]
##
## Coefficients:
## ar1 sar1
## -0.5068 -0.5394
## s.e. 0.0445 0.0511
##
## sigma^2 estimated as 36923: log likelihood=-2498.49
## AIC=5002.98 AICc=5003.04 BIC=5014.75Propuesta de modelo MA(1) estacional y no estacional
fit2<-Arima(ts.p,order = c(0,1,1),seasonal = c(0,1,1))
fit2## Series: ts.p
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.6521 -0.8892
## s.e. 0.0372 0.0314
##
## sigma^2 estimated as 24985: log likelihood=-2432.91
## AIC=4871.82 AICc=4871.88 BIC=4883.59Propuesta considerando un modelo de Autocorrelación extendida
eacf(diff(diff(ts.p,12)))## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o x x x o
## 1 x o o o o o o o o o o x x x
## 2 x x o o o o o o o o o x o o
## 3 o x o o o o o o o o o x x o
## 4 x x o o o o o o o o o x o x
## 5 x x o o o o o o o o o x x o
## 6 x x x o o o o o o o o x o o
## 7 o x o o o o o o o o o x o oGráfico de Criterio de Bayes
Bayes<-armasubsets(y=diff(diff(ts.p,12)),nar = 12,nma = 12,y.name = "test",ar.method = 'ols')
plot(Bayes)
Utilizando las herramientas ACF, PACF y EACF se sugiere un modelo MA (1) estacional y no estacional con un criterio de Akaike (AIC) de 4859.83 que es menor la AIC de un modelo AR(1), así mismo considerando el criterio de Bayes (BIC) se recomienda usar un MA(1).
8. Como parte del diagn´ostico del modelo, analizar los residuos y aplicar la prueba Ljung-Box
Modelo AR(1) estacional y no estacional
checkresiduals(fit1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(1,1,0)[12]
## Q* = 119.3, df = 22, p-value = 2.331e-15
##
## Model df: 2. Total lags used: 24El modelo AR(1) estacional y no estacional los rezagos sobrepasan las líneas de bondad de ajuste, aunado a esto, de acuerdo con la prueba Ljung-Box se obtiene un valor muy pequeño, próximo a cero, por lo cual se rechaza que este modelo no tenga correlación serial, por lo cual es preferible no utilizar este modelo, ya que al no cumplir con los criterios básicos como la no correlación y tener además un AIC y BIC mayor al MA(1), es un modelo mal especificado.
Modelo MA (1) estacional y no estacional
checkresiduals(fit2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 26.645, df = 22, p-value = 0.225
##
## Model df: 2. Total lags used: 24El modelo MA(1) estacional y no estacional recomendado usar de acuerdo a las herramientas ACF y EACF, mantiene sus residuos dentro de las bandas de ajuste, además su valor p en la prueba Ljung-Box arroja un valor de 0.23, con lo cual se acepta que el modelo no tiene correlación serial y es adecuado utilizar este modelo. Así mismo los residuos se distribuyen de forma similar a una distribución normal, aunque esto último puede ser mejorable.
9. Usar la función auto.arima() y compare sus resultados con el inciso anterior.
fit3<-auto.arima(ts.p)
fit3## Series: ts.p
## ARIMA(0,1,3)(0,0,2)[12]
##
## Coefficients:
## ma1 ma2 ma3 sma1 sma2
## -0.7152 0.0516 0.0792 0.1806 0.3667
## s.e. 0.0512 0.0653 0.0514 0.0544 0.0590
##
## sigma^2 estimated as 30282: log likelihood=-2538.76
## AIC=5089.53 AICc=5089.75 BIC=5113.26checkresiduals(fit3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,3)(0,0,2)[12]
## Q* = 25.86, df = 19, p-value = 0.1341
##
## Model df: 5. Total lags used: 24Introduciendo la función Auto.arima en R, se recomienda usar un modelo MA (3) no estacional y un MA(2), con diferencia únicamente típica. Los resultado de este modelo dan prueba que no existe correlación serial al tener un valor p de 0.13 y hay una distribución normal mejor que el modelo MA(1) sin embargo hay rezagos que salen de la banda de ajuste, y el criterio de AIC es mayor que el modelo propuesto con anterioridad
Auto.arima considerando D=1: AR(1)MA(2) no estacional y MA(2)I estacional
fit4<-auto.arima(ts.p,D=1,stepwise = FALSE,approximation = FALSE)
fit4## Series: ts.p
## ARIMA(1,0,2)(0,1,2)[12]
##
## Coefficients:
## ar1 ma1 ma2 sma1 sma2
## 0.9700 -0.6973 0.1053 -0.9628 0.0905
## s.e. 0.0162 0.0544 0.0569 0.0659 0.0619
##
## sigma^2 estimated as 24443: log likelihood=-2433.11
## AIC=4878.21 AICc=4878.44 BIC=4901.77checkresiduals(fit4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(0,1,2)[12]
## Q* = 21.412, df = 19, p-value = 0.3145
##
## Model df: 5. Total lags used: 24Finalmente con la función auto.arima con una diferención fija se recomienda utilizar un modelo ARIMA(1,0,2) no estacional y sin diferencia y un MA(2) estacional con diferencia estacional. El modelo recomendado cumple con la no correlación lineal al tener un valor p de 0.31 en la prueba Ljung-Box.Los rezagos se mantienen dentro de las bandas de ajuste del modelo en un ACF y se logra una distribución normal que mejora en comparación a los modelos anteriores.
De acuerdo a los criterios de Akaike y Bayes el modelo más conveniente es el MA (1) estacional y no estacional, sin embargo el modelo ARIMA (1,0,2)(0,1,2)12 también cumple con los criterios básicos para modelar una serie de tiempo, por lo anterior es conveniente realizar un gráfico cuantil cuantil en el que se compare la normalidad de ambos modelos para definir el mejor modelo.
Pruebas de normalidad
qqnorm(residuals(fit2));qqline(residuals(fit2))
qqnorm(residuals(fit4));qqline(residuals(fit4))
Con ambos modelos se logra una distribución normal, por una ligera diferencia el ultimo modelo resulta mas conveniente para la serie de Pasajeros del STCM.
10. Presente la ecuación final y describa.
La ecuación del modelo corresponde a un modelo estacional multiplicativo ARIMA (p,d,q)x(P,D,Q)12 , el cual se caracteriza por ser parsimonioso
11.Realizar el pronÓtico ARIMA para dos años.
getrmse <- function(x,h,...)
{
train.end <- time(x)[length(x)-h]
test.start <- time(x)[length(x)-h+1]
train <- window(x,end=train.end)
test <- window(x,start=test.start)
fit5 <- Arima(train,...)
fc <- forecast(fit5,h=h)
return(accuracy(fc,test)[2,"RMSE"])
}
getrmse(ts.p, h=24, order=c(1,1,0), seasonal=c(1,1,0))## [1] 516.0217getrmse(ts.p, h=24, order=c(1,1,1), seasonal=c(1,1,0))## [1] 438.4455getrmse(ts.p, h=24, order=c(1,1,0), seasonal=c(0,1,0))## [1] 577.1318getrmse(ts.p, h=24, order=c(2,1,0), seasonal=c(1,1,0))## [1] 453.0043getrmse(ts.p, h=24, order=c(2,1,1), seasonal=c(1,1,0))## [1] 447.0899getrmse(ts.p, h=24, order=c(2,1,1), seasonal=c(1,1,1))## [1] 356.3151getrmse(ts.p, h=24, order=c(2,1,1), seasonal=c(0,1,1))## [1] 361.6042getrmse(ts.p, h=24, order=c(1,0,0), seasonal=c(1,1,0))## [1] 366.5824getrmse(ts.p, h=24, order=c(1,0,1), seasonal=c(1,1,0))## [1] 381.4193getrmse(ts.p, h=24, order=c(1,0,1), seasonal=c(1,1,1))## [1] 353.3097getrmse(ts.p, h=24, order=c(2,0,1), seasonal=c(0,1,1))## [1] 361.2164getrmse(ts.p, h=24, order=c(1,1,0), seasonal=c(1,0,0))## [1] 478.8014getrmse(ts.p, h=24, order=c(1,1,1), seasonal=c(0,0,1))## [1] 366.6083Pronosticando varios modelos ARIMA se elige el que tenga el valor AICc más bajo dado que el modelo con el AICc más bajo tiende a tener ligeramente mejores resultados. Por lo tanto, se considera utilizar el modelo (2, 0, 1) (0, 1, 1) el cual tiene un valor de 309.8313. Estimando el modelo ARIMA propuesto se obtienen los siguientes resultados:
fit6<-Arima(ts.p,order = c(2,0,1),seasonal = c(0,1,1))
fit6## Series: ts.p
## ARIMA(2,0,1)(0,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 sma1
## 0.8281 0.1367 -0.5492 -0.8827
## s.e. 0.0819 0.0752 0.0690 0.0327
##
## sigma^2 estimated as 24579: log likelihood=-2434.55
## AIC=4879.11 AICc=4879.27 BIC=4898.74Gráfica del pronostico
autoplot(forecast(fit6))
El pronóstico permite predecir los valores del número de pasajeros en los próximos dos años, de acuerdo con el grafico la serie se mantendrá con el mismo ritmo de crecimiento que el de los últimos años.
13. Usando el análisis de transferencia, modelar observaciones cambios estructurales y/o observaciones atípicas en el contexto ARIMAX.
MODELO VAR
14. Analizar modelo VAR considerando:
a)Buscar otra variable y redactar la teoría que se encuentra detrás de al relación entre la variable que se propone.
En base al análisis realizado con modelos ARIMA respecto a la variable de número de pasajeros del Sistema de Transporte Colectivo Metro (STCM) de la Ciudad de México, dado que la serie se mide en miles de personas, a partir del cálculo del promedio de pasajeros transportados diariamente por dicho sistema de transporte, contempla en su cálculo a los pasajeros con boleto pagado y usuarios de acceso gratuito (personas de la tercera edad, discapacitados, menores de 5 años y derechohabientes del STCM), se busca demostrar que la serie no se explica solamente por el pasado de la misma variable.
Por consiguiente se pude intuir que la mejor variable independiente para explicar el número de pasajeros dentro del STCM es el precio de este sistema de transporte, sin embargo, debido a que a al inicio de operaciones, de 1969 a 1986 este sistema de transporte puso a la venta Planillas que contenían 5 boletos, que eran de un color distinto al boleto unitario, para diferenciar su costo.
De abril de 1986 a diciembre de 1995, el organismo sacó a la venta el tan cotizado Abono de transporte, el cual era un boleto especial único, que permitía el acceso al servicio cuantas veces fuera necesario por el mismo costo. Posteriormente, de agosto a diciembre de 1996 vendió otro Abono de transporte y de 1996 a 1998, dio paso nuevamente a las Planillas, las cuales contenían 25 boletos cada una, debido a esto no es posible unificar el precio del boleto.
Por lo cual es mejor optar por el número de kilómetros recorridos en miles de kilómetros el cual es un promedio diario de kilómetros recorridos que al igual que el número de pasajeros es calculado por el mismo sistema, cabe mencionar que esta variable cuenta con la misma temporalidad y frecuencia que la variable dependiente.
De manera que para sustentar esta relación se propone usar como base teórica el análisis de la oferta y demanda en el transporte, siendo el número de pasajeros la demanda y el número de kilómetros la oferta, debido a que es la capacidad que el STCM ha tenido y tiene para ofrecer en la actualidad. Estos son los dos lados de un mismo fenómeno que hemos reconocido como “mercado” de servicios de transporte.
En el caso del transporte una función de demanda muestra, por ejemplo, un número de pasajeros deseando utilizar un servicio, así mismo, para el caso de una empresa que ofrece un servicio de transporte de pasajeros, la función de servicio estará dada por la cantidad de trenes-kilómetro ofrecidos a determinada tarifa.
b) Evaluar la estacionalidad de la nueva variable.
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericggseasonplot(ts.km,polar = TRUE)
ts.ddP <- diff(diff(ts.p, 12))
ts.dkm <- diff(ts.km)En el grafico anterior se muestra que la variable de kilómetros recorridos en miles de kilómetros en los primeros años tiene mayor estacionalidad en los meses de abril, diciembre y enero, debido a que en estos años aún se ampliaban las líneas del sistema, mientras el tiempo avanza observamos que la serie mantiene el mismo acercamiento hacia todos los meses, siendo los últimos donde tiene una cierta tendencia en enero, abril, junio y octubre; a lo cual podemos decir que tiene que ver con la construcción de la línea 12 del STCM.
c) Proponga el orden de rezagos óptimo para el VAR y presentar la estimación.
base2 <- cbind.zoo(pasajeros = ts.ddP, km = ts.dkm)
View(base2)
base3 <- base2[-c(1:12),]
View(base3)
plot(base3)
library(vars)## Loading required package: MASS## Loading required package: strucchange## Loading required package: sandwich## Loading required package: lmtestVARselect(base3, lag.max = 11)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 11 2 2 11
##
## $criteria
## 1 2 3 4 5
## AIC(n) 12.21191 12.07611 12.05935 12.07026 12.07711
## HQ(n) 12.23749 12.11876 12.11905 12.14702 12.17093
## SC(n) 12.27628 12.18339 12.20954 12.26337 12.31314
## FPE(n) 201170.14791 175626.34130 172707.92036 174604.79569 175808.09347
## 6 7 8 9 10
## AIC(n) 12.07393 12.08878 12.10366 12.11046 12.09311
## HQ(n) 12.18480 12.21672 12.24865 12.27251 12.27222
## SC(n) 12.35286 12.41063 12.46842 12.51814 12.54370
## FPE(n) 175253.34126 177881.98386 180555.79119 181798.65432 178682.68963
## 11
## AIC(n) 11.98572
## HQ(n) 12.18189
## SC(n) 12.47923
## FPE(n) 160502.10491var11 <- VAR(base3, p=11)
var11##
## VAR Estimation Results:
## =======================
##
## Estimated coefficients for equation pasajeros:
## ==============================================
## Call:
## pasajeros = pasajeros.l1 + km.l1 + pasajeros.l2 + km.l2 + pasajeros.l3 + km.l3 + pasajeros.l4 + km.l4 + pasajeros.l5 + km.l5 + pasajeros.l6 + km.l6 + pasajeros.l7 + km.l7 + pasajeros.l8 + km.l8 + pasajeros.l9 + km.l9 + pasajeros.l10 + km.l10 + pasajeros.l11 + km.l11 + const
##
## pasajeros.l1 km.l1 pasajeros.l2 km.l2 pasajeros.l3
## -0.679596787 -5.492449621 -0.388957134 -1.693044420 -0.096500176
## km.l3 pasajeros.l4 km.l4 pasajeros.l5 km.l5
## 1.395374941 0.089555158 -6.633655151 0.067688818 7.067064418
## pasajeros.l6 km.l6 pasajeros.l7 km.l7 pasajeros.l8
## -0.022321590 11.832235496 -0.068187768 2.328850529 -0.124032594
## km.l8 pasajeros.l9 km.l9 pasajeros.l10 km.l10
## 0.977846741 -0.003282094 -2.754586344 0.107244806 0.322292693
## pasajeros.l11 km.l11 const
## 0.306066289 3.012885601 -4.319154654
##
##
## Estimated coefficients for equation km:
## =======================================
## Call:
## km = pasajeros.l1 + km.l1 + pasajeros.l2 + km.l2 + pasajeros.l3 + km.l3 + pasajeros.l4 + km.l4 + pasajeros.l5 + km.l5 + pasajeros.l6 + km.l6 + pasajeros.l7 + km.l7 + pasajeros.l8 + km.l8 + pasajeros.l9 + km.l9 + pasajeros.l10 + km.l10 + pasajeros.l11 + km.l11 + const
##
## pasajeros.l1 km.l1 pasajeros.l2 km.l2 pasajeros.l3
## 1.268200e-04 -4.279040e-01 -4.318689e-04 -2.188628e-01 3.888657e-04
## km.l3 pasajeros.l4 km.l4 pasajeros.l5 km.l5
## -5.972544e-02 -4.944571e-04 3.007343e-04 -5.144108e-04 1.984564e-02
## pasajeros.l6 km.l6 pasajeros.l7 km.l7 pasajeros.l8
## -8.423668e-05 1.291233e-01 2.044915e-04 -5.040390e-02 1.224410e-03
## km.l8 pasajeros.l9 km.l9 pasajeros.l10 km.l10
## -7.599650e-02 1.362645e-03 -1.294398e-01 5.481975e-04 -2.049955e-01
## pasajeros.l11 km.l11 const
## 1.792268e-04 -1.354260e-01 2.633780e-01Se propone un VAR 11 debido a que al realizar el análisis este se encuentra con el mejor AC, el cual es de 11.98572
d) Probar la correlación serial y normalidad en los residuos.
Prueba de normalidad
var.norm <- normality.test(var11, multivariate.only = T)
var.norm## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object var11
## Chi-squared = 816.11, df = 4, p-value < 2.2e-16
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object var11
## Chi-squared = 73.954, df = 2, p-value < 2.2e-16
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object var11
## Chi-squared = 742.16, df = 2, p-value < 2.2e-16Prueba de correlación serial
var.serial <- serial.test(var11, lags.pt = 11, type="PT.asymptotic")
var.serial##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var11
## Chi-squared = 28.82, df = 0, p-value < 2.2e-16La prueba de correlación serial arroja una p-value de 2.2e-16 la cual es mayor a 0.05, por lo tanto, se acepta que hay correlacion serial.
La prueba de normalidad arroja una p-value de 2.2e-16 la cual es mayor a 0.05, por lo tanto, se acepta que la serie tiene una distribución normal.
e) Resultados de la prueba de causalidad de Granger
causality(var11, cause='pasajeros')$Granger##
## Granger causality H0: pasajeros do not Granger-cause km
##
## data: VAR object var11
## F-Test = 0.87231, df1 = 11, df2 = 680, p-value = 0.5675causality(var11, cause='km')$Granger##
## Granger causality H0: km do not Granger-cause pasajeros
##
## data: VAR object var11
## F-Test = 1.2255, df1 = 11, df2 = 680, p-value = 0.2657De acuerdo a la prueba de causalidad de granger, se rechaza la hipotesis nula, por lo cual se confirma que la distancia promedio medida por km si causa al de número de pasajeros. En otras palabras se considera que la infraetructura de un sistema de transporte causa la demanda de pasajeros.
f ) Crear un pronóstico para dos años usando el VAR especificado.
PronosticoVAR<-predict(var11,n.ahead=10,ci=0.95)
class(PronosticoVAR)## [1] "varprd"plot(PronosticoVAR,names="pasajeros")
fanchart(PronosticoVAR,names = "km")
g) Interpretar las funciones impulso-respuesta.
plot(irf(var11, impulse='pasajeros', response = 'km'))
plot(irf(var11, impulse='km', response = 'pasajeros'))
COINTEGRACIÓN
15. Evaluar la potencial cointegración usando la metodología de Engle-Granger y Johansen
par(mfrow=c(1,2))
plot(ts.p,main="Pasajeros")
plot(ts.km,main="Distancia (KM)")
Prueba Dickey Fuller
ur.df(ts.p)##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -0.077ur.df(ts.km)##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: 1.3357Las series de tiempo no son estacionarias.
library(dynlm)
lr.regkd<-dynlm(ts.p~L(ts.p)+ts.km+L(ts.km))El modelo de regresión indica que la variables convergen a largo plazo en 0.7 meses.
error<-residuals(lr.regkd)
plot(error)
ur.df(error)##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -15.2078La prueba dickey fuller indica que existe cointegración entre el número de pasajeros y la distancia en km promedio del STCM.
library(lmtest)
data<-ts.intersect(ts.p,ts.km,error)
dyn<-dynlm(d(ts.p)~L(error)+L(d(ts.km))+L(d(ts.p)),data=data)
summary(dyn)##
## Time series regression with "ts" data:
## Start = 1986(4), End = 2018(3)
##
## Call:
## dynlm(formula = d(ts.p) ~ L(error) + L(d(ts.km)) + L(d(ts.p)),
## data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1232.92 -121.44 31.79 131.71 700.57
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.6481 10.0890 0.560 0.5759
## L(error) -0.7374 0.1300 -5.671 2.81e-08 ***
## L(d(ts.km)) -25.8463 5.9626 -4.335 1.87e-05 ***
## L(d(ts.p)) 0.1994 0.1194 1.669 0.0959 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 197.2 on 380 degrees of freedom
## Multiple R-squared: 0.2567, Adjusted R-squared: 0.2509
## F-statistic: 43.76 on 3 and 380 DF, p-value: < 2.2e-16dwtest(dyn)##
## Durbin-Watson test
##
## data: dyn
## DW = 2.1854, p-value = 0.9649
## alternative hypothesis: true autocorrelation is greater than 0El valor de la prueba Durbin Watson indica que no existe correlación entre las dos variables de análisis.