T12_mST_CP_Macro.R

Ver crudo
# ══════════════════════════════════════════════════
# T12 · Caso práctico Macro — IPC de España (INE)
# Abre primero mSeriesTemporales.Rproj en RStudio
# Datos: data/ipc_mensual_espana.RData (INE, serie IPC290751)
# ══════════════════════════════════════════════════

library(tidyverse)
library(timetk)
library(modeltime)
library(parsnip)
library(rsample)
library(workflows)
library(recipes)
library(prophet)

pausa <- function(msg = "\n  [Pulsa ENTER para continuar...]") {
  if (interactive()) { cat(msg); invisible(readline()) }
}

recm_f <- function(e) sqrt(mean(e^2, na.rm = TRUE))

# ── 1. CARGA Y PARTICION TEMPORAL ─────────────────
load("data/ipc_mensual_espana.RData")
ipc <- ipc_mensual_espana
ipc$fecha <- as.Date(ipc$fecha)
split_cmp <- time_series_split(ipc, date_var = fecha, assess = "12 months", cumulative = TRUE)
pausa("\n  [Datos cargados. Pulsa ENTER...]")

# ── 2. XGBOOST SOBRE EL NIVEL (SIN DIFERENCIAR): FALLA ─
ipc_feat_nivel <- ipc |> tk_augment_lags(ipc, .lags = c(1, 12)) |> drop_na()
split_nivel <- time_series_split(ipc_feat_nivel, date_var = fecha, assess = "12 months", cumulative = TRUE)
rec_nivel <- recipe(ipc ~ fecha + ipc_lag1 + ipc_lag12, data = training(split_nivel)) |>
  step_timeseries_signature(fecha) |>
  step_rm(matches("(iso$)|(xts$)|(hour)|(minute)|(second)|(am.pm)|(lbl)")) |>
  step_normalize(matches("(index.num)|(year)")) |> step_rm(fecha)
wflow_nivel <- workflow() |> add_model(boost_tree(mode = "regression") |> set_engine("xgboost")) |> add_recipe(rec_nivel)
m_xgb_nivel <- recursive(
  fit(wflow_nivel, data = training(split_nivel)),
  transform = function(data) data |> tk_augment_lags(ipc, .lags = c(1, 12)),
  train_tail = tail(training(split_nivel), 13)
)
fc_nivel <- predict(m_xgb_nivel, testing(split_nivel))$.pred
cat(sprintf("XGBoost (nivel): RECM = %.3f | rango train=[%.1f,%.1f] rango pred=[%.1f,%.1f]\n",
            recm_f(testing(split_nivel)$ipc - fc_nivel),
            min(training(split_nivel)$ipc), max(training(split_nivel)$ipc), min(fc_nivel), max(fc_nivel)))
pausa()

# ── 3. XGBOOST SOBRE LA SERIE DIFERENCIADA: MEJORA ─
ipc_d <- ipc |> mutate(dipc = c(NA, diff(ipc))) |> drop_na(dipc) |> select(fecha, dipc, ipc)
ipc_feat_diff <- ipc_d |> tk_augment_lags(dipc, .lags = c(1, 12)) |> drop_na()
split_diff <- time_series_split(ipc_feat_diff, date_var = fecha, assess = "12 months", cumulative = TRUE)
rec_diff <- recipe(dipc ~ fecha + dipc_lag1 + dipc_lag12, data = training(split_diff)) |>
  step_timeseries_signature(fecha) |>
  step_rm(matches("(iso$)|(xts$)|(hour)|(minute)|(second)|(am.pm)|(lbl)|(diff)")) |>
  step_normalize(matches("(index.num)|(year)")) |> step_rm(fecha)
wflow_diff <- workflow() |> add_model(boost_tree(mode = "regression") |> set_engine("xgboost")) |> add_recipe(rec_diff)
m_xgb_diff <- recursive(
  fit(wflow_diff, data = training(split_diff)),
  transform = function(data) data |> tk_augment_lags(dipc, .lags = c(1, 12)),
  train_tail = tail(training(split_diff), 13)
)
fc_diff <- predict(m_xgb_diff, testing(split_diff))$.pred
ultimo_nivel <- tail(training(split_diff)$ipc, 1)
fc_nivel_diff <- ultimo_nivel + cumsum(fc_diff)
cat(sprintf("XGBoost (diferencias, reconstruido): RECM = %.3f\n",
            recm_f(testing(split_diff)$ipc - fc_nivel_diff)))
pausa()

# ── 4. PROPHET ─────────────────────────────────────
df_prophet <- ipc |> transmute(ds = fecha, y = ipc)
m_prophet <- prophet(df_prophet, yearly.seasonality = TRUE,
                      weekly.seasonality = FALSE, daily.seasonality = FALSE, verbose = FALSE)
futuro <- make_future_dataframe(m_prophet, periods = 12, freq = "month")
fc_prophet <- predict(m_prophet, futuro)
cat("Componentes de Prophet calculados (tendencia + estacional anual).\n")
pausa()

# ── 5. COMPARACION FINAL CON modeltime ────────────
m_arima_mt   <- arima_reg() |> set_engine("auto_arima") |> fit(ipc ~ fecha, data = training(split_cmp))
m_prophet_mt <- prophet_reg(seasonality_yearly = TRUE, seasonality_weekly = FALSE, seasonality_daily = FALSE) |>
  set_engine("prophet") |> fit(ipc ~ fecha, data = training(split_cmp))

tabla <- modeltime_table(m_arima_mt, m_prophet_mt) |>
  modeltime_calibrate(new_data = testing(split_cmp)) |>
  modeltime_accuracy() |>
  select(.model_desc, mae, rmse, mape) |>
  add_row(.model_desc = "XGBOOST (diferencias)", mae = mean(abs(testing(split_diff)$ipc - fc_nivel_diff)),
          rmse = recm_f(testing(split_diff)$ipc - fc_nivel_diff),
          mape = mean(abs(testing(split_diff)$ipc - fc_nivel_diff) / testing(split_diff)$ipc) * 100) |>
  arrange(rmse)
print(tabla)
cat(sprintf("\nMejor modelo: %s\n", tabla$.model_desc[1]))

# === FIN Caso Práctico Macro — Tema 12 (IPC, aprendizaje automatico) ===