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
I check wether the seasonal difference of the log is stationary:
## 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
## [,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
##
## 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
## [,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
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
##
## 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
## [,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
##
## 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
## [,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
##
## 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
## [,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