# ══════════════════════════════════════════════════
# T04 · 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(urca)
library(strucchange)
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_ts  <- ts(log(pib$pib_indice), start = c(1995, 1), frequency = 4)
dlpib    <- diff(log(pib$pib_indice))
dlpib_ts <- ts(dlpib, start = c(1995, 2), frequency = 4)
pausa("\n  [Datos cargados. Pulsa ENTER...]")

# ── 2. BOX-COX: ¿ES EL LOGARITMO LA TRANSFORMACION OPTIMA? ─
lam_g <- forecast::BoxCox.lambda(pib$pib_indice, method = "guerrero")
lam_l <- forecast::BoxCox.lambda(pib$pib_indice, method = "loglik")
cat(sprintf("Lambda (Guerrero) = %.3f   Lambda (log-verosimilitud) = %.3f\n", lam_g, lam_l))
cat("Ninguno coincide exactamente con 0 (logaritmo); se usa log por interpretabilidad economica.\n")
pausa()

# ── 3. ADF: log(PIB) EN NIVELES Y EN DIFERENCIAS ──
adf_niv <- ur.df(lpib_ts, type = "trend", selectlags = "AIC")
adf_dif <- ur.df(dlpib_ts, type = "drift", selectlags = "AIC")
cat(sprintf("ADF log(PIB) niveles    : tau3 = %.3f (cv 5%% = %.2f) -> %s\n",
            adf_niv@teststat[1], adf_niv@cval["tau3", "5pct"],
            ifelse(adf_niv@teststat[1] > adf_niv@cval["tau3", "5pct"], "no se rechaza H0", "se rechaza H0")))
cat(sprintf("ADF diff(log(PIB))      : tau2 = %.3f (cv 5%% = %.2f) -> %s\n",
            adf_dif@teststat[1], adf_dif@cval["tau2", "5pct"],
            ifelse(adf_dif@teststat[1] > adf_dif@cval["tau2", "5pct"], "no se rechaza H0", "se rechaza H0")))
cat("Conclusion: log(PIB) ~ I(1).\n")
pausa()

# ── 4. PHILLIPS-PERRON Y KPSS (CONFIRMACION) ──────
pp_niv   <- ur.pp(lpib_ts, type = "Z-tau", model = "trend")
kpss_niv <- ur.kpss(lpib_ts, type = "tau")
cat(sprintf("PP log(PIB)   : Z-tau = %.3f (cv 5%% = %.2f)\n", pp_niv@teststat, pp_niv@cval[1, "5pct"]))
cat(sprintf("KPSS log(PIB) : %.3f (cv 5%% = %.3f) -> %s\n",
            kpss_niv@teststat, kpss_niv@cval[1, "5pct"],
            ifelse(kpss_niv@teststat > kpss_niv@cval[1, "5pct"], "se rechaza estacionariedad", "no se rechaza")))
pausa()

# ── 5. ZIVOT-ANDREWS: RUPTURA ENDOGENA ────────────
za <- ur.za(lpib_ts, model = "both", lag = NULL)
cat(sprintf("Zivot-Andrews: ruptura en %s | estadistico = %.3f (cv 5%% = %.2f)\n",
            pib$trimestre[za@bpoint], za@teststat, za@cval[2]))
cat("Ni siquiera con la ruptura mas favorable se rechaza la raiz unitaria.\n")
pausa()

# ── 6. MULTIPLES RUPTURAS (BAI-PERRON, strucchange) ─
t_idx <- 1:length(lpib_ts)
bp <- breakpoints(as.numeric(lpib_ts) ~ t_idx)
cat("Rupturas estimadas (BIC):\n")
print(pib$fecha[bp$breakpoints])

plot(lpib_ts, col = "steelblue", main = "log(PIB) con rupturas estimadas")
abline(v = time(lpib_ts)[bp$breakpoints], lty = 2, col = "grey40")

# === FIN Caso Práctico Macro — Tema 04 (PIB, raíces unitarias) ===