Analisis exploratorio de la tasa de desempleo en Colombia

Autores: Yeison Buitrago, Diego Robayo y Gilbert Rodriguez
Fecha: 13/04/2026

# --- 1.3 Inspeccionar la estructura ---
str(Tasa_Desempleo_raw)

## tibble [299 × 2] (S3: tbl_df/tbl/data.frame)
##  $ Fecha    : POSIXct[1:299], format: "2001-01-31" "2001-02-28" ...
##  $ Desempleo: num [1:299] 16.6 17.4 15.8 14.5 14 ...

head(Tasa_Desempleo_raw, 10)

## # A tibble: 10 × 2
##    Fecha               Desempleo
##    <dttm>                  <dbl>
##  1 2001-01-31 00:00:00      16.6
##  2 2001-02-28 00:00:00      17.4
##  3 2001-03-31 00:00:00      15.8
##  4 2001-04-30 00:00:00      14.5
##  5 2001-05-31 00:00:00      14.0
##  6 2001-06-30 00:00:00      15.3
##  7 2001-07-31 00:00:00      15.2
##  8 2001-08-31 00:00:00      14.9
##  9 2001-09-30 00:00:00      14.1
## 10 2001-10-31 00:00:00      14.5

tail(Tasa_Desempleo_raw, 5)

## # A tibble: 5 × 2
##   Fecha               Desempleo
##   <dttm>                  <dbl>
## 1 2025-07-31 00:00:00      8.84
## 2 2025-08-31 00:00:00      8.6 
## 3 2025-09-30 00:00:00      8.17
## 4 2025-10-31 00:00:00      8.2 
## 5 2025-11-30 00:00:00      7.02

dim(Tasa_Desempleo_raw)

## [1] 299   2

names(Tasa_Desempleo_raw)

## [1] "Fecha"     "Desempleo"

# --- 1.4 Convertir Fecha de POSIXct a Date ---
# R ya leyó las fechas bien desde Excel (POSIXct)
# Solo hay que quitarle la parte de la hora con as.Date()
Tasa_Desempleo_raw <- Tasa_Desempleo_raw %>%
  mutate(Fecha = as.Date(Fecha))

# Verificar que quedó bien
head(Tasa_Desempleo_raw)

## # A tibble: 6 × 2
##   Fecha      Desempleo
##   <date>         <dbl>
## 1 2001-01-31      16.6
## 2 2001-02-28      17.4
## 3 2001-03-31      15.8
## 4 2001-04-30      14.5
## 5 2001-05-31      14.0
## 6 2001-06-30      15.3

class(Tasa_Desempleo_raw$Fecha)   # Debe decir "Date"

## [1] "Date"

# --- 1.5 Crear objeto XTS ---
Fechas_yearmon <- as.yearmon(Tasa_Desempleo_raw$Fecha)

Desempleo_xts <- xts(
  x         = Tasa_Desempleo_raw$Desempleo,
  frequency = 12,
  order.by  = Fechas_yearmon
)

# --- 1.6 Exploración del objeto de serie de tiempo ---
ts_info(Desempleo_xts)    # Inicio, fin, frecuencia, Nº observaciones

##  The Desempleo_xts series is a xts object with 1 variable and 299 observations
##  Frequency: monthly 
##  Start time: ene. 2001 
##  End time: nov. 2025

class(Desempleo_xts)      # Debe decir "xts" "zoo"

## [1] "xts" "zoo"

# --- 1.8 Gráfico interactivo con TSstudio ---
ts_plot(Desempleo_xts,
        title  = "Tasa de Desempleo Mensual - Total Nacional",
        Ytitle = "Tasa de Desempleo (%)",
        Xtitle = "Año",
        Xgrid  = TRUE,
        Ygrid  = TRUE)

# ============================================================
# GRÁFICO 5: Comparación de tres niveles de suavidad LOESS
# Esto te permite ver cuál nivel captura mejor la tendencia
# ============================================================

# Suavidad baja (0.10) → sigue bien los ciclos
p1 <- tbl_desempleo %>%
  plot_time_series(
    .date_var    = Fecha,
    .value       = Desempleo,
    .title       = "LOESS span = 0.10 (detallado)",
    .interactive = FALSE,
    .smooth      = TRUE,
    .smooth_span = 0.10
  )

# Suavidad media (0.25) → balance entre detalle y suavidad
p2 <- tbl_desempleo %>%
  plot_time_series(
    .date_var    = Fecha,
    .value       = Desempleo,
    .title       = "LOESS span = 0.25 (medio)",
    .interactive = FALSE,
    .smooth      = TRUE,
    .smooth_span = 0.25
  )

# Suavidad alta (0.75) → solo muestra tendencia de largo plazo
p3 <- tbl_desempleo %>%
  plot_time_series(
    .date_var    = Fecha,
    .value       = Desempleo,
    .title       = "LOESS span = 0.75 (suave)",
    .interactive = FALSE,
    .smooth      = TRUE,
    .smooth_span = 0.75
  )

# Mostrar los tres juntos
library(patchwork)
p1 / p2 / p3

# --- Extraer valores y fechas de los objetos ya existentes ---
valores          <- as.numeric(Desempleo_xts)
fechas           <- zoo::index(Desempleo_xts)
tiempo_numerico  <- as.numeric(time(Desempleo_xts))

# --- Convertir a objeto ts (necesario para decompose y filter) ---
Desempleo_ts <- ts(valores,
                   start     = c(2001, 1),
                   frequency = 12)


# ============================================================
# MÉTODO 1: DESCOMPOSICIÓN CLÁSICA CON decompose()
# Promedio móvil simple con pesos iguales
# ============================================================

descomp_clasica    <- decompose(Desempleo_ts, type = "additive")
tendencia_decompose <- descomp_clasica$trend

# Gráfico de los cuatro componentes
# --- Gráfico corregido de descomposición clásica ---
plot(descomp_clasica)

# Nota: los primeros y últimos 6 valores serán NA
# porque decompose necesita una ventana completa de 12 meses


# ============================================================
# MÉTODO 2: PROMEDIO MÓVIL CON PESOS — FILTRO BOXCAR
# Pesos: 0.5 en extremos, 1 en el centro, dividido entre 12
# ============================================================

pesos_boxcar      <- c(0.5, rep(1, 11), 0.5) / 12
tendencia_boxcar  <- stats::filter(Desempleo_ts,
                                   filter = pesos_boxcar,
                                   sides  = 2)

# Gráfico
plot(Desempleo_ts,
     main = "Tendencia — Filtro Boxcar (PM con pesos)",
     ylab = "Desempleo (%)",
     col  = "gray50", lwd = 1)
lines(tendencia_boxcar, col = "orange", lwd = 2)
legend("topright",
       legend = c("Serie original", "Tendencia Boxcar"),
       col    = c("gray50", "orange"),
       lwd    = c(1, 2))

# ============================================================
# MÉTODO 3: SUAVIZAMIENTO KERNEL (NADARAYA-WATSON)
# Función ksmooth() con kernel gaussiano
# bandwidth en años: 0.5 = 6 meses, 1 = 1 año, 2 = 2 años
# ============================================================

tendencia_kernel_05 <- ksmooth(tiempo_numerico, valores,
                               kernel    = "normal",
                               bandwidth = 0.5)

tendencia_kernel_1  <- ksmooth(tiempo_numerico, valores,
                               kernel    = "normal",
                               bandwidth = 1)

tendencia_kernel_2  <- ksmooth(tiempo_numerico, valores,
                               kernel    = "normal",
                               bandwidth = 2)

# Gráfico
plot(Desempleo_ts,
     main = "Tendencia — Suavizamiento Kernel (Nadaraya-Watson)",
     ylab = "Desempleo (%)",
     col  = "gray50", lwd = 1)
lines(tendencia_kernel_05$x, tendencia_kernel_05$y, col = "red",       lwd = 2)
lines(tendencia_kernel_1$x,  tendencia_kernel_1$y,  col = "blue",      lwd = 2)
lines(tendencia_kernel_2$x,  tendencia_kernel_2$y,  col = "darkgreen", lwd = 2)
legend("topright",
       legend = c("Serie original",
                  "Kernel bw=0.5 (6 meses)",
                  "Kernel bw=1.0 (1 año)",
                  "Kernel bw=2.0 (2 años)"),
       col    = c("gray50", "red", "blue", "darkgreen"),
       lwd    = c(1, 2, 2, 2))

# ============================================================
# MÉTODO 4: SUAVIZAMIENTO LOWESS
# Función lowess() — regresión local ponderada por vecindad
# f = fracción de datos usados como vecindad en cada punto
# ============================================================

tendencia_lowess_010 <- lowess(tiempo_numerico, valores, f = 0.10)
tendencia_lowess_025 <- lowess(tiempo_numerico, valores, f = 0.25)
tendencia_lowess_075 <- lowess(tiempo_numerico, valores, f = 0.75)

# Gráfico
plot(Desempleo_ts,
     main = "Tendencia — Suavizamiento LOWESS",
     ylab = "Desempleo (%)",
     col  = "gray50", lwd = 1)
lines(fechas, tendencia_lowess_010$y, col = "red",       lwd = 2)
lines(fechas, tendencia_lowess_025$y, col = "blue",      lwd = 2)
lines(fechas, tendencia_lowess_075$y, col = "darkgreen", lwd = 2)
legend("topright",
       legend = c("Serie original",
                  "LOWESS span=0.10",
                  "LOWESS span=0.25",
                  "LOWESS span=0.75"),
       col    = c("gray50", "red", "blue", "darkgreen"),
       lwd    = c(1, 2, 2, 2))

# ============================================================
# MÉTODO 5: DESCOMPOSICIÓN STL
# Usa tsbl_desempleo ya creado en el Paso 3
# Más robusta que decompose() ante el pico del COVID (2020)
# robust = TRUE reduce el efecto de valores atípicos
# ============================================================

# --- Crear descomposición STL ---
descomp_STL <- tsbl_desempleo %>%
  model(
    STL(Desempleo ~ trend() + season(window = "periodic"),
        robust = TRUE)
  ) %>%
  components()

# --- Graficar ---
autoplot(descomp_STL) +
  labs(title = "Descomposición STL — Tasa de Desempleo (2001-2025)")

## Warning: `autoplot.dcmp_ts()` was deprecated in fabletools 0.6.0.
## ℹ Please use `ggtime::autoplot.dcmp_ts()` instead.
## ℹ Graphics functions have been moved to the {ggtime} package. Please use
##   `library(ggtime)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# --- Extraer tendencia y fechas STL ---

tendencia_STL       <- descomp_STL$trend
fechas_STL          <- as.Date(descomp_STL$Fecha)

# ============================================================
# GRÁFICO COMPARATIVO FINAL
# Los cinco métodos sobre la misma serie
# ============================================================

plot(Desempleo_ts,
     main = "Comparación de métodos de estimación de tendencia",
     ylab = "Desempleo (%)",
     xlab = "Año",
     col  = "gray70",
     lwd  = 1)

# M1: decompose
lines(tendencia_decompose,             col = "black",     lwd = 2)

# M2: Boxcar
lines(tendencia_boxcar,                col = "orange",    lwd = 2)

# M3: Kernel (bandwidth = 1 año)
lines(tendencia_kernel_1$x,
      tendencia_kernel_1$y,            col = "red",       lwd = 2)

# M4: LOWESS (span = 0.25)
lines(fechas, tendencia_lowess_025$y,  col = "blue",      lwd = 2)

# M5: STL
lines(fechas_STL, tendencia_STL,       col = "darkgreen", lwd = 2)

