1 Carga y Exploración Inicial de los Datos

¿Qué se hace aquí?
Se importa la serie de tiempo de la Tasa de Desempleo mensual de Colombia desde un archivo Excel. Luego se convierte la columna de fechas al formato Date (quitando la parte de hora que Excel agrega), y se construye un objeto xts (eXtensible Time Series), que es el formato base de trabajo para los análisis posteriores. El objeto xts permite indexar la serie por fechas tipo yearmon, lo que facilita las operaciones de ventanas temporales, subperiodos y gráficos interactivos.

# Inspeccionar la estructura del dataframe importado
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"
# Convertir Fecha de POSIXct a Date
# R lee las fechas de Excel como POSIXct (con hora). Se usa as.Date() para
# quedarse solo con la parte de la fecha (año-mes-día), sin la hora.
Tasa_Desempleo_raw <- Tasa_Desempleo_raw %>%
  mutate(Fecha = as.Date(Fecha))

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"
# Crear objeto XTS (eXtensible Time Series)
# xts es el formato más flexible para series de tiempo en R:
# permite indexar por fechas, hacer subconjuntos con sintaxis de fechas,
# y es compatible con la mayoría de paquetes de análisis temporal.
Fechas_yearmon <- as.yearmon(Tasa_Desempleo_raw$Fecha)

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

# 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"
# 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)

2 Construcción de Formatos de Objeto

¿Por qué se necesitan múltiples formatos?
Diferentes paquetes de R exigen diferentes formatos de entrada. En este análisis se usan tres formatos simultáneamente:

  • xts: formato base, compatible con TSstudio y funciones del ecosistema zoo.
  • tibble: formato rectangular del tidyverse, necesario para timetk y ggplot2.
  • tsibble: formato de series de tiempo moderno del ecosistema tidyverts (feasts, fable), que permite aplicar modelos STL, gráficos estacionales y subseries de forma directa.

3 Visualización Inicial: Tendencia, Ciclo y Estacionalidad con LOESS

¿Qué se hace aquí y por qué?
Antes de aplicar cualquier método formal de descomposición, se realiza una exploración visual de la tendencia usando suavizamiento LOESS (Locally Estimated Scatterplot Smoothing). El parámetro span controla qué fracción de los datos se usa como vecindad en cada punto:

  • span = 0.10 (bajo): la curva sigue de cerca los ciclos y la estacionalidad; útil para detectar variaciones rápidas.
  • span = 0.25 (medio): balance entre detalle y suavidad
  • span = 0.75 (alto): solo muestra la tendencia de largo plazo; ignora ciclos y ruido.

Esta comparación visual ayuda a decidir qué nivel de suavidad es apropiado antes de elegir un método formal.

p1 <- tbl_desempleo %>%
  plot_time_series(
    .date_var    = Fecha,
    .value       = Desempleo,
    .title       = "LOESS span = 0.10 (detallado — sigue ciclos y estacionalidad)",
    .interactive = FALSE,
    .smooth      = TRUE,
    .smooth_span = 0.10
  )

p2 <- tbl_desempleo %>%
  plot_time_series(
    .date_var    = Fecha,
    .value       = Desempleo,
    .title       = "LOESS span = 0.25 (medio — balance detalle/suavidad)",
    .interactive = FALSE,
    .smooth      = TRUE,
    .smooth_span = 0.25
  )

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

p1 / p2 / p3

4 Estimación Formal de la Tendencia: Cinco Métodos

¿Qué se hace aquí?
Se aplican cinco métodos distintos para estimar la componente de tendencia de la serie. El objetivo es comparar su comportamiento antes de elegir uno para remover la tendencia. Al final se superponen los cinco sobre la serie original.

Resumen de los cinco métodos:

Método Función Característica principal
1. Descomposición clásica decompose() Promedio móvil simple; sensible a atípicos; produce NA en los bordes
2. Filtro Boxcar stats::filter() Promedio Movil con pesos diferenciados en los extremos; variante más refinada del M1
3. Kernel Nadaraya-Watson ksmooth() Suaviza con kernel gaussiano; el bandwidth controla la escala temporal
4. LOWESS lowess() Regresión local ponderada; f controla la vecindad
5. STL (robusto) fable::STL() Descomposición basada en Loess; robust=TRUE reduce el efecto del COVID

¿Por qué se escogió STL para remover la tendencia?
La serie presenta un pico extremo en 2020 por el COVID-19. Los métodos 1 a 4 son sensibles a valores atípicos: el pico distorsiona la tendencia estimada en los meses cercanos. El método STL con robust=TRUE asigna menor peso a las observaciones atípicas durante la estimación, logrando una tendencia que refleja mejor el comportamiento estructural de largo plazo sin que el COVID “contamine” los períodos adyacentes.

# Extraer valores y fechas base para los métodos no-STL
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 stats::filter)
Desempleo_ts <- ts(valores,
                   start     = c(2001, 1),
                   frequency = 12)

