Importación y creación de objetos ts

library(readxl)
library(forecast)
library(TSstudio)
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(kableExtra)
library(magrittr)
library(dplyr)
library(tidyr)
library(ggplot2)


# --- PIB trimestral (2009 Q1 – 2025 Q4) ---
pib_raw <- read_excel("C:/Users/Paniagua/Downloads/PIB a precios corrientes  trimestral.xls",
                      sheet = "Datos", col_types = c("skip","numeric"), skip = 6)
names(pib_raw) <- "PIB"
PIB.ts <- ts(pib_raw$PIB, start = c(2009, 1), frequency = 4)

# --- IPC mensual (2009 enero – 2026 mayo) ---
ipc_raw <- read_excel("C:/Users/Paniagua/Downloads/Índice de precios al consumidor.xls",
                      sheet = "Datos", col_types = c("skip","numeric"), skip = 6)
names(ipc_raw) <- "IPC"
IPC.ts <- ts(ipc_raw$IPC, start = c(2009, 1), frequency = 12)

# --- Importación de Hidrocarburos mensual (2017 sep – 2026 abr) ---
hid_raw <- read_excel("C:/Users/Paniagua/Downloads/Importación de hidrocarburos.xls",
                      sheet = "Datos", col_types = c("skip","numeric"), skip = 6)
names(hid_raw) <- "HID"
HID.ts <- ts(hid_raw$HID, start = c(2017, 9), frequency = 12)

autoplot(PIB.ts, main = "PIB Trimestral, El Salvador 2009-2025", xlab = "Años/Trimestres", ylab = "Millones USD")

autoplot(IPC.ts, main = "IPC General, El Salvador 2009-2026", xlab = "Años/Meses", ylab = "Índice")

autoplot(HID.ts, main = "Importación de Hidrocarburos, El Salvador 2017-2026", xlab = "Años/Meses", ylab = "Millones USD")

# Serie 1: PIB Trimestral (2009–2025) {.tabset}

1A. SARIMA Manual

Órdenes de integración

d_pib <- ndiffs(PIB.ts)
D_pib <- nsdiffs(PIB.ts)
data.frame(d = d_pib, D = D_pib) %>%
  kable(caption = "Órdenes de Integración — PIB") %>% kable_material()
Órdenes de Integración — PIB
d D
1 1

Correlograma de la serie diferenciada

PIB.ts %>%
  diff(lag = 4, differences = D_pib) %>%
  diff(differences = d_pib) %>%
  ts_cor(lag.max = 4*4)

Estimación del modelo

pib.manual <- Arima(PIB.ts,
                    order    = c(0, 1, 0),
                    seasonal = c(1, 1, 1))
summary(pib.manual)
## Series: PIB.ts 
## ARIMA(0,1,0)(1,1,1)[4] 
## 
## Coefficients:
##         sar1     sma1
##       0.0879  -0.8758
## s.e.  0.1547   0.1413
## 
## sigma^2 = 64295:  log likelihood = -439.71
## AIC=885.41   AICc=885.82   BIC=891.84
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 15.24958 240.1594 127.0158 0.1125891 1.959043 0.3511589 -0.1762194

Diagnóstico de residuos

checkresiduals(pib.manual)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0)(1,1,1)[4]
## Q* = 4.365, df = 6, p-value = 0.6274
## 
## Model df: 2.   Total lags used: 8

Pronóstico 2026 (4 trimestres)

h_pib <- 4   # 4 trimestres = 2026 completo

pron.pib.manual <- forecast(pib.manual, h = h_pib, level = c(0.95))
pron.pib.manual %>% print()
##         Point Forecast    Lo 95     Hi 95
## 2026 Q1       9369.766 8871.988  9867.544
## 2026 Q2       9681.673 8977.842 10385.504
## 2026 Q3       9654.080 8792.121 10516.039
## 2026 Q4      10208.451 9213.178 11203.725
pron.pib.manual %>% autoplot() +
  labs(title = "Pronóstico PIB 2026 — SARIMA Manual", x = "Años/Trimestres", y = "Millones USD")

1B. SARIMA Automático

pib.auto <- auto.arima(PIB.ts)
summary(pib.auto)
## Series: PIB.ts 
## ARIMA(1,0,0)(0,1,1)[4] with drift 
## 
## Coefficients:
##          ar1     sma1    drift
##       0.8176  -0.7784  72.0990
## s.e.  0.0870   0.1507  10.9996
## 
## sigma^2 = 60944:  log likelihood = -443.83
## AIC=895.66   AICc=896.34   BIC=904.29
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -2.09545 233.8164 124.3587 -0.2416389 1.921945 0.3438129
##                     ACF1
## Training set -0.07426543
pron.pib.auto <- forecast(pib.auto, h = h_pib, level = c(0.95))
pron.pib.auto %>% print()
##         Point Forecast    Lo 95     Hi 95
## 2026 Q1       9230.008 8746.134  9713.882
## 2026 Q2       9468.117 8843.105 10093.129
## 2026 Q3       9368.868 8665.122 10072.615
## 2026 Q4       9905.815 9154.024 10657.606
pron.pib.auto %>% autoplot() +
  labs(title = "Pronóstico PIB 2026 — SARIMA Automático", x = "Años/Trimestres", y = "Millones USD")

