##Cargamos la base de datos

library(readxl)
## Warning: package 'readxl' was built under R version 4.1.3
library(forecast)
## Warning: package 'forecast' was built under R version 4.1.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
serie.ivae <- read_excel("C:/Users/Administrador/Desktop/ECO/Pronostico Modelo Sarima.xlsx",
      col_types = c("skip", "numeric"),
              skip = 0)

#Serie                      
serie.ivae.ts<- ts(data = serie.ivae, 
                   start = c(2009,1),
                   frequency = 12)



Yt<-serie.ivae.ts
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.1.3
library(forecast)

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

##Orden de Integración

library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.1.3
library(magrittr)
## Warning: package 'magrittr' was built under R version 4.1.3
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 = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "Yt estacionaria")

##Verificamos los valores para (p,q) & (P,Q)

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

##Estimación del modelo. Usando forecast

library(forecast)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.1.3
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)[12] 
## 
## Coefficients:
##          sar1     sma1
##       -0.0855  -0.7231
## s.e.   0.1307   0.1212
## 
## sigma^2 = 6.809:  log likelihood = -351.38
## AIC=708.76   AICc=708.93   BIC=717.71
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE     MASE
## Training set 0.04757629 2.483195 1.689155 0.02273075 1.663519 0.440321
##                    ACF1
## Training set 0.06586925
modelo_estimado %>% autoplot(type="both")+theme_solarized()

modelo_estimado %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations

##Estimación del modelo y generar el pronóstico

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.8540065
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   118.3463237
## b     0.1600306
## s1    0.9501653
## s2    1.0195813
## s3    1.0450993
## s4    0.9926690
## s5    1.0000193
## s6    0.9831838
## s7    0.9576173
## s8    1.0161329
## s9    1.0897384
## s10   0.9766794
## s11   0.9611902
## s12   1.0061994
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
##          Point Forecast     Lo 95    Hi 95
## Apr 2022       112.6006 107.11066 118.0906
## May 2022       120.9900 113.54350 128.4365
## Jun 2022       124.1854 115.25067 133.1201
## Jul 2022       118.1142 108.37761 127.8507
## Aug 2022       119.1488 108.28278 130.0148
## Sep 2022       117.3002 105.62196 128.9785
## Oct 2022       114.4032 102.08324 126.7232
## Nov 2022       121.5565 107.70542 135.4076
## Dec 2022       130.5361 114.99867 146.0734
## Jan 2023       117.1494 102.40141 131.8974
## Feb 2023       115.4453 100.18435 130.7063
## Mar 2023       121.0123  83.09391 158.9306
Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)

##Verificación de sobre ajuste/sub ajuste

library(tsibble)
## Warning: package 'tsibble' was built under R version 4.1.3
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(feasts)
## Warning: package 'feasts' was built under R version 4.1.3
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 4.1.3
## 
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
## 
##     accuracy, forecast
library(fable)
## Warning: package 'fable' was built under R version 4.1.3
library(fabletools)
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.3
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
## 
##     extract
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.3
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
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)[12]> <ARIMA(0,1,0)(0,1,1)[12]> <ARIMA(0,1,0)(1,1,0)[12]>
## # ... 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_automatico Orders   6.13   -347.  701.  702.  713. <cpl [1]>  <cpl [12]>
## 2 arima_010_011    Orders   6.74   -352.  707.  707.  713. <cpl [0]>  <cpl [12]>
## 3 arima_original   Orders   6.81   -351.  709.  709.  718. <cpl [12]> <cpl [12]>
## 4 arima_010_110    Orders   7.88   -360.  723.  723.  729. <cpl [12]> <cpl [0]>

##Validación Cruzada (Cross Validated)

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

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

TSCV<-data.cross.validation %>% 
model(ARIMA(IVAE ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt)
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
## 1 observation is missing at 2022 abr.
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(IVAE ~ pdq(0, ~ Test  0.0296  2.94  2.09 -0.0182  2.00 0.546 0.498 0.154