2.1. Modelo Aditivo.

Modelo Aditivo: La serie temporal es el resultado de la suma de los componentes teóricos. \[\color{teal}{Y_t=T_t+C_t+S_t+I_t}\]

ma2_12 <- ma(serie.ivae.ts, 12, centre = TRUE)
autoplot(serie.ivae.ts,main = "IVAE, El Salvador 2009-2021[marzo]",
           xlab = "Años/Meses",
           ylab = "Indice")+
  autolayer(ma2_12,series = "Tt")

library(magrittr)
Yt <- serie.ivae.ts #Serie original
Tt <- ma2_12 #Media móvil centrada (2x12-MA) como componente de Tendencia Ciclo
SI <- Yt - Tt #Diferencia que contiene componentes Estacional e Irregular

St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) #Promediando los resultados de cada mes
#Los factores estacionales deben sumar "0" en el modelo aditivo
St <- St - sum(St) / 12 
#Generar la serie de factores para cada valor de la serie original
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") 

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

descomposicion_aditiva<-decompose(serie.ivae.ts,type = "additive")
autoplot(descomposicion_aditiva,main="Descomposición Aditiva",xlab="Años/Meses")

library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, IVAE")+xlab("Años/Meses")

2.2. Modelo Multiplicativo.

Modelo Multiplicativo: La serie temporal es el resultado de la amplificación/atenuación de la tendencia a causa del resto de los componentes.

\[\color{teal}{Y_t=T_t· C_t·S_t·I_t}\]

Un modelo Multiplicativo puede expresarse como un un modelo aditivo, a través de una transformación logaritmica, (las minusculas indican el logaritmo natural de la variable/componente): \[\color{teal}{y_t=t_t+c_t+s_t+i_t}\]

Tt<- ma(serie.ivae.ts, 12, centre = TRUE)
autoplot(Tt,main = "Componente Tendencia [Ciclo]", xlab = "Años/Meses",ylab = "Tt")

SI<-Yt/Tt #Serie sin tendencia.
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) #Promediando los resultados de cada mes
#Los factores estacionales deben promediar "1" en el modelo multiplicativo
St <- St*12/sum(St) 
#Generar la serie de factores para cada valor de la serie original
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") 

It<-Yt/(Tt*St)
autoplot(It,
         main = "Componente Irregular",
         xlab = "Años/Meses",
         ylab = "It")

descomposicion_multiplicatica<-decompose(serie.ivae.ts,type = "multiplicative")
autoplot(descomposicion_multiplicatica,main="Descomposición Multiplicativa",xlab="Años/Meses")

library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
  model(classical_decomposition(value, type = "multiplicative")) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Multiplicativa, IVAE") + xlab("Años/Meses")

2.3. Descomposición usando la libreria TSstudio

library(TSstudio)
ts_decompose(Yt, type = "additive", showline = TRUE)

ts_seasonal(Yt,type = "box",title = "Análisis de Valores Estacionales")

3.1. Pronóstico Modelo de Holt Winters

Si la serie resulta tener un patrón estacional, está técnica permite realizar predicciones extrapolando los componentes de la serie temporal observada.

3.1.1 Usando Stats y forecast

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.8408163
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   117.0799442
## b     0.1600306
## s1    0.9502255
## s2    1.0233274
## s3    1.0518154
## s4    0.9900276
## s5    1.0007537
## s6    0.9807088
## s7    0.9600206
## s8    1.0149628
## s9    1.0915423
## s10   0.9796752
## s11   0.9584676
## s12   0.9994880

#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW

##          Point Forecast     Lo 95    Hi 95
## Apr 2021       111.4044 105.94952 116.8593
## May 2021       120.1386 112.77972 127.4976
## Jun 2021       123.6515 114.83356 132.4694
## Jul 2021       116.5461 107.01096 126.0813
## Aug 2021       117.9689 107.30373 128.6341
## Sep 2021       115.7630 104.33417 127.1918
## Oct 2021       113.4746 101.36814 125.5810
## Nov 2021       120.1312 106.57272 133.6897
## Dec 2021       129.3698 114.12877 144.6109
## Jan 2022       116.2681 101.78191 130.7543
## Feb 2022       113.9046  98.99575 128.8134
## Mar 2022       118.9394  81.61739 156.2614