1C. Holt-Winters

pib.hw <- HoltWinters(PIB.ts, seasonal = "multiplicative")
pron.pib.hw <- forecast(pib.hw, h = h_pib, level = c(0.95))
pron.pib.hw %>% print()
##         Point Forecast    Lo 95     Hi 95
## 2026 Q1       9256.911 8869.430  9644.392
## 2026 Q2       9606.445 9047.345 10165.545
## 2026 Q3       9409.430 8737.530 10081.330
## 2026 Q4      10074.125 9353.045 10795.204
pron.pib.hw %>% autoplot() +
  labs(title = "Pronóstico PIB 2026 — Holt-Winters", x = "Años/Trimestres", y = "Millones USD")

1D. Validación Cruzada

PIB.tsibble <- PIB.ts %>% as_tsibble() %>% rename(PIB = value)

cv_pib <- PIB.tsibble %>%
  stretch_tsibble(.init = 20, .step = 4) %>%   # saltar de 4 en 4 en vez de 1 en 1
  model(
    SARIMA_manual = ARIMA(PIB ~ pdq(0,1,0) + PDQ(1,1,1)),
    SARIMA_auto   = ARIMA(PIB, ic = "bic", stepwise = TRUE),  # TRUE es más rápido
    ETS_HW        = ETS(PIB)
  ) %>%
  forecast(h = 1) %>%
  accuracy(PIB.tsibble)

cv_pib %>%
  select(.model, RMSE, MAE, MAPE) %>%
  arrange(MAPE) %>%
  kable(caption = "Validación Cruzada — PIB Trimestral", digits = 4) %>% kable_material()
Validación Cruzada — PIB Trimestral
.model RMSE MAE MAPE
ETS_HW 171.5043 120.7483 1.8396
SARIMA_manual 184.0973 162.9993 2.3816
SARIMA_auto 182.1587 158.4847 2.3896

Serie 2: IPC General (2009–2026) {.tabset}

2A. SARIMA Manual

Órdenes de integración

d_ipc <- ndiffs(IPC.ts)
D_ipc <- nsdiffs(IPC.ts)
data.frame(d = d_ipc, D = D_ipc) %>%
  kable(caption = "Órdenes de Integración — IPC") %>% kable_material()
Órdenes de Integración — IPC
d D
1 0

Correlograma de la serie diferenciada

ipc_diff <- IPC.ts

if (D_ipc > 0) ipc_diff <- diff(ipc_diff, lag = 12, differences = D_ipc)
if (d_ipc > 0) ipc_diff <- diff(ipc_diff, differences = d_ipc)

ts_cor(ipc_diff, lag.max = 3 * 12)

Estimación del modelo

# Horizonte: meses restantes de 2026 (junio–diciembre = 7 meses)
h_ipc <- 7

ipc.manual <- Arima(IPC.ts,
                    order    = c(0, 1, 0),
                    seasonal = c(1, 1, 1))
summary(ipc.manual)
## Series: IPC.ts 
## ARIMA(0,1,0)(1,1,1)[12] 
## 
## Coefficients:
##         sar1     sma1
##       0.1259  -1.0000
## s.e.  0.0763   0.0786
## 
## sigma^2 = 0.2068:  log likelihood = -138.05
## AIC=282.1   AICc=282.23   BIC=291.94
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.02526395 0.4381596 0.3005031 0.02005916 0.2633305 0.1384969
##                   ACF1
## Training set 0.2946617

Diagnóstico de residuos

checkresiduals(ipc.manual)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0)(1,1,1)[12]
## Q* = 67.489, df = 22, p-value = 0.000001633
## 
## Model df: 2.   Total lags used: 24

Pronóstico hasta diciembre 2026

pron.ipc.manual <- forecast(ipc.manual, h = h_ipc, level = c(0.95))
pron.ipc.manual %>% print()
##          Point Forecast    Lo 95    Hi 95
## Jun 2026       134.0066 133.0898 134.9234
## Jul 2026       134.1504 132.8539 135.4469
## Aug 2026       134.0289 132.4411 135.6168
## Sep 2026       134.0387 132.2052 135.8722
## Oct 2026       134.3570 132.3070 136.4069
## Nov 2026       134.5747 132.3290 136.8203
## Dec 2026       134.3311 131.9056 136.7566
pron.ipc.manual %>% autoplot() +
  labs(title = "Pronóstico IPC hasta diciembre 2026 — SARIMA Manual", x = "Años/Meses", y = "Índice")

