1. IVAE SLV 2009 - Marzo 2022

1.1. Pronóstico para 2022 completo usando un modelo SARIMA

a) Usando la metodología paso a paso

Cargar serie

library(readxl)
library(tidyr)
library(dplyr)
library(forecast)
IVAE_03_22 <- read_excel("IVAE_03_22.xlsx", 
    col_names = FALSE, skip = 6, n_max = 10)

data.ivae<-pivot_longer(data = IVAE_03_22[1,],names_to = "vars",cols = 2:160,values_to = "indice") %>% select("indice")
serie.ivae.ts<- data.ivae %>% ts(start = c(2009,1),frequency = 12)


Descomposición Clásica Aditiva

#Componente TCt
ma2_12_1 <- ma(serie.ivae.ts, 12, centre = TRUE)

#Componente St
Yt_1 <- serie.ivae.ts
Tt_1 <- ma2_12_1 
SI_1 <- Yt_1 - Tt_1 
St_1 <- tapply(SI_1, cycle(SI_1), mean, na.rm = TRUE)
St_1 <- St_1 - sum(St_1) / 12 
St_1 <-
  rep(St_1, len = length(Yt_1)) %>% ts(start = c(2009, 1), frequency = 12) 

#Componente It
It_1<-Yt_1-Tt_1-St_1

#Descomposición aditiva
library(tsibble)
library(feasts)
library(ggplot2)
Yt_1 %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  )  %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, IVAE")+xlab("Años/Meses")


Serie IVAE

library(TSstudio)
library(forecast)

ts_plot(Yt_1,Xtitle = "Años/Meses")


Orden de integración

library(kableExtra)
library(magrittr)
d_1<-ndiffs(Yt_1)
D_1<-nsdiffs(Yt_1)
ordenes_integracion_1<-c(d_1,D_1)
names(ordenes_integracion_1)<-c("d","D")
ordenes_integracion_1 %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 1
D 1


Gráfico de la serie diferenciada

Yt_1 %>% 
  diff(lag = 12,diffences=D_1) %>% 
  diff(diffences=d_1) %>% 
  ts_plot(title = "Yt estacionaria")


Verificar los valores para (p,q) & (P,Q)

Yt_1 %>% 
  diff(lag = 12,diffences=D_1) %>% 
  diff(diffences=d_1) %>% 
  ts_cor(lag.max = 36)

Se concluye que:

  • “p” es 0
  • “P” es 1
  • “q” es 0
  • “Q” es 2

Con esto se estimará el modelo SARIMA:

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


b) Usando forecast.

library(forecast)
library(ggthemes)
modelo_estimado_1 <- Yt_1 %>% 
  Arima(order = c(0, 1, 0),
        seasonal = c(1, 1, 2))

summary(modelo_estimado_1)
## Series: . 
## ARIMA(0,1,0)(1,1,2)[12] 
## 
## Coefficients:
##          sar1     sma1     sma2
##       -0.1092  -0.7162  -0.0056
## s.e.   1.0571   1.0578   0.8278
## 
## sigma^2 = 6.826:  log likelihood = -351.22
## AIC=710.43   AICc=710.72   BIC=722.37
## 
## Training set error measures:
##                     ME   RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 0.0520504 2.4777 1.685366 0.0278499 1.663524 0.4366388 0.06841698


En la práctica nuestro modelo depende del valor del mismo mes en el año anterior y del error de medición que se cometió en el año anterior.

El MAPE es de 1.66%, quiere decir que de cada 100 valores pronosticados, la distancia o el error cometido es de 1.66%. 


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

Se puede comprobar que el modelo es estable, o se cumple el teorema de invertibilidad si los puntos se encuentran dentro del círculo. 


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


Yt_Sarima_1<-modelo_estimado_1$fitted
grafico_comparativo_1<-cbind(Yt_1,Yt_Sarima_1)
ts_plot(grafico_comparativo_1)


c) Verificación de sobre ajuste/sub ajuste

Se estimarán los modelos: Partiendo del modelo original