# ── Método 1: Descomposición clásica ───────────────────────────
# Usa un promedio móvil simple de ventana 12.
# Limitación: los primeros y últimos 6 valores serán NA porque
# necesita una ventana completa de 12 meses para operar en los bordes.
descomp_clasica     <- decompose(Desempleo_ts, type = "additive")
tendencia_decompose <- descomp_clasica$trend
plot(descomp_clasica)

# ── Método 2: Filtro Boxcar (Promedio Movil con pesos) ─────────────────────
# Los extremos de la ventana reciben peso 0.5 y el resto peso 1,
# todo dividido entre 12. Es más refinado que el PM simple porque
# reduce el salto artificial en los bordes de la ventana.
pesos_boxcar     <- c(0.5, rep(1, 11), 0.5) / 12
tendencia_boxcar <- stats::filter(Desempleo_ts,
                                  filter = pesos_boxcar,
                                  sides  = 2)

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: Kernel Nadaraya-Watson ───────────────────────────
# Estima la tendencia en cada punto como un promedio ponderado
# de los valores vecinos, donde los pesos siguen una distribución
# gaussiana. A mayor bandwidth, más suave es la curva resultante.
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)

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: LOWESS ───────────────────────────────────────────
# Regresión polinomial local ponderada. El parámetro f indica
# qué fracción de los datos se usa como vecindad en cada punto.
# f pequeño = sigue ciclos; f grande = solo tendencia de fondo.
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)

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: STL (MÉTODO SELECCIONADO) ────────────────────────
# STL = Seasonal-Trend decomposition using Loess.
# season(window="periodic") fuerza un patrón estacional constante.
# robust=TRUE: en cada iteración, las observaciones con residuos
# grandes (como el pico COVID de 2020) reciben menor peso, lo que
# evita que distorsionen la estimación de la tendencia adyacente.
descomp_STL <- tsbl_desempleo %>%
  model(
    STL(Desempleo ~ trend() + season(window = "periodic"),
        robust = TRUE)
  ) %>%
  components()

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 del objeto STL (fable)
# Nota: en fable los componentes son columnas directas del tsibble,
# no se accede con $time.series como en stl() de R base.
tendencia_STL <- descomp_STL$trend
fechas_STL    <- as.Date(descomp_STL$Fecha)

# ── Gráfico comparativo: los cinco métodos ─────────────────────
plot(Desempleo_ts,
     main = "Comparación de métodos de estimación de tendencia",
     ylab = "Desempleo (%)", xlab = "Año",
     col = "gray70", lwd = 1)
lines(tendencia_decompose,            col = "black",     lwd = 2)
lines(tendencia_boxcar,               col = "orange",    lwd = 2)
lines(tendencia_kernel_1$x,
      tendencia_kernel_1$y,           col = "red",       lwd = 2)
lines(fechas, tendencia_lowess_025$y, col = "blue",      lwd = 2)
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 — SELECCIONADO)"),
       col = c("gray70", "black", "orange", "red", "blue", "darkgreen"),
       lwd = c(1, 2, 2, 2, 2, 2), cex = 0.80)

5 Eliminación de la Tendencia con STL

¿Qué se hace aquí?
Se construye la serie sin tendencia restando, observación por observación, la tendencia estimada por STL a la serie original:

\[\text{Sin tendencia}_t = \text{Desempleo}_t - \hat{T}_t^{\text{STL}}\]

El resultado conserva la estacionalidad y el residuo, pero ya no tiene el componente de largo plazo. Esta es la serie sobre la cual se realizarán todos los análisis de estacionalidad en las secciones siguientes.

Se grafica el panel comparativo: arriba la serie original con la tendencia superpuesta, abajo la serie sin tendencia oscilando alrededor de cero.

# Construir data frame de trabajo desde el objeto STL (fable)
# descomp_STL contiene columnas: Fecha, Desempleo, trend, season_year, remainder
df_stl <- descomp_STL %>%
  as_tibble() %>%
  mutate(
    Fecha         = as.Date(Fecha),
    Sin_tendencia = Desempleo - trend      # Aquí ocurre la remoción
  )

# Las primeras filas tendrán NA heredados del STL en la columna trend
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 con tendencia STL superpuesta
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 original | Rojo = tendencia STL (método seleccionado)",
    x        = NULL,
    y        = "Tasa de desempleo (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

# Gráfico 2: Serie sin tendencia
# Esta serie oscila alrededor de cero y conserva estacionalidad + residuo
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 | Oscila alrededor de cero",
    x        = "Año",
    y        = "Desempleo − tendencia (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

p_original / p_sin_tend

6 Análisis de Autocorrelación (ACF y PACF)

¿Qué se hace aquí y por qué?
La Función de Autocorrelación (ACF) mide la correlación de la serie consigo misma a diferentes rezagos. Si hay picos significativos en los rezagos múltiplos de 12 (12, 24, 36…), eso es la evidencia estadística directa de que existe un ciclo anual en la serie sin tendencia.

