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)
#La funciĆ³n diff tiene como fin identificar el valor del registro de
hace doce meses. Lo que representa el cambio en porcentaje del cambio
anual, ya que es el periodo de tiempo que se esta considerando.
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)
#En este modelo, se entiende como valores actuales. Pero, utilizado los valores y variables anteriores. Diferente a un modelo de regresiĆ³n lineal.