1 Modelimg the Mexican Economy

1.1 Data Download

library(readxl)
library(xts)
library(zoo)
library(tseries)
library(forecast)
library(astsa)
library(lmtest)

igae <- read_xls("IGAE.xls", sheet = "Página 1") 

# I create the a ts dataset indicating the dates:
igae.ts<-ts(coredata(igae$IGAE),start=c(1993,1),frequency=12)

igae.ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1993  60.40769  61.02252  63.94325  61.86598  63.61290  62.88259  62.79795
## 1994  63.02927  62.73316  65.79027  65.89172  66.85880  66.63579  64.70206
## 1995  65.10234  60.29419  62.92987  59.23909  61.17677  61.08234  59.77917
## 1996  64.59736  63.25033  65.00435  63.82202  66.56896  65.57817  65.91223
## 1997  67.92226  66.01572  67.08871  70.18939  71.25363  70.71757  70.48390
## 1998  72.28125  71.26764  75.30751  72.95238  75.13802  74.64519  74.61446
## 1999  74.28381  72.70940  76.99211  74.32923  76.92304  76.55675  76.31262
## 2000  77.93003  77.36732  80.35795  77.00756  82.34331  81.54154  79.62973
## 2001  79.42346  76.41926  80.76787  77.53115  81.78409  80.60701  79.32531
## 2002  77.38629  75.10932  77.12227  80.83024  82.06994  79.72127  79.92718
## 2003  79.27239  77.17550  80.33976  79.95093  82.28842  81.44337  80.93570
## 2004  80.93007  79.19382  84.92785  82.93623  85.23905  85.62087  82.94023
## 2005  82.25837  80.89942  83.62640  86.59373  87.88889  86.35209  83.78420
## 2006  87.04556  84.08297  89.54008  87.10250  93.13935  91.36864  88.35889
## 2007  88.77729  85.76828  91.42205  89.32063  94.32237  93.53460  91.41102
## 2008  90.95435  88.98933  89.05326  94.71533  94.86436  94.45588  94.15931
## 2009  84.89588  81.84098  86.87355  84.31067  86.08159  88.18426  88.73169
## 2010  86.86702  85.22716  92.29229  90.58182  92.30005  93.23105  92.32869
## 2011  90.20546  88.56874  95.82899  91.60005  95.94631  96.48831  95.21792
## 2012  94.67463  94.07903  99.09222  95.47932 100.31525  99.68684  99.31091
## 2013  97.96858  94.67101  97.06079  99.91623 102.27529  99.51755 100.80982
## 2014  98.75513  96.95653 101.49481 100.31854 104.86081 102.89402 104.13869
## 2015 102.31977  99.87062 104.96625 103.68802 106.19806 107.31713 107.46008
## 2016 104.41881 104.57988 105.90821 107.29216 109.06156 109.96077 107.52441
## 2017 108.42006 105.28496 111.90569 106.20453 112.30643 112.80000 108.95594
## 2018 107.81001 111.19094 111.59738 115.82252 114.44374 112.81183 114.33744
## 2019 109.09518 112.46967 109.98886 115.45875 113.04393 113.54114 113.42805
## 2020 108.40448 109.91800  87.99004  89.47630  98.01972 102.22197 102.93410
##            Aug       Sep       Oct       Nov       Dec
## 1993  62.02789  62.01834  62.83480  63.62177  66.12119
## 1994  65.79634  65.76477  66.73068  67.84687  68.06744
## 1995  61.06820  60.84772  60.21824  63.06020  65.41346
## 1996  65.36457  65.03774  67.82464  68.55904  69.26845
## 1997  70.02845  70.75158  73.23728  72.75789  73.78044
## 1998  73.47925  73.76732  74.46418  74.44344  75.67735
## 1999  75.67001  76.24775  76.45465  77.67013  78.05347
## 2000  81.19074  80.16487  80.80958  80.91619  79.49830
## 2001  80.58703  78.26957  79.97355  80.09390  79.09247
## 2002  80.30576  78.38018  81.60774  80.09518  80.49051
## 2003  79.41995  79.18878  81.87911  80.63451  83.33037
## 2004  83.13244  82.19751  84.36268  85.94154  86.23057
## 2005  86.43120  84.64256  86.85310  88.65732  89.15110
## 2006  90.24619  87.86910  91.59749  91.14520  90.86138
## 2007  92.21733  89.47518  95.27999  93.80926  92.45720
## 2008  91.49309  90.57651  95.50683  92.04270  91.89505
## 2009  86.55360  86.50378  91.12057  91.24910  91.85630
## 2010  91.56219  90.55653  93.84529  95.96933  95.59592
## 2011  96.53903  94.10211  97.43771 101.08275  98.81575
## 2012  99.09917  95.18924 101.85991 104.38649 100.52817
## 2013 100.27357  96.80675 103.66719 104.23166 102.80157
## 2014 101.73624 100.01550 106.81109 106.47987 106.88143
## 2015 105.48271 105.02281 109.29460 109.09904 109.26256
## 2016 108.94810 106.06088 110.56483 114.14690 112.51169
## 2017 111.77091 106.32100 112.84684 116.17409 113.99218
## 2018 108.78172 115.99083 117.77276 113.17022 112.71950
## 2019 109.01316 115.46429 116.32435 113.47028 111.92756
## 2020 109.33333 111.49701 110.43280