La Función de Autocorrelación Parcial (PACF) mide la correlación en el rezago k eliminando el efecto de los rezagos intermedios. Complementa al ACF para identificar el orden de los procesos Auto-Regresivos que podrían modelar la dependencia temporal.

Se comparan las dos series (con y sin tendencia) para mostrar cómo la remoción de la tendencia elimina la autocorrelación de largo plazo y deja visible el patrón estacional.

# Construir Sin_tendencia_ts como objeto ts (necesario para acf/pacf de R base)
# Se eliminan los NA que el STL produce en los extremos antes de construir el ts.
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 hasta el rezago 48 (4 años)
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%: ±1.96/√n
n     <- length(Sin_tendencia_ts)
banda <- qnorm(0.975) / sqrt(n)

# Unir en un data frame para ggplot comparativo
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)")
)

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 = "Picos en rezagos 12, 24, 36 confirman ciclo anual | Líneas punteadas: banda 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))

# Calcular PACF de ambas series
# El PACF NO incluye el rezago 0 (a diferencia del ACF)
pacf_original <- pacf(Desempleo_ts,     lag.max = 48, plot = FALSE)
pacf_sin_tend <- pacf(Sin_tendencia_ts, lag.max = 48, plot = FALSE)

n     <- length(Sin_tendencia_ts)
banda <- qnorm(0.975) / sqrt(n)

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)")
)

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))

7 Gráficos de Dispersión de Retardos

¿Qué se hace aquí?
Se grafica \(y_t\) vs \(y_{t-k}\) para \(k = 1, 2, \ldots, 24\) sobre la serie sin tendencia. Cada panel corresponde a un rezago e incluye una curva loess que describe el patrón de dependencia. Esto complementa al ACF porque permite detectar dependencias no lineales que el ACF no captura.

¿Cómo se interpreta?

  • Nube con pendiente positiva → correlación positiva a ese rezago.
  • Nube con pendiente negativa → correlación negativa.
  • Nube circular / sin patrón → rezago no significativo.
  • Curva loess no lineal → dependencia no lineal.

Los rezagos más importantes para una serie mensual son: rezago 1 (memoria de corto plazo), rezago 12 (mismo mes del año anterior, confirma ciclo anual) y rezago 24 (segundo orden del ciclo anual).

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

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

8 Mapa de Calor Temporal

¿Qué se hace aquí?
El mapa de calor representa la intensidad de la serie en una matriz Año × Mes. Cada celda muestra el valor del desempleo codificado por color. Permite detectar visualmente si ciertos meses son sistemáticamente más altos o más bajos que otros a lo largo de todos los años, lo que es la definición intuitiva de estacionalidad.

Se comparan dos mapas: la serie original (donde el patrón estacional compite visualmente con la tendencia) y la serie sin tendencia (donde la estacionalidad queda expuesta en forma “pura”).

# Mapa de calor de la serie original (2001-2025)
ts_heatmap(Desempleo_ts,
           title = "Mapa de Calor — Tasa de Desempleo Mensual Colombia (2001-2025)")
# Mapa de calor de la serie SIN tendencia
# 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)")

9 Diagnóstico Estacional con Boxplots por Período

¿Qué se hace aquí?
Se usan boxplots para resumir la distribución de los valores de la serie sin tendencia agrupados por trimestre y por año. Si las medianas difieren sistemáticamente entre trimestres, eso confirma que la estacionalidad tiene una estructura trimestral. El boxplot por año ayuda a identificar si el patrón estacional ha cambiado a lo largo del tiempo (estacionalidad cambiante).

# Construir tibble desde df_stl_clean para plot_seasonal_diagnostics
tbl_sin_tendencia <- df_stl_clean %>%
  select(Fecha, Sin_tendencia) %>%
  rename(sin_tend = Sin_tendencia) %>%
  mutate(Fecha = as.Date(Fecha))

# Boxplot por trimestre
# Si las medianas de T1, T2, T3 y T4 difieren → hay estructura estacional trimestral
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)"
  )

# Boxplot por año
# Permite detectar si el patrón estacional cambia entre períodos (p.ej. antes/después COVID)
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)"
  ) +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

10 Subseries Mensuales y Trimestrales

¿Qué se hace aquí?
gg_subseries() divide la serie en paneles, uno por mes o trimestre, y muestra la evolución de ese período a lo largo de todos los años. La línea azul horizontal representa la media de ese mes/trimestre. Si las medias difieren entre paneles, la estacionalidad es evidente. Si los valores dentro de cada panel muestran una tendencia creciente o decreciente, la estacionalidad está cambiando en el tiempo.

# Construir tsibble desde df_stl_clean (feasts requiere índice yearmonth)
tsbl_sin_tendencia <- df_stl_clean %>%
  select(Fecha, Sin_tendencia) %>%
  mutate(Fecha = yearmonth(Fecha)) %>%
  as_tsibble(index = Fecha)

