Ejercicio de la clase

Determinación de Yt

datos_ivae<-pivot_longer(data = IVAE_03_22[1,],names_to = "vars", cols = 2:160,values_to = "indice") %>% select("indice")
data_ivae_ts<-datos_ivae %>% ts(start = c(2009,1),frequency = 12)

descomposicion_aditiva<-decompose(data_ivae_ts,type = "additive")
Yt<-descomposicion_aditiva$x
autoplot(Yt,main="Descomposición Aditiva del IVAE",xlab="Años/Meses",ylab = "indice")

Componente de tendencia Tt

ma2_12 <- ma(data_ivae_ts, 12, centre = TRUE)
autoplot(data_ivae_ts,main = "IVAE, El Salvador 2009-2021[marzo]",
           xlab = "Años/Meses",
           ylab = "Indice")+
  autolayer(ma2_12,series = "Tt")
## Warning: Removed 12 row(s) containing missing values (geom_path).

Factores estacionales

library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
## 
##     extract
Yt <- data_ivae_ts 
Tt <- ma2_12 
SI <- Yt - Tt 

St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) 
St <- St - sum(St) / 12 

St <-
  rep(St, len = length(Yt)) %>% ts(start = c(2009, 1), frequency = 12) 
autoplot(St,
         main = "Factores Estacionales",
         xlab = "Años/Meses",
         ylab = "Factor Estacional")

Cálculo del componente irregular

It<-Yt-Tt-St
autoplot(It,
         main = "Componente Irregular",
         xlab = "Aos/Meses",
         ylab = "It")

Metodología Box-Jenkins

library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.2.1
library(forecast)

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

Identificación

Orden de integración

library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.2.1
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
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 = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "Yt estacionaria")

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

##Se construyen las funciones de autocorrelación simple y parcial.

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

Estimación del modelo Sarima

Usando Forecast

library(forecast)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.2.1
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.1022  -0.7232
## s.e.   0.1266   0.1135
## 
## sigma^2 = 6.779:  log likelihood = -351.22
## AIC=708.43   AICc=708.6   BIC=717.38
## 
## Training set error measures:
##                      ME     RMSE     MAE        MPE     MAPE      MASE
## Training set 0.05202814 2.477766 1.68543 0.02782791 1.663587 0.4366552
##                    ACF1
## Training set 0.06837002
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
#Determinación del modelo HW 
modelo_HoWin<-HoltWinters(x=Yt,
                          seasonal = "multiplicative",
                          optim.start = c(0.9,0.9,0.9))

modelo_HoWin
## 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.8472467
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   118.3719784
## b     0.1600947
## s1    0.9529095
## s2    1.0192349
## s3    1.0440936
## s4    0.9962326
## s5    1.0009929
## s6    0.9838373
## s7    0.9554392
## s8    1.0145286
## s9    1.0872315
## s10   0.9748891
## s11   0.9626701
## s12   1.0059813
#Generar el pronóstico:
PronosticosHW<-forecast(object = modelo_HoWin,h = 12,level = c(0.95))
PronosticosHW
##          Point Forecast     Lo 95    Hi 95
## Apr 2022       112.9503 107.46492 118.4358
## May 2022       120.9752 113.57246 128.3779
## Jun 2022       124.0929 115.22236 132.9634
## Jul 2022       118.5640 108.86876 128.2592
## Aug 2022       119.2908 108.50116 130.0804
## Sep 2022       117.4038 105.81297 128.9947
## Oct 2022       114.1680 101.97016 126.3657
## Nov 2022       121.3911 107.66977 135.1125
## Dec 2022       130.2643 114.88360 145.6450
## Jan 2023       116.9603 102.34981 131.5708
## Feb 2023       115.6485 100.48405 130.8129
## Mar 2023       121.0126  83.19827 158.8270
Yt_Sarima<-modelo_estimado$fitted
Yt_modelo_HoWin<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_modelo_HoWin)
ts_plot(grafico_comparativo)

Verificación de sobre ajuste/sub ajuste

library(tsibble)
## Warning: package 'tsibble' was built under R version 4.2.1
## 
## 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.2.1
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 4.2.1
## 
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
## 
##     accuracy, forecast
library(fable)
## Warning: package 'fable' was built under R version 4.2.1
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)[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 × 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.  702.  702.  714. <cpl [1]>  <cpl [12]>
## 2 arima_010_011    Orders   6.72   -352.  707.  707.  713. <cpl [0]>  <cpl [12]>
## 3 arima_original   Orders   6.78   -351.  708.  709.  717. <cpl [12]> <cpl [12]>
## 4 arima_010_110    Orders   7.99   -361.  726.  726.  731. <cpl [12]> <cpl [0]>

Validación cruzada

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 × 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.0506  2.94  2.09 0.00243  2.00 0.542 0.494 0.155