Ejercicio de Pronóstico Modelos SARIMA

Serie temporal PIB trimestral

#carga del archivo fuente
       library(dplyr)
    library(tidyr)
    library(forecast)
  library(TSstudio)
  library(readxl)
PIB_trimestral<- read_excel("C:/Users/liizm/documents/econometria/PIB trimestral.xlsx",col_names = FALSE,col_types = c("skip", "numeric"), skip = 6)


PIB_trimestral.ts <- ts(data = PIB_trimestral,
                    start = c(2009, 1),
                    frequency = 4)
PIB_trimestral.ts %>% 
  autoplot(main = "PIB trimestral precios corrientes SLV 2009-2021",
                           xlab = "Años/Trimestres",
                           ylab = "Indice")

library(forecast)
  library(TSstudio)
   Yt <-PIB_trimestral.ts
ts_plot(Yt,Xtitle = "Años/Trimestres")

Estimacion de modelo HW

library(forecast)

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.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

Gráfico de la serie original y del pronóstico.

PronosticosHW %>% autoplot()

Orden de integración de la serie

library(kableExtra)
library(magrittr)
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

Gráfico de la serie diferenciada

Yt %>% 
  diff(lag = 4,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "Yt estacionaria")

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

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

Estimación del modelo Sarima usando la líbreria forecast

library(forecast)
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)[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)

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

Verificación de sobre ajuste/sub ajuste

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a<-Yt %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
        arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)),
        arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 0)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a)

## # A mable: 1 x 4
##             arima_original            arima_010_011            arima_010_110
##                    <model>                  <model>                  <model>
## 1 <ARIMA(0,1,0)(1,1,1)[4]> <ARIMA(0,1,0)(0,1,1)[4]> <ARIMA(0,1,0)(1,1,0)[4]>
## # ... with 1 more variable: arima_automatico <model>

a %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla
tabla

## # A tibble: 4 x 9
##   `Model name`     .model sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <chr>            <chr>   <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 arima_010_011    Orders 75930.   -332.  668.  669.  672. <cpl [0]> <cpl [4]>
## 2 arima_010_110    Orders 78594.   -333.  669.  670.  673. <cpl [4]> <cpl [0]>
## 3 arima_original   Orders 77377.   -332.  669.  670.  675. <cpl [4]> <cpl [4]>
## 4 arima_automatico Orders 62738.   -334.  676.  677.  683. <cpl [1]> <cpl [4]>

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt<-Yt %>% as_tsibble() %>% rename(PIB_trimestral=value)

data.cross.validation<-Yt %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 24,.step = 1)

