Workshop 3
2022-06-02
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
#1 Instalamos los paquetes y los activamos para poder ejecutar el codigo de manera correcta. Llamar al excel y colocarlo en la misma carpeta que correspone al proyecto de R, para que se pueda obtener la base de datos con la que se quiere trabajar.
library(readxl)
library(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(zoo)
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)
library(astsa)##
## Attaching package: 'astsa'## The following object is masked from 'package:forecast':
##
## gaslibrary(lmtest)
igae <- read_excel("igae.xls")
# I convert the igae data frame into a time series R dataset.
# This is convenient when working with ARIMA-SARIMA models:
igae<-ts(coredata(igae$IGAE),start=c(1993,1),frequency=12)
# The ts function transform the dataset from a data frame to a ts object
# I can see the content of IGAE in a table format:
igae## 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.17551 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.56875 95.82899 91.60005 95.94631 96.48831 95.21792
## 2012 94.67463 94.07904 99.09222 95.47932 100.31525 99.68684 99.31091
## 2013 97.96857 94.67100 97.06078 99.91621 102.27528 99.51753 100.80982
## 2014 98.75517 96.95658 101.49488 100.31859 104.86086 102.89406 104.13871
## 2015 102.31967 99.87044 104.96599 103.68780 106.19784 107.31698 107.46002
## 2016 104.41914 104.58039 105.90879 107.29295 109.06275 109.96190 107.52515
## 2017 108.41883 105.28225 111.90313 106.20035 112.30199 112.79610 108.95429
## 2018 111.01863 107.81545 111.20270 111.60594 115.83099 114.44769 112.80452
## 2019 112.50953 108.82637 112.27312 109.80338 115.22991 112.82079 113.37014
## 2020 112.36544 108.86540 110.03905 88.24539 88.99407 97.93032 102.37454
## 2021 105.36942 102.67215 111.51814 108.08017 112.06310 111.19094 109.89353
## 2022 107.31808 105.27670 111.95949
## 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.84530 95.96933 95.59592
## 2011 96.53903 94.10211 97.43771 101.08275 98.81575
## 2012 99.09917 95.18924 101.85991 104.38648 100.52816
## 2013 100.27357 96.80675 103.66720 104.23168 102.80161
## 2014 101.73625 100.01548 106.81104 106.47979 106.88127
## 2015 105.48268 105.02291 109.29472 109.09934 109.26326
## 2016 108.94857 106.06069 110.56410 114.14503 112.50874
## 2017 111.77004 106.32251 112.84965 116.18136 114.00214
## 2018 114.32657 108.76697 115.95791 117.71706 113.09200
## 2019 113.23208 108.85990 115.37576 116.29415 113.40311
## 2020 102.74109 103.36978 109.26091 111.07944 110.86607
## 2021 107.25443 104.26659 108.48678 112.96735 112.33393
## 2022plot(igae)
#La grafica representa las crisis que el paÃs ha sufrido a lo largo del
tiempo. Pero, en este caso el periodo de principios de 2020 empata con
covid-19, de esta manera, podemos ver que la caida es mucho mayor y más
pronunciada que en las dos crisis anteriores. Sin mencionar, que es más
prolongada.
logigae = log(igae)plot(logigae)
seasonaldiflog = diff(logigae,lag=12)
plot(seasonaldiflog)
adf.test(seasonaldiflog, k=1)## Warning in adf.test(seasonaldiflog, k = 1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: seasonaldiflog
## Dickey-Fuller = -5.737, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationarylibrary(astsa)
acf2(seasonaldiflog)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0.79 0.64 0.53 0.41 0.35 0.28 0.18 0.13 0.06 -0.05 -0.14 -0.23 -0.18
## PACF 0.79 0.04 0.03 -0.09 0.08 -0.03 -0.11 0.01 -0.07 -0.17 -0.09 -0.11 0.30
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF -0.11 -0.09 -0.06 -0.03 -0.05 -0.06 -0.03 -0.05 -0.05 -0.03 -0.09 -0.07
## PACF 0.08 -0.01 0.00 0.06 -0.12 -0.07 0.09 -0.09 -0.13 0.03 -0.20 0.26
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF -0.05 -0.08 -0.09 -0.08 -0.09 -0.08 -0.07 -0.09 -0.08 -0.09 -0.05 -0.06
## PACF 0.07 -0.06 -0.04 0.07 -0.09 -0.07 0.09 -0.08 -0.14 0.04 0.00 0.04
## [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF -0.06 -0.05 -0.03 -0.04 -0.03 0.00 -0.01 0.00 0.02 0.01 -0.03
## PACF 0.11 -0.04 0.02 -0.01 -0.06 0.02 0.03 -0.07 -0.06 0.00 -0.15# Load package forecast
library(forecast)
# Use function Arima
m1a <- Arima(seasonaldiflog, order = c(2,0,0),
include.constant = TRUE)
# we can display the results of the coefficients and their p-values
# with the coeftest function:
coeftest(m1a)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.7546609 0.0541815 13.9284 < 2.2e-16 ***
## ar2 0.0431344 0.0541729 0.7962 0.425896
## intercept 0.0197113 0.0074153 2.6582 0.007856 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Use function Arima
m1b <- Arima(igae, order = c(2,0,0),
seasonal = list(order=c(0,1,0), period=12),
include.constant = TRUE,
lambda = 0)
# we can display the results of the coefficients and their p-values
# with the coeftest function:
coeftest(m1b)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.75466220 0.05418155 13.9284 <2e-16 ***
## ar2 0.04313591 0.05417288 0.7963 0.4259
## drift 0.00164252 0.00062032 2.6479 0.0081 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forecast_igae <- forecast(m1b, h = 48)
autoplot(forecast_igae)