knitr::opts_chunk$set(echo = TRUE)

#IVAE SLV 2009 - Marzo 2022.

#Pronóstico 2022 completo usando modelo SARIMA.

#Usando metodología paso a paso.

library(tidyr)
library(dplyr)
library(forecast)
library(readxl)
IVAE_03_22 <- read_excel("C:/Users/isaac/OneDrive/Escritorio/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]

#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 el 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%, se entiende 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)

#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 

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

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

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

#PIB trimestral precios corrientes SLV 2009-2021. #Pronóstico 2022 completo usando modelo SARIMA. #Usando la metodología paso a paso.

library(readxl)
library(forecast)
serie.pib <- 
  read_excel("C:/Users/isaac/OneDrive/Escritorio/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")

#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.3454 -0.5083 0.0665 s.e. 0.5413 0.5373 0.3901

sigma^2 = 79135: 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.557001 258.7677 120.9477 -0.004254091 2.054141 0.3680989 ACF1 Training set -0.1014549 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)

#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 

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
1 arima_010_110 Orders 69639. -332. 669. 670. 675. <cpl [4]> <cpl [4]> 2 arima_010_011 Orders 69840. -332. 670. 670. 675. <cpl [0]> <cpl [8]> 3 arima_original Orders 71060. -332. 671. 672. 679. <cpl [4]> <cpl [8]> 4 arima_automatico Orders 56463. -334. 676. 677. 684. <cpl [1]> <cpl [4]>

#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 1 ARIMA(PIB ~ pdq(0, 1, … Test 39.4 377. 197. 0.368 3.20 0.600 0.794 0.0347 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.

#IPC general SLV 2009-2022 (mayo). #Pronóstico 2022 completo usando modelo SARIMA. #Usando la metodología paso a paso.

library(readxl)
library(forecast)
serie.ipc <-
  read_excel("C:/Users/isaac/OneDrive/Escritorio/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")

#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.9437 -0.7979 -0.9390 1.0743 0.0586 -0.0962 s.e. 0.0775 0.1300 0.2112 0.2465 0.1290 0.0877

sigma^2 = 0.2226: log likelihood = -104.7 AIC=223.41 AICc=224.14 BIC=244.93

Training set error measures: ME RMSE MAE MPE MAPE MASE Training set 0.0559556 0.4614454 0.3030184 0.04969345 0.2778146 0.1642047 ACF1 Training set 0.07746987 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)

#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 arima_010_110 arima_automatico 1 <ARIMA(3,1,0)(2,0,0)[12]>

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

A tibble: 1 × 9

Model name .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
1 arima_automatico Orders 0.218 -105. 223. 224. 242. <cpl [27]> <cpl [0]>

#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 1 ARIMA(IPC ~ pdq(1, … Test -0.0230 0.275 0.208 -0.0207 0.185 0.113 0.101 0.264 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.