Ejercicio de Pronóstico Modelos SARIMA

1. Presentacion de la serie temporal

    #carga del archivo fuente
       library(dplyr)
    library(tidyr)
    library(forecast)
  library(TSstudio)
  library(readxl)
PIB_trimestral<- read_excel("C:/Users/Familia/Downloads/PIB trimestral precios corrientes SLV 2009-2021.xlsx",col_names = FALSE,col_types = c("skip", "numeric"), skip = 6)


PIB_trimestral.ts <- ts(data = PIB_trimestral,
                    start = c(2009, 1),
                    frequency = 4)
PIB_trimestral.ts %>% 
  autoplot(main = "PIB trimestral precios corrientes SLV 2009-2021",
                           xlab = "Años/Trimestres",
                           ylab = "Indice")

library(forecast)
  library(TSstudio)
   Yt <-PIB_trimestral.ts
ts_plot(Yt,Xtitle = "Años/Trimestres")

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

PronosticosHW<-forecast(object = ModeloHW,h = 4,level = c(0.95))
PronosticosHW

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

2.Orden de integracion de la serie

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 Integración") %>% kable_material()

Ordenes de Integración

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

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

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

4.Estimación del modelo (Usando forecast)

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

summary(modelo_estimado)

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

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

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

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

5.Verificación de sobre ajuste/sub ajuste

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

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

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

TSCV<-data.cross.validation %>% 
model(ARIMA(PIB_trimestral ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
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_trimestral ~~ Test   43.4  394.  215. 0.422  3.48 0.654 0.830 0.0397

# en la verificacion de sobre ajuste el MAPE presento una variacion superior a un punto porcentual (3.48) sobre la estimacion del modelo (2.05). no obstante no determina que las predicciones de nuestro modelo dejen de tener un buen poder predictivo.