\(SARIMA(0,1,0)(1,1,2)_{[12]}\), se estima un nuevo modelo con P-1

\(SARIMA(0,1,0)(0,1,2)_{[12]}\), y otro con Q-1:

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


library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a_1<-Yt_1 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 2)),
        arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 2)),
        arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a_1)
## # A mable: 1 x 4
##              arima_original             arima_010_011             arima_010_110
##                     <model>                   <model>                   <model>
## 1 <ARIMA(0,1,0)(1,1,2)[12]> <ARIMA(0,1,0)(0,1,2)[12]> <ARIMA(0,1,0)(1,1,1)[12]>
## # … with 1 more variable: arima_automatico <model>


a_1 %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) -> tabla_1
tabla_1
## # A tibble: 4 × 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   6.13   -347.  702.  702.  714. <cpl [1]>  <cpl [12]>
## 2 arima_010_110    Orders   6.78   -351.  708.  709.  717. <cpl [12]> <cpl [12]>
## 3 arima_010_011    Orders   6.78   -351.  708.  709.  717. <cpl [0]>  <cpl [24]>
## 4 arima_original   Orders   6.83   -351.  710.  711.  722. <cpl [12]> <cpl [24]>


1.2. Validación cruzada para el modelo estimado en 1

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

Yt_1<-Yt_1 %>% as_tsibble() %>% rename(IVAE = value)

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

TSCV_1<-data.cross.validation_1 %>% 
model(ARIMA(IVAE ~ pdq(0, 1, 0) + PDQ(1, 1, 2))) %>% 
forecast(h=1) %>% accuracy(Yt_1)

print(TSCV_1)
## # A tibble: 1 × 10
##   .model                  .type    ME  RMSE   MAE    MPE  MAPE  MASE RMSSE  ACF1
##   <chr>                   <chr> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA(IVAE ~ pdq(0, 1,… Test  0.109  3.10  2.19 0.0524  2.09 0.566 0.521 0.185

Podemos observar que el MAPE es de 2.08%, a diferencia del anterior que era de 1.66% en la estimación del modelo, podríamos decir entonces que nuestro modelo sigue teniendo un gran poder predictivo.


2. PIB trimestral precios corrientes SLV 2009-2021

2.1. Pronóstico para 2022 completo usando un modelo SARIMA

a) Usando la metodología paso a paso

Cargar serie

library(readxl)
library(forecast)
serie.pib <-
  read_excel("PIB_trimestral_SLV_2021.xlsx",
             col_types = c("skip", "numeric"))

serie.pib.ts <- ts(data = serie.pib,
                    start = c(2009, 1),
                    frequency = 4)


Descomposición Clásica Aditiva

#Componente TCt
ma2_12_2 <- ma(serie.pib.ts, 4, centre = TRUE)

#Componente St
Yt_2 <- serie.pib.ts
Tt_2 <- ma2_12_2 
SI_2 <- Yt_2 - Tt_2 
St_2 <- tapply(SI_2, cycle(SI_2), mean, na.rm = TRUE)
St_2 <- St_2 - sum(St_2) / 4 
St_2 <-
  rep(St_2, len = length(Yt_2)) %>% ts(start = c(2009, 1), frequency = 4) 

#Componente It
It_2<-Yt_2-Tt_2-St_2

#Descomposición aditiva
library(tsibble)
library(feasts)
library(ggplot2)
Yt_2 %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, PIB")+xlab("Años/Meses")


Serie PIB trimestral

library(TSstudio)
library(forecast)

ts_plot(Yt_2,Xtitle = "Años/Meses")


Orden de integración

library(kableExtra)
library(magrittr)
d_2<-ndiffs(Yt_2)
D_2<-nsdiffs(Yt_2)
ordenes_integracion_2<-c(d_2,D_2)
names(ordenes_integracion_2)<-c("d","D")
ordenes_integracion_2 %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 1
D 1


Gráfico de la serie diferenciada

Yt_2 %>% 
  diff(lag = 4,diffences=D_2) %>% 
  diff(diffences=d_2) %>% 
  ts_plot(title = "Yt estacionaria")


