# ══════════════════════════════════════════════════
# T11 · Series Temporales — Evaluación y Combinación de Predicciones
# Abre primero mSeriesTemporales.Rproj en RStudio
# ══════════════════════════════════════════════════

library(tidyverse)
library(forecast)

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

set.seed(2026)

# ── 1. MASE: EJEMPLO SOBRE UNA SERIE SIMULADA ─────
y <- ts(50 + cumsum(rnorm(150, mean = 0.2, sd = 1)) + 3 * sin(2 * pi * (1:150) / 12), frequency = 12)
naive_err <- diff(y, lag = 12)
mae_naive <- mean(abs(naive_err))
m_ets_sim <- ets(y)
mase_ets  <- mean(abs(residuals(m_ets_sim))) / mae_naive
cat(sprintf("MAE naive estacional = %.3f | MASE del ETS ajustado = %.3f\n", mae_naive, mase_ets))
pausa()

# ── 2. VALIDACION CRUZADA TEMPORAL (tsCV) ─────────
f_ets  <- function(x, h) forecast(ets(x), h = h)
f_naive_est <- function(x, h) snaive(x, h = h)
e_ets   <- tsCV(y, f_ets, h = 1, window = 60)
e_naive <- tsCV(y, f_naive_est, h = 1, window = 60)
rmse <- function(e) sqrt(mean(e^2, na.rm = TRUE))
cat(sprintf("RMSE (rolling-origin): ETS=%.3f  naive estacional=%.3f\n", rmse(e_ets), rmse(e_naive)))
pausa()

# ── 3. CONTRASTE DE DIEBOLD-MARIANO ───────────────
comun <- !is.na(e_ets) & !is.na(e_naive)
dm <- dm.test(e_ets[comun], e_naive[comun], h = 1, power = 2)
cat(sprintf("DM ETS vs naive: estadistico=%.3f  p-valor=%.4f\n", dm$statistic, dm$p.value))
cat(ifelse(dm$p.value < 0.05,
           "-> El ETS predice significativamente mejor que el naive estacional.\n",
           "-> No hay evidencia suficiente de que el ETS supere al naive estacional.\n"))
pausa()

# ── 4. COMBINACION DE PREDICCIONES: PESO OPTIMO BATES-GRANGER ─
f_holt <- function(x, h) forecast::holt(x, h = h)
e_holt <- tsCV(y, f_holt, h = 1, window = 60)
comun2 <- !is.na(e_ets) & !is.na(e_holt)
idx <- which(comun2)
mitad <- floor(length(idx) / 2)
idx1 <- idx[1:mitad]; idx2 <- idx[(mitad + 1):length(idx)]

s1 <- sd(e_ets[idx1]); s2 <- sd(e_holt[idx1]); rho <- cor(e_ets[idx1], e_holt[idx1])
w_star <- (s2^2 - rho * s1 * s2) / (s1^2 + s2^2 - 2 * rho * s1 * s2)
cat(sprintf("Peso optimo (ETS) = %.3f | sigma_ets=%.3f sigma_holt=%.3f rho=%.3f\n", w_star, s1, s2, rho))

actual <- as.numeric(y)
fc_ets  <- actual - as.numeric(e_ets)
fc_holt <- actual - as.numeric(e_holt)
e_opt <- actual[idx2] - (w_star * fc_ets[idx2] + (1 - w_star) * fc_holt[idx2])
e_eq  <- actual[idx2] - (0.5 * fc_ets[idx2] + 0.5 * fc_holt[idx2])
cat(sprintf("RMSE fuera de muestra: ETS solo=%.3f Holt solo=%.3f Combo optimo=%.3f Combo igual=%.3f\n",
            rmse(e_ets[idx2]), rmse(e_holt[idx2]), rmse(e_opt), rmse(e_eq)))

# === FIN Script T11 — Evaluación y Combinación de Predicciones ===