Ejercicio de Pronóstico Modelos SARIMA

PARTE I: PIB trimestral precios corrientes SLV 2009-2021.

Cargamos la base de datos

library(readxl)
library(forecast)

PIB_SLV<- read_excel("F:/Practica Pronostico Modelos SARIMA/TAREA/PIB trimestral SLV 2009-2021.xlsx",  col_types = c("skip","numeric"),
                            skip = 5)

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

Yt_1 <- serie.PIB.ts
serie.PIB.ts %>% 
  autoplot(main = "PIB trimestral precios corrientes SLV 2009-2021",
           xlab = "años/ trimestres",
           ylab = "Indice")

1. Gráfico de la Serie (Metodología Box-Jenkins)

library(TSstudio)
library(forecast)

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

2. Identificacion

2.1 Orden de Integracion

library(kableExtra)
library(magrittr)
d_1<-ndiffs(Yt_1)
D_1<-nsdiffs(Yt_1)
ordenes_integracion<-c(d_1,D_1)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integracion") %>% kable_material()

Ordenes de Integracion
	x
d	1
D	1

#Gráfico de la serie diferenciada
Yt_1 %>%
  diff(lag=4,diffences=D_1) %>%
  diff(diffences=d_1) %>%
  ts_plot(title= "Yt estacionaria")

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

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

2.3 Estimacion del modelo

2.3.1 Usando forecast

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

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

#grafico de las raices (estabilidad/convertibilidad)
modelo_estimado_PIB %>% autoplot(type="both")+theme_solarized()

#verificacion
modelo_estimado_PIB %>% check_res(lag.max = 4)

library(forecast)

#Estimar el modelo
modeloHW_PIB<-HoltWinters(x = Yt_1,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
modeloHW_PIB

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

#Generar el pronostico:
pronosticosMW_PIB<-forecast(object = modeloHW_PIB,h=4, level=c(0.95))
pronosticosMW_PIB

##         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 comparativo
Yt_Sarima_PIB<-modelo_estimado_PIB$fitted
Yt_HW_PIB<-pronosticosMW_PIB$fitted
grafico_comparativo_1<-cbind(Yt_1, Yt_Sarima_PIB,Yt_HW_PIB)
ts_plot(grafico_comparativo_1)

3. Verificacion de sobre ajuste/sub ajuste

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a<-Yt_1 %>% 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 × 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]>

Validacion Cruzada (Cross Validated)

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

#convertimos
Yt<-Yt_1 %>% 
  as_tsibble() %>% 
  rename(PIB=value)

#generador de la dta (roja)
data.cross.validation.PIB<-Yt %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 20,.step = 1)