legend("topright",
       legend = c("Serie original",
                  "M1: decompose",
                  "M2: Boxcar ",
                  "M3: Kernel (bw=1 año)",
                  "M4: LOWESS (span=0.25)",
                  "M5: STL (robusto)"),
       col    = c("gray70", "black", "orange", "red", "blue", "darkgreen"),
       lwd    = c(1, 2, 2, 2, 2, 2),
       cex    = 0.80)

# --- Extraer tendencia STL y construir data frame de trabajo ---
# descomp_STL ya contiene las columnas: Fecha, Desempleo, trend, season_year, remainder

df_stl <- descomp_STL %>%
  as_tibble() %>%
  mutate(
    Fecha         = as.Date(Fecha),
    Sin_tendencia = Desempleo - trend      # Serie original menos la tendencia STL
  )

# Verificar (las primeras filas tendrán NA en trend → heredados del STL)
head(df_stl %>% select(Fecha, Desempleo, trend, Sin_tendencia), 12)

## # A tibble: 12 × 4
##    Fecha      Desempleo trend Sin_tendencia
##    <date>         <dbl> <dbl>         <dbl>
##  1 2001-01-01      16.6  14.5      2.10    
##  2 2001-02-01      17.4  14.6      2.86    
##  3 2001-03-01      15.8  14.6      1.18    
##  4 2001-04-01      14.5  14.7     -0.178   
##  5 2001-05-01      14.0  14.7     -0.714   
##  6 2001-06-01      15.3  14.8      0.508   
##  7 2001-07-01      15.2  14.9      0.325   
##  8 2001-08-01      14.9  14.9      0.000777
##  9 2001-09-01      14.1  15.0     -0.892   
## 10 2001-10-01      14.5  15.0     -0.490   
## 11 2001-11-01      13.7  15.1     -1.38    
## 12 2001-12-01      13.6  15.1     -1.51

# ============================================================
# GRÁFICO 1: Serie original
# ============================================================