2B. SARIMA Automático

ipc.auto <- auto.arima(IPC.ts)
summary(ipc.auto)
## Series: IPC.ts 
## ARIMA(0,1,1)(1,0,0)[12] with drift 
## 
## Coefficients:
##          ma1    sar1   drift
##       0.2663  0.1957  0.1639
## s.e.  0.0650  0.0719  0.0477
## 
## sigma^2 = 0.1996:  log likelihood = -126.26
## AIC=260.52   AICc=260.72   BIC=273.87
## 
## Training set error measures:
##                      ME      RMSE       MAE           MPE     MAPE      MASE
## Training set 0.00190199 0.4424161 0.2993996 -0.0005750672 0.264985 0.1379883
##                    ACF1
## Training set 0.01670461
pron.ipc.auto <- forecast(ipc.auto, h = h_ipc, level = c(0.95))
pron.ipc.auto %>% print()
##          Point Forecast    Lo 95    Hi 95
## Jun 2026       134.0254 133.1499 134.9010
## Jul 2026       134.2414 132.8287 135.6541
## Aug 2026       134.2793 132.4835 136.0751
## Sep 2026       134.4268 132.3164 136.5373
## Oct 2026       134.6311 132.2472 137.0150
## Nov 2026       134.7923 132.1632 137.4214
## Dec 2026       134.8224 131.9691 137.6757
pron.ipc.auto %>% autoplot() +
  labs(title = "Pronóstico IPC hasta diciembre 2026 — SARIMA Automático", x = "Años/Meses", y = "Índice")

2C. Holt-Winters

ipc.hw <- HoltWinters(IPC.ts, seasonal = "multiplicative")
pron.ipc.hw <- forecast(ipc.hw, h = h_ipc, level = c(0.95))
pron.ipc.hw %>% print()
##          Point Forecast    Lo 95    Hi 95
## Jun 2026       134.0374 132.9560 135.1189
## Jul 2026       134.0454 132.5886 135.5023
## Aug 2026       133.6077 131.8140 135.4014
## Sep 2026       133.5216 131.4034 135.6397
## Oct 2026       133.8208 131.3797 136.2619
## Nov 2026       134.0343 131.2734 136.7952
## Dec 2026       133.7209 130.6500 136.7918
pron.ipc.hw %>% autoplot() +
  labs(title = "Pronóstico IPC hasta diciembre 2026 — Holt-Winters", x = "Años/Meses", y = "Índice")

2D. Validación Cruzada

IPC.tsibble <- IPC.ts %>% as_tsibble() %>% rename(IPC = value)

cv_ipc <- IPC.tsibble %>%
  stretch_tsibble(.init = 60, .step = 12) %>%   # saltar de 12 en 12 (un año)
  model(
    SARIMA_manual = ARIMA(IPC ~ pdq(0,1,0) + PDQ(1,1,1)),
    SARIMA_auto   = ARIMA(IPC, ic = "bic", stepwise = TRUE),  # TRUE más rápido
    ETS_HW        = ETS(IPC)
  ) %>%
  forecast(h = 1) %>%
  accuracy(IPC.tsibble)

cv_ipc %>%
  select(.model, RMSE, MAE, MAPE) %>%
  arrange(MAPE) %>%
  kable(caption = "Validación Cruzada — IPC General", digits = 4) %>% kable_material()
Validación Cruzada — IPC General
.model RMSE MAE MAPE
SARIMA_manual 0.4564 0.3140 0.2762
ETS_HW 0.4372 0.3415 0.2951
SARIMA_auto 0.4439 0.3441 0.2984

Serie 3: Importación de Hidrocarburos (2017–2026) {.tabset}

3A. SARIMA Manual

Órdenes de integración

d_hid <- ndiffs(HID.ts)
D_hid <- nsdiffs(HID.ts)
data.frame(d = d_hid, D = D_hid) %>%
  kable(caption = "Órdenes de Integración — Hidrocarburos") %>% kable_material()
Órdenes de Integración — Hidrocarburos
d D
1 0

Correlograma de la serie diferenciada

hid_diff <- HID.ts
if (D_hid > 0) hid_diff <- diff(hid_diff, lag = 12, differences = D_hid)
if (d_hid > 0) hid_diff <- diff(hid_diff, differences = d_hid)
ts_cor(hid_diff, lag.max = 3 * 12)

Estimación del modelo

# Horizonte: meses restantes de 2026 (mayo–diciembre = 8 meses)
h_hid <- 8

hid.manual <- Arima(HID.ts,
                    order    = c(0, 1, 0),
                    seasonal = c(1, 1, 1))