#generador de pronostico (verde) y simulacion 
TSCV<-data.cross.validation.PIB %>% 
model(ARIMA(PIB ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt)

print(TSCV)

## # 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   30.1  370.  197. 0.228  3.21 0.600 0.780 0.0422

PARTE II: IPC general SLV 2009-2022 (mayo).

Cargamos la base de datos

library(readxl)
library(forecast)

IPC_SLV<- read_excel("F:/Practica Pronostico Modelos SARIMA/TAREA/IPC general SLV 2009-2022 (mayo).xlsx",  col_types = c("skip","numeric"),
                            skip = 6)

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

Yt_2 <- serie.IPC.ts
serie.IPC.ts %>% 
  autoplot(main = "IPC general SLV 2009-2022 (mayo)",
           xlab = "años/ trimestres",
           ylab = "Indice")

1. Gráfico de la Serie (Metodología Box-Jenkins)

library(TSstudio)
library(forecast)

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

2. Identificacion

2.1 Orden de Integracion

library(kableExtra)
library(magrittr)
d_2<-ndiffs(Yt_2)
D_2<-nsdiffs(Yt_2)
ordenes_integracion<-c(d_2,D_2)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integracion") %>% kable_material()

Ordenes de Integracion
	x
d	1
D	0

#Gráfico de la serie diferenciada
Yt_2 %>%
  diff(lag=12,diffences=D_2) %>%
  diff(diffences=d_2) %>%
  ts_plot(title= "Yt estacionaria")

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

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

2.3 Estimacion del modelo

2.3.1 Usando forecast

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

## Series: . 
## ARIMA(1,1,1)(1,0,2)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1    sma1     sma2
##       0.9406  -0.7934  0.0571  0.0615  -0.0683
## s.e.  0.0859   0.1447  0.6621  0.6595   0.1106
## 
## sigma^2 = 0.2251:  log likelihood = -104.87
## AIC=221.74   AICc=222.29   BIC=240.15
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.05582688 0.4654425 0.3091819 0.04952956 0.2834401 0.1668998
##                    ACF1
## Training set 0.08746749

#grafico de las raices (estabilidad/convertibilidad)
modelo_estimado_IPC %>% autoplot(type="both")+theme_solarized()

#verificacion
modelo_estimado_IPC %>% check_res(lag.max = 36)

library(forecast)

#Estimar el modelo
modeloHW_IPC<-HoltWinters(x = Yt_2,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
modeloHW_IPC

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Yt_2, seasonal = "multiplicative", optim.start = c(0.9,     0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.8419073
##  beta : 0.0743462
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   122.7484331
## b     0.4728397
## s1    1.0043101
## s2    1.0017791
## s3    0.9980880
## s4    0.9970592
## s5    0.9984875
## s6    0.9981783
## s7    0.9959355
## s8    0.9977953
## s9    1.0010455
## s10   1.0045211
## s11   1.0052312
## s12   1.0055526

#Generar el pronostico:
pronosticosMW_IPC<-forecast(object = modeloHW_IPC,h=12, level=c(0.95))
pronosticosMW_IPC

##          Point Forecast    Lo 95    Hi 95
## May 2022       123.7524 122.6054 124.8994
## Jun 2022       123.9142 122.3694 125.4590
## Jul 2022       123.9295 122.0324 125.8267
## Aug 2022       124.2733 122.0409 126.5056
## Sep 2022       124.9234 122.3615 127.4853
## Oct 2022       125.3567 122.4738 128.2396
## Nov 2022       125.5459 122.3495 128.7424
## Dec 2022       126.2522 122.7300 129.7744
## Jan 2023       127.1368 123.2817 130.9918
## Feb 2023       128.0532 123.8606 132.2457
## Mar 2023       128.6190 124.0961 133.1419
## Apr 2023       129.1356 123.8244 134.4468

#Gráfico comparativo
Yt_Sarima_IPC<-modelo_estimado_IPC$fitted
Yt_HW_IPC<-pronosticosMW_IPC$fitted
grafico_comparativo_2<-cbind(Yt_2, Yt_Sarima_IPC,Yt_HW_IPC)
ts_plot(grafico_comparativo_2)

3. Verificacion de sobre ajuste/sub ajuste

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
#generando modelo automatico
a_2<-Yt_2 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 2)),
        arima_010_011 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(0, 0, 2)),
        arima_010_110 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 1)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a_2)

## # A mable: 1 x 4
##                       arima_original                      arima_010_011
##                              <model>                            <model>
## 1 <ARIMA(1,1,1)(1,0,2)[12] w/ drift> <ARIMA(1,1,1)(0,0,2)[12] w/ drift>
## # … with 2 more variables: arima_010_110 <model>, arima_automatico <model>

#ordenar  
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_011    Orders  0.219   -102.  217.  218.  235. <cpl [1]>  <cpl [25]>
## 2 arima_automatico Orders  0.218   -102.  217.  218.  235. <cpl [26]> <cpl [0]> 
## 3 arima_010_110    Orders  0.219   -103.  217.  218.  236. <cpl [13]> <cpl [13]>
## 4 arima_original   Orders  0.220   -102.  219.  220.  240. <cpl [13]> <cpl [25]>

Validacion Cruzada (Cross Validated)

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

#convertimos
Yt<-Yt_2 %>% as_tsibble() %>% rename(IPC=value)

#generador de la dta (roja)
data.cross.validation.IPC<-Yt %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 60,.step = 1)