1.2 Plotting the index

plot(igae.ts)

I check wether the seasonal difference of the log is stationary:

logigae.ts=log(igae.ts)
plot(diff(log(igae.ts), lag=12))

adf.test(diff(log(igae.ts), lag =12),k=0)
## Warning in adf.test(diff(log(igae.ts), lag = 12), k = 0): p-value smaller than
## printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(log(igae.ts), lag = 12)
## Dickey-Fuller = -6.6059, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary
acf2(diff(log(igae.ts), lag = 12), max.lag = 24)

##      [,1] [,2] [,3]  [,4] [,5]  [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  0.76 0.68 0.62  0.47 0.42  0.34  0.22  0.15 0.09  0.01  0.00 -0.08 -0.06
## PACF 0.76 0.25 0.10 -0.18 0.04 -0.03 -0.13 -0.08 0.02 -0.04  0.07 -0.14  0.14
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.06 -0.07 -0.04 -0.05 -0.09 -0.07 -0.04 -0.05 -0.04 -0.03  -0.1
## PACF  0.02  0.02 -0.01 -0.03 -0.11  0.00  0.10 -0.03 -0.04  0.06  -0.2

I see that this is a typical AR signature with p=2

I run the Arima

model1 <- Arima(igae.ts, order = c(2,0,0), 
            seasonal = list(order=c(0,1,0),period=12),
            include.constant = TRUE,
            lambda = 0) 

model1
## Series: igae.ts 
## ARIMA(2,0,0)(0,1,0)[12] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1     ar2   drift
##       0.5735  0.2518  0.0015
## s.e.  0.0539  0.0540  0.0007
## 
## sigma^2 estimated as 0.0007434:  log likelihood=704.12
## AIC=-1400.24   AICc=-1400.12   BIC=-1385.14
coeftest(model1)
## 
## z test of coefficients:
## 
##         Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.57351410 0.05385938 10.6484 < 2.2e-16 ***
## ar2   0.25175152 0.05398221  4.6636 3.107e-06 ***
## drift 0.00149951 0.00071005  2.1118    0.0347 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check residuals

acf2(model1$residuals,max.lag = 24)

##       [,1]  [,2] [,3]  [,4] [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.02 -0.02 0.15 -0.11 0.06 0.10 -0.05  0.01 0.02 -0.10  0.07 -0.18 -0.02
## PACF -0.02 -0.02 0.15 -0.11 0.06 0.07 -0.01 -0.02 0.00 -0.08  0.06 -0.21  0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.01 -0.06  0.08  0.04 -0.12 -0.06  0.10 -0.05  0.03  0.13 -0.17
## PACF -0.04  0.02  0.05  0.06 -0.09 -0.09  0.09  0.00 -0.02  0.13 -0.18

2 I add AR(3) term

model2 <- Arima(igae.ts, order = c(3,0,0), 
            seasonal = list(order=c(0,1,0),period=12),
            include.constant = TRUE,
            lambda = 0) 
model2
## Series: igae.ts 
## ARIMA(3,0,0)(0,1,0)[12] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1     ar2     ar3   drift
##       0.5500  0.1994  0.0916  0.0015
## s.e.  0.0555  0.0624  0.0556  0.0008
## 
## sigma^2 estimated as 0.0007395:  log likelihood=705.47
## AIC=-1400.94   AICc=-1400.75   BIC=-1382.07
coeftest(model2)
## 
## z test of coefficients:
## 
##         Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.55003097 0.05550205  9.9101 < 2.2e-16 ***
## ar2   0.19935928 0.06243262  3.1932  0.001407 ** 
## ar3   0.09162553 0.05561315  1.6476  0.099445 .  
## drift 0.00149636 0.00077318  1.9353  0.052949 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acf2(model2$residuals,max.lag = 24)