# Subseries mensuales (period=12 → ciclo anual)
# La línea azul = media de cada mes a lo largo de todos los años
tsbl_sin_tendencia %>%
  gg_subseries(Sin_tendencia, period = 12) +
  labs(
    title    = "Subseries mensuales — Desempleo sin tendencia (STL)",
    subtitle = "Línea azul = media de cada mes | Variación dentro de panel = cambio en el tiempo",
    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"),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

# Subseries trimestrales
# Se agrega primero a trimestre para evitar múltiples valores por período
tsbl_sin_tend_trim <- df_stl_clean %>%
  select(Fecha, Sin_tendencia) %>%
  mutate(trimestre = yearquarter(Fecha)) %>%
  group_by(trimestre) %>%
  summarise(sin_tend_trim = mean(Sin_tendencia, na.rm = TRUE), .groups = "drop") %>%
  as_tsibble(index = trimestre)

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"))

11 Regresión Trigonométrica para Estimación de Ciclos

¿Qué se hace aquí?
Se ajustan modelos de regresión donde las variables explicativas son funciones coseno y seno de diferentes períodos. El modelo es:

\[y_t = \beta_1 \cos\!\left(\frac{2\pi t}{P}\right) + \beta_2 \sin\!\left(\frac{2\pi t}{P}\right) + \varepsilon_t\]

La amplitud estimada del ciclo es \(A = \sqrt{\beta_1^2 + \beta_2^2}\) y el R² indica qué proporción de la varianza explica ese ciclo. Se evalúan tres períodos candidatos: 12 meses (ciclo anual), 6 meses (ciclo semestral) y 60 meses (ciclo económico de 5 años), además de un modelo combinado 12+6.

n_obs  <- length(Sin_tendencia_ts)
tiempo <- 1:n_obs
y      <- as.numeric(Sin_tendencia_ts)

# ── Visualización inicial ───────────────────────────────────────
# Se superpone cada onda candidata sobre la serie para evaluar
# visualmente cuál período se ajusta mejor a las oscilaciones observadas.
par(mfrow = c(3, 1), mar = c(3, 4, 3, 2))

plot.ts(y, main = "Serie sin tendencia — con ciclo de 12 meses superpuesto",
        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)

plot.ts(y, main = "Serie sin tendencia — con ciclo de 6 meses superpuesto",
        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)

plot.ts(y, main = "Serie sin tendencia — con ciclo de 60 meses superpuesto",
        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))

# ── Modelos de regresión trigonométrica ────────────────────────
# Se ajustan SIN intercepto (~ 0 +) porque la serie ya está centrada en cero.

# Ciclo anual (P=12)
z1_12 <- cos(2 * pi * tiempo / 12); z2_12 <- sin(2 * pi * tiempo / 12)
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
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 ciclo anual (12 meses):", round(A_12, 4), "\n")
## Amplitud ciclo anual (12 meses): 0.757
cat("Fase estimada (radianes):", round(atan2(-beta2_12, beta1_12), 4), "\n")
## Fase estimada (radianes): -1.6838
# Ciclo semestral (P=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
A_6 <- sqrt(coef(fit_6)["z1_6"]^2 + coef(fit_6)["z2_6"]^2)
cat("Amplitud ciclo semestral (6 meses):", round(A_6, 4), "\n")
## Amplitud ciclo semestral (6 meses): 0.7735
# Ciclo económico largo (P=60)
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
A_60 <- sqrt(coef(fit_60)["z1_60"]^2 + coef(fit_60)["z2_60"]^2)
cat("Amplitud ciclo largo (60 meses):", round(A_60, 4), "\n")
## Amplitud ciclo largo (60 meses): 0.328
# Modelo combinado (P=12 + P=6)
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
# ── Gráficos comparativos: serie vs ajuste ─────────────────────
par(mfrow = c(3, 1), mar = c(3, 4, 3, 2))
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", "Ajuste trigonométrico"),
       col = c("gray60", "firebrick"), lwd = c(1, 2), cex = 0.8)

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", "Ajuste trigonométrico"),
       col = c("gray60", "darkgreen"), lwd = c(1, 2), cex = 0.8)

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", "Ajuste combinado"),
       col = c("gray60", "darkorange"), lwd = c(1, 2), cex = 0.8)

par(mfrow = c(1, 1))

# ── Tabla resumen de amplitudes y R² ───────────────────────────
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

12 Periodograma y Periodograma Suavizado

¿Qué se hace aquí?
El periodograma es la representación de la serie en el dominio de la frecuencia. Descompone la varianza total de la serie en la contribución de cada frecuencia. Un pico en la frecuencia \(f = 1/12\) confirma que existe un ciclo anual; picos en \(f = 2/12\) y \(f = 3/12\) son sus armónicos.

Se calculan tres versiones:

  1. Periodograma crudo con spectrum(): altamente ruidoso, pero identifica frecuencias dominantes.
  2. Periodograma vía FFT: calcula directamente la transformada de Fourier.
  3. mvspec() de astsa: incluye bandas de confianza al 95%.
  4. Periodograma suavizado: el argumento span promedia frecuencias vecinas para reducir el ruido y hacer más visibles los picos reales.
