## Cargando Datos
library(readxl)
library(forecast)
library(readxl)
PIB_trimestral_precios_corrientes_SLV_2009_20211 <- read_excel("PIB trimestral precios corrientes SLV 2009-20211.xlsx", 
    range = "B7:B58", col_names = FALSE)
## Convirtiendo datos en serie temporal

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

Metodología paso a paso

library(TSstudio)
library(forecast)

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

Identificación

Orden de Integracion

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 Integracion") %>% kable_material()

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

Verificando los valores para (p, q) y (P, Q)

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

En la parte del ACF debido a que los primeros 3 ordenes están dentro de las bandas (no son significativos), el componente autoregresivo de la parte no ordinaria el valor de p es igual a 0

PACF es para identificar valor de q y debido a que los coeficientes de autocorrelacion parcial se encuentran dentro de las bandas se puede garantizar que no hay un componente de media movil ordinaria (q)

En la parte ordinaria, el PIB trimestral para los años 2009-2021 no depende del valor inmediato anterior ni de los errores de medición.

El componente autoregresivo estacional (P) es 1

En el caso de media movil (Q), se observa que para el primer y segundo periodo resulta ser muy significativa (quiere decir que el valor actual de la serie depende de la fluctuacion del año anterior en este periodo). El componente de media movil es estacional con un valor de 2.

Estimación del modelo

Usando librería forecast

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

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

## Raices de los polinomios de la parte autoregresiva y parte media movil.
modelo_estimado %>% autoplot(type="both")+theme_solarized()

## Construye el correlograma
modelo_estimado %>% check_res(lag.max = 52)

El error porcentual absoluto medio (MAPE) es de 2.05%, quiere decir que de cada 100 valores pronosticados la distancia o error cometido es de 2.05%.

library(forecast)

#Estimar el modelo
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

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

##         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
## 2023 Q1       7481.263 6598.506 8364.020
## 2023 Q2       7765.084 6751.622 8778.546
## 2023 Q3       7600.414 6518.808 8682.020
## 2023 Q4       8134.887 6984.348 9285.426
## 2024 Q1       7670.040 6503.665 8836.415
## 2024 Q2       7959.795 6677.070 9242.520
## 2024 Q3       7789.809 6465.102 9114.515
## 2024 Q4       8336.345 6927.229 9745.462

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

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

## # 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 %>% 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_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]>

Validación Cruzada (Cross Validated)

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

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

TSCV <- data.cross.validation %>%
  model(ARIMA(PIB ~ pdq(0,1,0)+ PDQ(1,1,2))) %>%
  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 ~ pdq(0, 1, ~ Test   52.6  416.  230. 0.547  3.71 0.698 0.877 0.0492

El MAPE original era 2.05%, en este caso al hacer el contraste con la data desconocida, el MAPE es de 3.71 por lo que nuestro modelo tiene un buen poder predictivo.

IPC General SLV 2009-2022 mayo

Cargando Data

library(readxl)
library(readxl)
IPC_general_SLV_2009_2022_mayo_ <- read_excel("IPC general SLV 2009-2022 (mayo).xltx",range = "B7:B167",col_names = FALSE)

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

Metodología paso a paso

library(TSstudio)
library(forecast)
ts_plot(Yt2, Xtitle = "Años/Meses")

Identificación

Orden de Integración

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

Ordenes de Integracion
	x
d	1
D	0

# El valor de d es 1 y el de D es 0

## Graficando la serie diferenciada
Yt2 %>%
  diff(lag=12, diffences= D) %>%
  diff(diffences=d) %>%
  ts_plot(title= "Yt estacionaria")

## Verificando los valores para (p,q) y (P,Q)
Yt2 %>%
  diff(lag=12, diffences = D) %>%
  diff(diffences=d) %>%
  ts_cor(lag.max = 36)

En la parte del ACF debido a que uno de los ordenes está fuera de las bandas (es significativo), el valor del componente autoregresivo de la parte no ordinaria es igual a 1

PACF es para identificar valor de q y debido a que uno de los coeficientes de autocorrelacion parcial se encuentra fuera de las bandas se puede garantizar que hay un componente de media movil ordinaria (q), el valor es 0

El componente autoregresivo estacional (P) es 1

En el caso de media movil (Q) es 2

Estimación del modelo usando forecast

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

## Series: . 
## ARIMA(1,1,0)(1,0,2)[12] 
## 
## Coefficients:
##          ar1     sar1    sma1     sma2
##       0.3150  -0.0171  0.2325  -0.0576
## s.e.  0.0766   0.8200  0.8189   0.1958
## 
## sigma^2 = 0.2269:  log likelihood = -106.73
## AIC=223.47   AICc=223.86   BIC=238.84
## 
## Training set error measures:
##                      ME      RMSE      MAE        MPE      MAPE      MASE
## Training set 0.09418575 0.4688496 0.307792 0.08409974 0.2815712 0.1667916
##                     ACF1
## Training set -0.06000596

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

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

Comparativa

library(forecast)

#Estimar el modelo
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.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

#Generando el pronostico

pronosticosMW2<-forecast(object = modeloHW2,h=12, level=c(0.95))
pronosticosMW2

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

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

Con el grafico comparativo vemos lo bien que nuestra serie original se reproduce ya que logra replicar bien la serie Yt_Sarima2 y el HW

Validación Cruzada (Cross Validated)

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

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

TSCV2<-data.cross.validation %>% 
model(ARIMA(IPC ~ pdq(1, 1, 0) + PDQ(1, 0, 2))) %>% 
forecast(h=1) %>% 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(IPC ~ pdq(1, 1~ Test  0.0538 0.439 0.292 0.0434 0.258 0.158 0.161 0.0951

al hacer el contraste con la data desconocida, el MAPE es de 0.257 por lo que nuestro modelo tiene un buen poder predictivo

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

Cargando Data

library(readxl)
library(readxl)
Importacion_de_Hidrocarburos1 <- read_excel("Importacion de Hidrocarburos1.xlsx", 
    range = "B7:B71", col_names = FALSE)

## New names:
## * `` -> `...1`

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

Metodología paso a paso

library(TSstudio)
library(forecast)
ts_plot(Yt3, Xtitle = "Años/Meses")