p_original <- ggplot(df_stl, aes(x = Fecha)) +
  geom_line(aes(y = Desempleo), color = "steelblue", linewidth = 0.7) +
  geom_line(aes(y = trend),     color = "firebrick", linewidth = 1.2, na.rm = TRUE) +
  labs(
    title    = "Serie original con tendencia STL",
    subtitle = "Azul = serie | Rojo = tendencia STL",
    x        = NULL,
    y        = "Tasa de desempleo (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

# ============================================================
# GRÁFICO 2: Serie sin tendencia (residuo + estacionalidad)
# ============================================================

p_sin_tend <- ggplot(df_stl %>% filter(!is.na(Sin_tendencia)),
                     aes(x = Fecha, y = Sin_tendencia)) +
  geom_line(color = "darkorange", linewidth = 0.7) +
  geom_hline(yintercept = 0, linetype = "dashed",
             color = "gray50", linewidth = 0.5) +
  labs(
    title    = "Serie sin tendencia (Desempleo − tendencia STL)",
    subtitle = "Conserva la estacionalidad y el residuo",
    x        = "Año",
    y        = "Desempleo − tendencia (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

# ============================================================
# PANEL FINAL: los dos gráficos apilados
# ============================================================

p_original / p_sin_tend

# --- Construir Sin_tendencia_ts (faltaba en el código anterior) ---
df_stl_clean <- df_stl %>% filter(!is.na(Sin_tendencia))

Sin_tendencia_ts <- ts(
  df_stl_clean$Sin_tendencia,
  start     = c(as.integer(format(min(df_stl_clean$Fecha), "%Y")),
                as.integer(format(min(df_stl_clean$Fecha), "%m"))),
  frequency = 12
)

# --- Calcular ACF de ambas series ---
acf_original <- acf(Desempleo_ts,     lag.max = 48, plot = FALSE)
acf_sin_tend <- acf(Sin_tendencia_ts, lag.max = 48, plot = FALSE)

# --- Banda de confianza al 95% ---
n     <- length(Sin_tendencia_ts)
banda <- qnorm(0.975) / sqrt(n)

# --- Data frame unificado ---
df_acf <- bind_rows(
  data.frame(
    rezago = as.numeric(acf_original$lag),
    acf    = as.numeric(acf_original$acf),
    serie  = "Con tendencia"
  ),
  data.frame(
    rezago = as.numeric(acf_sin_tend$lag),
    acf    = as.numeric(acf_sin_tend$acf),
    serie  = "Sin tendencia (STL)"
  )
)

# --- Gráfico comparativo en un solo panel ---
ggplot(df_acf, aes(x = rezago, y = acf, color = serie)) +
  geom_hline(yintercept = 0,      color = "gray40",    linewidth = 0.4) +
  geom_hline(yintercept =  banda, color = "steelblue",
             linetype = "dashed", linewidth = 0.5) +
  geom_hline(yintercept = -banda, color = "steelblue",
             linetype = "dashed", linewidth = 0.5) +
  geom_segment(aes(xend = rezago, yend = 0),
               linewidth = 0.7, alpha = 0.8) +
  geom_point(size = 1.5, alpha = 0.9) +
  facet_wrap(~ serie, ncol = 1) +
  scale_color_manual(values = c("Con tendencia"       = "steelblue",
                                "Sin tendencia (STL)" = "darkorange")) +
  scale_x_continuous(breaks = seq(0, 48, by = 6)) +
  labs(
    title    = "Comparación ACF — Serie original vs sin tendencia (STL)",
    subtitle = "Líneas punteadas azules: bandas de confianza al 95%",
    x        = "Rezago (meses)",
    y        = "Autocorrelación"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    plot.title      = element_text(face = "bold"),
    strip.text      = element_text(face = "bold", size = 11)
  )

# ============================================================
# OPCIÓN 2: PACF en un solo gráfico comparativo con ggplot2
# (mismo enfoque que usamos para el ACF)
# ============================================================

# --- Calcular PACF de ambas series ---
pacf_original <- pacf(Desempleo_ts,     lag.max = 48, plot = FALSE)
pacf_sin_tend <- pacf(Sin_tendencia_ts, lag.max = 48, plot = FALSE)

# --- Banda de confianza al 95% (±1.96/√n) ---
n     <- length(Sin_tendencia_ts)
banda <- qnorm(0.975) / sqrt(n)

# --- Data frame unificado ---
# Nota: el PACF no incluye el rezago 0 (a diferencia del ACF)
df_pacf <- bind_rows(
  data.frame(
    rezago = as.numeric(pacf_original$lag),
    pacf   = as.numeric(pacf_original$acf),
    serie  = "Con tendencia"
  ),
  data.frame(
    rezago = as.numeric(pacf_sin_tend$lag),
    pacf   = as.numeric(pacf_sin_tend$acf),
    serie  = "Sin tendencia (STL)"
  )
)

# --- Gráfico comparativo en un solo panel ---
ggplot(df_pacf, aes(x = rezago, y = pacf, color = serie)) +
  geom_hline(yintercept = 0,      color = "gray40",    linewidth = 0.4) +
  geom_hline(yintercept =  banda, color = "steelblue",
             linetype = "dashed", linewidth = 0.5) +
  geom_hline(yintercept = -banda, color = "steelblue",
             linetype = "dashed", linewidth = 0.5) +
  geom_segment(aes(xend = rezago, yend = 0),
               linewidth = 0.7, alpha = 0.8) +
  geom_point(size = 1.5, alpha = 0.9) +
  facet_wrap(~ serie, ncol = 1) +
  scale_color_manual(values = c("Con tendencia"       = "steelblue",
                                "Sin tendencia (STL)" = "darkorange")) +
  scale_x_continuous(breaks = seq(0, 48, by = 6)) +
  labs(
    title    = "Comparación PACF — Serie original vs sin tendencia (STL)",
    subtitle = "Líneas punteadas azules: bandas de confianza al 95%",
    x        = "Rezago (meses)",
    y        = "Autocorrelación Parcial"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    plot.title      = element_text(face = "bold"),
    strip.text      = element_text(face = "bold", size = 11)
  )

# ============================================================
# MÉTODO 18: GRÁFICOS DE DISPERSIÓN DE RETARDOS
# Aplicado a la serie SIN tendencia (STL)
# Réplica de los chunks "dispersión retardos_1" y
# "dispersion con TSstudio" del profesor
# ============================================================

# ============================================================
# OPCIÓN 4: Ampliar a 24 rezagos (2 años)
# Útil para detectar si hay patrón estacional en los rezagos
# múltiplos de 12 (rezago 12 y 24 son los más importantes
# para una serie mensual con posible ciclo anual)
# ============================================================

par(mar = c(2, 2, 2, 2))

astsa::lag1.plot(Sin_tendencia_ts,
                 max.lag = 24,
                 corr    = FALSE)

# ============================================================
# LECTURA DE RESULTADOS — qué buscar en cada panel
# ============================================================

# RELACIÓN LINEAL POSITIVA: nube de puntos con pendiente
# ascendente → la serie en t está positivamente correlacionada
# con su valor k períodos atrás.

# RELACIÓN LINEAL NEGATIVA: nube con pendiente descendente →
# correlación negativa a ese rezago.

# SIN RELACIÓN: nube circular o sin patrón → rezago no
# significativo.

# RELACIÓN NO LINEAL: la curva loess (línea roja) se curva
# pero la nube no sigue una línea recta → dependencia no
# lineal, que el ACF no detecta pero este gráfico sí muestra.

# REZAGOS CLAVE PARA DESEMPLEO MENSUAL:
# - Rezago 1:  memoria de corto plazo (mes anterior)
# - Rezago 12: ciclo anual (mismo mes del año anterior)
# - Rezago 24: ciclo anual de segundo orden

# ============================================================
# MÉTODO 23: MAPA DE CALOR TEMPORAL
# Réplica de los chunks "ciclos1", "ciclos2" y "mapa de calor"
# del profesor, adaptada a la serie de Desempleo Colombia
# ============================================================

# ============================================================
# MAPA DE CALOR 1: Serie completa (2001-2025)
# Réplica del chunk "ciclos1":
# ts_heatmap(USgas, title = "...")
# ============================================================

ts_heatmap(Desempleo_ts,
           title = "Mapa de Calor — Tasa de Desempleo Mensual Colombia (2001-2025)")

# ============================================================
# MAPA DE CALOR 3: Serie sin tendencia (STL)
# Permite ver la estacionalidad pura sin el efecto
# de la tendencia de largo plazo
# ============================================================

ts_heatmap(Sin_tendencia_ts,
           title = "Mapa de Calor — Desempleo sin tendencia STL (2001-2025)")

# --- Construir tibble directamente desde df_stl_clean
# que ya tiene fechas bien formadas como columna Date ---
tbl_sin_tendencia <- df_stl_clean %>%
  select(Fecha, Sin_tendencia) %>%
  rename(sin_tend = Sin_tendencia) %>%
  mutate(Fecha = as.Date(Fecha))

head(tbl_sin_tendencia)

## # A tibble: 6 × 2
##   Fecha      sin_tend
##   <date>        <dbl>
## 1 2001-01-01    2.10 
## 2 2001-02-01    2.86 
## 3 2001-03-01    1.18 
## 4 2001-04-01   -0.178
## 5 2001-05-01   -0.714
## 6 2001-06-01    0.508

class(tbl_sin_tendencia$Fecha)   # debe decir "Date"

## [1] "Date"

# ============================================================
# BLOQUE 2: BOXPLOT POR TRIMESTRE
# ============================================================

tbl_sin_tendencia %>%
  plot_seasonal_diagnostics(
    .date_var    = Fecha,
    .value       = sin_tend,
    .feature_set = c("quarter"),
    .geom        = "boxplot",
    .interactive = FALSE,
    .title       = "Diagnóstico estacional por trimestre — Desempleo sin tendencia (STL)"
  )

# ============================================================
# BLOQUE 3: BOXPLOT POR AÑO
# ============================================================

tbl_sin_tendencia %>%
  plot_seasonal_diagnostics(
    .date_var    = Fecha,
    .value       = sin_tend,
    .feature_set = c("year"),
    .geom        = "boxplot",
    .interactive = FALSE,
    .title       = "Diagnóstico por año — Desempleo sin tendencia (STL)"
  )

# ============================================================
# MÉTODO 27: SUBSERIES Y GRÁFICOS DE TEMPORADA
# Réplica del chunk "subseries" del profesor:
# UKgrid_tstbl %>% gg_subseries(ND, period=24)
# UKgrid_tstbl %>% gg_season(ND, period=24)
# Adaptada a Desempleo SIN tendencia (STL)

# ============================================================
# PASO PREVIO: construir tsibble desde df_stl_clean
# feasts necesita un objeto tsibble con índice yearmonth
# ============================================================

tsbl_sin_tendencia <- df_stl_clean %>%
  select(Fecha, Sin_tendencia) %>%
  mutate(Fecha = yearmonth(Fecha)) %>%
  as_tsibble(index = Fecha)

# Verificar
tsbl_sin_tendencia

## # A tsibble: 299 x 2 [1M]
##        Fecha Sin_tendencia
##        <mth>         <dbl>
##  1 2001 ene.      2.10    
##  2 2001 feb.      2.86    
##  3 2001 mar.      1.18    
##  4 2001 abr.     -0.178   
##  5 2001 may.     -0.714   
##  6 2001 jun.      0.508   
##  7 2001 jul.      0.325   
##  8 2001 ago.      0.000777
##  9 2001 sep.     -0.892   
## 10 2001 oct.     -0.490   
## # ℹ 289 more rows

class(tsbl_sin_tendencia)

## [1] "tbl_ts"     "tbl_df"     "tbl"        "data.frame"

# ============================================================
# BLOQUE 1: GRÁFICO DE SUBSERIES POR MES
# Réplica de:
# UKgrid_tstbl %>% gg_subseries(ND, period=24)
# Muestra la evolución de cada mes a lo largo de todos
# los años con la media horizontal de ese mes
# period=12 porque la serie es mensual (ciclo anual)
# ============================================================

tsbl_sin_tendencia %>%
  gg_subseries(Sin_tendencia, period = 12) +
  labs(
    title    = "Subseries mensuales — Desempleo sin tendencia (STL)",
    subtitle = "Línea azul horizontal = media de cada mes a lo largo de todos los años",
    y        = "Desempleo − tendencia (%)",
    x        = "Año"
  ) +
  geom_hline(yintercept = 0,
             color      = "firebrick",
             linetype   = "dashed",
             linewidth  = 0.5) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

# ============================================================
# BLOQUE 4 CORREGIDO: SUBSERIES POR TRIMESTRE
# ============================================================

# --- Construir el resumen trimestral desde df_stl_clean
# directamente, sin pasar por el tsibble mensual ---

tsbl_sin_tend_trim <- df_stl_clean %>%
  select(Fecha, Sin_tendencia) %>%
  mutate(
    trimestre = yearquarter(Fecha)   # asignar trimestre a cada fila
  ) %>%
  group_by(trimestre) %>%
  summarise(
    sin_tend_trim = mean(Sin_tendencia, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  as_tsibble(index = trimestre)      # ahora sí hay un valor único por trimestre

# Verificar que no hay duplicados
duplicates(tsbl_sin_tend_trim, index = trimestre)

## # A tibble: 0 × 2
## # ℹ 2 variables: trimestre <qtr>, sin_tend_trim <dbl>

# Verificar estructura
tsbl_sin_tend_trim

## # A tsibble: 100 x 2 [1Q]
##    trimestre sin_tend_trim
##        <qtr>         <dbl>
##  1   2001 Q1        2.05  
##  2   2001 Q2       -0.128 
##  3   2001 Q3       -0.189 
##  4   2001 Q4       -1.13  
##  5   2002 Q1        1.18  
##  6   2002 Q2        0.256 
##  7   2002 Q3       -0.255 
##  8   2002 Q4       -0.192 
##  9   2003 Q1       -0.0369
## 10   2003 Q2       -0.798 
## # ℹ 90 more rows

# --- Gráfico de subseries trimestrales ---
tsbl_sin_tend_trim %>%
  gg_subseries(sin_tend_trim) +
  labs(
    title    = "Subseries trimestrales — Desempleo sin tendencia (STL)",
    subtitle = "Línea azul = media de cada trimestre a lo largo de todos los años",
    y        = "Desempleo − tendencia (%)",
    x        = "Año"
  ) +
  geom_hline(yintercept = 0,
             color      = "firebrick",
             linetype   = "dashed",
             linewidth  = 0.5) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

# ============================================================
# MÉTODO 28: REGRESIÓN TRIGONOMÉTRICA PARA ESTIMACIÓN DE CICLOS
# Réplica de los chunks "regresión estimar ciclo_1",
# "regresión estimar ciclo_2" y "regresión estimar ciclo_3"
# del profesor, adaptada a Desempleo SIN tendencia (STL)
# ============================================================

# ============================================================
# PASO PREVIO: definir el vector de tiempo y la serie
# ============================================================

n_obs  <- length(Sin_tendencia_ts)
tiempo <- 1:n_obs          # índice entero 1, 2, ..., n
y      <- as.numeric(Sin_tendencia_ts)

# ============================================================
# BLOQUE 1: VISUALIZACIÓN INICIAL
# Réplica del chunk "regresión estimar ciclo_1":
# plot.ts(cs), plot.ts(cs+w), plot.ts(cs+5*w)
# Permite decidir visualmente qué frecuencias explorar
# Para desempleo mensual los candidatos naturales son:
#   ω = 1/12  → ciclo anual      (periodo = 12 meses)
#   ω = 1/6   → ciclo semestral  (periodo = 6 meses)
#   ω = 1/60  → ciclo de 5 años  (periodo = 60 meses)
# ============================================================

par(mfrow = c(3, 1), mar = c(3, 4, 3, 2))

# --- Gráfico 1: Ciclo 12 meses ---
plot.ts(y,
        main = "Serie sin tendencia — Desempleo Colombia",
        ylab = "Desempleo − tendencia (%)",
        col  = "steelblue")
lines(2 * cos(2 * pi * tiempo / 12), col = "firebrick", lwd = 2)
legend("topright",
       legend = c("Serie", "Ciclo 12 meses"),
       col    = c("steelblue", "firebrick"),
       lwd    = c(1, 2), cex = 0.8)

# --- Gráfico 2: Ciclo 6 meses ---
plot.ts(y,
        main = "Serie sin tendencia — Desempleo Colombia",
        ylab = "Desempleo − tendencia (%)",
        col  = "steelblue")
lines(2 * cos(2 * pi * tiempo / 6), col = "darkgreen", lwd = 2)
legend("topright",
       legend = c("Serie", "Ciclo 6 meses"),
       col    = c("steelblue", "darkgreen"),
       lwd    = c(1, 2), cex = 0.8)

# --- Gráfico 3: Ciclo 60 meses ---
plot.ts(y,
        main = "Serie sin tendencia — Desempleo Colombia",
        ylab = "Desempleo − tendencia (%)",
        col  = "steelblue")
lines(2 * cos(2 * pi * tiempo / 60), col = "darkorange", lwd = 2)
legend("topright",
       legend = c("Serie", "Ciclo 60 meses"),
       col    = c("steelblue", "darkorange"),
       lwd    = c(1, 2), cex = 0.8)

par(mfrow = c(1, 1))

# ============================================================
# BLOQUE 2: REGRESIÓN TRIGONOMÉTRICA — CICLO ANUAL (ω = 1/12)
# Réplica del chunk "regresión estimar ciclo_2":
# z1 = cos(2*pi*1:500/50)
# z2 = sin(2*pi*1:500/50)
# summary(fit <- lm(x~0+z1+z2))
# ============================================================

# --- Construir las covariables trigonométricas ---
# Para ω = 1/12 (ciclo anual, periodo = 12 meses)
z1_12 <- cos(2 * pi * tiempo / 12)
z2_12 <- sin(2 * pi * tiempo / 12)

# --- Ajustar el modelo sin intercepto ---
fit_12 <- lm(y ~ 0 + z1_12 + z2_12)
summary(fit_12)

## 
## Call:
## lm(formula = y ~ 0 + z1_12 + z2_12)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6610 -0.5885 -0.1394  0.4750  8.8002 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## z1_12 -0.08534    0.10832  -0.788    0.431    
## z2_12  0.75220    0.10796   6.967 2.09e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.322 on 297 degrees of freedom
## Multiple R-squared:  0.142,  Adjusted R-squared:  0.1362 
## F-statistic: 24.58 on 2 and 297 DF,  p-value: 1.321e-10

# --- Amplitud estimada del ciclo anual ---
# A² = β1² + β2²  →  A = √(β1² + β2²)
beta1_12 <- coef(fit_12)["z1_12"]
beta2_12 <- coef(fit_12)["z2_12"]
A_12     <- sqrt(beta1_12^2 + beta2_12^2)
cat("Amplitud estimada ciclo anual (12 meses):", round(A_12, 4), "\n")

## Amplitud estimada ciclo anual (12 meses): 0.757

# --- Fase estimada ---
fase_12 <- atan2(-beta2_12, beta1_12)
cat("Fase estimada (radianes):", round(fase_12, 4), "\n")

## Fase estimada (radianes): -1.6838

# ============================================================
# BLOQUE 3: REGRESIÓN TRIGONOMÉTRICA — CICLO SEMESTRAL (ω = 1/6)
# ============================================================

z1_6 <- cos(2 * pi * tiempo / 6)
z2_6 <- sin(2 * pi * tiempo / 6)

fit_6 <- lm(y ~ 0 + z1_6 + z2_6)
summary(fit_6)

## 
## Call:
## lm(formula = y ~ 0 + z1_6 + z2_6)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8966 -0.6578 -0.1504  0.6539  9.9139 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## z1_6  0.01215    0.10793   0.113     0.91    
## z2_6  0.77337    0.10757   7.190 5.29e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.317 on 297 degrees of freedom
## Multiple R-squared:  0.1483, Adjusted R-squared:  0.1425 
## F-statistic: 25.85 on 2 and 297 DF,  p-value: 4.463e-11

beta1_6 <- coef(fit_6)["z1_6"]
beta2_6 <- coef(fit_6)["z2_6"]
A_6     <- sqrt(beta1_6^2 + beta2_6^2)
cat("Amplitud estimada ciclo semestral (6 meses):", round(A_6, 4), "\n")

## Amplitud estimada ciclo semestral (6 meses): 0.7735

# ============================================================
# BLOQUE 4: REGRESIÓN TRIGONOMÉTRICA — CICLO LARGO (ω = 1/60)
# Ciclo de 5 años, posible ciclo económico
# ============================================================

z1_60 <- cos(2 * pi * tiempo / 60)
z2_60 <- sin(2 * pi * tiempo / 60)

fit_60 <- lm(y ~ 0 + z1_60 + z2_60)
summary(fit_60)

## 
## Call:
## lm(formula = y ~ 0 + z1_60 + z2_60)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1353 -0.6912 -0.1154  0.5592  8.9356 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)  
## z1_60   0.2958     0.1154   2.564   0.0108 *
## z2_60  -0.1416     0.1150  -1.231   0.2192  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.408 on 297 degrees of freedom
## Multiple R-squared:  0.02652,    Adjusted R-squared:  0.01996 
## F-statistic: 4.045 on 2 and 297 DF,  p-value: 0.01849

beta1_60 <- coef(fit_60)["z1_60"]
beta2_60 <- coef(fit_60)["z2_60"]
A_60     <- sqrt(beta1_60^2 + beta2_60^2)
cat("Amplitud estimada ciclo largo (60 meses):", round(A_60, 4), "\n")

## Amplitud estimada ciclo largo (60 meses): 0.328

# ============================================================
# BLOQUE 5: MODELO COMBINADO — ciclo anual + semestral
# Réplica de la idea del profesor de combinar frecuencias:
# x_t = β1·cos(2πt/12) + β2·sin(2πt/12)
#      + β3·cos(2πt/6)  + β4·sin(2πt/6) + w_t
# ============================================================

fit_combinado <- lm(y ~ 0 + z1_12 + z2_12 + z1_6 + z2_6)
summary(fit_combinado)

## 
## Call:
## lm(formula = y ~ 0 + z1_12 + z2_12 + z1_6 + z2_6)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6494 -0.3898 -0.0828  0.3997  9.4643 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## z1_12 -0.08527    0.09886  -0.863    0.389    
## z2_12  0.75220    0.09852   7.635 3.16e-13 ***
## z1_6   0.01157    0.09886   0.117    0.907    
## z2_6   0.77337    0.09852   7.850 7.77e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.207 on 295 degrees of freedom
## Multiple R-squared:  0.2903, Adjusted R-squared:  0.2807 
## F-statistic: 30.17 on 4 and 295 DF,  p-value: < 2.2e-16

cat("R² modelo combinado (anual + semestral):",
    round(summary(fit_combinado)$r.squared, 4), "\n")

## R² modelo combinado (anual + semestral): 0.2903

# ============================================================
# BLOQUE 6: GRÁFICO COMPARATIVO — serie vs ajuste
# Réplica del chunk "regresión estimar ciclo_3":
# plot.ts(x, col=8)
# lines(fitted(fit), col=2)
# ============================================================

par(mfrow = c(3, 1), mar = c(3, 4, 3, 2))

# --- Ciclo anual ---
plot.ts(y, col = "gray60",
        main = "Serie vs ajuste — Ciclo anual (12 meses)",
        ylab = "Desempleo − tendencia (%)")
lines(fitted(fit_12), col = "firebrick", lwd = 2)
legend("topright",
       legend = c("Serie sin tendencia", "Ajuste trigonométrico"),
       col    = c("gray60", "firebrick"), lwd = c(1, 2), cex = 0.8)

# --- Ciclo semestral ---
plot.ts(y, col = "gray60",
        main = "Serie vs ajuste — Ciclo semestral (6 meses)",
        ylab = "Desempleo − tendencia (%)")
lines(fitted(fit_6), col = "darkgreen", lwd = 2)
legend("topright",
       legend = c("Serie sin tendencia", "Ajuste trigonométrico"),
       col    = c("gray60", "darkgreen"), lwd = c(1, 2), cex = 0.8)

# --- Modelo combinado ---
plot.ts(y, col = "gray60",
        main = "Serie vs ajuste — Modelo combinado (12 + 6 meses)",
        ylab = "Desempleo − tendencia (%)")
lines(fitted(fit_combinado), col = "darkorange", lwd = 2)
legend("topright",
       legend = c("Serie sin tendencia", "Ajuste combinado"),
       col    = c("gray60", "darkorange"), lwd = c(1, 2), cex = 0.8)

par(mfrow = c(1, 1))

# ============================================================
# BLOQUE 7: TABLA RESUMEN DE AMPLITUDES Y R²
# Permite comparar cuál ciclo explica más varianza
# ============================================================

resumen_ciclos <- data.frame(
  Ciclo      = c("Anual (12 meses)",
                 "Semestral (6 meses)",
                 "Largo (60 meses)",
                 "Combinado (12+6)"),
  Periodo    = c(12, 6, 60, NA),
  Amplitud   = c(round(A_12, 4),
                 round(A_6,  4),
                 round(A_60, 4),
                 NA),
  R_cuadrado = c(round(summary(fit_12)$r.squared,        4),
                 round(summary(fit_6)$r.squared,         4),
                 round(summary(fit_60)$r.squared,        4),
                 round(summary(fit_combinado)$r.squared, 4))
)

print(resumen_ciclos)

##                     Ciclo Periodo Amplitud R_cuadrado
## z1_12    Anual (12 meses)      12   0.7570     0.1420
## z1_6  Semestral (6 meses)       6   0.7735     0.1483
## z1_60    Largo (60 meses)      60   0.3280     0.0265
##          Combinado (12+6)      NA       NA     0.2903

# ============================================================
# MÉTODOS 29 Y 30: PERIODOGRAMA Y PERIODOGRAMA SUAVIZADO
# Réplica de los chunks "Periodograma_1", "Periodograma_2",
# "espectral_soi" y "espectral_rec" del profesor,
# adaptada a Desempleo SIN tendencia (STL)
# ============================================================

# ============================================================
# BLOQUE 1: PERIODOGRAMA CRUDO — función spectrum()
# Réplica del chunk "Periodograma_1":
# spectrum(x, log='no')
# ============================================================

par(mfrow = c(1, 1), mar = c(4, 4, 3, 2))

per_crudo <- spectrum(y,
                      log  = "no",
                      main = "Periodograma — Desempleo sin tendencia (STL)",
                      ylab = "Densidad espectral",
                      xlab = "Frecuencia (ciclos por mes)")

# --- Identificar la frecuencia dominante ---
ubicacion_max <- which.max(per_crudo$spec)

frec_max    <- per_crudo$freq[ubicacion_max]
periodo_max <- 1 / frec_max

cat("Frecuencia donde se maximiza el periodograma:",
    round(frec_max, 6), "\n")

## Frecuencia donde se maximiza el periodograma: 0.166667

cat("Periodo correspondiente (meses):",
    round(periodo_max, 2), "\n")

## Periodo correspondiente (meses): 6

cat("Periodo correspondiente (años):",
    round(periodo_max / 12, 2), "\n")

## Periodo correspondiente (años): 0.5

# Marcar la frecuencia dominante en el gráfico
abline(v   = frec_max,
       lty = 2,
       col = "firebrick",
       lwd = 2)

# --- Identificar la segunda frecuencia más importante ---
n_spec          <- length(per_crudo$spec)
valor_seg       <- sort(per_crudo$spec,
                        partial = n_spec - 1)[n_spec - 1]
ubicacion_seg   <- which(per_crudo$spec == valor_seg)
frec_seg        <- per_crudo$freq[ubicacion_seg]
periodo_seg     <- 1 / frec_seg

cat("Segunda frecuencia dominante:",
    round(frec_seg, 6), "\n")

## Segunda frecuencia dominante: 0.083333

cat("Periodo correspondiente (meses):",
    round(periodo_seg, 2), "\n")

## Periodo correspondiente (meses): 12

abline(v   = frec_seg,
       lty = 2,
       col = "darkorange",
       lwd = 2)

legend("topright",
       legend = c(
         paste0("Frec. principal: ", round(frec_max, 4),
                " (", round(periodo_max, 1), " meses)"),
         paste0("Frec. secundaria: ", round(frec_seg, 4),
                " (", round(periodo_seg, 1), " meses)")
       ),
       col    = c("firebrick", "darkorange"),
       lty    = 2, lwd = 2, cex = 0.85)

# ============================================================
# BLOQUE 2: PERIODOGRAMA VÍA TRANSFORMADA DE FOURIER (FFT)
# Réplica del chunk "Periodograma_1":
# Px = Mod(fft(x))^2
# ============================================================

Px   <- Mod(fft(y))^2
Freq <- 0:(n_obs - 1) / n_obs

# Solo graficar la mitad útil (hasta la frecuencia Nyquist)
mitad <- 1:(n_obs %/% 2)

plot(Freq[mitad], Px[mitad],
     type = "l",
     main = "Periodograma vía FFT — Desempleo sin tendencia",
     ylab = "Potencia espectral",
     xlab = "Frecuencia (ciclos por mes)",
     col  = "steelblue")

u <- which.max(Px[mitad])
cat("Frecuencia máxima vía FFT:",
    round(Freq[u], 6), "\n")

## Frecuencia máxima vía FFT: 0.167224

cat("Periodo correspondiente (meses):",
    round(1 / Freq[u], 2), "\n")

## Periodo correspondiente (meses): 5.98

abline(v   = Freq[u],
       lty = 2,
       col = "firebrick",
       lwd = 2)

# ============================================================
# BLOQUE 3: PERIODOGRAMA CON mvspec() de astsa
# Réplica del chunk "espectral_soi":
# soi.per = astsa::mvspec(soi, log="no")
# Incluye bandas de confianza al 95%
# ============================================================

par(mar = c(4, 4, 3, 2))

per_mvspec <- astsa::mvspec(y,
                            log  = "no",
                            main = "Periodograma mvspec — Desempleo sin tendencia (STL)")

# --- Frecuencia dominante con mvspec ---
ubicacion_mv <- which.max(per_mvspec$spec)
frec_mv      <- per_mvspec$freq[ubicacion_mv]
periodo_mv   <- 1 / frec_mv

cat("Frecuencia dominante (mvspec):",
    round(frec_mv, 6), "\n")

## Frecuencia dominante (mvspec): 0.166667

cat("Periodo correspondiente (meses):",
    round(periodo_mv, 2), "\n")

## Periodo correspondiente (meses): 6

# --- Segunda frecuencia (mvspec) ---
n_mv          <- length(per_mvspec$spec)
valor_seg_mv  <- sort(per_mvspec$spec,
                      partial = n_mv - 1)[n_mv - 1]
ubicacion_seg_mv <- which(per_mvspec$spec == valor_seg_mv)
frec_seg_mv      <- per_mvspec$freq[ubicacion_seg_mv]
periodo_seg_mv   <- 1 / frec_seg_mv

cat("Segunda frecuencia dominante (mvspec):",
    round(frec_seg_mv, 6), "\n")

## Segunda frecuencia dominante (mvspec): 0.083333

cat("Periodo correspondiente (meses):",
    round(periodo_seg_mv, 2), "\n")

## Periodo correspondiente (meses): 12

abline(v   = frec_mv,
       lty = 2, col = "firebrick", lwd = 2)
abline(v   = frec_seg_mv,
       lty = 2, col = "darkorange", lwd = 2)

# ============================================================
# MÉTODO 30: PERIODOGRAMA SUAVIZADO
# Réplica del chunk "Periodograma_2":
# spectrum(x, log='no', span=5)
# spectrum(x, log='no', span=c(5,5))
# El argumento span promedia frecuencias vecinas
# para reducir el ruido del periodograma crudo
# ============================================================

par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))

# --- span = 5 (suavizado leve) ---
spectrum(y,
         log  = "no",
         span = 5,
         main = "Periodograma suavizado — span = 5",
         ylab = "Densidad espectral",
         xlab = "Frecuencia",
         col  = "steelblue")
abline(v   = frec_max,
       lty = 2, col = "firebrick", lwd = 1.5)

# --- span = c(5,5) (dos pasadas de suavizado) ---
spectrum(y,
         log  = "no",
         span = c(5, 5),
         main = "Periodograma suavizado — span = c(5,5)",
         ylab = "Densidad espectral",
         xlab = "Frecuencia",
         col  = "steelblue")
abline(v   = frec_max,
       lty = 2, col = "firebrick", lwd = 1.5)

# --- span = c(3,3) (suavizado moderado) ---
spectrum(y,
         log  = "no",
         span = c(3, 3),
         main = "Periodograma suavizado — span = c(3,3)",
         ylab = "Densidad espectral",
         xlab = "Frecuencia",
         col  = "steelblue")
abline(v   = frec_max,
       lty = 2, col = "firebrick", lwd = 1.5)

# --- span = c(9,9) (suavizado fuerte) ---
spectrum(y,
         log  = "no",
         span = c(9, 9),
         main = "Periodograma suavizado — span = c(9,9)",
         ylab = "Densidad espectral",
         xlab = "Frecuencia",
         col  = "steelblue")
abline(v   = frec_max,
       lty = 2, col = "firebrick", lwd = 1.5)

par(mfrow = c(1, 1))

# ============================================================
# BLOQUE 5: TABLA RESUMEN DE FRECUENCIAS DETECTADAS
# Réplica del estilo de interpretación del profesor
# en los chunks "espectral_soi" y "espectral_rec"
# ============================================================

# Extraer las 5 frecuencias más importantes
top5_idx  <- order(per_crudo$spec, decreasing = TRUE)[1:5]

resumen_frecuencias <- data.frame(
  Rango        = 1:5,
  Frecuencia   = round(per_crudo$freq[top5_idx], 6),
  Periodo_mes  = round(1 / per_crudo$freq[top5_idx], 2),
  Periodo_anio = round(1 / per_crudo$freq[top5_idx] / 12, 2),
  Potencia     = round(per_crudo$spec[top5_idx], 4)
)

print(resumen_frecuencias)

##   Rango Frecuencia Periodo_mes Periodo_anio Potencia
## 1     1   0.166667           6         0.50  40.2163
## 2     2   0.083333          12         1.00  39.4192
## 3     3   0.250000           4         0.33  16.1453
## 4     4   0.333333           3         0.25  15.8235
## 5     5   0.013333          75         6.25   9.5824

# ============================================================
# MÉTODO 31: ANÁLISIS DE ESTACIONALIDAD CON VARIABLES DUMMY
# Réplica del chunk "Dummy Est" del profesor,
# adaptada a Desempleo SIN tendencia (STL)
# ============================================================


# ============================================================
# BLOQUE 1: CONSTRUCCIÓN DE LA MATRIZ DE DUMMIES MENSUALES
# Réplica del chunk "Dummy Est":
# s=4, l=diag(s), X=matrix(rep(t(l),T/s),...)
# Aquí s=12 porque la serie es mensual
# ============================================================

s     <- 12                          # número de estaciones (meses)
T_obs <- length(Sin_tendencia_ts)    # número total de observaciones

# --- Construir matriz de dummies estacionales ---
# Cada columna indica si la observación pertenece a ese mes
# (1=enero, 2=febrero, ..., 12=diciembre)
identidad <- diag(s)
X_dummy   <- matrix(
  rep(t(identidad), ceiling(T_obs / s)),
  ncol   = ncol(identidad),
  byrow  = TRUE
)[1:T_obs, ]

# Asignar nombres de columnas
colnames(X_dummy) <- month.abb

# Verificar estructura
cat("Dimensión de la matriz de dummies:", dim(X_dummy), "\n")

## Dimensión de la matriz de dummies: 299 12

head(X_dummy, 14)

##       Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
##  [1,]   1   0   0   0   0   0   0   0   0   0   0   0
##  [2,]   0   1   0   0   0   0   0   0   0   0   0   0
##  [3,]   0   0   1   0   0   0   0   0   0   0   0   0
##  [4,]   0   0   0   1   0   0   0   0   0   0   0   0
##  [5,]   0   0   0   0   1   0   0   0   0   0   0   0
##  [6,]   0   0   0   0   0   1   0   0   0   0   0   0
##  [7,]   0   0   0   0   0   0   1   0   0   0   0   0
##  [8,]   0   0   0   0   0   0   0   1   0   0   0   0
##  [9,]   0   0   0   0   0   0   0   0   1   0   0   0
## [10,]   0   0   0   0   0   0   0   0   0   1   0   0
## [11,]   0   0   0   0   0   0   0   0   0   0   1   0
## [12,]   0   0   0   0   0   0   0   0   0   0   0   1
## [13,]   1   0   0   0   0   0   0   0   0   0   0   0
## [14,]   0   1   0   0   0   0   0   0   0   0   0   0

# ============================================================
# BLOQUE 2: AJUSTE DEL MODELO CON DUMMIES
# Réplica del chunk "Dummy Est":
# salida_dummy_ruido = stan_glm(resp~.-1, data=simul)
# Aquí usamos lm() sin intercepto (igual que el profesor)
# El modelo es:
# Y_t = δ1·D1t + δ2·D2t + ... + δ12·D12t + εt
# donde Dit = 1 si la observación t es del mes i
# ============================================================

df_dummy <- data.frame(
  Y      = as.numeric(Sin_tendencia_ts),
  X_dummy
)

# Ajuste sin intercepto (las 12 dummies lo reemplazan)
fit_dummy <- lm(Y ~ . - 1, data = df_dummy)
summary(fit_dummy)

## 
## Call:
## lm(formula = Y ~ . - 1, data = df_dummy)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2452 -0.4339 -0.1157  0.2014  8.9598 
## 
## Coefficients:
##     Estimate Std. Error t value Pr(>|t|)    
## Jan  2.27104    0.22031  10.308  < 2e-16 ***
## Feb  1.26849    0.22031   5.758 2.19e-08 ***
## Mar  0.16881    0.22031   0.766 0.444163    
## Apr  0.43371    0.22031   1.969 0.049957 *  
## May  0.29042    0.22031   1.318 0.188467    
## Jun  0.03254    0.22031   0.148 0.882673    
## Jul  0.47047    0.22031   2.135 0.033567 *  
## Aug -0.04492    0.22031  -0.204 0.838575    
## Sep -0.39547    0.22031  -1.795 0.073689 .  
## Oct -0.90654    0.22031  -4.115 5.07e-05 ***
## Nov -1.21057    0.22031  -5.495 8.63e-08 ***
## Dec -0.75222    0.22485  -3.345 0.000931 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.102 on 287 degrees of freedom
## Multiple R-squared:  0.4246, Adjusted R-squared:  0.4005 
## F-statistic: 17.65 on 12 and 287 DF,  p-value: < 2.2e-16

# --- Extraer coeficientes (efectos estacionales estimados) ---
efectos_estacionales <- coef(fit_dummy)
cat("\nEfectos estacionales estimados por mes:\n")

## 
## Efectos estacionales estimados por mes:

print(round(efectos_estacionales, 4))

##     Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep     Oct 
##  2.2710  1.2685  0.1688  0.4337  0.2904  0.0325  0.4705 -0.0449 -0.3955 -0.9065 
##     Nov     Dec 
## -1.2106 -0.7522

# ============================================================
# BLOQUE 3: GRÁFICO DE LOS EFECTOS ESTACIONALES ESTIMADOS
# Permite ver qué meses tienen desempleo sistemáticamente
# por encima o por debajo de la tendencia
# ============================================================

df_efectos <- data.frame(
  mes    = factor(month.abb, levels = month.abb),
  efecto = as.numeric(efectos_estacionales)
)

ggplot(df_efectos, aes(x = mes, y = efecto,
                       fill = efecto > 0)) +
  geom_col(alpha = 0.85, width = 0.7) +
  geom_hline(yintercept = 0,
             color      = "black",
             linetype   = "dashed",
             linewidth  = 0.6) +
  geom_text(aes(label = round(efecto, 2),
                vjust = ifelse(efecto > 0, -0.4, 1.2)),
            size = 3.2) +
  scale_fill_manual(values = c("TRUE"  = "steelblue",
                               "FALSE" = "firebrick"),
                    guide  = "none") +
  labs(
    title    = "Efectos estacionales estimados — Modelo Dummy",
    subtitle = "Azul = meses con desempleo sobre la tendencia | Rojo = por debajo",
    x        = "Mes",
    y        = "Efecto estacional estimado (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

# ============================================================
# BLOQUE 4 CORREGIDO: SERIE OBSERVADA vs AJUSTE POR DUMMIES
# ============================================================

df_ajuste <- data.frame(
  tiempo    = as.numeric(time(Sin_tendencia_ts)),
  observado = as.numeric(Sin_tendencia_ts),
  ajustado  = as.numeric(fitted(fit_dummy)),
  residuo   = as.numeric(residuals(fit_dummy))
)

par(mfrow = c(2, 1), mar = c(3, 4, 3, 2))

# --- Panel 1: Serie observada vs ajuste ---
plot(df_ajuste$tiempo,
     df_ajuste$observado,
     type = "l",
     col  = "gray60",
     lwd  = 1,
     main = "Serie sin tendencia vs ajuste por dummies estacionales",
     ylab = "Desempleo − tendencia (%)",
     xlab = "Año",
     ylim = range(c(df_ajuste$observado,
                    df_ajuste$ajustado)))

lines(df_ajuste$tiempo,
      df_ajuste$ajustado,
      col = "firebrick",
      lwd = 2)

abline(h   = 0,
       lty = 2,
       col = "black",
       lwd = 0.8)

legend("topright",
       legend = c("Serie sin tendencia", "Ajuste dummies"),
       col    = c("gray60", "firebrick"),
       lwd    = c(1, 2),
       cex    = 0.85)

# --- Panel 2: Residuos ---
plot(df_ajuste$tiempo,
     df_ajuste$residuo,
     type = "l",
     col  = "steelblue",
     lwd  = 1,
     main = "Residuos — Modelo dummy estacional",
     ylab = "Residuo (%)",
     xlab = "Año")

abline(h   = 0,
       lty = 2,
       col = "firebrick",
       lwd = 1)

par(mfrow = c(1, 1))

# ============================================================
# BLOQUE 5: PERIODOGRAMA — serie y residuos del modelo dummy
# ============================================================

par(mfrow = c(2, 1), mar = c(4, 4, 3, 2))

per_dummy <- spectrum(df_ajuste$observado,
                      log  = "no",
                      main = "Periodograma — Serie sin tendencia",
                      ylab = "Densidad espectral",
                      xlab = "Frecuencia")

ubicacion_dummy <- which.max(per_dummy$spec)
frec_dom        <- per_dummy$freq[ubicacion_dummy]

cat("Frecuencia dominante:",       round(frec_dom, 6), "\n")

## Frecuencia dominante: 0.166667

cat("Periodo (meses):",            round(1 / frec_dom, 2), "\n")

## Periodo (meses): 6

abline(v   = frec_dom,
       lty = 2, col = "firebrick", lwd = 2)

per_residuos <- spectrum(df_ajuste$residuo,
                         log  = "no",
                         main = "Periodograma — Residuos del modelo dummy",
                         ylab = "Densidad espectral",
                         xlab = "Frecuencia")

abline(v   = 1/12,
       lty = 2, col = "firebrick", lwd = 1.5)

text(x      = 1/12,
     y      = max(per_residuos$spec) * 0.9,
     labels = "f = 1/12",
     col    = "firebrick",
     cex    = 0.85,
     pos    = 4)

par(mfrow = c(1, 1))

# ============================================================
# BLOQUE 6: RESUMEN FINAL
# ============================================================

cat("\n=== RESUMEN DEL MODELO DUMMY ESTACIONAL ===\n")

## 
## === RESUMEN DEL MODELO DUMMY ESTACIONAL ===

cat("R²:          ", round(summary(fit_dummy)$r.squared,     4), "\n")

## R²:           0.4246

cat("R² ajustado: ", round(summary(fit_dummy)$adj.r.squared, 4), "\n")

## R² ajustado:  0.4005

cat("Frecuencia principal:", round(frec_dom, 6), "\n")

## Frecuencia principal: 0.166667

cat("Periodo:",              round(1 / frec_dom, 2), "meses\n")

## Periodo: 6 meses

coef_tabla  <- summary(fit_dummy)$coefficients
dummies_sig <- coef_tabla[coef_tabla[, 4] < 0.05, ]

cat("\nDummies significativas al 5%:\n")

## 
## Dummies significativas al 5%:

print(round(dummies_sig, 4))

##     Estimate Std. Error t value Pr(>|t|)
## Jan   2.2710     0.2203 10.3085   0.0000
## Feb   1.2685     0.2203  5.7578   0.0000
## Apr   0.4337     0.2203  1.9686   0.0500
## Jul   0.4705     0.2203  2.1355   0.0336
## Oct  -0.9065     0.2203 -4.1149   0.0001
## Nov  -1.2106     0.2203 -5.4949   0.0000
## Dec  -0.7522     0.2249 -3.3454   0.0009

# ============================================================
# PREPARACIÓN: Crear la serie sin tendencia usando STL
# Serie sin tendencia = Serie original − Tendencia STL
# ============================================================

# --- Verificar que los objetos necesarios existen ---
# tendencia_STL  → extraída en el paso anterior
# fechas_STL     → fechas correspondientes a la tendencia STL
# valores        → valores originales de Desempleo_xts

# --- Crear serie sin tendencia ---
# Nota: tendencia_STL tiene la misma longitud que valores (sin NA)
Desempleo_sintend <- valores - tendencia_STL

# --- Crear tibble de la serie sin tendencia ---
tbl_sintend <- tibble(
  Fecha     = fechas_STL,
  Desempleo = Desempleo_sintend
)

# --- Agregar variables de calendario necesarias para Fourier y GAM ---
tbl_sintend <- tbl_sintend %>%
  mutate(
    Mes      = month(Fecha),              # Mes (1 a 12)
    Anio     = year(Fecha),               # Año
    t        = as.numeric(Fecha),         # Índice numérico de fecha
    t_scaled = scale(t)[, 1]             # Índice estandarizado
  )

# --- Verificar ---
str(tbl_sintend)

## tibble [299 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Fecha    : Date[1:299], format: "2001-01-01" "2001-02-01" ...
##  $ Desempleo: num [1:299] 2.104 2.857 1.176 -0.178 -0.714 ...
##  $ Mes      : num [1:299] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Anio     : num [1:299] 2001 2001 2001 2001 2001 ...
##  $ t        : num [1:299] 11323 11354 11382 11413 11443 ...
##  $ t_scaled : num [1:299] -1.72 -1.71 -1.7 -1.69 -1.68 ...

head(tbl_sintend, 12)

## # A tibble: 12 × 6
##    Fecha      Desempleo   Mes  Anio     t t_scaled
##    <date>         <dbl> <dbl> <dbl> <dbl>    <dbl>
##  1 2001-01-01  2.10         1  2001 11323    -1.72
##  2 2001-02-01  2.86         2  2001 11354    -1.71
##  3 2001-03-01  1.18         3  2001 11382    -1.70
##  4 2001-04-01 -0.178        4  2001 11413    -1.69
##  5 2001-05-01 -0.714        5  2001 11443    -1.68
##  6 2001-06-01  0.508        6  2001 11474    -1.67
##  7 2001-07-01  0.325        7  2001 11504    -1.65
##  8 2001-08-01  0.000777     8  2001 11535    -1.64
##  9 2001-09-01 -0.892        9  2001 11566    -1.63
## 10 2001-10-01 -0.490       10  2001 11596    -1.62
## 11 2001-11-01 -1.38        11  2001 11627    -1.61
## 12 2001-12-01 -1.51        12  2001 11657    -1.60

# --- Gráfico de la serie sin tendencia ---
ggplot(tbl_sintend, aes(x = Fecha, y = Desempleo)) +
  geom_line(color = "gray40", linewidth = 0.8) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(
    title    = "Serie sin tendencia (STL) — Tasa de Desempleo",
    subtitle = "Lo que queda son ciclos + estacionalidad + ruido",
    y        = "Desviación respecto a la tendencia (%)",
    x        = "Fecha"
  ) +
  theme_minimal()

# ============================================================
# MÉTODO 34: MODELAMIENTO ESTACIONAL CON SERIES DE FOURIER
# Función: lm() con fourier_vec() de timetk
# Captura la estacionalidad como combinación de senos y cosenos
# K = número de armónicos (mayor K = más detalle estacional)
# ============================================================

# --- Crear variables de Fourier con diferentes armónicos ---
# Para serie mensual el período es 12

tbl_fourier <- tbl_sintend %>%
  mutate(
    # Armónico 1 (K=1): captura el ciclo anual principal
    C1 = fourier_vec(Fecha, period = 12, K = 1, type = "cos"),
    S1 = fourier_vec(Fecha, period = 12, K = 1, type = "sin"),
    # Armónico 2 (K=2): captura variaciones semestrales
    C2 = fourier_vec(Fecha, period = 12, K = 2, type = "cos"),
    S2 = fourier_vec(Fecha, period = 12, K = 2, type = "sin"),
    # Armónico 3 (K=3): captura variaciones trimestrales
    C3 = fourier_vec(Fecha, period = 12, K = 3, type = "cos"),
    S3 = fourier_vec(Fecha, period = 12, K = 3, type = "sin")
  )

# --- Modelo con 1 armónico (K=1) ---
# Captura solo el ciclo anual
modelo_fourier_K1 <- lm(Desempleo ~ C1 + S1,
                        data = tbl_fourier)
summary(modelo_fourier_K1)

## 
## Call:
## lm(formula = Desempleo ~ C1 + S1, data = tbl_fourier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3997 -0.7507 -0.2358  0.4110  8.4778 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.1361     0.0778   1.749   0.0813 .  
## C1           -0.4739     0.1100  -4.308 2.24e-05 ***
## S1            0.4519     0.1100   4.108 5.18e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.345 on 296 degrees of freedom
## Multiple R-squared:  0.1068, Adjusted R-squared:  0.1008 
## F-statistic:  17.7 on 2 and 296 DF,  p-value: 5.489e-08

# --- Modelo con 2 armónicos (K=2) ---
# Captura ciclo anual + variaciones semestrales
modelo_fourier_K2 <- lm(Desempleo ~ C1 + S1 + C2 + S2,
                        data = tbl_fourier)
summary(modelo_fourier_K2)

## 
## Call:
## lm(formula = Desempleo ~ C1 + S1 + C2 + S2, data = tbl_fourier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4985 -0.7626 -0.2494  0.4031  8.4758 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.13612    0.07795   1.746   0.0818 .  
## C1          -0.47551    0.11025  -4.313 2.20e-05 ***
## S1           0.45236    0.11024   4.103 5.28e-05 ***
## C2          -0.08488    0.11011  -0.771   0.4414    
## S2          -0.04934    0.11037  -0.447   0.6551    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.348 on 294 degrees of freedom
## Multiple R-squared:  0.1092, Adjusted R-squared:  0.0971 
## F-statistic: 9.012 on 4 and 294 DF,  p-value: 7.051e-07

# --- Modelo con 3 armónicos (K=3) ---
# Captura ciclo anual + semestral + trimestral
modelo_fourier_K3 <- lm(Desempleo ~ C1 + S1 + C2 + S2 + C3 + S3,
                        data = tbl_fourier)
summary(modelo_fourier_K3)

## 
## Call:
## lm(formula = Desempleo ~ C1 + S1 + C2 + S2 + C3 + S3, data = tbl_fourier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4532 -0.7636 -0.2731  0.4099  8.3442 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.13600    0.07801   1.743   0.0823 .  
## C1          -0.47572    0.11032  -4.312 2.22e-05 ***
## S1           0.45274    0.11032   4.104 5.27e-05 ***
## C2          -0.08371    0.11019  -0.760   0.4481    
## S2          -0.04818    0.11045  -0.436   0.6630    
## C3           0.04643    0.11025   0.421   0.6739    
## S3           0.13184    0.11038   1.194   0.2333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.349 on 292 degrees of freedom
## Multiple R-squared:  0.1141, Adjusted R-squared:  0.09588 
## F-statistic: 6.267 on 6 and 292 DF,  p-value: 3.284e-06

# --- Comparar los tres modelos por R² ajustado ---
cat("R² ajustado K=1:", round(summary(modelo_fourier_K1)$adj.r.squared, 4), "\n")

## R² ajustado K=1: 0.1008

cat("R² ajustado K=2:", round(summary(modelo_fourier_K2)$adj.r.squared, 4), "\n")

## R² ajustado K=2: 0.0971

cat("R² ajustado K=3:", round(summary(modelo_fourier_K3)$adj.r.squared, 4), "\n")

## R² ajustado K=3: 0.0959

# --- Gráfico comparativo de los tres modelos Fourier ---
tbl_fourier <- tbl_fourier %>%
  mutate(
    ajuste_K1 = fitted(modelo_fourier_K1),
    ajuste_K2 = fitted(modelo_fourier_K2),
    ajuste_K3 = fitted(modelo_fourier_K3)
  )

ggplot(tbl_fourier, aes(x = Fecha)) +
  geom_line(aes(y = Desempleo),  color = "gray50",    linewidth = 0.7,
            alpha = 0.8) +
  geom_line(aes(y = ajuste_K1), color = "red",        linewidth = 1.2) +
  geom_line(aes(y = ajuste_K2), color = "blue",       linewidth = 1.2) +
  geom_line(aes(y = ajuste_K3), color = "darkgreen",  linewidth = 1.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title    = "Modelamiento estacional con Series de Fourier",
    subtitle = "Rojo=K1 (anual), Azul=K2 (+semestral), Verde=K3 (+trimestral)",
    y        = "Desviación respecto a la tendencia (%)",
    x        = "Fecha"
  ) +
  theme_minimal()

# ============================================================
# MÉTODO 35: MODELAMIENTO ESTACIONAL CON GAM
# Función: gam() del paquete mgcv
# Usa splines cúbicos cíclicos — más flexible que Fourier
# bs="cc" = spline cúbico cíclico (ideal para estacionalidad)
# ============================================================

# --- Agregar variable de mes numérico ---
tbl_gam <- tbl_sintend %>%
  mutate(Mes_num = as.numeric(Mes))

# --- Modelo GAM con spline cíclico sobre el mes ---
# k = número de nodos (máximo 12 para datos mensuales)
# bs = "cc" = cúbico cíclico (el inicio y el fin del año se conectan)
modelo_gam <- gam(Desempleo ~ s(Mes_num, bs = "cc", k = 12),
                  data = tbl_gam)

summary(modelo_gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Desempleo ~ s(Mes_num, bs = "cc", k = 12)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.13845    0.07307   1.895   0.0591 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df     F p-value    
## s(Mes_num) 7.712     10 7.955  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.206   Deviance explained = 22.7%
## GCV = 1.6445  Scale est. = 1.5965    n = 299

# --- Ver la curva del efecto estacional estimado ---
plot(modelo_gam,
     main   = "Efecto estacional estimado por GAM",
     xlab   = "Mes",
     ylab   = "Efecto sobre el desempleo (%)",
     shade  = TRUE,        # Muestra banda de confianza
     col    = "blue",
     lwd    = 2)
abline(h = 0, lty = 2, col = "red")

# --- Extraer valores ajustados del GAM ---
tbl_gam <- tbl_gam %>%
  mutate(ajuste_gam = fitted(modelo_gam))

# --- Gráfico GAM sobre la serie sin tendencia ---
ggplot(tbl_gam, aes(x = Fecha)) +
  geom_line(aes(y = Desempleo),  color = "gray50",  linewidth = 0.7) +
  geom_line(aes(y = ajuste_gam), color = "purple",  linewidth = 1.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title    = "Modelamiento estacional con GAM (spline cíclico)",
    subtitle = "Púrpura = componente estacional estimada por GAM",
    y        = "Desviación respecto a la tendencia (%)",
    x        = "Fecha"
  ) +
  theme_minimal()

# ============================================================
# GRÁFICO COMPARATIVO FINAL: FOURIER vs GAM
# ============================================================

# Unir los ajustes de Fourier K3 y GAM en un solo tibble
tbl_comparacion <- tbl_sintend %>%
  mutate(
    ajuste_fourier = fitted(modelo_fourier_K3),
    ajuste_gam     = fitted(modelo_gam)
  )

ggplot(tbl_comparacion, aes(x = Fecha)) +
  geom_line(aes(y = Desempleo),
            color = "gray60", linewidth = 0.7, alpha = 0.8) +
  geom_line(aes(y = ajuste_fourier),
            color = "darkgreen", linewidth = 1.3) +
  geom_line(aes(y = ajuste_gam),
            color = "purple",    linewidth = 1.3,
            linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +
  labs(
    title    = "Comparación: Fourier (K=3) vs GAM — Estacionalidad",
    subtitle = "Verde sólido = Fourier K3 | Púrpura punteado = GAM",
    y        = "Desviación respecto a la tendencia (%)",
    x        = "Fecha"
  ) +
  theme_minimal()

# ============================================================
# TABLA RESUMEN DE BONDAD DE AJUSTE
# ============================================================

cat("=== BONDAD DE AJUSTE — MODELOS ESTACIONALES ===\n\n")

## === BONDAD DE AJUSTE — MODELOS ESTACIONALES ===

cat("Fourier K=1  R² ajustado:", 
    round(summary(modelo_fourier_K1)$adj.r.squared, 4), "\n")

## Fourier K=1  R² ajustado: 0.1008

cat("Fourier K=2  R² ajustado:", 
    round(summary(modelo_fourier_K2)$adj.r.squared, 4), "\n")

## Fourier K=2  R² ajustado: 0.0971

cat("Fourier K=3  R² ajustado:", 
    round(summary(modelo_fourier_K3)$adj.r.squared, 4), "\n")

## Fourier K=3  R² ajustado: 0.0959

cat("GAM          R² ajustado:", 
    round(summary(modelo_gam)$r.sq,                 4), "\n")

## GAM          R² ajustado: 0.2064

# ============================================================
# MÉTODO 36: GRÁFICOS DE DENSIDAD POR MES
# Serie: Tasa de Desempleo SIN TENDENCIA (STL)
# Objetivo: confirmar presencia de estacionalidad
# ============================================================

# ============================================================
# PREPARACIÓN: Crear data frame con variable de mes
# ============================================================

df_sintend <- data.frame(
  Fecha     = fechas_STL,
  Desempleo = Desempleo_sintend,
  Mes       = factor(month.abb[month(fechas_STL)],
                     levels = month.abb)
)

# Verificar
head(df_sintend, 12)

##         Fecha     Desempleo Mes
## 1  2001-01-01  2.1042055689 Jan
## 2  2001-02-01  2.8570385644 Feb
## 3  2001-03-01  1.1757190264 Mar
## 4  2001-04-01 -0.1778309758 Apr
## 5  2001-05-01 -0.7137696451 May
## 6  2001-06-01  0.5079885642 Jun
## 7  2001-07-01  0.3247080855 Jul
## 8  2001-08-01  0.0007774362 Aug
## 9  2001-09-01 -0.8922146811 Sep
## 10 2001-10-01 -0.4897514609 Oct
## 11 2001-11-01 -1.3807930091 Nov
## 12 2001-12-01 -1.5081431878 Dec

str(df_sintend)

## 'data.frame':    299 obs. of  3 variables:
##  $ Fecha    : Date, format: "2001-01-01" "2001-02-01" ...
##  $ Desempleo: num  2.104 2.857 1.176 -0.178 -0.714 ...
##  $ Mes      : Factor w/ 12 levels "Jan","Feb","Mar",..: 1 2 3 4 5 6 7 8 9 10 ...

# ============================================================
# GRÁFICO 1: Densidad por mes en paneles separados
# Replica exactamente lo que hizo el profesor con USgas
# ============================================================

ggplot(df_sintend, aes(x = Desempleo)) +
  geom_density(aes(fill = Mes), alpha = 0.6) +
  facet_grid(rows = vars(Mes)) +
  geom_vline(xintercept = 0,
             linetype   = "dashed",
             color      = "black",
             linewidth  = 0.5) +
  labs(
    title    = "Densidad por mes — Serie sin tendencia (STL)",
    subtitle = "Distribuciones distintas entre meses indican estacionalidad",
    x        = "Desviación respecto a la tendencia (%)",
    y        = "Densidad"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# ============================================================
# GRÁFICO 2: Densidad superpuesta de todos los meses
# Permite ver de un vistazo si los meses se separan entre sí
# ============================================================

ggplot(df_sintend, aes(x = Desempleo,
                       fill  = Mes,
                       color = Mes)) +
  geom_density(alpha = 0.20, linewidth = 0.8) +
  geom_vline(xintercept = 0,
             linetype   = "dashed",
             color      = "black",
             linewidth  = 0.8) +
  labs(
    title    = "Densidad superpuesta por mes — Serie sin tendencia",
    subtitle = "Curvas separadas = estacionalidad fuerte | Curvas mezcladas = sin estacionalidad",
    x        = "Desviación respecto a la tendencia (%)",
    y        = "Densidad",
    fill     = "Mes",
    color    = "Mes"
  ) +
  theme_minimal()

# ============================================================
# GRÁFICO 3: Boxplot por mes
# Complementa las densidades mostrando mediana y dispersión
# ============================================================

ggplot(df_sintend, aes(x = Mes, y = Desempleo, fill = Mes)) +
  geom_boxplot(alpha         = 0.7,
               outlier.color = "red",
               outlier.shape = 16,
               outlier.size  = 2) +
  geom_hline(yintercept = 0,
             linetype   = "dashed",
             color      = "black",
             linewidth  = 0.8) +
  labs(
    title    = "Boxplot por mes — Serie sin tendencia (STL)",
    subtitle = "Medianas distintas entre meses confirman estacionalidad",
    x        = "Mes",
    y        = "Desviación respecto a la tendencia (%)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# ============================================================
# GRÁFICO 4: Medias por mes con barras de error
# Si las medias varían sistemáticamente → estacionalidad presente
# ============================================================

df_medias <- df_sintend %>%
  group_by(Mes) %>%
  summarise(
    Media = mean(Desempleo, na.rm = TRUE),
    DE    = sd(Desempleo,   na.rm = TRUE),
    .groups = "drop"
  )

ggplot(df_medias, aes(x = Mes, y = Media, fill = Mes)) +
  geom_bar(stat  = "identity", alpha = 0.8) +
  geom_errorbar(aes(ymin = Media - DE,
                    ymax = Media + DE),
                width     = 0.3,
                linewidth = 0.8) +
  geom_hline(yintercept = 0,
             linetype   = "dashed",
             color      = "black",
             linewidth  = 0.8) +
  labs(
    title    = "Media mensual — Serie sin tendencia (STL)",
    subtitle = "Barras de error = ±1 desviación estándar",
    x        = "Mes",
    y        = "Desviación promedio respecto a la tendencia (%)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# --- Tabla de medias por mes ---
cat("=== MEDIAS POR MES — SERIE SIN TENDENCIA ===\n\n")

## === MEDIAS POR MES — SERIE SIN TENDENCIA ===

print(df_medias)

## # A tibble: 12 × 3
##    Mes     Media    DE
##    <fct>   <dbl> <dbl>
##  1 Jan    2.27   0.558
##  2 Feb    1.27   0.462
##  3 Mar    0.169  0.665
##  4 Apr    0.434  1.70 
##  5 May    0.290  1.95 
##  6 Jun    0.0325 1.62 
##  7 Jul    0.470  1.54 
##  8 Aug   -0.0449 0.878
##  9 Sep   -0.395  0.698
## 10 Oct   -0.907  0.501
## 11 Nov   -1.21   0.358
## 12 Dec   -0.752  0.464

# ============================================================
# MÉTODO 37: GRÁFICOS DE TEMPORADA Y ESTACIONALES
# Replica chunks "mas acerca estacionalidad_1" y "_2"
# Usa tanto la serie original como la serie sin tendencia
# ============================================================

# ============================================================
# PREPARACIÓN: Convertir serie sin tendencia a objeto ts
# Necesario para ggseasonplot() y ts_seasonal()
# ============================================================

# --- Serie sin tendencia como objeto ts ---
Desempleo_sintend_ts <- ts(Desempleo_sintend,
                           start     = c(2001, 1),
                           frequency = 12)

# --- Tsibble de la serie sin tendencia ---
# Necesario para gg_season() y gg_subseries() de feasts
tsbl_sintend <- tsibble(
  Fecha     = yearmonth(fechas_STL),
  Desempleo = Desempleo_sintend,
  index     = Fecha
)

# Verificar
ts_info(Desempleo_sintend_ts)

##  The Desempleo_sintend_ts series is a ts object with 1 variable and 299 observations
##  Frequency: 12 
##  Start time: 2001 1 
##  End time: 2025 11

tsbl_sintend

## # A tsibble: 299 x 2 [1M]
##        Fecha Desempleo
##        <mth>     <dbl>
##  1 2001 ene.  2.10    
##  2 2001 feb.  2.86    
##  3 2001 mar.  1.18    
##  4 2001 abr. -0.178   
##  5 2001 may. -0.714   
##  6 2001 jun.  0.508   
##  7 2001 jul.  0.325   
##  8 2001 ago.  0.000777
##  9 2001 sep. -0.892   
## 10 2001 oct. -0.490   
## # ℹ 289 more rows

# ============================================================
# BLOQUE A: GRÁFICOS ESTACIONALES — SIN forecast
# Usando feasts y ggplot2 que ya están cargados
# ============================================================

# --- A1: Gráfico estacional — Serie ORIGINAL ---
# Replica ggseasonplot usando gg_season() de feasts
tsbl_desempleo %>%
  gg_season(Desempleo, period = "year") +
  labs(
    title    = "Gráfico estacional — Serie ORIGINAL",
    subtitle = "Cada línea de color representa un año distinto",
    x        = "Mes",
    y        = "Tasa de Desempleo (%)"
  ) +
  theme_minimal()

# --- A2: Gráfico estacional — Serie SIN TENDENCIA ---
tsbl_sintend %>%
  gg_season(Desempleo, period = "year") +
  labs(
    title    = "Gráfico estacional — Serie SIN TENDENCIA (STL)",
    subtitle = "Si las líneas se superponen → patrón estacional estable en el tiempo",
    x        = "Mes",
    y        = "Desviación respecto a la tendencia (%)"
  ) +
  theme_minimal()

# --- A4: Subseries por mes (equivalente al polar de forecast) ---
tsbl_sintend %>%
  gg_subseries(Desempleo, period = "year") +
  labs(
    title    = "Subseries por mes — Serie sin tendencia",
    subtitle = "Línea azul = media del mes a lo largo de todos los años",
    x        = "Año",
    y        = "Desviación respecto a la tendencia (%)"
  ) +
  theme_minimal()

# ============================================================
# BLOQUE B: GRÁFICOS CON TSstudio::ts_seasonal()
# Replica chunk "mas acerca estacionalidad_2"
# Cuatro vistas del patrón estacional
# ============================================================

# --- B1: Vista normal — líneas por año ---
ts_seasonal(Desempleo_sintend_ts, type = "normal")

# --- B2: Vista de ciclo — orden cronológico por mes ---
# Permite ver si el patrón cambia a lo largo de los años
ts_seasonal(Desempleo_sintend_ts, type = "cycle")

# --- B3: Vista de caja — boxplot por mes ---
ts_seasonal(Desempleo_sintend_ts, type = "box")

# --- B4: Todas las vistas juntas ---
ts_seasonal(Desempleo_sintend_ts, type = "all")

# ============================================================
# BLOQUE C: GRÁFICOS CON feasts
# gg_season() y gg_subseries()
# ============================================================

# --- C1: gg_season — Superpone años en el mismo gráfico ---
tsbl_sintend %>%
  gg_season(Desempleo, period = "year") +
  labs(
    title    = "Gráfico de temporada — Serie sin tendencia (feasts)",
    subtitle = "Cada línea de color representa un año distinto",
    x        = "Mes",
    y        = "Desviación respecto a la tendencia (%)"
  ) +
  theme_minimal()

# --- C2: gg_subseries — Subseries por mes a través del tiempo ---
# Muestra la evolución de cada mes a lo largo de todos los años
# La línea horizontal azul es la media del mes
tsbl_sintend %>%
  gg_subseries(Desempleo, period = "year") +
  labs(
    title    = "Subseries por mes — Serie sin tendencia (feasts)",
    subtitle = "Línea azul = media del mes | Cambios en la línea = estacionalidad cambiante",
    x        = "Año",
    y        = "Desviación respecto a la tendencia (%)"
  ) +
  theme_minimal()

# ============================================================
# BLOQUE D: COMPARACIÓN SERIE ORIGINAL vs SIN TENDENCIA
# Para ver el impacto de remover la tendencia en la lectura
# del patrón estacional
# ============================================================

par(mfrow = c(1, 2),
    mar   = c(4, 4, 3, 2))

# Serie original
monthplot(Desempleo_ts,
          main = "Monthplot — Serie ORIGINAL",
          ylab = "Desempleo (%)",
          xlab = "Mes")

# Serie sin tendencia
monthplot(Desempleo_sintend_ts,
          main = "Monthplot — Serie SIN TENDENCIA",
          ylab = "Desviación (%)",
          xlab = "Mes")

par(mfrow = c(1, 1))

# ============================================================
# MÉTODOS 39 Y 40: WAVELET Y SUAVIZAMIENTO SPLINES
# Serie: Tasa de Desempleo SIN TENDENCIA (STL)
# ============================================================

# MÉTODO 39: ANÁLISIS WAVELET
# Función: biwavelet::wt()
# Detecta ciclos cuya frecuencia o amplitud cambian en el tiempo
# Es un análisis TIEMPO-FRECUENCIA simultáneo:
# → el eje X muestra el tiempo
# → el eje Y muestra el período (en meses)
# → el color muestra la intensidad del ciclo
# ============================================================

# --- Preparar el input para biwavelet ---
# Requiere una matriz de 2 columnas: [índice numérico, valores]
indice_wavelet <- 1:length(Desempleo_sintend)

serie_wavelet <- cbind(
  indice_wavelet,
  as.numeric(Desempleo_sintend)
)

# Verificar
head(serie_wavelet)

##      indice_wavelet           
## [1,]              1  2.1042056
## [2,]              2  2.8570386
## [3,]              3  1.1757190
## [4,]              4 -0.1778310
## [5,]              5 -0.7137696
## [6,]              6  0.5079886

dim(serie_wavelet)

## [1] 299   2

# --- Calcular el wavelet ---
# Puede tardar unos segundos
wt_resultado <- biwavelet::wt(serie_wavelet)

# --- Gráfico 1: Potencia corregida por sesgo (recomendada) ---
par(oma = c(0, 0, 0, 1),
    mar = c(5, 4, 4, 5) + 0.1)

plot(wt_resultado,
     type    = "power.corr.norm",
     main    = "Wavelet — Potencia corregida (serie sin tendencia)",
     plot.cb = TRUE)         # Muestra barra de color

# --- Gráfico 2: Potencia con sesgo ---
plot(wt_resultado,
     type    = "power.norm",
     main    = "Wavelet — Potencia con sesgo (serie sin tendencia)",
     plot.cb = TRUE)

# --- Gráfico 3: Los dos juntos para comparar ---
par(mfrow = c(1, 2),
    oma   = c(0, 0, 0, 1),
    mar   = c(5, 4, 4, 5) + 0.1)

plot(wt_resultado,
     type    = "power.corr.norm",
     main    = "Potencia corregida",
     plot.cb = TRUE)

plot(wt_resultado,
     type    = "power.norm",
     main    = "Potencia con sesgo",
     plot.cb = TRUE)

par(mfrow = c(1, 1))

# --- Cómo leer los gráficos wavelet ---
# Eje X   → tiempo (meses desde enero 2001)
# Eje Y   → período en meses (escala logarítmica)
#            período=12 → ciclo anual
#            período=6  → ciclo semestral
#            período=60 → ciclo de 5 años
# Color   → intensidad del ciclo en ese momento
#            rojo = ciclo muy fuerte
#            azul = ciclo débil o ausente
# Cono    → región fuera del cono = resultados menos confiables
#            (efecto de borde por los extremos de la serie)


# ============================================================
# MÉTODO 40: SUAVIZAMIENTO SPLINES
# Función: smooth.spline() de R base
# Minimiza el balance entre ajuste y suavidad
# spar controla la suavidad:
#      spar cercano a 0 → muy detallado, sigue bien la serie
#      spar cercano a 1 → muy suave, solo muestra tendencia residual
# ============================================================

# --- Preparar índice numérico ---
tiempo_spline <- as.numeric(fechas_STL)

# --- Ajustar splines con diferentes niveles de suavidad ---
spline_05  <- smooth.spline(tiempo_spline, Desempleo_sintend, spar = 0.5)
spline_07  <- smooth.spline(tiempo_spline, Desempleo_sintend, spar = 0.7)
spline_10  <- smooth.spline(tiempo_spline, Desempleo_sintend, spar = 1.0)

# --- Gráfico 1: Comparación de tres niveles de suavidad ---
plot(fechas_STL, Desempleo_sintend,
     type = "l",
     main = "Suavizamiento Splines — Serie sin tendencia",
     ylab = "Desviación respecto a la tendencia (%)",
     xlab = "Fecha",
     col  = "gray60",
     lwd  = 1)

lines(as.Date(spline_05$x), spline_05$y,
      col = "red",       lwd = 2)
lines(as.Date(spline_07$x), spline_07$y,
      col = "blue",      lwd = 2)
lines(as.Date(spline_10$x), spline_10$y,
      col = "darkgreen", lwd = 2)

abline(h   = 0,
       lty = 2,
       col = "black")

legend("topright",
       legend = c("Serie sin tendencia",
                  "Spline spar=0.5 (detallado)",
                  "Spline spar=0.7 (medio)",
                  "Spline spar=1.0 (suave)"),
       col    = c("gray60", "red", "blue", "darkgreen"),
       lwd    = c(1, 2, 2, 2),
       cex    = 0.85)

# --- Gráfico 2: Tres paneles separados para ver mejor cada nivel ---
par(mfrow = c(3, 1),
    mar   = c(3, 4, 3, 2))

# Panel 1: spar = 0.5 (detallado)
plot(fechas_STL, Desempleo_sintend,
     type = "l", col = "gray60", lwd = 1,
     main = "Spline spar=0.5 (detallado — sigue bien los ciclos)",
     ylab = "Desviación (%)")
lines(as.Date(spline_05$x), spline_05$y,
      col = "red", lwd = 2)
abline(h = 0, lty = 2)

# Panel 2: spar = 0.7 (medio)
plot(fechas_STL, Desempleo_sintend,
     type = "l", col = "gray60", lwd = 1,
     main = "Spline spar=0.7 (medio — balance entre detalle y suavidad)",
     ylab = "Desviación (%)")
lines(as.Date(spline_07$x), spline_07$y,
      col = "blue", lwd = 2)
abline(h = 0, lty = 2)

# Panel 3: spar = 1.0 (suave)
plot(fechas_STL, Desempleo_sintend,
     type = "l", col = "gray60", lwd = 1,
     main = "Spline spar=1.0 (suave — solo estructura general)",
     ylab = "Desviación (%)")
lines(as.Date(spline_10$x), spline_10$y,
      col = "darkgreen", lwd = 2)
abline(h = 0, lty = 2)

par(mfrow = c(1, 1))

# --- Gráfico 3: versión ggplot2 del spline elegido ---
# Usamos spar=0.7 como balance recomendado
df_spline <- data.frame(
  Fecha       = fechas_STL,
  Desempleo   = as.numeric(Desempleo_sintend),
  Spline_07   = spline_07$y
)

ggplot(df_spline, aes(x = Fecha)) +
  geom_line(aes(y = Desempleo),
            color = "gray60", linewidth = 0.7) +
  geom_line(aes(y = Spline_07),
            color = "blue",   linewidth = 1.5) +
  geom_hline(yintercept = 0,
             linetype = "dashed", color = "black") +
  labs(
    title    = "Suavizamiento Spline (spar=0.7) — Serie sin tendencia",
    subtitle = "Azul = curva spline | Captura la estructura residual tras remover tendencia STL",
    x        = "Fecha",
    y        = "Desviación respecto a la tendencia (%)"
  ) +
  theme_minimal()

# ============================================================
# SERIE SIN TENDENCIA Y SIN ESTACIONALIDAD
# Residuo puro = Serie original − Tendencia STL − Estacionalidad STL
# ============================================================

# ── Extraer los tres componentes del objeto fable ──────────────
tendencia_STL      <- descomp_STL$trend
estacionalidad_STL <- descomp_STL$season_year   # así se llama en fable, NO "seasonal"
residuo_STL        <- descomp_STL$remainder      # así se llama en fable, NO "remainder" de stl()
fechas_STL         <- as.Date(descomp_STL$Fecha)

# ── Verificación cruzada ────────────────────────────────────────
residuo_manual <- valores - tendencia_STL - estacionalidad_STL

# ── Data frame para ggplot ──────────────────────────────────────
df_residuo <- data.frame(
  Fecha   = fechas_STL,
  Residuo = as.numeric(residuo_STL)
)

# ── Gráfico ────────────────────────────────────────────────────
ggplot(df_residuo, aes(x = Fecha, y = Residuo)) +
  geom_line(color = "gray30", linewidth = 0.7) +
  geom_hline(yintercept = 0,
             linetype   = "dashed",
             color      = "red",
             linewidth  = 0.8) +
  geom_smooth(method = "loess", span = 0.15,
              color  = "steelblue",
              se     = FALSE,
              linewidth = 1.2) +
  labs(
    title    = "Residuo STL — Serie sin tendencia y sin estacionalidad",
    subtitle = "Lo que queda es ruido + ciclos irregulares no capturados por STL",
    x        = "Fecha",
    y        = "Residuo (%)"
  ) +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'