par(mfrow = c(1, 1), mar = c(4, 4, 3, 2))

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

ubicacion_max <- which.max(per_crudo$spec)
frec_max      <- per_crudo$freq[ubicacion_max]
periodo_max   <- 1 / frec_max

cat("Frecuencia dominante:", round(frec_max, 6), "\n")
## Frecuencia dominante: 0.166667
cat("Periodo (meses):",      round(periodo_max, 2), "\n")
## Periodo (meses): 6
abline(v = frec_max, lty = 2, col = "firebrick", lwd = 2)

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 (meses):", round(periodo_seg, 2), "\n")
## Periodo (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)

# ── 2. Periodograma vía FFT ────────────────────────────────────
Px   <- Mod(fft(y))^2
Freq <- 0:(n_obs - 1) / n_obs
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 (FFT):", round(Freq[u], 6), "\n")
## Frecuencia máxima (FFT): 0.167224
cat("Periodo (meses):", round(1 / Freq[u], 2), "\n")
## Periodo (meses): 5.98
abline(v = Freq[u], lty = 2, col = "firebrick", lwd = 2)

# ── 3. mvspec() con bandas de confianza ───────────────────────
par(mar = c(4, 4, 3, 2))
per_mvspec <- astsa::mvspec(y, log = "no",
                             main = "Periodograma mvspec — Desempleo sin tendencia")
ubicacion_mv <- which.max(per_mvspec$spec)
frec_mv      <- per_mvspec$freq[ubicacion_mv]
cat("Frecuencia dominante (mvspec):", round(frec_mv, 6), "\n")
## Frecuencia dominante (mvspec): 0.166667
cat("Periodo (meses):", round(1 / frec_mv, 2), "\n")
## Periodo (meses): 6
abline(v = frec_mv, lty = 2, col = "firebrick", lwd = 2)

# ── Tabla resumen de 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
# El argumento span promedia las estimaciones de frecuencias vecinas,
# reduciendo el ruido y haciendo los picos más visibles.
# Más span = más suavizado = picos más anchos pero más claros.
par(mfrow = c(2, 2), mar = c(4, 4, 3, 2), cex.main = 0.85)

# --- span = 5 (leve) ---
sp5 <- spectrum(y, log = "no", span = 5,
                main = "Periodograma suavizado\nspan = 5 (leve)",
                ylab = "Densidad espectral", xlab = "Frecuencia", col = "steelblue")
abline(v = frec_max, lty = 2, col = "firebrick", lwd = 1.5)
frec_sp5   <- sp5$freq[which.max(sp5$spec)]
periodo_sp5 <- 1 / frec_sp5
cat("span=5        → Frec. principal:", round(frec_sp5, 4),
    "| Periodo:", round(periodo_sp5, 2), "meses\n")
## span=5        → Frec. principal: 0.08 | Periodo: 12.5 meses
# --- span = c(5,5) (dos pasadas) ---
sp55 <- spectrum(y, log = "no", span = c(5, 5),
                 main = "Periodograma suavizado\nspan = c(5,5) — dos pasadas",
                 ylab = "Densidad espectral", xlab = "Frecuencia", col = "steelblue")
abline(v = frec_max, lty = 2, col = "firebrick", lwd = 1.5)
frec_sp55    <- sp55$freq[which.max(sp55$spec)]
periodo_sp55 <- 1 / frec_sp55
cat("span=c(5,5)   → Frec. principal:", round(frec_sp55, 4),
    "| Periodo:", round(periodo_sp55, 2), "meses\n")
## span=c(5,5)   → Frec. principal: 0.0833 | Periodo: 12 meses
# --- span = c(3,3) (moderado) ---
sp33 <- spectrum(y, log = "no", span = c(3, 3),
                 main = "Periodograma suavizado\nspan = c(3,3) — moderado",
                 ylab = "Densidad espectral", xlab = "Frecuencia", col = "steelblue")
abline(v = frec_max, lty = 2, col = "firebrick", lwd = 1.5)
frec_sp33    <- sp33$freq[which.max(sp33$spec)]
periodo_sp33 <- 1 / frec_sp33
cat("span=c(3,3)   → Frec. principal:", round(frec_sp33, 4),
    "| Periodo:", round(periodo_sp33, 2), "meses\n")
## span=c(3,3)   → Frec. principal: 0.0833 | Periodo: 12 meses
# --- span = c(9,9) (fuerte) ---
sp99 <- spectrum(y, log = "no", span = c(9, 9),
                 main = "Periodograma suavizado\nspan = c(9,9) — fuerte",
                 ylab = "Densidad espectral", xlab = "Frecuencia", col = "steelblue")
abline(v = frec_max, lty = 2, col = "firebrick", lwd = 1.5)

