R Notebook

library(readr)

Attaching package: <U+393C><U+3E31>readr<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:TSA<U+393C><U+3E32>:

    spec

MODELOS DE REGRESIÓN / ARIMA / VAR / COINTEGRACIÓN.

1. BREVE REFERENCIA DE LA SERIE DE TIEMPO ELEGIDA, ASI COMO FECUENCIA, UNIDAD DE MEDIDA, FUENTE ETC.

 - La base a elegir sera la cantidad de pasajeros que viajan en el metrobus de
   manera mensual por cada año, esto en miles de pasajeros, desde (2008-2018).

 - informacion proporcionada por:http://www.beta.inegi.org.mx/app/tabulados/?nc=100100045

 - la serie fue escogida para analizar el comportamiento de la cantidad de pasajeros que 
   utilizan el metrobus, y ver que cambios ha tenido con el paso del tiempo.

2. graficar la serie de tiempo, analizar patrones y observaciones atipicas. apoye su analisis con una descomposición clasica.

• Descomposición clasica.

   En las graficas anteriores podemos notar que la serie tiene una tendencia positiva 
   ya que cada año aumentan la cantidad de los pasajeros,En la última grafica (random)                   podemos apreciar que las variables que construyen la serie son independientes                         aleatorias e identicamente distribuidas.

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

  Notamos que al aplicar la diferencia de los logaritmos se estabiliza la variaza y la           media

cor(diff(log(ST))[-1],tasa[-1])

[1] 0.9937032

BoxCox.ar(ST,lambda = seq(-5,1.8,5))

possible convergence problem: optim gave code = 1possible convergence problem: optim gave code = 1possible convergence problem: optim gave code = 1possible convergence problem: optim gave code = 1possible convergence problem: optim gave code = 1possible convergence problem: optim gave code = 1

4. modelar la serie de tiempo elegida como función de una tendencia, tendencia cuadratica y/o medias estacionales. presente el dianostico del modelo elegido.

mes. <- season(ST)
modelo1 <- lm(ST ~ mes. + time(ST) + I(time(ST)^2))
stargazer(modelo1, type = "text")

length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changed

===============================================
                        Dependent variable:    
                    ---------------------------
                                ST             
-----------------------------------------------
mes.February                 -174.627          
                             (673.551)         
                                               
mes.March                     331.556          
                             (673.585)         
                                               
mes.April                     95.186           
                             (673.641)         
                                               
mes.May                     1,233.312*         
                             (673.718)         
                                               
mes.June                      812.601          
                             (673.819)         
                                               
mes.July                      392.164          
                             (690.831)         
                                               
mes.August                  1,431.934**        
                             (690.863)         
                                               
mes.September                -146.235          
                             (690.906)         
                                               
mes.October                 1,516.936**        
                             (690.961)         
                                               
mes.November                  291.548          
                             (691.026)         
                                               
mes.December                 -877.608          
                             (691.103)         
                                               
time(ST)                    13,535.250         
                           (69,354.810)        
                                               
I(time(ST)2)                  -2.894           
                             (17.225)          
                                               
Constant                  -15,504,975.000      
                         (69,812,714.000)      
                                               
-----------------------------------------------
Observations                    126            
R2                             0.937           
Adjusted R2                    0.930           
Residual Std. Error    1,579.592 (df = 112)    
F Statistic          128.519*** (df = 13; 112) 
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

notamos que el intercepto es negativo, y que solo dos meses son estadisticamente                       significativos por lo cual quiza no sea una serie que se esplique con el tiempo.

Y el valor de R cuadrado explica en un 93% la variacion o varianza de la variables                    explicativas de manera conjunta respecto al tiempo.

5. En caso de estacionalidad, aplicar diferencias estacionales.

En este caso como la serie de número de miles de pasajeros en el metrobus de la cdmx, no presenta comportamiento estacional, por lo tanto no se podria aplicar diferencias estacionales.

6. Usar prueba Dickey-Fuller para evaluar el orden de integracion de la serie. Difereniar hasta que la serie sea estacionaria.