Verificar los valores para (p,q) & (P,Q)

Yt_2 %>% 
  diff(lag = 4,diffences=D_2) %>% 
  diff(diffences=d_2) %>% 
  ts_cor(lag.max = 12)

Se concluye que:

  • “p” es 0
  • “P” es 1
  • “q” es 0
  • “Q” es 2

Con esto se estimará el modelo SARIMA:

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


b) Usando forecast.

library(forecast)
library(ggthemes)
modelo_estimado_2 <- Yt_2 %>% 
  Arima(order = c(0, 1, 0),
        seasonal = c(1, 1, 2))

summary(modelo_estimado_2)
## Series: . 
## ARIMA(0,1,0)(1,1,2)[4] 
## 
## Coefficients:
##          sar1     sma1    sma2
##       -0.3494  -0.5042  0.0640
## s.e.   0.5359   0.5323  0.3901
## 
## sigma^2 = 79136:  log likelihood = -331.65
## AIC=671.31   AICc=672.26   BIC=678.71
## 
## Training set error measures:
##                    ME     RMSE      MAE          MPE     MAPE      MASE
## Training set 6.615822 258.7679 121.0037 -0.002908518 2.055323 0.3682692
##                    ACF1
## Training set -0.1013843


En la práctica nuestro modelo depende del valor del mismo mes en el año anterior y del error de medición que se cometió en el año anterior.

El MAPE es de 2.05%, quiere decir que de cada 100 valores pronosticados, la distancia o el error cometido es de 2.05%. 


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

Se puede comprobar que el modelo es estable, o se cumple el teorema de invertibilidad si los puntos se encuentran dentro del círculo. 


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


Yt_Sarima_2<-modelo_estimado_2$fitted
grafico_comparativo_2<-cbind(Yt_2,Yt_Sarima_2)
ts_plot(grafico_comparativo_2)


c) Verificación de sobre ajuste/sub ajuste

Se estimarán los modelos: Partiendo del modelo original

\(SARIMA(0,1,0)(1,1,2)_{[4]}\), se estima un nuevo modelo con P-1

\(SARIMA(0,1,0)(0,1,2)_{[4]}\), y otro con Q-1:

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


library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a_2<-Yt_2 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 2)),
        arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 2)),
        arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a_2)
## # A mable: 1 x 4
##             arima_original            arima_010_011            arima_010_110
##                    <model>                  <model>                  <model>
## 1 <ARIMA(0,1,0)(1,1,2)[4]> <ARIMA(0,1,0)(0,1,2)[4]> <ARIMA(0,1,0)(1,1,1)[4]>
## # … with 1 more variable: arima_automatico <model>


a_2 %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) -> tabla_2
tabla_2
## # A tibble: 4 × 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_110    Orders 77377.   -332.  669.  670.  675. <cpl [4]> <cpl [4]>
## 2 arima_010_011    Orders 77738.   -332.  670.  670.  675. <cpl [0]> <cpl [8]>
## 3 arima_original   Orders 79136.   -332.  671.  672.  679. <cpl [4]> <cpl [8]>
## 4 arima_automatico Orders 62738.   -334.  676.  677.  683. <cpl [1]> <cpl [4]>


2.2. Validación cruzada para el modelo estimado en 1

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

Yt_2<-Yt_2 %>% as_tsibble() %>% rename(PIB = value)

data.cross.validation_2<-Yt_2 %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 20,.step = 1)

TSCV_2<-data.cross.validation_2 %>% 
model(ARIMA(PIB ~ pdq(0, 1, 0) + PDQ(1, 1, 2))) %>% 
forecast(h=1) %>% accuracy(Yt_2)

print(TSCV_2)
## # A tibble: 1 × 10
##   .model                  .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <chr>                   <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 ARIMA(PIB ~ pdq(0, 1, … Test   52.6  416.  230. 0.547  3.71 0.698 0.877 0.0492