frec_sp99    <- sp99$freq[which.max(sp99$spec)]
periodo_sp99 <- 1 / frec_sp99
cat("span=c(9,9)   → Frec. principal:", round(frec_sp99, 4),
    "| Periodo:", round(periodo_sp99, 2), "meses\n")
## span=c(9,9)   → Frec. principal: 0.0833 | Periodo: 12 meses
par(mfrow = c(1, 1))

# Tabla resumen de los cuatro niveles de suavizado
cat("\n=== RESUMEN — FRECUENCIA PRINCIPAL POR NIVEL DE SUAVIZADO ===\n")
## 
## === RESUMEN — FRECUENCIA PRINCIPAL POR NIVEL DE SUAVIZADO ===
resumen_span <- data.frame(
  Span        = c("5 (leve)", "c(5,5) dos pasadas",
                  "c(3,3) moderado", "c(9,9) fuerte"),
  Frecuencia  = round(c(frec_sp5, frec_sp55, frec_sp33, frec_sp99), 6),
  Periodo_mes = round(c(periodo_sp5, periodo_sp55, periodo_sp33, periodo_sp99), 2)
)
print(resumen_span)
##                 Span Frecuencia Periodo_mes
## 1           5 (leve)   0.080000        12.5
## 2 c(5,5) dos pasadas   0.083333        12.0
## 3    c(3,3) moderado   0.083333        12.0
## 4      c(9,9) fuerte   0.083333        12.0

13 Estacionalidad con Variables Dummy

¿Qué se hace aquí?
Se construye una matriz de 12 variables indicadoras (dummy), una por mes, y se ajusta una regresión sin intercepto:

\[y_t = \delta_1 D_{1t} + \delta_2 D_{2t} + \cdots + \delta_{12} D_{12t} + \varepsilon_t\]

Cada coeficiente \(\delta_i\) estima el efecto promedio del mes \(i\) sobre el desempleo en relación a la tendencia. Si los coeficientes difieren significativamente entre meses, la estacionalidad está bien definida. Los residuos del modelo deberían tener un periodograma más plano que la serie original, señal de que las dummies capturaron el ciclo estacional.

s     <- 12
T_obs <- length(Sin_tendencia_ts)

# Construir matriz de dummies: cada columna = indicador de un mes
identidad <- diag(s)
X_dummy   <- matrix(
  rep(t(identidad), ceiling(T_obs / s)),
  ncol  = ncol(identidad), byrow = TRUE
)[1:T_obs, ]
colnames(X_dummy) <- month.abb

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
# Ajuste sin intercepto (~ . - 1): las 12 dummies lo reemplazan
df_dummy <- data.frame(Y = as.numeric(Sin_tendencia_ts), X_dummy)
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
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
# Gráfico de efectos: azul = meses sobre la tendencia, rojo = por debajo
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"))

# Panel: serie observada vs ajuste por dummies, y residuos
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))
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)

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))