summary(ur.df(ST))

############################################### 
# 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 
-7859.0  -602.6    54.7   969.8  4742.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
z.lag.1     0.012805   0.008354   1.533    0.128    
z.diff.lag -0.560181   0.076373  -7.335 2.69e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1615 on 122 degrees of freedom
Multiple R-squared:  0.3071,    Adjusted R-squared:  0.2958 
F-statistic: 27.04 on 2 and 122 DF,  p-value: 1.904e-10


Value of test-statistic is: 1.5329 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

adfTest(diff(log(ST)),lags=0,type=c("c"))

p-value smaller than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 0
  STATISTIC:
    Dickey-Fuller: -18.4442
  P VALUE:
    0.01 

Description:
 Fri Nov 30 12:21:21 2018 by user: dullu

7. usar herramientas ACF,PACF, EACF y/o criterios de akaike/bayes para construir propuestas de modelos.

Se hace un diagnostico de residuos de la serie.

8. como parte del diagnostico del modelo, analizar los residuos y aplicar prueba Ljung-box y ## 9. Usar la funcion auto.arima() y compare sus resultados con el inciso anterior.

Autoarima

Series: st 
ARIMA(0,1,1) with drift 

Coefficients:
          ma1   drift
      -0.6085  0.0117
s.e.   0.0758  0.0031

sigma^2 estimated as 0.00782:  log likelihood=126.6
AIC=-247.21   AICc=-247.01   BIC=-238.72

checkresiduals(Autoarima)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1) with drift
Q* = 13.787, df = 22, p-value = 0.9089

Model df: 2.   Total lags used: 24

El autoarima nos arroja un modelo en el que solamente sobresale un rezago, ahora se propondran mas modelos par saber si aun podemos corregir ms nuestra serie, de acuerdo al criterio de AIC.

PROPUESTAS DE MODELO

Se propone:

m1

Call:
arima(x = st, order = c(0, 1, 3), seasonal = c(1, 0, 1))

Coefficients:
          ma1     ma2      ma3    sar1     sma1
      -0.5538  0.0410  -0.0269  0.9990  -0.9747
s.e.   0.0898  0.1001   0.0947  0.0065   0.0836

sigma^2 estimated as 0.006439:  log likelihood = 130.24,  aic = -248.49

checkresiduals(m1)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,3)(1,0,1)[12]
Q* = 10.728, df = 19, p-value = 0.9326

Model df: 5.   Total lags used: 24

m2

Call:
arima(x = st, order = c(0, 1, 2), seasonal = c(1, 0, 1))

Coefficients:
          ma1     ma2    sar1     sma1
      -0.5530  0.0250  0.9987  -0.9709
s.e.   0.0887  0.0857  0.0085   0.0942

sigma^2 estimated as 0.006472:  log likelihood = 130.2,  aic = -250.4

checkresiduals(m2)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,2)(1,0,1)[12]
Q* = 10.409, df = 20, p-value = 0.9601

Model df: 4.   Total lags used: 24

Para escoger el modelo adecuado para nuestra serie nos basamos en el criterio de AIC:Para autoarima nos presenta: AIC=2196.47, para el modelo 1 nos presenta: aic = -248.49, para el segundo modelo nos presenta: aic = -250.4.

Segun el criterio de AIC el mejor modelo es el más pequeño y positivo por lo que el modelo elegido es el autoarima.

10. Presente la ecuacion final y describa.

\[ Y_t = -0.5530 \phi_{et-1}0.250\phi_{et-3}0.9987\theta_{et-2}-0.9709\theta\Theta_{et-12} \]

       esta es la ecuacion que nos representa el modelo.
       
       
       

11. realizar el pronostico arima para dos años.

pronostico

$`mean`
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep
2018                                                       10.23730 10.24905 10.26080
2019 10.30778 10.31953 10.33128 10.34302 10.35477 10.36651 10.37826 10.39001 10.40175
2020 10.44874 10.46049 10.47223 10.48398 10.49572 10.50747                           
          Oct      Nov      Dec