Podemos observar que el MAPE es de 3.71%, a diferencia del anterior que era de 2.05% en la estimación del modelo, podríamos decir entonces que nuestro modelo sigue teniendo un gran poder predictivo.


3. IPC general SLV 2009-2022 (mayo)

3.1. Pronóstico para 2022 completo usando un modelo SARIMA

a) Usando la metodología paso a paso

Cargar serie

library(readxl)
library(forecast)
serie.ipc <-
  read_excel("IPC_general_SLV_2022.xlsx",
             col_types = c("skip", "numeric"))

serie.ipc.ts <- ts(data = serie.ipc,
                    start = c(2009, 1),
                    frequency = 12)


Descomposición Clásica Aditiva

#Componente TCt
ma2_12_3 <- ma(serie.ipc.ts, 12, centre = TRUE)

#Componente St
Yt_3 <- serie.ipc.ts
Tt_3 <- ma2_12_3 
SI_3 <- Yt_3 - Tt_3 
St_3 <- tapply(SI_3, cycle(SI_3), mean, na.rm = TRUE)
St_3 <- St_3 - sum(St_3) / 12 
St_3 <-
  rep(St_3, len = length(Yt_3)) %>% ts(start = c(2009, 1), frequency = 12) 

#Componente It
It_3<-Yt_3-Tt_3-St_3

#Descomposición aditiva
library(tsibble)
library(feasts)
library(ggplot2)
Yt_3 %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, IPC")+xlab("Años/Meses")


Serie IPC

library(TSstudio)
library(forecast)

ts_plot(Yt_3,Xtitle = "Años/Meses")


Orden de integración

library(kableExtra)
library(magrittr)
d_3<-ndiffs(Yt_3)
D_3<-nsdiffs(Yt_3)
ordenes_integracion_3<-c(d_3,D_3)
names(ordenes_integracion_3)<-c("d","D")
ordenes_integracion_3 %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 1
D 0


Gráfico de la serie diferenciada

Yt_3 %>% 
  diff(lag = 12,diffences=D_3) %>% 
  diff(diffences=d_3) %>% 
  ts_plot(title = "Yt estacionaria")


Verificar los valores para (p,q) & (P,Q)

Yt_3 %>% 
  diff(lag = 12,diffences=D_3) %>% 
  diff(diffences=d_3) %>% 
  ts_cor(lag.max = 36)

Se concluye que:

  • “p” es 1
  • “P” es 1
  • “q” es 1
  • “Q” es 3

Con esto se estimará el modelo SARIMA:

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


b) Usando forecast.

library(forecast)
library(ggthemes)
modelo_estimado_3 <- Yt_3 %>% 
  Arima(order = c(1, 1, 1),
        seasonal = c(1, 0, 3))

summary(modelo_estimado_3)
## Series: . 
## ARIMA(1,1,1)(1,0,3)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1     sma2    sma3
##       0.9289  -0.7894  0.9905  -0.9135  -0.1630  0.1131
## s.e.  0.0840   0.1307  0.0416   0.1244   0.1065  0.0880
## 
## sigma^2 = 0.2175:  log likelihood = -104.18
## AIC=222.37   AICc=223.1   BIC=243.89
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.05115671 0.4561007 0.3004523 0.04545039 0.2756045 0.1628142
##                    ACF1
## Training set 0.07731549


En la práctica nuestro modelo depende del valor del mismo mes en el año anterior y del error de medición que se cometió en el año anterior.

El MAPE es de 0.27%, quiere decir que de cada 100 valores pronosticados, la distancia o el error cometido es de 0.27%. 


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

Se puede comprobar que el modelo es estable, o se cumple el teorema de invertibilidad si los puntos se encuentran dentro del círculo. 


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


Yt_Sarima_3<-modelo_estimado_3$fitted
grafico_comparativo_3<-cbind(Yt_3,Yt_Sarima_3)
ts_plot(grafico_comparativo_3)


c) Verificación de sobre ajuste/sub ajuste

Se estimarán los modelos: Partiendo del modelo original

\(SARIMA(1,1,1)(1,0,3)_{[12]}\), se estima un nuevo modelo con P-1