summary(hid.manual)
## Series: HID.ts 
## ARIMA(0,1,0)(1,1,1)[12] 
## 
## Coefficients:
##          sar1     sma1
##       -0.1901  -1.0000
## s.e.   0.1149   0.1746
## 
## sigma^2 = 1196:  log likelihood = -465.53
## AIC=937.07   AICc=937.34   BIC=944.6
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 0.00404975 31.99418 24.77386 -2.001188 16.5362 0.5366904
##                    ACF1
## Training set -0.4283658

Diagnóstico de residuos

checkresiduals(hid.manual)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0)(1,1,1)[12]
## Q* = 44.344, df = 19, p-value = 0.000846
## 
## Model df: 2.   Total lags used: 21

Pronóstico hasta diciembre 2026

pron.hid.manual <- forecast(hid.manual, h = h_hid, level = c(0.95))
pron.hid.manual %>% print()
##          Point Forecast     Lo 95    Hi 95
## May 2026       238.9895 166.92585 311.0531
## Jun 2026       250.1643 148.25098 352.0776
## Jul 2026       246.5473 121.72952 371.3651
## Aug 2026       242.9220  98.79482 387.0492
## Sep 2026       238.9709  77.83188 400.1100
## Oct 2026       247.3777  71.05398 423.7013
## Nov 2026       249.0047  58.70418 439.3052
## Dec 2026       228.4295  25.11073 431.7483
pron.hid.manual %>% autoplot() +
  labs(title = "Pronóstico Hidrocarburos hasta diciembre 2026 — SARIMA Manual",
       x = "Años/Meses", y = "Millones USD")

3B. SARIMA Automático

hid.auto <- auto.arima(HID.ts)
summary(hid.auto)
## Series: HID.ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.5891
## s.e.   0.0747
## 
## sigma^2 = 960.8:  log likelihood = -499.55
## AIC=1003.1   AICc=1003.22   BIC=1008.37
## 
## Training set error measures:
##                    ME    RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 2.738386 30.6967 24.01432 -2.163135 16.67212 0.520236 -0.03305953
pron.hid.auto <- forecast(hid.auto, h = h_hid, level = c(0.95))
pron.hid.auto %>% print()
##          Point Forecast    Lo 95    Hi 95
## May 2026       220.0131 159.2616 280.7645
## Jun 2026       220.0131 154.3338 285.6923
## Jul 2026       220.0131 149.7507 290.2754
## Aug 2026       220.0131 145.4488 294.5773
## Sep 2026       220.0131 141.3819 298.6442
## Oct 2026       220.0131 137.5153 302.5109
## Nov 2026       220.0131 133.8219 306.2042
## Dec 2026       220.0131 130.2804 309.7457
pron.hid.auto %>% autoplot() +
  labs(title = "Pronóstico Hidrocarburos hasta diciembre 2026 — SARIMA Automático",
       x = "Años/Meses", y = "Millones USD")

3C. Holt-Winters

hid.hw <- HoltWinters(HID.ts, seasonal = "multiplicative")
pron.hid.hw <- forecast(hid.hw, h = h_hid, level = c(0.95))
pron.hid.hw %>% print()
##          Point Forecast    Lo 95    Hi 95
## May 2026       207.1278 167.6315 246.6241
## Jun 2026       214.8071 166.0649 263.5494
## Jul 2026       222.2817 165.3113 279.2521
## Aug 2026       202.3469 142.4713 262.2226
## Sep 2026       205.0296 138.5657 271.4935
## Oct 2026       237.3487 157.3179 317.3794
## Nov 2026       222.1036 141.6939 302.5133
## Dec 2026       215.5426 132.5158 298.5694
pron.hid.hw %>% autoplot() +
  labs(title = "Pronóstico Hidrocarburos hasta diciembre 2026 — Holt-Winters",
       x = "Años/Meses", y = "Millones USD")

3D. Validación Cruzada

HID.tsibble <- HID.ts %>% as_tsibble() %>% rename(HID = value)

cv_hid <- HID.tsibble %>%
  stretch_tsibble(.init = 36, .step = 12) %>% 
  model(
    SARIMA_manual = ARIMA(HID ~ pdq(0,1,0) + PDQ(1,1,1)),
    SARIMA_auto   = ARIMA(HID, ic = "bic", stepwise = TRUE),
    ETS_HW        = ETS(HID)
  ) %>%
  forecast(h = 1) %>%
  accuracy(HID.tsibble)

cv_hid %>%
  select(.model, RMSE, MAE, MAPE) %>%
  arrange(MAPE) %>%
  kable(caption = "Validación Cruzada — Importación de Hidrocarburos", digits = 4) %>% kable_material()
Validación Cruzada — Importación de Hidrocarburos
.model RMSE MAE MAPE
ETS_HW 27.2430 21.7247 13.9994
SARIMA_auto 27.7536 22.7408 14.6048
SARIMA_manual 31.2287 22.3734 16.0565