##      [,1] [,2] [,3]  [,4] [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  0.02 0.04  0.1 -0.11 0.05 0.07 -0.04 -0.01 0.01 -0.11  0.06 -0.19 -0.03
## PACF 0.02 0.04  0.1 -0.12 0.05 0.07 -0.03 -0.04 0.01 -0.09  0.05 -0.21  0.01
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.00 -0.05  0.08  0.03 -0.10 -0.07  0.08 -0.03  0.03  0.12 -0.16
## PACF -0.02  0.01  0.05  0.04 -0.09 -0.08  0.09  0.01 -0.04  0.13 -0.18

phi3 is not significantly bigger than zero, then I try MA(3)

model3 <- Arima(igae.ts, order = c(2,0,3), 
            seasonal = list(order=c(0,1,0),period=12),
            include.constant = TRUE,
            lambda = 0) 
model3
## Series: igae.ts 
## ARIMA(2,0,3)(0,1,0)[12] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1     ar2     ma1      ma2     ma3   drift
##       0.0037  0.7115  0.5896  -0.1756  0.1873  0.0015
## s.e.  0.0666  0.0662  0.0812   0.0937  0.0676  0.0007
## 
## sigma^2 estimated as 0.0007117:  log likelihood=712.33
## AIC=-1410.65   AICc=-1410.3   BIC=-1384.23
coeftest(model3)
## 
## z test of coefficients:
## 
##          Estimate  Std. Error z value  Pr(>|z|)    
## ar1    0.00372040  0.06657151  0.0559   0.95543    
## ar2    0.71148542  0.06622464 10.7435 < 2.2e-16 ***
## ma1    0.58962888  0.08120348  7.2611 3.839e-13 ***
## ma2   -0.17555623  0.09372034 -1.8732   0.06104 .  
## ma3    0.18729138  0.06760578  2.7703   0.00560 ** 
## drift  0.00152152  0.00068076  2.2350   0.02542 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acf2(model3$residuals,max.lag = 24)

##      [,1] [,2] [,3]  [,4] [,5] [,6]  [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF     0    0 0.02 -0.04 0.03 0.13 -0.06 0.02    0 -0.06  0.04 -0.16 -0.04
## PACF    0    0 0.02 -0.04 0.03 0.13 -0.06 0.01    0 -0.04  0.03 -0.17 -0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.01 -0.04  0.09  0.03 -0.12 -0.07  0.07 -0.02  0.04  0.10 -0.16
## PACF  0.00 -0.03  0.10  0.02 -0.07 -0.09  0.08  0.00 -0.01  0.12 -0.17
model4 <- Arima(igae.ts, order = c(2,0,0), 
            seasonal = list(order=c(0,1,1),period=12),
            include.constant = TRUE,
            lambda = 0) 
model4
## Series: igae.ts 
## ARIMA(2,0,0)(0,1,1)[12] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1     ar2     sma1   drift
##       0.6117  0.2496  -0.6647  0.0018
## s.e.  0.0543  0.0545   0.0716  0.0003
## 
## sigma^2 estimated as 0.0006354:  log likelihood=726.51
## AIC=-1443.03   AICc=-1442.84   BIC=-1424.15
coeftest(model4)
## 
## z test of coefficients:
## 
##          Estimate  Std. Error z value  Pr(>|z|)    
## ar1    0.61168638  0.05432471 11.2598 < 2.2e-16 ***
## ar2    0.24962675  0.05450779  4.5797 4.657e-06 ***
## sma1  -0.66473926  0.07158761 -9.2857 < 2.2e-16 ***
## drift  0.00175145  0.00030246  5.7907 7.009e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acf2(model4$residuals,max.lag = 24)

##       [,1]  [,2] [,3]  [,4] [,5]  [,6]  [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.04 -0.08 0.13 -0.06 0.05  0.00 -0.01 0.04 0.05 -0.08  0.06  0.08 -0.10
## PACF -0.04 -0.09 0.12 -0.06 0.07 -0.02  0.02 0.02 0.06 -0.08  0.06  0.05 -0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.04 -0.02  0.06  0.12 -0.13 -0.07  0.12 -0.05 -0.01  0.13 -0.14
## PACF  0.02 -0.03  0.09  0.09 -0.10 -0.09  0.08 -0.03  0.02  0.09 -0.13

2.1 Forecasting the Model:

forecast_igae <- forecast(model4, h=48)
# I stored the information of the forecast in forecast_air

# I plot the forecast
autoplot(forecast_igae)