2018 10.27254 10.28429 10.29604
2019 10.41350 10.42525 10.43699
2020                           

$lower
              80%      95%
Jul 2018 10.12397 10.06398
Aug 2018 10.12734 10.06291
Sep 2018 10.13125 10.06267
Oct 2018 10.13561 10.06312
Nov 2018 10.14034 10.06414
Dec 2018 10.14541 10.06567
Jan 2019 10.15075 10.06763
Feb 2019 10.15635 10.06997
Mar 2019 10.16217 10.07265
Apr 2019 10.16819 10.07564
May 2019 10.17439 10.07891
Jun 2019 10.18076 10.08243
Jul 2019 10.18728 10.08618
Aug 2019 10.19394 10.09015
Sep 2019 10.20073 10.09431
Oct 2019 10.20764 10.09866
Nov 2019 10.21465 10.10317
Dec 2019 10.22178 10.10785
Jan 2020 10.22899 10.11267
Feb 2020 10.23631 10.11763
Mar 2020 10.24370 10.12273
Apr 2020 10.25118 10.12794
May 2020 10.25874 10.13328
Jun 2020 10.26636 10.13873

$upper
              80%      95%
Jul 2018 10.35063 10.41063
Aug 2018 10.37076 10.43519
Sep 2018 10.39034 10.45892
Oct 2018 10.40948 10.48197
Nov 2018 10.42823 10.50443
Dec 2018 10.44667 10.52640
Jan 2019 10.46481 10.54794
Feb 2019 10.48271 10.56909
Mar 2019 10.50038 10.58990
Apr 2019 10.51785 10.61040
May 2019 10.53514 10.63063
Jun 2019 10.55227 10.65060
Jul 2019 10.56924 10.67034
Aug 2019 10.58607 10.68987
Sep 2019 10.60278 10.70920
Oct 2019 10.61936 10.72834
Nov 2019 10.63584 10.74732
Dec 2019 10.65221 10.76614
Jan 2020 10.66848 10.78481
Feb 2020 10.68467 10.80334
Mar 2020 10.70076 10.82174
Apr 2020 10.71678 10.84001
May 2020 10.73271 10.85817
Jun 2020 10.74858 10.87621

      con este pronostico notamos que que la serie mantendrá un comportamiento constante    por lo que es probable (95%) que en los proximos 2 años el número de pasajeros del    metrobus no aumente exponencialmente.
 
       
      

12. utilizando los metodos de pronosticos simples y/o suavizamientoexponencial, generar un pronostico para dos años. elegir el metodo que presenta la raiz de error medio al cuadrado más pequeño

   aqui se esta creando la ventana en la cual generaremos nuetros pronosticos
   

pronosticos.

ST3 <- window(ST, start=2008)
  accuracy (STfit1,ST3)#mean

                        ME     RMSE      MAE       MPE     MAPE     MASE       ACF1
Training set -3.524095e-13 4160.820 3716.157 -11.93560 32.82269 1.737421  0.9300643
Test set      6.099221e+03 6304.002 6099.221  30.33576 30.33576 2.851578 -0.3683302
             Theil's U
Training set        NA
Test set      2.290417

accuracy (STfit2,ST3)#naive

                   ME     RMSE       MAE       MPE     MAPE      MASE       ACF1
Training set 144.0601 1222.182  966.5532 0.8917718 6.912104 0.4518941 -0.4096249
Test set     871.9492 1816.658 1567.5052 3.7528064 7.835335 0.7328580 -0.3683302
             Theil's U
Training set        NA
Test set     0.6850087

accuracy (STfit3,ST3)#seasonal naive

                    ME     RMSE      MAE      MPE      MAPE      MASE       ACF1
Training set 1844.7206 2573.506 2138.893 13.24892 15.177397 1.0000000  0.6727412
Test set      843.2108 1573.798 1322.181  3.91727  6.686061 0.6181611 -0.4843001
             Theil's U
Training set        NA
Test set     0.6205445

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

