# ══════════════════════════════════════════════════
# T03 · 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)

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  <- log(pib$pib_indice)
dlpib <- diff(lpib)                                  # tasa de crecimiento trimestral
dlpib_ts <- ts(dlpib, start = c(1995, 2), frequency = 4)

cat(sprintf("diff(log(PIB)): %d obs, media=%.4f, sd=%.4f\n",
            length(dlpib_ts), mean(dlpib_ts), sd(dlpib_ts)))
pausa()

# ── 2. IDENTIFICACION: FAC Y FACP ─────────────────
op <- par(mfrow = c(1, 2))
acf(dlpib_ts,  lag.max = 12, main = "FAC diff(log(PIB))")
pacf(dlpib_ts, lag.max = 12, main = "FACP diff(log(PIB))")
par(op)
cat("Un unico coeficiente significativo en el retardo 1 (negativo): firma de MA(1).\n")
pausa()

# ── 3. ESTIMACION DEL MA(1) ───────────────────────
m_ma1 <- Arima(dlpib_ts, order = c(0, 0, 1), include.mean = TRUE, method = "ML")
cat(sprintf("MA(1): theta1=%.4f (se=%.4f)  mu=%.4f (se=%.4f)\n",
            coef(m_ma1)["ma1"], sqrt(diag(vcov(m_ma1)))["ma1"],
            coef(m_ma1)["intercept"], sqrt(diag(vcov(m_ma1)))["intercept"]))
pausa()

# ── 4. DIAGNOSTICO DE RESIDUOS ────────────────────
res <- residuals(m_ma1)
op <- par(mfrow = c(1, 2))
plot(res, type = "l", col = "grey30", main = "Residuos MA(1)", ylab = "residuo")
abline(h = 0, lty = 2)
acf(res, main = "FAC de los residuos")
par(op)
lb <- Box.test(res, lag = 8, type = "Ljung-Box", fitdf = 1)
cat(sprintf("Ljung-Box(8) residuos: Q=%.3f  p=%.4f -> %s\n", lb$statistic, lb$p.value,
            ifelse(lb$p.value > 0.05, "no se rechaza incorrelacion (OK)", "quedan restos de estructura")))
pausa()

# ── 5. COMPARACION POR AIC Y BIC ──────────────────
especificaciones <- list("AR(1)" = c(1,0,0), "MA(1)" = c(0,0,1), "AR(2)" = c(2,0,0),
                          "MA(2)" = c(0,0,2), "ARMA(1,1)" = c(1,0,1))
comp <- map_dfr(names(especificaciones), function(nm) {
  m <- Arima(dlpib_ts, order = especificaciones[[nm]], include.mean = TRUE, method = "ML")
  tibble(modelo = nm, AIC = AIC(m), BIC = BIC(m))
}) |> arrange(AIC)
print(comp)
cat(sprintf("Mejor modelo por AIC y BIC: %s\n", comp$modelo[1]))

# === FIN Caso Práctico Macro — Tema 03 (PIB, modelos ARMA) ===