# Periodograma antes y después de aplicar las dummies
# Si las dummies capturaron bien la estacionalidad, el periodograma
# de los residuos no debería tener pico dominante en f=1/12.
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")
frec_dom <- per_dummy$freq[which.max(per_dummy$spec)]
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))

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("Dummies significativas al 5%:\n")
## Dummies significativas al 5%:
coef_tabla <- summary(fit_dummy)$coefficients
print(round(coef_tabla[coef_tabla[, 4] < 0.05, ], 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

14 Modelamiento Estacional: Series de Fourier y GAM

¿Qué se hace aquí?
Se aplican dos métodos más flexibles para modelar la estacionalidad sobre la serie sin tendencia:

Series de Fourier (timetk::fourier_vec()): Regresión con armónicos seno/coseno de período 12. El número de armónicos K controla la complejidad: - K=1: solo el ciclo anual principal. - K=2: ciclo anual + variaciones semestrales. - K=3: ciclo anual + semestral + trimestral.

GAM con spline cíclico (mgcv::gam()): Estima la estacionalidad de forma no paramétrica usando un spline cúbico cíclico (bs="cc"), que garantiza que el valor de enero siguiente sea continuo con el diciembre anterior. Es más flexible que Fourier porque no impone una forma funcional fija.

# Crear serie sin tendencia y tibble base para Fourier y GAM
Desempleo_sintend <- valores - tendencia_STL

tbl_sintend <- tibble(Fecha = fechas_STL, Desempleo = Desempleo_sintend) %>%
  mutate(
    Mes      = month(Fecha),
    Anio     = year(Fecha),
    t        = as.numeric(Fecha),
    t_scaled = scale(t)[, 1]
  )

# ── Series de Fourier ──────────────────────────────────────────
tbl_fourier <- tbl_sintend %>%
  mutate(
    C1 = fourier_vec(Fecha, period = 12, K = 1, type = "cos"),
    S1 = fourier_vec(Fecha, period = 12, K = 1, type = "sin"),
    C2 = fourier_vec(Fecha, period = 12, K = 2, type = "cos"),
    S2 = fourier_vec(Fecha, period = 12, K = 2, type = "sin"),
    C3 = fourier_vec(Fecha, period = 12, K = 3, type = "cos"),
    S3 = fourier_vec(Fecha, period = 12, K = 3, type = "sin")
  )

modelo_fourier_K1 <- lm(Desempleo ~ C1 + S1,             data = tbl_fourier)
modelo_fourier_K2 <- lm(Desempleo ~ C1 + S1 + C2 + S2,   data = tbl_fourier)
modelo_fourier_K3 <- lm(Desempleo ~ C1 + S1 + C2 + S2 + C3 + S3, data = tbl_fourier)

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
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()

# ── GAM con spline cúbico cíclico ─────────────────────────────
# bs="cc" garantiza continuidad entre diciembre y enero del año siguiente
tbl_gam <- tbl_sintend %>% mutate(Mes_num = as.numeric(Mes))
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
plot(modelo_gam, main = "Efecto estacional estimado por GAM",
     xlab = "Mes", ylab = "Efecto sobre el desempleo (%)",
     shade = TRUE, col = "blue", lwd = 2)
abline(h = 0, lty = 2, col = "red")

tbl_gam <- tbl_gam %>% mutate(ajuste_gam = fitted(modelo_gam))

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()

# ── Comparativo final Fourier K3 vs GAM ───────────────────────
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_bondad <- data.frame(
  Modelo          = c("Fourier K=1 (ciclo anual)",
                      "Fourier K=2 (anual + semestral)",
                      "Fourier K=3 (anual + semestral + trimestral)",
                      "GAM (spline cíclico)"),
  R2_ajustado     = c(round(summary(modelo_fourier_K1)$adj.r.squared, 4),
                      round(summary(modelo_fourier_K2)$adj.r.squared, 4),
                      round(summary(modelo_fourier_K3)$adj.r.squared, 4),
                      round(summary(modelo_gam)$r.sq,                 4))
)

knitr::kable(tabla_bondad,
             col.names = c("Modelo", "R² Ajustado"),
             caption   = "Bondad de ajuste — Modelos estacionales",
             align     = c("l", "c"))
Bondad de ajuste — Modelos estacionales
Modelo R² Ajustado
Fourier K=1 (ciclo anual) 0.1008
Fourier K=2 (anual + semestral) 0.0971
Fourier K=3 (anual + semestral + trimestral) 0.0959
GAM (spline cíclico) 0.2064

15 Análisis Wavelet y Suavizamiento Splines

Wavelet (biwavelet): Análisis tiempo-frecuencia simultáneo. A diferencia del periodograma (que asume que los ciclos son estacionarios), el wavelet detecta ciclos cuya frecuencia o amplitud cambian en el tiempo. El eje X es el tiempo, el eje Y es el período en meses (escala logarítmica), y el color indica la intensidad del ciclo (rojo = fuerte, azul = débil). La región fuera del cono de influencia tiene menor confiabilidad.

Splines (smooth.spline): Suaviza la serie sin tendencia para revelar la estructura residual. El parámetro spar controla el balance entre ajuste y suavidad: valores cercanos a 0 siguen bien los ciclos, valores cercanos a 1 solo muestran la tendencia residual.

# ── Wavelet ────────────────────────────────────────────────────
indice_wavelet <- 1:length(Desempleo_sintend)
serie_wavelet  <- cbind(indice_wavelet, as.numeric(Desempleo_sintend))
wt_resultado   <- biwavelet::wt(serie_wavelet)

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)

plot(wt_resultado, type = "power.norm",
     main = "Wavelet — Potencia con sesgo (serie sin tendencia)", plot.cb = TRUE)

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))

# ── Suavizamiento Splines ──────────────────────────────────────
tiempo_spline <- as.numeric(fechas_STL)
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)

# Comparación de los tres niveles en un solo gráfico
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)

# Versión ggplot con spline recomendado (spar=0.7)
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 estructura residual tras remover tendencia STL",
       x = "Fecha", y = "Desviación respecto a la tendencia (%)") +
  theme_minimal()

16 Eliminación de la Estacionalidad: Residuo del Modelo Dummy

¿Qué se hace aquí y por qué se eligió este método?
Se extrae el residuo puro del modelo de variables dummy restando, observación por observación, el efecto estacional mensual estimado a la serie sin tendencia:

\[\text{Residuo}_t = \text{Sin tendencia}_t - \hat{\delta}_{mes(t)}\]

donde \(\hat{\delta}_{mes(t)}\) es el coeficiente dummy del mes al que pertenece la observación \(t\), estimado en la regresión sin intercepto \(y_t = \delta_1 D_{1t} + \cdots + \delta_{12} D_{12t} + \varepsilon_t\). Este residuo es equivalente a residuals(fit_dummy) y representa la variación de la tasa de desempleo que no es explicada ni por la tendencia de largo plazo ni por el patrón estacional mensual promedio.