Analisis VAR

library(vars)

Error in library(vars) : there is no package called ‘vars’

stf<- read.csv(file.choose())

Error in file.choose() : file choice cancelled

Una vez declaradas ambas variables como serie de tiempo se grafican.Esta funcionnos permite visualizar el comportamiento de nuestras variables.

A continuacion se graficara usando la funcion ggseasonplot para corroborar si la variable presenta estacionalidad.

Ahora con la segunda variable kilometros recorridos.

Con esa grafica corroboramos que la segunda variable no tiene estacionalidad.

summary(reg.w)

Call:
lm(formula = km ~ time(km) + mes.)

Residuals:
    Min      1Q  Median      3Q     Max 
-472.04 -110.91  -16.47   92.66  413.90 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -5.646e+05  1.026e+04 -55.031   <2e-16 ***
time(km)       2.816e+02  5.097e+00  55.246   <2e-16 ***
mes.February  -1.380e+02  7.386e+01  -1.868   0.0643 .  
mes.March      3.434e+01  7.386e+01   0.465   0.6429    
mes.April      1.729e+01  7.387e+01   0.234   0.8154    
mes.May        1.243e+02  7.388e+01   1.682   0.0953 .  
mes.June       1.090e+02  7.389e+01   1.475   0.1430    
mes.July       1.345e+02  7.568e+01   1.777   0.0782 .  
mes.August     9.111e+01  7.568e+01   1.204   0.2312    
mes.September -7.247e+00  7.569e+01  -0.096   0.9239    
mes.October    9.110e+01  7.569e+01   1.204   0.2313    
mes.November  -3.428e+01  7.570e+01  -0.453   0.6515    
mes.December  -4.612e+01  7.571e+01  -0.609   0.5436    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 173.2 on 113 degrees of freedom
Multiple R-squared:  0.9647,    Adjusted R-squared:  0.9609 
F-statistic: 257.3 on 12 and 113 DF,  p-value: < 2.2e-16

La tabla nos muestra que los meses de febrero, mayo y julio existe un 10% de significancia.      

summary(reg.p)

Call:
lm(formula = pass ~ time(pass) + mes.)

Residuals:
    Min      1Q  Median      3Q     Max 
-6273.4  -726.4  -147.7  1036.6  3479.2 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -3.777e+06  9.316e+04 -40.542   <2e-16 ***
time(pass)     1.884e+03  4.628e+01  40.711   <2e-16 ***
mes.February  -1.746e+02  6.706e+02  -0.260   0.7951    
mes.March      3.317e+02  6.707e+02   0.495   0.6219    
mes.April      9.531e+01  6.707e+02   0.142   0.8873    
mes.May        1.233e+03  6.708e+02   1.839   0.0686 .  
mes.June       8.126e+02  6.709e+02   1.211   0.2284    
mes.July       3.972e+02  6.872e+02   0.578   0.5644    
mes.August     1.437e+03  6.872e+02   2.091   0.0388 *  
mes.September -1.411e+02  6.872e+02  -0.205   0.8378    
mes.October    1.522e+03  6.873e+02   2.215   0.0288 *  
mes.November   2.967e+02  6.874e+02   0.432   0.6668    
mes.December  -8.725e+02  6.875e+02  -1.269   0.2070    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1573 on 113 degrees of freedom
Multiple R-squared:  0.9372,    Adjusted R-squared:  0.9305 
F-statistic: 140.4 on 12 and 113 DF,  p-value: < 2.2e-16

La grafica nos dice que la variable pasajeros es significativa al 5% en los meses de agosto y octubre respecto a enero.      

El modelo indicado para modelar esta variable es el que nos marco el autoarima en el punto pasado y es: AR(0,1,1).

Ahora haremos lo mismo con la variable kilometros recorridos.

Usaremos la funcion autoarima para buscar el modelo indicado para corregir esta variable.

Autoarima2<- auto.arima(stf, d=1,stepwise=FALSE, approximation=FALSE)
Autoarima

checkresiduals()