Descomposición de Series Temporales (Enfoque Tradicional)

library(readxl)
library(forecast)

# Leemos el archivo completo omitiendo las 5 filas de títulos institucionales
serie.completa <- read_excel("C:/Users/Rebeca/Downloads/El Salvador IMAE por actividades.xlsx", skip = 5)
## New names:
## • `` -> `...1`
# Seleccionamos únicamente la segunda columna (donde están los valores del IMAE)
# Esto evita tener que adivinar los 'col_types' manualmente
serie.imae <- serie.completa[, 2]

# Creamos el objeto de serie de tiempo desde enero de 2018
serie.ivae.ts <- ts(data = serie.imae, 
                    start = c(2018, 1), 
                    frequency = 12)

# Graficamos la serie original
serie.ivae.ts %>% 
  autoplot(main = "IMAE, El Salvador 2018-2026[marzo]", 
           xlab = "Años/Meses", 
           ylab = "Indice")

Modelo Aditivo

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)
## Warning: package 'ggplot2' was built under R version 4.5.3
Yt <- serie.ivae.ts 

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 every 8 hours.
## 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()`).

### Componente de Tendencia Tt

library(forecast)
ma2_12 <- ma(serie.ivae.ts, 12, centre = TRUE)
autoplot(serie.ivae.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

library(magrittr)
## Warning: package 'magrittr' was built under R version 4.5.3
library(forecast) # Asegurar que autoplot funcione correctamente

Yt <- serie.ivae.ts # Serie original (2018 - 2026[marzo])
Tt <- ma2_12        # Media móvil centrada (2x12-MA) como tendencia

# 1. Obtener la diferencia (Modelo Aditivo)
SI <- Yt - Tt 

# 2. Promediar los resultados de cada mes
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) 

# 3. Ajustar para que los factores sumen exactamente "0"
St <- St - sum(St) / 12 

# 4. CORRECCIÓN: Replicar los factores usando el inicio correcto (2018)
St <- rep(St, len = length(Yt)) %>% 
  ts(start = start(Yt), frequency = 12) # Copia automáticamente el inicio de tu serie original (2018, 1)

# 5. Graficar los factores estacionales ajustados
autoplot(St,
         main = "Factores Estacionales - IMAE El Salvador",
         xlab = "Años/Meses",
         ylab = "Factor Estacional")

Cálculo del Componente Irregular It

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

### Descomposición Aditiva (usando la libreria stats)

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

Modelo Multiplicativo

Componente Tendencia Ciclo

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

Cálculo de Factores Estacionales

library(magrittr)
library(forecast)

# 1. Obtener la serie sin tendencia (Modelo Multiplicativo)
SI <- Yt / Tt 

# 2. Promediar los resultados de cada mes
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) 

# 3. Ajustar para que los factores promedien exactamente "1"
St <- St * 12 / sum(St) 

# 4. CORRECCIÓN: Replicar los factores usando el inicio real de tu serie (2018)
St <- rep(St, len = length(Yt)) %>% 
  ts(start = start(Yt), frequency = 12) # Detecta automáticamente el inicio correcto

# 5. Graficar los factores estacionales multiplicativos
autoplot(St,
         main = "Factores Estacionales (Modelo Multiplicativo) - IMAE",
         xlab = "Años/Meses",
         ylab = "Factor Estacional (Porcentaje)")

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)

# Realizar la descomposición multiplicativa automática
descomposicion_multiplicativa <- decompose(serie.ivae.ts, type = "multiplicative")

# Graficar los 4 componentes multiplicativos
autoplot(descomposicion_multiplicativa, 
         main = "Descomposición Multiplicativa - IMAE El Salvador",
         xlab = "Años/Meses")

Descomposición Multiplicativa usando libreria feasts

library(tsibble)
library(feasts)
library(ggplot2)
library(magrittr) # Por si acaso para asegurar el operador %>%