TSCV<-data.cross.validation %>% 
model(ARIMA(PIB_trimestral ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt)

print(TSCV)

## # A tibble: 1 x 10
##   .model                  .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <chr>                   <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 ARIMA(PIB_trimestral ~~ Test   43.4  394.  215. 0.422  3.48 0.654 0.830 0.0397

Serie temporal hidrocarburos

#carga del archivo fuente
       library(dplyr)
    library(tidyr)
    library(forecast)
  library(TSstudio)
  library(readxl)
hidrocarburos<- read_excel("C:/Users/liizm/documents/econometria/hidrocarburos.xlsx",col_names = FALSE,col_types = c("skip", "numeric"), skip = 6)


hidroc.ts <- ts(data = hidrocarburos,
                    start = c(2017, 1),
                    frequency = 12)
hidroc.ts %>% 
  autoplot(main = "Hidrocarburos 2017 - 2022",
                           xlab = "Años/Meses",
                           ylab = "Volumen")

library(forecast)
  library(TSstudio)
   Yt2 <-hidroc.ts
ts_plot(Yt2,Xtitle = "Años/Meses")

Estimacion de modelo HW

library(forecast)

ModeloHW2<-HoltWinters(x = Yt2,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
ModeloHW2

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Yt2, 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

PronosticosHW2<-forecast(object = ModeloHW2,h = 12,level = c(0.95))
PronosticosHW2

##          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

Gráfico de la serie original y del pronóstico.

PronosticosHW2 %>% autoplot()

Orden de integración de la serie

library(kableExtra)
library(magrittr)
d2<-ndiffs(Yt2)
D2<-nsdiffs(Yt2)
ordenes_integracion2<-c(d2,D2)
names(ordenes_integracion2)<-c("d2","D2")
ordenes_integracion2 %>% kable(caption = "Ordenes de Integración 2") %>% kable_material()

Ordenes de Integración 2

	x
d2	0
D2	0

Gráfico de la serie diferenciada

Yt2 %>% 
  diff(lag = 12,diffences=D2) %>% 
  diff(diffences=d2) %>% 
  ts_plot(title = "Yt2 estacionaria")

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

Yt2 %>% 
  diff(lag = 12,diffences=D2) %>% 
  diff(diffences=d2) %>% 
  ts_cor(lag.max = 64)

Estimación del modelo Arima usando la líbreria forecast

library(forecast)
library(ggthemes)
modelo_estimado2 <- Yt2 %>% 
  Arima(order = c(0, 1, 0),
        seasonal = c(1, 1, 1))

summary(modelo_estimado2)

## Series: . 
## ARIMA(0,1,0)(1,1,1)[12] 
## 
## Coefficients:
##          sar1     sma1
##       -0.2376  -0.6941
## s.e.   0.2483   0.3848
## 
## sigma^2 = 1.936e+09:  log likelihood = -622.57
## AIC=1251.14   AICc=1251.65   BIC=1256.94
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 1657.922 38496.49 27289.55 0.09035022 12.73477 0.6830186
##                    ACF1
## Training set -0.4642396

modelo_estimado2 %>% autoplot(type="both")+theme_solarized()

modelo_estimado2 %>% check_res(lag.max = 64)

Yt_Sarima2<-modelo_estimado2$fitted
Yt_HW2<-PronosticosHW2$fitted
grafico_comparativo2<-cbind(Yt2,Yt_Sarima2,Yt_HW2)
ts_plot(grafico_comparativo2)

Verificación de sobre ajuste/sub ajuste

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a2<-Yt2 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
        arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)),
        arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 0)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a2)

## # A mable: 1 x 4
##              arima_original             arima_010_011             arima_010_110
##                     <model>                   <model>                   <model>
## 1 <ARIMA(0,1,0)(1,1,1)[12]> <ARIMA(0,1,0)(0,1,1)[12]> <ARIMA(0,1,0)(1,1,0)[12]>
## # ... with 1 more variable: arima_automatico <model>

a2 %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla2
tabla2

## # A tibble: 4 x 9
##   `Model name`     .model     sigma2 log_lik   AIC  AICc   BIC ar_roots ma_roots
##   <chr>            <chr>       <dbl>   <dbl> <dbl> <dbl> <dbl> <list>   <list>  
## 1 arima_010_011    Orders     1.65e9   -623. 1250. 1250. 1254. <cpl>    <cpl>   
## 2 arima_original   Orders     1.94e9   -623. 1251. 1252. 1257. <cpl>    <cpl>   
## 3 arima_010_110    Orders     2.34e9   -624. 1253. 1253. 1257. <cpl>    <cpl>   
## 4 arima_automatico Orders     1.10e9   -757. 1517. 1517. 1521. <cpl>    <cpl>

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt2<-Yt2 %>% as_tsibble() %>% rename(Hidrocarburos=value)

data.cross.validation2<-Yt2 %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 60,.step = 1)

