## Cargando Datos
library(readxl)
library(forecast)
PIB_trimestral_precios_corrientes_SLV_2009_2021 <- read_excel("C:/Users/50379/Downloads/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_2021, start = c(2009,1),
                   frequency = 4)
Yt <- serie.pib.ts

##METODOLOGIA PASO A PASO

library(TSstudio)
library(forecast)

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

#IDENTIFICACION

##ORDEN DE INTEGRACION

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

##ESTIMACION DEL MODELO ###USANDO LIBRERYA FORECATS

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

##VALIDACION CRUZADA

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

##CARGANDO DATA

library(readxl)
IPC_general_SLV_2009_2022_mayo_ <- read_excel("C:/Users/50379/Downloads/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

#METODOLOGIA PPASO A PASO

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

#IDENTIFICACION

##ORDEN DE INTEGRACION

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

##ESTIMACION DEL MODELO

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)