Clase

Carga de datos

library(readxl)
library(forecast)
library(magrittr)
SERIEIVAE <- read_excel("~/EMA1182022/IVAE/IVAE_SLV.xlsx",col_types = c("skip", "numeric"), skip = 5)

SERIEIVAE.ts <- ts(data = SERIEIVAE, start = c(2009, 1), frequency = 12)
SERIEIVAE.ts %>% autoplot(main = "IVAE, El Salvador 2009-2021[marzo]", xlab = "Años/Meses", ylab = "Indice")

Yt<- SERIEIVAE.ts

Estimación de modelo HW

ModeloHW<-HoltWinters(x = Yt,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
ModeloHW
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Yt, seasonal = "multiplicative", optim.start = c(0.9,     0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.8408163
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   117.0799442
## b     0.1600306
## s1    0.9502255
## s2    1.0233274
## s3    1.0518154
## s4    0.9900276
## s5    1.0007537
## s6    0.9807088
## s7    0.9600206
## s8    1.0149628
## s9    1.0915423
## s10   0.9796752
## s11   0.9584676
## s12   0.9994880
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
##          Point Forecast     Lo 95    Hi 95
## Apr 2021       111.4044 105.94952 116.8593
## May 2021       120.1386 112.77972 127.4976
## Jun 2021       123.6515 114.83356 132.4694
## Jul 2021       116.5461 107.01096 126.0813
## Aug 2021       117.9689 107.30373 128.6341
## Sep 2021       115.7630 104.33417 127.1918
## Oct 2021       113.4746 101.36814 125.5810
## Nov 2021       120.1312 106.57272 133.6897
## Dec 2021       129.3698 114.12877 144.6109
## Jan 2022       116.2681 101.78191 130.7543
## Feb 2022       113.9046  98.99575 128.8134
## Mar 2022       118.9394  81.61739 156.2614
PronosticosHW %>% autoplot()

Aproximación de espacio de los estados

PronosticosHW2<-hw(y = Yt,
                   h = 12,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
PronosticosHW2
##          Point Forecast     Lo 95    Hi 95
## Apr 2021       107.4928  99.59377 115.3918
## May 2021       115.3370 106.85846 123.8155
## Jun 2021       115.6946 107.18632 124.2029
## Jul 2021       109.6301 101.56399 117.6962
## Aug 2021       112.2047 103.94477 120.4646
## Sep 2021       110.2284 102.10923 118.3476
## Oct 2021       107.9393  99.98348 115.8951
## Nov 2021       114.4306 105.99022 122.8710
## Dec 2021       122.6459 113.59234 131.6994
## Jan 2022       108.5830 100.56060 116.6054
## Feb 2022       107.1239  99.20182 115.0460
## Mar 2022       113.1678 104.79016 121.5455
PronosticosHW2 %>% autoplot()

identificando orden de integración

library(kableExtra)
d<-ndiffs(Yt)
D<-nsdiffs(Yt)
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 1
D 1

Grafico de la serie estacionaria

library(TSstudio)
Yt %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "IVAE estacionaria")

Valores para (p,q) y (P,Q)

Yt %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_cor(lag.max = 36)

En caso observando el PACF los Componentes autoregresivos el valor de p=0 mientras que P=1 mientras que los Componentes de media movil los valores de Q=1 y q=0

Estimando SARIMA

\[SARIMA(0,1,0)(1,1,1)_{[12]}\]

library(ggthemes)
modelo_estimado <- Yt %>% 
  Arima(order = c(0, 1, 0),
        seasonal = c(1, 1, 1))

summary(modelo_estimado)
## Series: . 
## ARIMA(0,1,0)(1,1,1)[12] 
## 
## Coefficients:
##          sar1     sma1
##       -0.2040  -0.6422
## s.e.   0.1586   0.1554
## 
## sigma^2 = 6.897:  log likelihood = -323.43
## AIC=652.85   AICc=653.04   BIC=661.55
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.04712358 2.488576 1.666188 0.02203353 1.658867 0.4981404
##                    ACF1
## Training set 0.08970896
modelo_estimado %>% autoplot(type="both")+theme_solarized()

modelo_estimado %>% check_res(lag.max = 36)

Para el valor de MAPE es de 1.65% lo que quiere decir que de cada 100 valores pronosticados el error cometido sera de 1.65% y ante los residuos estamos en presencia de ruido blanco

Grafico comparativo

Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)

PIB trimestral

Carga de datos

library(readxl)
library(forecast)
PIB_Trim <- read_excel("~/EMA1182022/SARIMA/PIB_Trim.xlsx", 
    col_names = F, skip = 6)

Serie_PIBT<-ts(data =PIB_Trim,start = c(2009,1), frequency = 4 )

Serie_PIBT[,2] %>%  autoplot(main = "PIB trimestral a precios corrientes, El Salvador 2009-2021", xlab = "Años/Meses", ylab = "Indice")

Estimación de modelo HW

ModeloHW<-HoltWinters(x = Serie_PIBT[,2],
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
ModeloHW
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Serie_PIBT[, 2], seasonal = "multiplicative",     optim.start = c(0.9, 0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.8742302
##  beta : 0
##  gamma: 0
## 
## Coefficients:
##            [,1]
## a  7429.8128051
## b    48.3963750
## s1    0.9751646
## s2    1.0058150
## s3    0.9783521
## s4    1.0406683
PronosticosHW<-forecast(object = ModeloHW,h = 4,level = c(0.95))
PronosticosHW
##         Point Forecast    Lo 95    Hi 95
## 2022 Q1       7292.485 6847.413 7737.557
## 2022 Q2       7570.373 6930.977 8209.768
## 2022 Q3       7411.019 6646.235 8175.804
## 2022 Q4       7933.429 7119.875 8746.982
PronosticosHW %>% autoplot()

Aproximación de espacio de los estados

PronosticosHW2<-hw(y = Serie_PIBT[,2],
                   h = 4,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
PronosticosHW2
##         Point Forecast    Lo 95    Hi 95
## 2022 Q1       7284.036 6676.405 7891.668
## 2022 Q2       7541.782 6703.701 8379.863
## 2022 Q3       7515.116 6517.016 8513.215
## 2022 Q4       8063.144 6844.733 9281.555
PronosticosHW2 %>% autoplot()

identificando orden de integración

library(kableExtra)
d<-ndiffs(Serie_PIBT[,2])
D<-nsdiffs(Serie_PIBT[,2])
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 1
D 1

Grafico de la serie estacionaria

library(TSstudio)
Serie_PIBT[,2] %>% 
  diff(lag = 4,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "PIB trimestral estacionaria")

Valores para (p,q) y (P,Q)

Serie_PIBT[,2] %>% 
  diff(lag = 4,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_cor(lag.max = 12)

En caso observando el PACF los Componentes autoregresivos el valor de p=0 mientra que P=1 ya que el retardo 8 no es tan significativo mientras que los Componentes de media movil los valores de Q=1 y q=0

Estimando SARIMA

\[SARIMA(0,1,0)(1,1,1)_{[4]}\]

library(ggthemes)
modelo_estimado <- Serie_PIBT[,2] %>% 
  Arima(order = c(0, 1, 0),
        seasonal = c(1, 1, 1))

summary(modelo_estimado)
## Series: . 
## ARIMA(0,1,0)(1,1,1)[4] 
## 
## Coefficients:
##          sar1     sma1
##       -0.3980  -0.4550
## s.e.   0.3732   0.3688
## 
## sigma^2 = 77377:  log likelihood = -331.67
## AIC=669.34   AICc=669.9   BIC=674.89
## 
## Training set error measures:
##                    ME     RMSE      MAE          MPE     MAPE      MASE
## Training set 6.530071 258.7683 120.9372 -0.004468918 2.054645 0.3680668
##                     ACF1
## Training set -0.09806694
modelo_estimado %>% autoplot(type="both")+theme_solarized()

modelo_estimado %>% check_res(lag.max = 12)

Para el valor de MAPE es de 2.05% lo que quiere decir que de cada 100 valores pronosticados el error cometido sera de 2.05% y ante los residuos estamos en presencia de ruido blanco

Grafico comparativo

PIB_Sarima<-modelo_estimado$fitted
PIB_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Serie_PIBT[,2],PIB_Sarima,PIB_HW)
ts_plot(grafico_comparativo)

IPC

Carga de datos

library(readxl)
library(forecast)
IPC <- read_excel("~/EMA1182022/SARIMA/IPC.xlsx", 
    col_names = F, skip = 6)

Serie_IPCT<-ts(data =IPC,start = c(2009,1), frequency = 12 )

Serie_IPCT[,2] %>%  autoplot(main = "Indice de precios al consumidor general, El Salvador 2009-2022[mayo]", xlab = "Años/Meses", ylab = "Indice")

Estimación de modelo HW

ModeloHW<-HoltWinters(x = Serie_IPCT[,2],
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
ModeloHW
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Serie_IPCT[, 2], seasonal = "multiplicative",     optim.start = c(0.9, 0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.8353873
##  beta : 0.07612657
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   122.7690242
## b     0.4801476
## s1    1.0043291
## s2    1.0021040
## s3    0.9986014
## s4    0.9975846
## s5    0.9989586
## s6    0.9984650
## s7    0.9959074
## s8    0.9975598
## s9    1.0007205
## s10   1.0041558
## s11   1.0048925
## s12   1.0053839
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       123.7827 122.6365 124.9289
## Jul 2022       123.9896 122.4499 125.5294
## Aug 2022       124.0357 122.1459 125.9256
## Sep 2022       124.3884 122.1648 126.6120
## Oct 2022       125.0394 122.4870 127.5918
## Nov 2022       125.4570 122.5842 128.3298
## Dec 2022       125.6139 122.4280 128.7997
## Jan 2023       126.3012 122.7897 129.8128
## Feb 2023       127.1819 123.3369 131.0269
## Mar 2023       128.1007 123.9171 132.2842
## Apr 2023       128.6771 124.1616 133.1927
## May 2023       129.2228 123.9338 134.5118
PronosticosHW %>% autoplot()

Aproximación de espacio de los estados

PronosticosHW2<-hw(y = Serie_IPCT[,2],
                   h = 12,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
PronosticosHW2
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       123.8302 121.9984 125.6620
## Jul 2022       124.5773 122.6371 126.5174
## Aug 2022       125.1928 123.0924 127.2932
## Sep 2022       125.6452 123.3312 127.9591
## Oct 2022       126.4964 123.9071 129.0858
## Nov 2022       127.6980 124.7742 130.6219
## Dec 2022       127.9518 124.6694 131.2342
## Jan 2023       128.7913 125.0953 132.4873
## Feb 2023       129.3087 125.1716 133.4457
## Mar 2023       130.6440 126.0045 135.2834
## Apr 2023       131.4546 126.2980 136.6113
## May 2023       132.0442 126.3504 137.7380
PronosticosHW2 %>% autoplot()

identificando orden de integración

library(kableExtra)
d<-ndiffs(Serie_IPCT[,2])
D<-nsdiffs(Serie_IPCT[,2])
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 1
D 0

Grafico de la serie estacionaria

library(TSstudio)
Serie_IPCT[,2] %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "IPC estacionaria")

Valores para (p,q) y (P,Q)

Serie_IPCT[,2] %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_cor(lag.max = 36)

En caso observando el PACF los Componentes autoregresivos el valor de p=4 mientras que P=2 ya que el retardo 36 no es tan significativo mientras que los Componentes de media movil los valores de Q=1 y q=2

Estimando SARIMA

\[SARIMA(2,1,4)(1,1,2)_{[12]}\]

library(ggthemes)
modelo_estimado <- Serie_IPCT[,2] %>% 
  Arima(order = c(2, 1, 4),
        seasonal = c(1, 1, 2))

summary(modelo_estimado)
## Series: . 
## ARIMA(2,1,4)(1,1,2)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     ma3     ma4     sar1     sma1
##       1.7824  -0.9301  -1.5629  0.5888  0.1283  0.0660  -0.2522  -0.5454
## s.e.  0.0820   0.0819   0.1312  0.1897  0.1683  0.0867   0.3401   0.3397
##          sma2
##       -0.4545
## s.e.   0.3160
## 
## sigma^2 = 0.2179:  log likelihood = -107.11
## AIC=234.21   AICc=235.82   BIC=264.18
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.01898683 0.4337828 0.2809931 0.01539364 0.2556043 0.1522693
##                     ACF1
## Training set 0.001775823
modelo_estimado %>% autoplot(type="both")+theme_solarized()

modelo_estimado %>% check_res(lag.max = 36)

Para el valor de MAPE es de 0.26% lo que quiere decir que de cada 100 valores pronosticados el error cometido sera de 0.26% y ante los residuos estamos en presencia de ruido blanco

Grafico comparativo

IPC_Sarima<-modelo_estimado$fitted
IPC_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Serie_IPCT[,2],IPC_Sarima,IPC_HW)
ts_plot(grafico_comparativo)

Hidrocarburos

Carga de datos

library(readxl)
library(forecast)
library(readxl)
Importacion_Hidrocarburos <- read_excel("~/EMA1182022/SARIMA/Importacion_Hidrocarburos.xlsx", 
    col_names = F, skip = 6)

Serie_HCT<-ts(data =Importacion_Hidrocarburos,start = c(2017,1), frequency = 12 )

Serie_HCT[,2] %>%  autoplot(main = "Importación de hidrocarburos, El Salvador 2017-2022[mayo]", xlab = "Años/Meses", ylab = "Indice")

Estimación de modelo HW

ModeloHW<-HoltWinters(x = Serie_HCT[,2],
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
ModeloHW
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Serie_HCT[, 2], seasonal = "multiplicative",     optim.start = c(0.9, 0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.1356726
##  beta : 0
##  gamma: 0.4908596
## 
## Coefficients:
##             [,1]
## a   2.550794e+05
## b   1.239453e+03
## s1  9.784035e-01
## s2  9.249193e-01
## s3  8.986508e-01
## s4  8.789500e-01
## s5  8.949436e-01
## s6  9.477674e-01
## s7  1.050046e+00
## s8  9.488853e-01
## s9  9.671890e-01
## s10 9.777809e-01
## s11 1.117223e+00
## s12 9.485809e-01
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
##          Point Forecast    Lo 95    Hi 95
## May 2022       250783.2 206916.8 294649.7
## Jun 2022       238220.6 193218.1 283223.1
## Jul 2022       232568.8 186457.9 278679.7
## Aug 2022       228559.7 181384.1 275735.2
## Sep 2022       233827.9 185312.7 282343.0
## Oct 2022       248804.2 198479.1 299129.3
## Nov 2022       276955.4 223949.1 329961.7
## Dec 2022       251449.9 199109.2 303790.5
## Jan 2023       257499.0 203752.7 311245.4
## Feb 2023       261530.9 206527.5 316534.2
## Mar 2023       300212.8 240973.3 359452.3
## Apr 2023       256072.1 221085.1 291059.1
PronosticosHW %>% autoplot()

Aproximación de espacio de los estados

PronosticosHW2<-hw(y = Serie_HCT[,2],
                   h = 12,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
PronosticosHW2
##          Point Forecast    Lo 95    Hi 95
## May 2022       244622.8 170166.0 319079.6
## Jun 2022       244390.9 168759.5 320022.3
## Jul 2022       221197.8 151640.5 290755.0
## Aug 2022       228978.4 155855.4 302101.5
## Sep 2022       222672.1 150496.1 294848.1
## Oct 2022       242786.0 162949.5 322622.5
## Nov 2022       270194.0 180098.8 360289.2
## Dec 2022       238763.4 158067.5 319459.2
## Jan 2023       257765.2 169500.7 346029.7
## Feb 2023       235361.5 153739.8 316983.3
## Mar 2023       287089.4 186295.0 387883.7
## Apr 2023       264957.2 170813.3 359101.1
PronosticosHW2 %>% autoplot()

identificando orden de integración

library(kableExtra)
d<-ndiffs(Serie_HCT[,2])
D<-nsdiffs(Serie_HCT[,2])
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 0
D 0

Grafico de la serie estacionaria

library(TSstudio)
Serie_HCT[,2] %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "IPC estacionaria")

Valores para (p,q) y (P,Q)

Serie_HCT[,2] %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_cor(lag.max = 36)

En caso observando el PACF los Componentes autoregresivos el valor de p=2 ya que el retardo 4 no es significativo mientras que P=0 ya que el retardo 12 no es tan significativo mientras que los Componentes de media movil los valores de Q=1 y q=1 ya que el retardo 5 no es tan significativo

Estimando SARIMA

\[SARIMA(1,1,1)(2,1,0)_{[12]}\]

library(ggthemes)
modelo_estimado <- Serie_HCT[,2] %>% 
  Arima(order = c(1, 1, 1),
        seasonal = c(2, 1, 0))

summary(modelo_estimado)
## Series: . 
## ARIMA(1,1,1)(2,1,0)[12] 
## 
## Coefficients:
##          ar1      ma1     sar1     sar2
##       0.0555  -0.7634  -0.7556  -0.2893
## s.e.  0.0159   0.0828   0.0133   0.0151
## 
## sigma^2 = 1.499e+09:  log likelihood = -612.98
## AIC=1235.95   AICc=1237.29   BIC=1245.61
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE       ACF1
## Training set 3549.924 33174.07 24583.42 0.3676404 11.56442 0.615288 -0.0229271
modelo_estimado %>% autoplot(type="both")+theme_solarized()

modelo_estimado %>% check_res(lag.max = 36)

Para el valor de MAPE es de 11.56% lo que quiere decir que de cada 100 valores pronosticados el error cometido sera de 11.56% y ante los residuos estamos en presencia de ruido blanco

Grafico comparativo

HC_Sarima<-modelo_estimado$fitted
HC_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Serie_HCT[,2],HC_Sarima,HC_HW)
ts_plot(grafico_comparativo)