#generador de pronostico (verde) y simulacion 
TSCV_2<-data.cross.validation.IPC %>% 
model(ARIMA(IPC ~ pdq(1, 1, 1) + PDQ(1, 0, 2))) %>% 
forecast(h=1) %>% accuracy(Yt)

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(IPC ~ pdq(1, 1… Test  0.0695 0.437 0.295 0.0579 0.261 0.159 0.160 0.0775

PARTE III: Importación de Hidrocarburos SLV 2017-2022 (mayo).

Cargamos la base de datos

library(readxl)
library(forecast)

M_Hidrocarburos<- read_excel("F:/Practica Pronostico Modelos SARIMA/TAREA/Hidrocarburos SLV 2017-2022 (mayo).xlsx",  col_types = c("skip","numeric"),
                            skip = 7)

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

Yt_3 <- serie.M_HC.ts
serie.M_HC.ts %>% 
  autoplot(main = "Importación de Hidrocarburos SLV 2017-2022 (mayo)",
           xlab = "años/ trimestres",
           ylab = "Indice")

1. Gráfico de la Serie (Metodología Box-Jenkins)

library(TSstudio)
library(forecast)

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

2. Identificacion

2.1 Orden de Integracion

library(kableExtra)
library(magrittr)
d_3<-ndiffs(Yt_3)
D_3<-nsdiffs(Yt_3)
ordenes_integracion<-c(d_3,D_3)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integracion") %>% kable_material()

Ordenes de Integracion
	x
d	0
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")

2.2 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)

2.3 Estimacion del modelo

2.3.1 Usando forecast

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

summary(modelo_estimado_M_HC)

## Series: . 
## ARIMA(1,0,1)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     sar1        mean
##       0.6652  -0.4851  -0.2022  215790.730
## s.e.  0.2952   0.3398   0.1519    5147.504
## 
## sigma^2 = 1.042e+09:  log likelihood = -741.69
## AIC=1493.39   AICc=1494.44   BIC=1504.1
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -428.6553 31240.12 25005.59 -2.421527 12.06669 0.6217076
##                      ACF1
## Training set -0.007961808

#grafico de las raices (estabilidad/convertibilidad)
modelo_estimado_M_HC %>% autoplot(type="both")+theme_solarized()

#verificacion
modelo_estimado_M_HC %>% check_res(lag.max = 36)

library(forecast)

#Estimar el modelo
modeloHW_M_HC<-HoltWinters(x = Yt_3,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
modeloHW_M_HC

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = Yt_3, seasonal = "multiplicative", optim.start = c(0.9,     0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.1580326
##  beta : 0
##  gamma: 0.4946113
## 
## Coefficients:
##             [,1]
## a   2.501763e+05
## b   1.683628e+03
## s1  9.111372e-01
## s2  8.894296e-01
## s3  8.823914e-01
## s4  8.896954e-01
## s5  9.398533e-01
## s6  1.039846e+00
## s7  9.375516e-01
## s8  9.546773e-01
## s9  9.615739e-01
## s10 1.097096e+00
## s11 9.326123e-01
## s12 8.393040e-01

#Generar el pronostico:
pronosticosMW_M_HC<-forecast(object = modeloHW_M_HC,h=12, level=c(0.95))
pronosticosMW_M_HC

##          Point Forecast    Lo 95    Hi 95
## Apr 2022       229478.9 182755.6 276202.2
## May 2022       225509.1 177130.5 273887.7
## Jun 2022       225210.3 175204.2 275216.3
## Jul 2022       228572.3 176863.3 280281.4
## Aug 2022       243040.8 189049.4 297032.2
## Sep 2022       270649.0 213366.1 327932.0
## Oct 2022       245602.6 188944.1 302261.0
## Nov 2022       251696.2 193221.1 310171.3
## Dec 2022       255133.4 195092.9 315173.9
## Jan 2023       292938.5 227754.8 358122.3
## Feb 2023       250589.4 188983.5 312195.2
## Mar 2023       226930.8 188989.5 264872.2

#Gráfico comparativo
Yt_Sarima_M_HC<-modelo_estimado_M_HC$fitted
Yt_HW_M_HC<-pronosticosMW_M_HC$fitted
grafico_comparativo_3<-cbind(Yt_3, Yt_Sarima_M_HC,Yt_HW_M_HC)
ts_plot(grafico_comparativo_3)

3. Verificacion de sobre ajuste/sub ajuste

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)

#generando modelo automatico
a_3<-Yt_3 %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(1, 0, 1) + PDQ(1, 0, 0)),
        arima_010_011 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 0)),
        arima_010_110 = ARIMA(value ~ pdq(1, 0, 1) + PDQ(1, 1, 0)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE) #modelo automatico 
  )
print(a_3)

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

Validacion Cruzada (Cross Validated)

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

#convertimos
Yt<-Yt_3 %>% 
  as_tsibble() %>% 
  rename(M_HC= value)

#generador de la dta (roja)
data.cross.validation.M_HC<-Yt %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 60,.step = 1)

#generador de pronostico (verde) y simulacion 
TSCV_3<-data.cross.validation.M_HC %>% 
model(ARIMA(M_HC ~ pdq(1, 0, 1) + PDQ(1, 0, 1))) %>% 
forecast(h=1) %>% accuracy(Yt)

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(M_HC ~ pdq(1, … Test  82270. 82270. 82270.  27.5  27.5  2.05  1.67    NA


h1,h2,h3,h4{
    font-size:14px;
    font-weight: bold;
    font-family: 'Times New Roman', Times, serif;
    line-height:2;
}
h1 {
    text-align: center;
}
h4{
    text-align:center;
    font-weight: bold
}
p, li{
  text-indent: 40px;
    font-size:12 px;
    font-family: 'Times New Roman', Times, serif;
      line-height: 2;
}