#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()

3.1.12 Usando Forecast (Aproximación por Espacios de los Estados ETS )

library(forecast)

#Generar el pronóstico:
PronosticosHW2<-hw(y = Yt,
                   h = 12,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
PronosticosHW2

##          Point Forecast     Lo 95    Hi 95
## Apr 2021       107.4928  99.59377 115.3918
## May 2021       115.3370 106.85846 123.8155
## Jun 2021       115.6946 107.18632 124.2029
## Jul 2021       109.6301 101.56399 117.6962
## Aug 2021       112.2047 103.94477 120.4646
## Sep 2021       110.2284 102.10923 118.3476
## Oct 2021       107.9393  99.98348 115.8951
## Nov 2021       114.4306 105.99022 122.8710
## Dec 2021       122.6459 113.59234 131.6994
## Jan 2022       108.5830 100.56060 116.6054
## Feb 2022       107.1239  99.20182 115.0460
## Mar 2022       113.1678 104.79016 121.5455

#Gráfico de la serie original y del pronóstico.
PronosticosHW2 %>% autoplot()

Metodología Box-Jenkins

library(TSstudio)
library(forecast)

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

Una serie se considera estacionaria, en sentido aproximado, si su media es constante a lo largo del tiempo, y su varianza es finita.

Normalmente una serie no es estacionaria, y debe ser sometida a transformaciones para lograr que lo sea. Dichas transformaciones pretenden estabilizar la media y la varianza de la serie, de tal manera que sea relativamente fácil identificar los componentes del proceso generador de los datos. > Las transformaciones para estabilizar la varianza están representadas por la familia de transformaciones de Box & Cox, la más popular es la transformación de logaritmo natural.

Las trasformaciones para estabilizar la media se consiguen a través de lo que se denomina diferencias regulares \(\Delta^d{Y_t}\) y de las diferencias estacionales \(\Delta_m^D{Y_t}\), puede ser necesario aplicar ambas d>0,D>0, o ninguna dependiendo lo que indiquen las pruebas de raices unitarias, apropiadas.

La serie se observa con poca volatilidad, a excepción del periodo de pandemia, lo cual indicaría que habrá que intervenir el modelo, para captar el efecto de la cuarentena sobre la actividad económica. Por lo que no habrá necesidad de aplicar ninguna transformación para estabilizar la varianza.

Identificación

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

#Gráfico de la serie diferenciada
Yt %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "Yt estacionaria")

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

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)[12] 
## 
## Coefficients:
##          sar1     sma1
##       -0.2040  -0.6422
## s.e.   0.1586   0.1554
## 
## sigma^2 = 6.897:  log likelihood = -323.43
## AIC=652.85   AICc=653.04   BIC=661.55
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.04712358 2.488576 1.666188 0.02203353 1.658867 0.4981404
##                    ACF1
## Training set 0.08970896

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

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

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

Se estimarán los modelos: Partiendo del modelo original \(SARIMA(0,1,0)(1,1,1)_{[12]}\), se estima un nuevo modelo con P-1 \(SARIMA(0,1,0)(0,1,1)_{[12]}\), y otro con Q-1: \(SARIMA(0,1,0)(1,1,0)_{[12]}\)

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)[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.30   -320.  649.  649.  660. <cpl [1]>  <cpl [12]>
## 2 arima_010_011    Orders   6.88   -324.  652.  653.  658. <cpl [0]>  <cpl [12]>
## 3 arima_original   Orders   6.90   -323.  653.  653.  662. <cpl [12]> <cpl [12]>
## 4 arima_010_110    Orders   7.60   -328.  660.  660.  666. <cpl [12]> <cpl [0]>