# ══════════════════════════════════════════════════
# T18 · Caso práctico Financiero — IBEX 35 (Decision Desk)
# Abre primero mSeriesTemporales.Rproj en RStudio
# ══════════════════════════════════════════════════

library(tidyverse)
library(evd)
library(ismev)
library(POT)

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

# ── 1. CARGA Y MAXIMOS POR BLOQUES MENSUALES ──────
source("scripts/deskR.R")
ret_ibex <- desk_returns("_IBEX", type = "log")
ret_ibex$fecha <- as.Date(ret_ibex$date)
perdidas <- -ret_ibex$return
n <- length(perdidas)

ret_ibex$mes <- format(ret_ibex$fecha, "%Y-%m")
bloques <- ret_ibex |> mutate(perdida = -return) |> group_by(mes) |> summarise(max_perdida = max(perdida))
cat(sprintf("Bloques mensuales: %d\n", nrow(bloques)))
pausa()

# ── 2. AJUSTE GEV ──────────────────────────────────
fit_gev <- gev.fit(bloques$max_perdida, show = FALSE)
mu_gev <- fit_gev$mle[1]; sigma_gev <- fit_gev$mle[2]; xi_gev <- fit_gev$mle[3]
cat(sprintf("GEV: mu=%.4f sigma=%.4f xi=%.4f (se=%.4f)\n", mu_gev, sigma_gev, xi_gev, fit_gev$se[3]))

return_level <- function(mu, sigma, xi, T) mu + (sigma / xi) * ((-log(1 - 1 / T))^(-xi) - 1)
cat(sprintf("Nivel de retorno a 1 anio: %.2f%% | a 10 anios: %.2f%%\n",
            return_level(mu_gev, sigma_gev, xi_gev, 12) * 100,
            return_level(mu_gev, sigma_gev, xi_gev, 120) * 100))
cat(sprintf("Peor mes observado en toda la muestra: %.2f%%\n", max(bloques$max_perdida) * 100))
pausa()

# ── 3. PEAKS-OVER-THRESHOLD Y GPD ──────────────────
u <- quantile(perdidas, 0.95)
Nu <- sum(perdidas > u)
fit_gpd <- fitgpd(perdidas, threshold = u, est = "mle")
sigma_gpd <- fit_gpd$fitted.values["scale"]; xi_gpd <- fit_gpd$fitted.values["shape"]
cat(sprintf("GPD: umbral=%.2f%% N_u=%d xi=%.4f (se=%.4f)\n", u * 100, Nu, xi_gpd, fit_gpd$std.err["shape"]))
pausa()

# ── 4. VaR Y ES EXTREMOS AL 99.9% ─────────────────
alpha <- 0.001
VaR_evt <- u + (sigma_gpd / xi_gpd) * (((n / Nu) * alpha)^(-xi_gpd) - 1)
ES_evt  <- (VaR_evt + sigma_gpd - xi_gpd * u) / (1 - xi_gpd)
VaR_hist <- as.numeric(quantile(perdidas, 1 - alpha))
mu_r <- mean(ret_ibex$return); sd_r <- sd(ret_ibex$return)
VaR_norm <- -(mu_r + qnorm(alpha) * sd_r)

cat(sprintf("VaR normal (99.9%%):    %.2f%%\n", VaR_norm * 100))
cat(sprintf("VaR historico (99.9%%): %.2f%%\n", VaR_hist * 100))
cat(sprintf("VaR EVT (99.9%%):       %.2f%%\n", VaR_evt * 100))
cat(sprintf("ES EVT (99.9%%):        %.2f%%\n", ES_evt * 100))

# === FIN Caso Práctico Financiero — Tema 18 (IBEX 35, valores extremos) ===