# Descomposición clásica multiplicativa usando el enfoque moderno (Tidyverts)
Yt %>% 
  as_tsibble() %>%  
  model(classical_decomposition(value, type = "multiplicative")) %>%  
  components() %>%  
  autoplot() +  
  labs(title = "Descomposición Clásica Multiplicativa - IVAE El Salvador",
       xlab = "Años/Meses")
## Ignoring unknown labels:
## • 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)
library(TSstudio)
ts_seasonal(Yt, 
            type = "box", 
            title = "Análisis de Valores Estacionales - IVAE El Salvador")

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

Pronóstico Modelo de Holt Winters

Usando Stats y forecast

library(forecast)

# 1. Estimar el modelo usando tus objetos definidos (Yt)
ModeloHW <- HoltWinters(x = Yt,                      
                        seasonal = "multiplicative",                      
                        optim.start = c(0.9, 0.9, 0.9))

# 2. Desplegar el objeto en consola
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.1133145
##  beta : 0
##  gamma: 0.5581176
## 
## Coefficients:
##           [,1]
## a   88.6094932
## b    0.4723193
## s1   0.7400912
## s2   1.0366972
## s3   0.9419354
## s4   0.8557054
## s5   0.8889486
## s6   0.7175855
## s7   0.7582607
## s8   0.9668398
## s9   1.0821564
## s10  0.7806829
## s11  0.6526374
## s12  0.8292333
#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
##          Point Forecast    Lo 95     Hi 95
## May 2026       65.92866 45.46336  86.39397
## Jun 2026       92.84052 71.68972 113.99131
## Jul 2026       84.79910 63.48287 106.11534
## Aug 2026       77.44029 55.98874  98.89184
## Sep 2026       80.86862 58.97950 102.75775
## Oct 2026       65.61846 43.99444  87.24249
## Nov 2026       69.69608 47.57104  91.82113
## Dec 2026       89.32445 65.71526 112.93363
## Jan 2027      100.48944 75.77846 125.20042
## Feb 2027       72.86323 49.92380  95.80266
## Mar 2027       61.22066 38.76982  83.67149
## Apr 2027       78.17790 40.07261 116.28318
#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
## May 2026       70.06399 -2.680502 142.8085
## Jun 2026      103.54552 -4.487947 211.5790
## Jul 2026       97.33869 -4.705973 199.3834
## Aug 2026       85.99014 -4.580852 176.5611
## Sep 2026       75.62609 -4.395537 155.6477
## Oct 2026       64.52669 -4.058697 133.1121
## Nov 2026       61.29085 -4.143693 126.7254
## Dec 2026       90.63573 -6.548164 187.8196
## Jan 2027      119.29628 -9.164570 247.7571
## Feb 2027       64.95123 -5.282738 135.1852
## Mar 2027       55.96343 -4.800833 116.7277
## Apr 2027       77.59913 -6.997720 162.1960
#Gráfico de la serie original y del pronóstico.
PronosticosHW2 %>% autoplot()

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

Metodología Box-Jenkins

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

Estimación del Modelo

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.0386  -0.8560
## s.e.  0.1457   0.2286
## 
## sigma^2 = 293:  log likelihood = -376.57
## AIC=759.14   AICc=759.42   BIC=766.53
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 1.138996 15.78034 9.310496 -40.59759 73.79015 0.6028872 -0.4090379
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

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)
## 
## 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   293.   -377.  757.  757.  762. <cpl [0]>  <cpl [12]>
## 2 arima_original   Orders   293.   -377.  759.  759.  767. <cpl [12]> <cpl [12]>
## 3 arima_010_110    Orders   401.   -385.  774.  774.  779. <cpl [12]> <cpl [0]> 
## 4 arima_automatico Orders   223.   -409.  833.  834.  851. <cpl [24]> <cpl [4]>

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 2026 may
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, 1,… Test  0.929  12.6  10.1 0.538  18.1 0.651 0.586 -0.569