TSCV2<-data.cross.validation2 %>% 
model(ARIMA(hidrocarburos ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=12) %>% accuracy(Yt2)

print(TSCV2)

## # A tibble: 1 x 10
##   .model                  .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <chr>                   <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 ARIMA(hidrocarburos ~ ~ Test     NA    NA    NA    NA    NA    NA    NA -0.344

Serie temporal IPC

#carga del archivo fuente
       library(dplyr)
    library(tidyr)
    library(forecast)
  library(TSstudio)
  library(readxl)
ipc<- read_excel("C:/Users/liizm/documents/econometria/IPC.xlsx",col_names = FALSE,col_types = c("skip", "numeric"), skip = 6)


ipc.ts <- ts(data = ipc,
                    start = c(2009, 1),
                    frequency = 12)
ipc.ts %>% 
  autoplot(main = "IPC 2009 - 2022",
                           xlab = "Años/Meses",
                           ylab = "Indice")

library(forecast)
  library(TSstudio)
   Yt3 <-ipc.ts
ts_plot(Yt3,Xtitle = "Años/Meses")

Estimacion de modelo HW

library(forecast)

ModeloHW3<-HoltWinters(x = Yt3,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
ModeloHW3

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Yt3, 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

PronosticosHW3<-forecast(object = ModeloHW3,h = 12,level = c(0.95))
PronosticosHW3

##          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

Gráfico de la serie original y del pronóstico.

PronosticosHW3 %>% autoplot()

Orden de integración de la serie

library(kableExtra)
library(magrittr)
d3<-ndiffs(Yt3)
D3<-nsdiffs(Yt3)
ordenes_integracion3<-c(d3,D3)
names(ordenes_integracion2)<-c("d3","D3")
ordenes_integracion3 %>% kable(caption = "Ordenes de Integración 3") %>% kable_material()

Ordenes de Integración 3

x
1
0

Gráfico de la serie diferenciada

Yt3 %>% 
  diff(lag = 12,diffences=D3) %>% 
  diff(diffences=d3) %>% 
  ts_plot(title = "Yt3 estacionaria")

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

Yt3 %>% 
  diff(lag = 12,diffences=D3) %>% 
  diff(diffences=d3) %>% 
  ts_cor(lag.max = 64)

Estimación del modelo Arima usando la líbreria forecast

library(forecast)
library(ggthemes)
modelo_estimado3 <- Yt3 %>% 
  Arima(order = c(0, 1, 0),
        seasonal = c(1, 1, 1))

summary(modelo_estimado3)

## Series: . 
## ARIMA(0,1,0)(1,1,1)[12] 
## 
## Coefficients:
##         sar1     sma1
##       0.1112  -1.0000
## s.e.  0.0950   0.1139
## 
## sigma^2 = 0.2312:  log likelihood = -114.66
## AIC=235.33   AICc=235.49   BIC=244.32
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.02246671 0.4579215 0.3010507 0.01772013 0.2729747 0.1631385
##                   ACF1
## Training set 0.2438271

modelo_estimado3 %>% autoplot(type="both")+theme_solarized()

modelo_estimado3 %>% check_res(lag.max = 64)

Yt_Sarima3<-modelo_estimado3$fitted
Yt_HW3<-PronosticosHW3$fitted
grafico_comparativo3<-cbind(Yt3,Yt_Sarima3,Yt_HW3)
ts_plot(grafico_comparativo3)

Verificación de sobre ajuste/sub ajuste

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a3<-Yt3 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
        arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)),
        arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 0)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a3)

## # A mable: 1 x 4
##              arima_original             arima_010_011             arima_010_110
##                     <model>                   <model>                   <model>
## 1 <ARIMA(0,1,0)(1,1,1)[12]> <ARIMA(0,1,0)(0,1,1)[12]> <ARIMA(0,1,0)(1,1,0)[12]>
## # ... with 1 more variable: arima_automatico <model>

a3 %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla3
tabla3

## # A tibble: 4 x 9
##   `Model name`     .model sigma2 log_lik   AIC  AICc   BIC ar_roots   ma_roots  
##   <chr>            <chr>   <dbl>   <dbl> <dbl> <dbl> <dbl> <list>     <list>    
## 1 arima_automatico Orders  0.225   -105.  223.  224.  241. <cpl [27]> <cpl [0]> 
## 2 arima_010_011    Orders  0.228   -115.  235.  235.  241. <cpl [0]>  <cpl [12]>
## 3 arima_original   Orders  0.231   -115.  235.  235.  244. <cpl [12]> <cpl [12]>
## 4 arima_010_110    Orders  0.348   -132.  268.  268.  274. <cpl [12]> <cpl [0]>

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt3<-Yt3 %>% as_tsibble() %>% rename(IPC=value)

data.cross.validation3<-Yt3 %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 24,.step = 1)

TSCV3<-data.cross.validation3 %>% 
model(ARIMA(ipc ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt3)

print(TSCV3)

## # A tibble: 1 x 10
##   .model                   .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>                    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA(ipc ~ pdq(0, 1, 0~ Test     NA    NA    NA    NA    NA    NA    NA 0.922