¿Por qué variables dummy y no STL, Fourier o GAM para remover la estacionalidad?
El modelo dummy fue el método que permitió la interpretación más directa y transparente del patrón estacional: cada coeficiente \(\hat{\delta}_i\) es el efecto promedio del mes \(i\) expresado en puntos porcentuales de desempleo, sin suponer ninguna forma funcional. A diferencia de Fourier o GAM, que modelan la estacionalidad como una curva suave continua, las dummies no imponen ninguna restricción sobre la forma del patrón, lo que las hace más apropiadas cuando los efectos mensuales son irregulares entre sí. Frente a STL, las dummies ofrecen la ventaja de que cada efecto estacional tiene una prueba de significancia estadística asociada (valor-p), lo que permite identificar cuáles meses contribuyen significativamente al patrón y cuáles no.

# ============================================================
# ELIMINACIÓN DE LA ESTACIONALIDAD CON VARIABLES DUMMY
# Serie sin tendencia y sin estacionalidad =
#   Serie sin tendencia − efectos estacionales estimados por dummy
#
# Lógica: fit_dummy ya estimó el efecto promedio de cada mes
# (δ1...δ12). Al restar esos efectos a la serie sin tendencia,
# lo que queda es el residuo puro: variación no explicada ni
# por la tendencia ni por el patrón estacional mensual.
#
# Esto equivale directamente a residuals(fit_dummy), que es
# la forma más directa de extraer ese residuo.
# ============================================================

# --- Paso 1: extraer el residuo del modelo dummy ---
# fit_dummy fue ajustado en el chunk estacionalidad_Dummy
residuo_dummy <- as.numeric(residuals(fit_dummy))

# --- Paso 2: recuperar el vector de tiempo para graficar ---
fechas_residuo <- df_stl_clean$Fecha   # fechas sin los NA del STL

# --- Paso 3: data frame para ggplot ---
df_residuo_dummy <- data.frame(
  Fecha   = fechas_residuo,
  Residuo = residuo_dummy
)

# --- Paso 4: estadísticas descriptivas del residuo ---
tabla_residuo <- data.frame(
  Estadístico = c("Media", "Desviación estándar", "Mínimo", "Máximo"),
  Valor       = c(round(mean(residuo_dummy, na.rm = TRUE), 6),
                  round(sd(residuo_dummy,   na.rm = TRUE), 4),
                  round(min(residuo_dummy,  na.rm = TRUE), 4),
                  round(max(residuo_dummy,  na.rm = TRUE), 4))
)

knitr::kable(tabla_residuo,
             col.names = c("Estadístico", "Valor"),
             caption   = "Estadísticas descriptivas — Residuo sin tendencia y sin estacionalidad (Dummy)",
             align     = c("l", "c"))
Estadísticas descriptivas — Residuo sin tendencia y sin estacionalidad (Dummy)
Estadístico Valor
Media 0.0000
Desviación estándar 1.0810
Mínimo -2.2452
Máximo 8.9598 
# --- Paso 5: gráfico del residuo puro ---
# Si el residuo es ruido blanco → dummy capturó toda la estacionalidad.
# Si aún hay patrones → queda estructura no capturada por los 12 efectos mensuales.
# La curva loess (azul) ayuda a visualizar si permanece alguna tendencia residual.
ggplot(df_residuo_dummy, 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 — Serie sin tendencia (STL) y sin estacionalidad (Dummy)",
    subtitle = "Gris = residuo puro | Azul (loess) = estructura remanente | Rojo = referencia cero",
    x        = "Fecha",
    y        = "Residuo (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
## `geom_smooth()` using formula = 'y ~ x'

# --- Paso 6: ACF del residuo para verificar que no queda estructura ---
# Si las barras del ACF caen dentro de las bandas punteadas → ruido blanco.
# Si hay picos significativos en rezago 12 o 24 → la dummy no eliminó
# completamente la estacionalidad.
acf_residuo <- acf(residuo_dummy, lag.max = 48, plot = FALSE)
n_res       <- length(residuo_dummy)
banda_res   <- qnorm(0.975) / sqrt(n_res)

df_acf_res <- data.frame(
  rezago = as.numeric(acf_residuo$lag),
  acf    = as.numeric(acf_residuo$acf)
)

ggplot(df_acf_res, aes(x = rezago, y = acf)) +
  geom_hline(yintercept = 0,          color = "gray40",    linewidth = 0.4) +
  geom_hline(yintercept =  banda_res, color = "steelblue",
             linetype = "dashed",     linewidth = 0.6) +
  geom_hline(yintercept = -banda_res, color = "steelblue",
             linetype = "dashed",     linewidth = 0.6) +
  geom_segment(aes(xend = rezago, yend = 0),
               color = "gray30", linewidth = 0.7) +
  geom_point(color = "gray30", size = 1.5) +
  scale_x_continuous(breaks = seq(0, 48, by = 6)) +
  labs(
    title    = "ACF del residuo — Verificación de ruido blanco",
    subtitle = "Barras dentro de las bandas punteadas → estacionalidad eliminada correctamente",
    x        = "Rezago (meses)",
    y        = "Autocorrelación"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))