# ══════════════════════════════════════════════════
# T05 · Caso práctico Macro — PIB de España (INE)
# Abre primero mSeriesTemporales.Rproj en RStudio
# Datos: data/pib_trimestral_espana.RData (INE CNTR4893 -> nivel)
# ══════════════════════════════════════════════════

library(tidyverse)
library(forecast)
library(tseries)

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

# ── 1. CARGA Y PREPARACION ────────────────────────
load("data/pib_trimestral_espana.RData")
pib <- pib_trimestral_espana
pib$fecha <- as.Date(pib$fecha)
lpib_ts <- ts(log(pib$pib_indice), start = c(1995, 1), frequency = 4)
pausa("\n  [Datos cargados. Pulsa ENTER...]")

# ── 2. LA TRAMPA DE LA DERIVA (con datos reales) ──
m_sin <- Arima(lpib_ts, order = c(0, 1, 1), include.drift = FALSE, method = "ML")
m_con <- Arima(lpib_ts, order = c(0, 1, 1), include.drift = TRUE,  method = "ML")
cat(sprintf("Sin deriva: theta1=%.4f  AIC=%.2f\n", coef(m_sin)["ma1"], AIC(m_sin)))
cat(sprintf("Con deriva: theta1=%.4f  AIC=%.2f (coincide con el MA(1) del Cap.3)\n",
            coef(m_con)["ma1"], AIC(m_con)))
pausa()

# ── 3. auto.arima() COMO VALIDACION CRUZADA ───────
auto_pib <- auto.arima(lpib_ts, seasonal = FALSE, stepwise = FALSE, approximation = FALSE)
cat("Modelo de auto.arima():\n")
print(auto_pib)
pausa()

# ── 4. DIAGNOSTICO: LJUNG-BOX Y JARQUE-BERA ───────
res <- residuals(m_con)
lb <- Box.test(res, lag = 12, type = "Ljung-Box", fitdf = 1)
jb <- jarque.bera.test(res)
cat(sprintf("Ljung-Box(12) : Q=%.3f  p=%.4f -> %s\n", lb$statistic, lb$p.value,
            ifelse(lb$p.value > 0.05, "sin autocorrelacion (OK)", "queda estructura")))
cat(sprintf("Jarque-Bera   : stat=%.1f  p=%.4f -> %s\n", jb$statistic, jb$p.value,
            ifelse(jb$p.value > 0.05, "normal", "NO normal (revisar valores atipicos)")))
idx_outlier <- which.max(abs(res))
cat(sprintf("Mayor residuo: %s (%.4f)\n", pib$trimestre[idx_outlier], res[idx_outlier]))
pausa()

# ── 5. PREDICCION: GAUSSIANA VS BOOTSTRAP ─────────
h <- 8
fc_normal <- forecast(m_con, h = h, level = 95)
set.seed(2026)
fc_boot   <- forecast(m_con, h = h, level = 95, bootstrap = TRUE, npaths = 2000)

pib_fc_normal <- exp(fc_normal$mean)
cat(sprintf("Prediccion PIB a %d trimestres: %.1f -> %.1f (ultimo dato: %.1f)\n",
            h, tail(pib$pib_indice, 1), tail(pib_fc_normal, 1), tail(pib$pib_indice, 1)))

amp_normal <- exp(fc_normal$upper) - exp(fc_normal$lower)
amp_boot   <- exp(fc_boot$upper)   - exp(fc_boot$lower)
cat(sprintf("Amplitud IC95%% h=1 : gaussiano=%.1f  bootstrap=%.1f\n", amp_normal[1], amp_boot[1]))
cat(sprintf("Amplitud IC95%% h=%d: gaussiano=%.1f  bootstrap=%.1f\n", h, amp_normal[h], amp_boot[h]))

plot(fc_normal, main = "Prediccion PIB (log), IC gaussiano al 95%")

# === FIN Caso Práctico Macro — Tema 05 (PIB, ARIMA y predicción) ===