Descomposición de Series Temporales (Enfoque Tradicional)

library(readxl)
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
serie.imae <-  read_excel("C:/Users/laura/Downloads/El Salvador IMAE por actividades.xlsx",
             col_types = c("skip", "numeric"),
             skip = 6)

serie.imae.ts <- ts(data = serie.imae,
                    start = c(2018, 1),
                    frequency = 12)
serie.imae.ts %>% 
  autoplot(main = "IMAE, El Salvador 2018-2026[marzo]",
                           xlab = "Años/Meses",
                           ylab = "Indice")

Modelo Aditivo

Componente de Tendencia Tt [Componente TCt]

ma2_12 <- ma(serie.imae.ts, 12, centre = TRUE)
autoplot(serie.imae.ts,main = "IMAE, El Salvador 2018-2026[marzo]",
           xlab = "Años/Meses",
           ylab = "Indice")+
  autolayer(ma2_12,series = "Tt")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

Cálculo de los Factores Estacionales [Componente St]

library(magrittr)
Yt <- serie.imae.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(2018, 1), frequency = 12) 
autoplot(St,
         main = "Factores Estacionales",
         xlab = "Años/Meses",
         ylab = "Factor Estacional") 

Cálculo del Componente Irregular It.

descomposicion <- decompose(serie.imae.ts, type = "additive")

It <- descomposicion$random

autoplot(It,
         main = "Componente Irregular 2018-2026",
         xlab = "Años/Meses",
         ylab = "It")

Descomposición Aditiva (usando la libreria stats):

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

Descomposición Aditiva usando libreria feasts

library(tsibble)
## Warning: package 'tsibble' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(feasts)
## Warning: package 'feasts' was built under R version 4.5.3
## Cargando paquete requerido: fabletools
## Warning: package 'fabletools' was built under R version 4.5.3
library(ggplot2)
Yt %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, IMAE")+xlab("Años/Meses")
## Warning: `autoplot.dcmp_ts()` was deprecated in fabletools 0.6.0.
## ℹ Please use `ggtime::autoplot.dcmp_ts()` instead.
## ℹ Graphics functions have been moved to the {ggtime} package. Please use
##   `library(ggtime)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).

Modelo Multiplicativo.

Componente Tendencia Ciclo [Tt=TCt]

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

Cálculo de Factores Estacionales [St]

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(2018, 1), frequency = 12) 
autoplot(St,
         main = "Factores Estacionales",
         xlab = "Años/Meses",
         ylab = "Factor Estacional") 

Cálculo del Componente Irregular [It]

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

Descomposición Multiplicativa (usando la libreria stats):

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

Descomposición Multiplicativa usando libreria feasts

library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
  model(classical_decomposition(value, type = "multiplicative")) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Multiplicativa, IMAE") + xlab("Años/Meses")
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).

Descomposición usando la libreria TSstudio

library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.5.3
ts_decompose(Yt, type = "additive", showline = TRUE)
ts_seasonal(Yt,type = "box",title = "Análisis de Valores Estacionales")

Pronóstico de series temporales, enfoque deterministico (clásico)

Pronóstico Modelo de Holt Winters

Usando Stats y forecast

library(forecast)

#Estimar el modelo
ModeloHW<-HoltWinters(x = Yt,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
## Warning in HoltWinters(x = Yt, seasonal = "multiplicative", optim.start =
## c(0.9, : dificultades de optimización: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
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.9633383
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   132.2274689
## b     0.2152098
## s1    0.9770831
## s2    1.0264055
## s3    1.0200485
## s4    0.9932670
## s5    1.0166883
## s6    0.9794076
## s7    0.9611332
## s8    1.0135212
## s9    1.0809612
## s10   0.9696077
## s11   0.9553065
## s12   1.0058424
#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
##          Point Forecast     Lo 95    Hi 95
## Apr 2026       129.4075 122.91843 135.8966
## May 2026       136.1608 126.92877 145.3928
## Jun 2026       135.5370 124.43331 146.6407
## Jul 2026       132.1922 119.69673 144.6877
## Aug 2026       135.5281 121.29717 149.7591
## Sep 2026       130.7693 115.69492 145.8436
## Oct 2026       128.5361 112.47300 144.5992
## Nov 2026       135.7603 117.71435 153.8062
## Dec 2026       145.0265 124.80033 165.2526
## Jan 2027       130.2955 111.09070 149.5002
## Feb 2027       128.5793 108.64966 148.5089
## Mar 2027       135.5976  98.76312 172.4321
#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()

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 2026       128.4312 119.6066 137.2558
## May 2026       136.2651 125.1483 147.3819
## Jun 2026       134.3411 121.8859 146.7963
## Jul 2026       132.4598 118.8600 146.0596
## Aug 2026       134.1745 119.1778 149.1712
## Sep 2026       130.7962 115.0737 146.5186
## Oct 2026       130.3138 113.6196 147.0081
## Nov 2026       136.5318 118.0220 155.0416
## Dec 2026       145.2733 124.5480 165.9986
## Jan 2027       130.3703 110.8878 149.8528
## Feb 2027       128.9350 108.8290 149.0409
## Mar 2027       135.4347 113.4679 157.4015
#Gráfico de la serie original y del pronóstico.
PronosticosHW2 %>% autoplot()

Pronóstico de series temporales, enfoque moderno (estocástico)

library(TSstudio)
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.5.3
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)

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