\(SARIMA(1,1,1)(0,0,3)_{[12]}\), y otro con Q-1:

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


library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a_3<-Yt_3 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 3)),
        arima_010_011 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(0, 0, 3)),
        arima_010_110 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 2)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a_3)
## # A mable: 1 x 4
##                       arima_original                      arima_010_011
##                              <model>                            <model>
## 1 <ARIMA(1,1,1)(1,0,3)[12] w/ drift> <ARIMA(1,1,1)(0,0,3)[12] w/ drift>
## # … with 2 more variables: arima_010_110 <model>, arima_automatico <model>


a_3 %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) -> tabla_3
tabla_3
## # A tibble: 4 × 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  0.219   -103.  220.  221.  241. <cpl [1]>  <cpl [37]>
## 2 arima_010_110    Orders  0.219   -103.  220.  221.  241. <cpl [13]> <cpl [25]>
## 3 arima_original   Orders  0.221   -103.  222.  223.  246. <cpl [13]> <cpl [37]>
## 4 arima_automatico Orders  0.225   -105.  223.  224.  241. <cpl [27]> <cpl [0]>


3.2. Validación cruzada para el modelo estimado en 1

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

Yt_3<-Yt_3 %>% as_tsibble() %>% rename(IPC = value)

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

TSCV_3<-data.cross.validation_3 %>% 
model(ARIMA(IPC ~ pdq(1, 1, 1) + PDQ(1, 0, 3))) %>% 
forecast(h=1) %>% accuracy(Yt_3)

print(TSCV_3)
## # A tibble: 1 × 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(1, 1,… Test  0.0411 0.349 0.265 0.0330 0.232 0.143 0.128 0.398

Podemos observar que el MAPE es de 0.23%, a diferencia del anterior que era de 0.27% en la estimación del modelo, podríamos decir entonces que nuestro modelo sigue teniendo un gran poder predictivo.


4. Importación de Hidrocarburos SLV 2017-2022 (mayo)

4.1. Pronóstico para 2022 completo usando un modelo SARIMA

a) Usando la metodología paso a paso

Cargar serie

library(readxl)
library(forecast)
serie.hid <-
  read_excel("Hidrocarburos_SLV_2022.xlsx",
             col_types = c("skip", "numeric"))

serie.hid.ts <- ts(data = serie.hid,
                    start = c(2017, 1),
                    frequency = 12)


Descomposición Clásica Aditiva

#Componente TCt
ma2_12_4 <- ma(serie.hid.ts, 12, centre = TRUE)

#Componente St
Yt_4 <- serie.hid.ts
Tt_4 <- ma2_12_4 
SI_4 <- Yt_4 - Tt_4 
St_4 <- tapply(SI_4, cycle(SI_4), mean, na.rm = TRUE)
St_4 <- St_4 - sum(St_4) / 12 
St_4 <-
  rep(St_4, len = length(Yt_4)) %>% ts(start = c(2017, 1), frequency = 12) 

#Componente It
It_4<-Yt_4-Tt_4-St_4

#Descomposición aditiva
library(tsibble)
library(feasts)
library(ggplot2)
Yt_4 %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, Importación de Hidrocarburos")+xlab("Años/Meses")


Serie de importación de Hidrocarburos

library(TSstudio)
library(forecast)

ts_plot(Yt_4,Xtitle = "Años/Meses")


Orden de integración

library(kableExtra)
library(magrittr)
d_4<-ndiffs(Yt_4)
D_4<-nsdiffs(Yt_4)
ordenes_integracion_4<-c(d_4,D_4)
names(ordenes_integracion_4)<-c("d","D")
ordenes_integracion_4 %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d 0
D 0


Gráfico de la serie diferenciada

Yt_4 %>% 
  diff(lag = 12,diffences=D_4) %>% 
  diff(diffences=d_4) %>% 
  ts_plot(title = "Yt estacionaria")


Verificar los valores para (p,q) & (P,Q)

Yt_4 %>% 
  diff(lag = 12,diffences=D_4) %>% 
  diff(diffences=d_4) %>% 
  ts_cor(lag.max = 36)

Se concluye que:

  • “p” es 1
  • “P” es 1
  • “q” es 2
  • “Q” es 1

Con esto se estimará el modelo SARIMA:

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


b) Usando forecast.

library(forecast)
library(ggthemes)
modelo_estimado_4 <- Yt_4 %>% 
  Arima(order = c(1, 0, 2),
        seasonal = c(1, 0, 1))

summary(modelo_estimado_4)
## Series: . 
## ARIMA(1,0,2)(1,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     ma2     sar1   sma1        mean
##       0.7427  -0.6214  0.1503  -0.9949  0.953  216942.159
## s.e.  0.1939   0.2269  0.1332   0.0497  0.231    7263.433
## 
## sigma^2 = 941483707:  log likelihood = -752.39
## AIC=1518.78   AICc=1520.78   BIC=1533.9
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 105.1099 29209.92 23410.96 -1.982901 11.27784 0.5859431
##                     ACF1
## Training set 0.006624694


En la práctica nuestro modelo depende del valor del mismo mes en el año anterior y del error de medición que se cometió en el año anterior.

El MAPE es de 11.27%, quiere decir que de cada 100 valores pronosticados, la distancia o el error cometido es de 11.27%. 


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

Se puede comprobar que el modelo es estable, o se cumple el teorema de invertibilidad si los puntos se encuentran dentro del círculo. 


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


Yt_Sarima_4<-modelo_estimado_4$fitted
grafico_comparativo_4<-cbind(Yt_4,Yt_Sarima_4)
ts_plot(grafico_comparativo_4)


c) Verificación de sobre ajuste/sub ajuste

Se estimarán los modelos: Partiendo del modelo original

\(SARIMA(1,0,2)(1,0,1)_{[12]}\), se estima un nuevo modelo con P-1

\(SARIMA(1,0,2)(0,0,1)_{[12]}\), y otro con Q-1:

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


library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a_4<-Yt_4 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(1, 0, 2) + PDQ(1, 0, 1)),
        arima_010_011 = ARIMA(value ~ pdq(1, 0, 2) + PDQ(0, 0, 1)),
        arima_010_110 = ARIMA(value ~ pdq(1, 0, 2) + PDQ(1, 0, 0)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a_4)
## # A mable: 1 x 4
##   arima_original                     arima_010_011
##          <model>                           <model>
## 1   <NULL model> <ARIMA(1,0,2)(0,0,1)[12] w/ mean>
## # … with 2 more variables: arima_010_110 <model>, arima_automatico <model>


a_4 %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) -> tabla_4
tabla_4
## # A tibble: 3 × 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     1.10e9   -757. 1517. 1517. 1521. <cpl>    <cpl>   
## 2 arima_010_110    Orders     1.05e9   -753. 1519. 1520. 1531. <cpl>    <cpl>   
## 3 arima_010_011    Orders     1.06e9   -753. 1519. 1520. 1532. <cpl>    <cpl>


4.2. Validación cruzada para el modelo estimado en 1

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

Yt_4<-Yt_4 %>% as_tsibble() %>% rename(Hidrocarburo = value)

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

TSCV_4<-data.cross.validation_4 %>% 
model(ARIMA(Hidrocarburo ~ pdq(1, 0, 2) + PDQ(1, 0, 1))) %>% 
forecast(h=1) %>% accuracy(Yt_4)

print(TSCV_4)
## # A tibble: 1 × 10
##   .model             .type     ME   RMSE    MAE   MPE  MAPE  MASE RMSSE     ACF1
##   <chr>              <chr>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>
## 1 ARIMA(Hidrocarbur… Test  44092. 58252. 47296.  15.3  16.8  1.18  1.19 -0.00562

Podemos observar que el MAPE es de 16.80%, a diferencia del anterior que era de 11.27% en la estimación del modelo, podríamos decir entonces que nuestro modelo sigue teniendo un gran poder predictivo.