Componentes Autorregresivos

Para elegir “p”, la parte autorregresiva regular, verificamos las barras Non-Seasonal que sobrepasan los intervalos de confianza en el gráfico PACF al inicio del mismo. Se observa que los lags no estacionales no sobrepasan claramente el intervalo de confianza en los primeros retardos, y las barras significativas se encuentran en órdenes superiores (lag 3), lo cual no constituye un patrón AR claro al inicio.

Se concluye entonces que el valor de “p” es 0

Para la elección de “P”, la parte autorregresiva estacional, verificamos las barras Seasonal Lag que sobrepasan los intervalos de confianza en el gráfico PACF. Se evidencia que la primera barra estacional (retardo 12) es altamente significativa, mientras que los retardos 24 y 36 no lo son.

Se concluye que el valor de “P” es 1

Componentes de Media Móvil

Para elegir “q”, la parte de media móvil regular, verificamos las barras Non-Seasonal que sobrepasan los intervalos de confianza en el gráfico ACF al inicio del mismo. Se observa que los lags no estacionales no sobrepasan el intervalo de confianza en los primeros retardos de forma consistente, y las barras apenas significativas se encuentran en órdenes superiores.

Se concluye entonces que el valor de “q” es 0

Para la elección de “Q”, la parte de media móvil estacional, verificamos las barras Seasonal Lag que sobrepasan los intervalos de confianza en el gráfico ACF. Se evidencia que la primera barra estacional (retardo 12) es altamente significativa, mientras que los retardos 24 y 36 no lo son.

Se concluye que el valor de “Q” es 1

Con esto se estimará el modelo SARIMA:

                                  SARIMA(0,1,0)(1,1,1)[12]

Estimación del modelo.

Usando forecast

library(forecast)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.5.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.1817  -0.9979
## s.e.  0.1186   0.6130
## 
## sigma^2 = 10.27:  log likelihood = -231.8
## AIC=469.6   AICc=469.89   BIC=476.96
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.06265305 2.952445 2.011228 0.01413767 1.754354 0.3573143
##                    ACF1
## Training set 0.01975159
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
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)
## Warning: package 'fable' was built under R version 4.5.3
library(fabletools)
library(tidyr)
## 
## Adjuntando el paquete: 'tidyr'
## The following object is masked from 'package:magrittr':
## 
##     extract
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Adjuntando el paquete: '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]>
## # ℹ 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_010_011    Orders   11.9   -232.  469.  469.  473. <cpl [0]>  <cpl [12]>
## 2 arima_original   Orders   10.3   -232.  470.  470.  477. <cpl [12]> <cpl [12]>
## 3 arima_automatico Orders   10.9   -231.  470.  471.  480. <cpl [1]>  <cpl [12]>
## 4 arima_010_110    Orders   14.4   -237.  479.  479.  484. <cpl [12]> <cpl [0]>

Validación Cruzada (Cross Validated)

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

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

TSCV<-data.cross.validation %>% 
model(ARIMA(IMAE ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt)
## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs
## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs
## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
## 1 observation is missing at 2026 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(IMAE ~ pdq(0, 1,… Test  0.164  3.40  2.75 0.101  2.21 0.489 0.437 -0.398