Del PIB original al ajuste estacional con X-13

Una guía paso a paso para entender regARIMA, SEATS y la serie desestacionalizada del PIB chileno.
R
time-series
economics
spanish
Autor/a

Joshua Kunst

Fecha de publicación

septiembre 1, 2026

Comparar el PIB de un trimestre con el anterior parece sencillo, pero la serie observada mezcla movimientos económicos, patrones estacionales, efectos de calendario y episodios excepcionales. En este post seguimos el recorrido completo del PIB trimestral de Chile dentro de X-13ARIMA-SEATS. La idea no es tratar seas() como una caja negra, sino entender qué ocurre en cada etapa:

  1. observar la serie oficial;
  2. transformarla a logaritmos;
  3. identificar la regresión y su matriz \(X\);
  4. construir la serie linealizada;
  5. entender el ARIMA y sus innovaciones;
  6. separar tendencia, estacionalidad e irregular mediante SEATS y obtener la serie desestacionalizada final;
  7. usar la serie desestacionalizada para medir la variación trimestral y contrastarla con la publicada por el Banco Central;
  8. usar el regARIMA para producir un pronóstico como extensión opcional.

Aunque seas() ejecuta regARIMA y SEATS en una sola llamada, aquí sus resultados se presentan en ese orden conceptual.

El mismo recorrido, como historia visual

Este ejercicio también puede recorrerse como una historia visual interactiva. Allí la serie ocupa el centro de la escena y se transforma paso a paso para mostrar, sin repetir el código, cómo se pasa del PIB observado a la serie desestacionalizada final.

Recorrer la historia visual del PIB en X-13 →

¿Qué es X-13ARIMA-SEATS y por qué usarlo?

X-13ARIMA-SEATS es un programa de ajuste estacional desarrollado y mantenido por el U.S. Census Bureau. No es solamente un modelo ARIMA ni una herramienta de pronóstico: es un sistema que reúne en un procedimiento reproducible varias etapas necesarias para hacer comparables periodos consecutivos de una serie económica.

Para ello combina el preajuste regARIMA con un método de descomposición y puede extender los extremos de la serie mediante pronósticos. Esto es especialmente útil para el PIB trimestral: X-13 intenta evitar que calendario, estacionalidad y episodios excepcionales se confundan entre sí, y entrega como resultado principal una serie apropiada para comparar un trimestre con el anterior:

\[ \text{PIB observado} \xrightarrow{\text{regARIMA + SEATS}} \text{PIB desestacionalizado}. \]

Entre sus ventajas prácticas se encuentran los regresores de calendario predefinidos, la detección automática de AO, LS y TC, la selección ARIMA, los métodos X-11 y SEATS y los diagnósticos de estacionalidad residual. Además, la especificación puede conservarse y volver a ejecutarse cuando se actualizan los datos.

X-13 no convierte el resultado en una verdad libre de supuestos. La serie desestacionalizada es una estimación latente, puede revisarse al incorporar nuevas observaciones y depende de la transformación, el modelo y las intervenciones seleccionadas. Tampoco identifica las causas económicas de los outliers ni garantiza ser el mejor método de pronóstico. Su principal aporte es ofrecer un marco especializado, diagnosticable y replicable para el ajuste estacional. Véase la descripción oficial del U.S. Census Bureau y el manual de X-13ARIMA-SEATS.

¿Y entonces qué son X-11 y X-12-ARIMA?

Los números identifican generaciones de programas y métodos; no significan once o doce meses y también pueden utilizarse con series trimestrales.

Nombre Qué es Función principal
X-11 Método de ajuste estacional basado en filtros lineales iterativos Usa medias móviles para estimar sucesivamente tendencia, estacionalidad e irregular.
X-12-ARIMA Programa que amplió X-11 Agrega el preajuste regARIMA, la detección de outliers y extensiones mediante pronósticos antes de aplicar X-11.
X-13ARIMA-SEATS Programa actual que amplía X-12-ARIMA Permite realizar el preajuste regARIMA y escoger entre X-11 o SEATS para la descomposición.
SEATS Método de descomposición incluido en X-13 Extrae tendencia, estacionalidad e irregular utilizando la estructura del modelo ARIMA.

La linealización pertenece al regARIMA y ocurre antes de cualquiera de los dos caminos de descomposición. En este post utilizaremos SEATS; X-11 aplicaría en su lugar filtros de medias móviles sobre la serie ya preajustada.

1. Configuración

La credencial del Banco Central se lee desde el entorno. No debe escribirse directamente en este documento. El paquete bcchr se obtiene desde su repositorio de GitHub, ya que no está publicado en CRAN.

Para configurar la credencial localmente, puede abrir ~/.Renviron con usethis::edit_r_environ(), agregar una línea con BCCH_TOKEN=su_token, guardar el archivo y reiniciar R. .Renviron permanece fuera del repositorio: el valor real del token nunca debe copiarse al .qmd ni a otro archivo versionado.

library(bcchr)
library(dplyr)
library(ggforce)
library(ggfortify)
library(ggplot2)
library(seasonal)
library(tibble)
library(tidyr)

x13_navy <- "#0B1D31"
for (geom in c("line", "path", "step")) {
  ggplot2::update_geom_defaults(geom, list(colour = x13_navy))
}
rm(geom)

invisible(checkX13())
options(bcchr.verbose = FALSE)

codigo_pib <- "F032.PIB.FLU.R.CLP.EP18.Z.Z.0.T"
codigo_pib_desestacionalizado_bcch <- "F032.PIB.FLU.R.CLP.EP18.Z.Z.1.T"
fecha_corte <- as.Date("2026-06-30")

# TRUE permite estudiar una descomposición aditiva en logaritmos.
usar_log <- TRUE

# Convierte una tabla con date y series en columnas al formato largo de ggplot.
series_a_formato_largo <- function(datos) {
  pivot_longer(
    datos,
    cols = -all_of("date"),
    names_to = "serie",
    values_to = "value"
  )
}

2. Serie oficial observada

Trabajaremos con la serie oficial del Banco Central de Chile PIB, volumen a precios del año anterior encadenado, referencia 2018 (miles de millones de pesos encadenados), identificada por el código F032.PIB.FLU.R.CLP.EP18.Z.Z.0.T. Para que el modelo explicado sea reproducible, la descarga se limita al segundo trimestre de 2026; las cifras históricas aún pueden estar sujetas a revisiones.

Esta es la única serie directamente observada del ejercicio. La tendencia, la estacionalidad, el irregular y la serie desestacionalizada serán estimaciones latentes de X-13.

pib_original <- get_series(
  codigo_pib,
  from = as.Date("1996-01-01"),
  to = fecha_corte
) |>
  transmute(
    date = as.Date(date),
    value = as.numeric(value)
  ) |>
  filter(!is.na(value)) |>
  arrange(date)
pib_original
# A tibble: 122 × 2
   date        value
   <date>      <dbl>
 1 1996-01-01 20265.
 2 1996-04-01 20369.
 3 1996-07-01 19581.
 4 1996-10-01 21421.
 5 1997-01-01 21401.
 6 1997-04-01 21756.
 7 1997-07-01 21178.
 8 1997-10-01 23335.
 9 1998-01-01 23040.
10 1998-04-01 23266.
# ℹ 112 more rows
if (anyDuplicated(pib_original$date)) {
  stop("La serie contiene fechas duplicadas.", call. = FALSE)
}

if (nrow(pib_original) != length(seq.Date(
  min(pib_original$date),
  max(pib_original$date),
  by = "3 months"
))) {
  stop("La serie contiene trimestres faltantes.", call. = FALSE)
}
Código del gráfico
ggplot(
  pib_original,
  aes(x = date, y = value)
) +
  geom_line(linewidth = 0.7) +
  labs(title = "PIB trimestral original de Chile", x = NULL,
    y = "Miles de millones de pesos encadenados")

3. ¿Por qué aplicar logaritmos?

Sea \(Y_t\) el PIB en niveles (escala original): la serie observada antes de aplicar logaritmos o diferencias. Definimos:

\[ y_t = \log(Y_t). \]

El logaritmo convierte relaciones multiplicativas en sumas. Si en niveles la serie se entiende como el producto de tendencia, estacionalidad e irregular, entonces en logaritmos puede escribirse como una suma. Además, una diferencia logarítmica pequeña puede interpretarse aproximadamente como un cambio porcentual.

La transformación se realiza antes de crear el objeto ts. Después indicaremos a X-13 que no aplique otra transformación.

if (any(pib_original$value <= 0)) {
  stop("El logaritmo requiere valores positivos.", call. = FALSE)
}

pib_modelo_valor <- if (usar_log) {
  log(pib_original$value)
} else {
  pib_original$value
}

escala_modelo <- if (usar_log) {
  "Logaritmo del PIB"
} else {
  "PIB en niveles (escala original)"
}

pib_modelo_ts <- stats::ts(
  pib_modelo_valor,
  start = c(
    as.integer(format(pib_original$date[[1]], "%Y")),
    (as.integer(format(pib_original$date[[1]], "%m")) - 1L) %/% 3L + 1L
  ),
  frequency = 4
)
pib_modelo_ts
          Qtr1      Qtr2      Qtr3      Qtr4
1996  9.916639  9.921753  9.882328  9.972118
1997  9.971207  9.987632  9.960708 10.057709
1998 10.044996 10.054765  9.999273 10.044036
1999 10.019358 10.021937  9.992257 10.096446
2000 10.075702 10.078717 10.041042 10.129639
2001 10.109002 10.119951 10.068302 10.152266
2002 10.120870 10.145501 10.109694 10.198970
2003 10.172832 10.195689 10.154587 10.237013
2004 10.218908 10.250583 10.226352 10.321313
2005 10.278974 10.308309 10.282011 10.375069
2006 10.336800 10.369723 10.336745 10.435776
2007 10.390599 10.423025 10.384255 10.482948
2008 10.449824 10.471879 10.419169 10.490409
2009 10.426971 10.443233 10.411969 10.503122
2010 10.446113 10.505386 10.484289 10.574981
2011 10.530146 10.566075 10.528336 10.628802
2012 10.588614 10.629019 10.592518 10.682611
2013 10.629133 10.662787 10.621391 10.710063
2014 10.653385 10.679690 10.634358 10.727037
2015 10.678065 10.699933 10.656768 10.745149
2016 10.707475 10.712496 10.676474 10.753627
2017 10.703990 10.719604 10.695184 10.784285
2018 10.748841 10.771672 10.718829 10.819943
2019 10.761294 10.786701 10.747587 10.791477
2020 10.757320 10.623902 10.650567 10.792505
2021 10.761593 10.792625 10.801064 10.902840
2022 10.819463 10.839152 10.802431 10.882706
2023 10.818989 10.835542 10.816122 10.899589
2024 10.853197 10.847212 10.839716 10.939742
2025 10.882049 10.883877 10.856618 10.955126
2026 10.878804 10.882012                    

Para comparar ambas escalas, organizamos la serie original y su transformación logarítmica en formato largo y las mostramos en dos facetas.

Código del gráfico
comparacion_transformacion_ancho <- tibble(
  date = pib_original$date,
  `PIB en niveles (escala original)` = pib_original$value,
  `Logaritmo del PIB` = log(pib_original$value)
)

comparacion_transformacion <- series_a_formato_largo(
  comparacion_transformacion_ancho
)

ggplot(
  comparacion_transformacion,
  aes(x = date, y = value)
) +
  geom_line(linewidth = 0.6) +
  facet_wrap(vars(serie),
    ncol = 1, scales = "free_y") +
  labs(title = "PIB antes y después de aplicar logaritmos",
    subtitle = "El logaritmo comprime las variaciones y puede estabilizar su varianza",
    x = NULL, y = NULL)

Las dos facetas representan las mismas observaciones. Las escalas verticales se separan porque sus unidades no son comparables: la primera está en niveles y la segunda en logaritmos. Lo importante es observar cómo el logaritmo comprime las variaciones absolutas a medida que aumenta el nivel del PIB. Esto suele ayudar cuando existe heterocedasticidad asociada al nivel: si las fluctuaciones en pesos crecen junto con el tamaño de la economía, el logaritmo las expresa como cambios aproximadamente proporcionales y puede estabilizar la varianza.

Esta transformación no garantiza que desaparezca toda heterocedasticidad, que debe evaluarse posteriormente sobre las innovaciones del modelo.

4. Un modelo, dos etapas

La llamada siguiente estima el modelo regARIMA y luego ejecuta SEATS. Se guardan también la matriz de regresores y sus efectos para poder inspeccionar la etapa intermedia.

modelo <- seas(
  pib_modelo_ts,
  transform.function = "none",
  outlier.types = c("ao", "ls"),
  estimate.save = "ref",
  regression.save = "rmx"
)
modelo

Call:
seas(x = pib_modelo_ts, transform.function = "none", outlier.types = c("ao", 
    "ls"), estimate.save = "ref", regression.save = "rmx")

Coefficients:
     Leap Year         Weekday        LS1998.4        LS2019.4        LS2020.2  
     0.0102198       0.0006426      -0.0492743      -0.0487984      -0.1450470  
      AO2020.3        LS2020.4        LS2021.3  MA-Seasonal-04  
     0.0603011       0.1092473       0.0396587       0.5886613  

transformfunction(modelo) devuelve "none" porque el logaritmo ya fue aplicado antes de entregar la serie a X-13.

summary(modelo)

Call:
seas(x = pib_modelo_ts, transform.function = "none", outlier.types = c("ao", 
    "ls"), estimate.save = "ref", regression.save = "rmx")

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
Leap Year       0.0102198  0.0027840   3.671 0.000242 ***
Weekday         0.0006426  0.0004471   1.437 0.150604    
LS1998.4       -0.0492743  0.0100152  -4.920 8.66e-07 ***
LS2019.4       -0.0487984  0.0101100  -4.827 1.39e-06 ***
LS2020.2       -0.1450470  0.0102093 -14.207  < 2e-16 ***
AO2020.3        0.0603011  0.0101302   5.953 2.64e-09 ***
LS2020.4        0.1092473  0.0143159   7.631 2.33e-14 ***
LS2021.3        0.0396587  0.0101382   3.912 9.16e-05 ***
MA-Seasonal-04  0.5886613  0.0752607   7.822 5.21e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

SEATS adj.  ARIMA: (0 1 0)(0 1 1)  Obs.: 122  Transform: none
AICc: -697.4, BIC: -671.9  QS (no seasonality in final):0.01927  
Box-Ljung (no autocorr.): 20.77   Shapiro (normality): 0.9883  
[1] "none"

No repasaremos aquí cómo interpretar cada prueba de significancia de la salida. Supondremos que esos diagnósticos fueron revisados y nos concentraremos en el papel de los regresores seleccionados y de las innovaciones, que permiten evaluar si el regARIMA dejó dependencia temporal relevante sin explicar.

4.1 ¿Qué muestra outlier()?

outliers_detectados <- outlier(modelo, full = TRUE)
outliers_detectados
         Qtr1     Qtr2     Qtr3     Qtr4
1996 NA       NA       NA       NA      
1997 NA       NA       NA       NA      
1998 NA       NA       NA       ls1998.4
1999 NA       NA       NA       NA      
2000 NA       NA       NA       NA      
2001 NA       NA       NA       NA      
2002 NA       NA       NA       NA      
2003 NA       NA       NA       NA      
2004 NA       NA       NA       NA      
2005 NA       NA       NA       NA      
2006 NA       NA       NA       NA      
2007 NA       NA       NA       NA      
2008 NA       NA       NA       NA      
2009 NA       NA       NA       NA      
2010 NA       NA       NA       NA      
2011 NA       NA       NA       NA      
2012 NA       NA       NA       NA      
2013 NA       NA       NA       NA      
2014 NA       NA       NA       NA      
2015 NA       NA       NA       NA      
2016 NA       NA       NA       NA      
2017 NA       NA       NA       NA      
2018 NA       NA       NA       NA      
2019 NA       NA       NA       ls2019.4
2020 NA       ls2020.2 ao2020.3 ls2020.4
2021 NA       NA       ls2021.3 NA      
2022 NA       NA       NA       NA      
2023 NA       NA       NA       NA      
2024 NA       NA       NA       NA      
2025 NA       NA       NA       NA      
2026 NA       NA                        

El resultado es una serie de etiquetas: deja NA en los periodos normales y muestra el outlier en el trimestre correspondiente. full = TRUE conserva el tipo y la fecha completa, por ejemplo LS2020.2 o AO2020.3; con FALSE mostraría solamente LS o AO.

  • AO, additive outlier, representa un impacto concentrado en un solo trimestre.
  • LS, level shift, representa un cambio persistente en el nivel desde una fecha.
  • TC, temporary change, representa un impacto transitorio que comienza en una fecha y disminuye gradualmente.

X-13 admite los tres tipos, pero este ejercicio busca solamente AO y LS porque se especificó outlier.types = c("ao", "ls"). Para incluir TC habría que agregar "tc" a ese argumento. La codificación exacta de los AO y LS seleccionados y el efecto \(X_{tj}\beta_j\) se muestran una sola vez en la sección 5.3, al pasar de la matriz \(X\) al efecto de regresión.

Esta función permite localizar las intervenciones, pero no muestra por sí sola sus coeficientes ni su significancia. Esos resultados aparecen en summary(modelo) y posteriormente en coeficientes_regresion.

Preparación del gráfico
serie_observada_outliers <- tibble(
  date = pib_original$date,
  value = as.numeric(original(modelo)),
  evento = toupper(as.character(outliers_detectados))
)

outliers_df <- serie_observada_outliers |>
  filter(!is.na(evento)) |>
  mutate(
    tipo_evento = case_when(
      startsWith(evento, "AO") ~ "Outlier aditivo",
      startsWith(evento, "LS") ~ "Cambio de nivel",
      startsWith(evento, "TC") ~ "Cambio temporal",
      .default = "Intervención"
    )
  )
outliers_df
# A tibble: 6 × 4
  date       value evento   tipo_evento    
  <date>     <dbl> <chr>    <chr>          
1 1998-10-01  10.0 LS1998.4 Cambio de nivel
2 2019-10-01  10.8 LS2019.4 Cambio de nivel
3 2020-04-01  10.6 LS2020.2 Cambio de nivel
4 2020-07-01  10.7 AO2020.3 Outlier aditivo
5 2020-10-01  10.8 LS2020.4 Cambio de nivel
6 2021-07-01  10.8 LS2021.3 Cambio de nivel
Código del gráfico
ggplot(
  serie_observada_outliers,
  aes(x = date, y = value)
) +
  geom_line(linewidth = 0.6) +
  geom_point(data = outliers_df, colour = "#D55E00", size = 1.8) +
  geom_mark_circle(
    data = outliers_df,
    aes(
      group = evento,
      label = evento,
      description = tipo_evento
    ),
    colour = "#D55E00",
    fill = "transparent",
    expand = grid::unit(0, "mm"),
    radius = grid::unit(0, "mm"),
    label.margin = margin(1.5, 1.5, 1.5, 1.5, "mm"),
    label.minwidth = grid::unit(20, "mm"),
    label.fontsize = c(8, 7),
    label.family = plot_font_family,
    label.fill = "#F9F9F9D9",
    label.colour = c("#D55E00", "#5F6873"),
    label.buffer = grid::unit(12, "mm"),
    con.colour = "#D55E00",
    con.size = 0.35,
    con.cap = grid::unit(0.5, "mm"),
    show.legend = FALSE
  ) +
  scale_x_date(
    limits = c(
      min(serie_observada_outliers$date),
      max(serie_observada_outliers$date) + round(365.25 * 3)
    )
  ) +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.24))) +
  labs(title = "PIB observado y outliers detectados por X-13",
    subtitle = "Los puntos localizan AO y LS; las etiquetas no prueban causalidad",
    x = NULL, y = escala_modelo)

Los puntos están ubicados sobre la serie efectivamente observada. Un AO marca el trimestre puntual en que ocurre el impacto; un LS marca la fecha desde la cual cambia el nivel, aunque visualmente se señale un único punto de inicio.

En esta estimación, los LS cercanos a 2020 son compatibles con una caída y una recuperación escalonada durante la crisis del COVID-19. El AO identifica un movimiento excepcional concentrado en un trimestre. Esta correspondencia es una interpretación económica plausible, no una relación causal demostrada por X-13: el programa detecta quiebres estadísticos y desconoce qué acontecimiento los produjo.

X-13 evalúa candidatos AO y LS en las fechas admisibles, incorpora los que superan su valor crítico y vuelve a estimar el regARIMA. La matriz \(X\) final contiene solamente las intervenciones que sobreviven ese proceso. Su selección depende del ARIMA y debe interpretarse como evidencia estadística condicional al modelo, no como una explicación causal del episodio. Véase la sección outlier del manual oficial.

5. La regresión con errores ARIMA

5.1 El modelo regARIMA

La ecuación básica es:

\[ y_t = x_t^\top\beta + z_t, \]

donde \(x_t^\top\beta\) reúne calendario y outliers, mientras que \(z_t\) sigue un proceso ARIMA.

ImportanteNo son dos modelos estimados por separado

La escritura anterior ayuda a interpretar el resultado, pero regARIMA estima conjuntamente los coeficientes de regresión \(\beta\) y los parámetros ARIMA. No realiza primero una regresión OLS y luego ajusta un ARIMA independiente a los residuos obtenidos.

5.2 La matriz de diseño \(X\)

Cada fila de \(X\) corresponde a un trimestre y cada columna a un regresor. Sus valores no siempre son indicadores 0/1: también existen escalones y contrastes centrados.

matriz_x <- series(
  modelo,
  "regression.regressionmatrix",
  reeval = FALSE
) |>
  stats::window(
    start = stats::start(pib_modelo_ts),
    end = stats::end(pib_modelo_ts)
  )

matriz_x_df <- as_tibble(as.data.frame(matriz_x)) |>
  mutate(date = pib_original$date, .before = 1)
matriz_x_df
# A tibble: 122 × 9
   date       Leap.Year Weekday LS1998.4 LS2019.4 LS2020.2 AO2020.3 LS2020.4
   <date>         <dbl>   <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
 1 1996-01-01      0.75       0       -1       -1       -1        0       -1
 2 1996-04-01      0          0       -1       -1       -1        0       -1
 3 1996-07-01      0          1       -1       -1       -1        0       -1
 4 1996-10-01      0          1       -1       -1       -1        0       -1
 5 1997-01-01     -0.25      -1       -1       -1       -1        0       -1
 6 1997-04-01      0          0       -1       -1       -1        0       -1
 7 1997-07-01      0          1       -1       -1       -1        0       -1
 8 1997-10-01      0          1       -1       -1       -1        0       -1
 9 1998-01-01     -0.25      -1       -1       -1       -1        0       -1
10 1998-04-01      0          0       -1       -1       -1        0       -1
# ℹ 112 more rows
# ℹ 1 more variable: LS2021.3 <dbl>

Las principales codificaciones se interpretan así:

  • AO es un outlier aditivo puntual: \(AO_t=\mathbb{1}(t=t_0)\). Vale 1 únicamente en el trimestre afectado.
  • LS es un cambio permanente de nivel: \(LS_t=-\mathbb{1}(t<t_0)\). Vale -1 antes del cambio y 0 desde la fecha.
  • Leap.Year compara febrero con su duración promedio de 28,25 días. Toma 0.75 en el primer trimestre de un año bisiesto, -0.25 en los demás primeros trimestres y 0 fuera del primer trimestre.
  • Weekday es un contraste entre días hábiles y fines de semana dentro del trimestre. No significa que la serie sea semanal.

Para Leap.Year, el contraste queda centrado porque:

\[ 0.75 - 0.25 - 0.25 - 0.25 = 0. \]

De esta manera no introduce por sí solo un cambio en el nivel promedio de largo plazo. Las definiciones completas se encuentran en el manual oficial de X-13.

Código del gráfico
matriz_x_larga <- series_a_formato_largo(matriz_x_df)

ggplot(
  matriz_x_larga,
  aes(x = date, y = value)
) +
  geom_hline(yintercept = 0, colour = "grey80") +
  geom_step(linewidth = 0.5) +
  facet_wrap(vars(serie),
    ncol = 2, scales = "free_y") +
  labs(title = "Variables de la matriz X del regARIMA",
    subtitle = "LS = cambio permanente; AO = pulso 0/1; calendario = contraste",
    x = NULL, y = "Valor de X")

5.3 De la matriz \(X\) al efecto de regresión

Cada columna de \(X\) tiene un coeficiente \(\beta_j\) estimado conjuntamente con el ARIMA. Para cada trimestre, X-13 suma las contribuciones de todos los regresores:

\[ (X\beta)_t = \sum_{j=1}^{k}X_{tj}\beta_j. \]

El resultado es una sola serie: el efecto de regresión que retiraremos en la sección siguiente para construir la serie linealizada. Primero extraemos los coeficientes estimados por X-13.

coeficientes_regresion <- as_tibble(modelo$est$reg) |>
  transmute(
    variable,
    beta = as.numeric(estimate),
    error_estandar = as.numeric(standard.error)
  )
coeficientes_regresion
# A tibble: 8 × 3
  variable       beta error_estandar
  <chr>         <dbl>          <dbl>
1 Leap Year  0.0102         0.00278 
2 Weekday    0.000643       0.000447
3 LS1998.4  -0.0493         0.0100  
4 LS2019.4  -0.0488         0.0101  
5 LS2020.2  -0.145          0.0102  
6 AO2020.3   0.0603         0.0101  
7 LS2020.4   0.109          0.0143  
8 LS2021.3   0.0397         0.0101  
nombres_matriz_x <- make.names(colnames(matriz_x))
beta_por_variable <- setNames(
  coeficientes_regresion$beta,
  make.names(coeficientes_regresion$variable)
)
beta <- unname(beta_por_variable[nombres_matriz_x])

if (anyNA(beta)) {
  stop("No fue posible alinear los coeficientes con las columnas de X.")
}

names(beta) <- colnames(matriz_x)

Para obtener el efecto en cada trimestre se multiplica el valor de cada regresor por su coeficiente y se suman las contribuciones. Los coeficientes se alinean por nombre con las columnas de \(X\), y el resultado se compara con la serie que informa X-13. Las tablas se recortan al periodo observado porque el programa también construye valores para sus extensiones de pronóstico.

Comprobar el efecto X beta
efectos_regarima <- series(
  modelo,
  "estimate.regressioneffects",
  reeval = FALSE
) |>
  stats::window(
    start = stats::start(pib_modelo_ts),
    end = stats::end(pib_modelo_ts)
  )

efecto_x_beta_reconstruido <- as.numeric(as.matrix(matriz_x) %*% beta)
efecto_x_beta_x13 <- as.numeric(efectos_regarima[, "Total.Reg"])

if (!isTRUE(all.equal(
  efecto_x_beta_reconstruido,
  efecto_x_beta_x13,
  tolerance = 1e-8
))) {
  stop("El efecto X beta reconstruido no coincide con el informado por X-13.")
}

comprobacion_x_beta <- tibble(
  date = pib_original$date,
  `Efecto X beta reconstruido` = efecto_x_beta_reconstruido,
  `X beta informado por X-13` = efecto_x_beta_x13
)
comprobacion_x_beta
# A tibble: 122 × 3
   date       `Efecto X beta reconstruido` `X beta informado por X-13`
   <date>                            <dbl>                       <dbl>
 1 1996-01-01                       0.102                       0.102 
 2 1996-04-01                       0.0942                      0.0942
 3 1996-07-01                       0.0949                      0.0949
 4 1996-10-01                       0.0949                      0.0949
 5 1997-01-01                       0.0910                      0.0910
 6 1997-04-01                       0.0942                      0.0942
 7 1997-07-01                       0.0949                      0.0949
 8 1997-10-01                       0.0949                      0.0949
 9 1998-01-01                       0.0910                      0.0910
10 1998-04-01                       0.0942                      0.0942
# ℹ 112 more rows

La coincidencia confirma que las contribuciones por periodo reconstruyen el efecto total informado por X-13. En esta especificación también podemos agruparlas en calendario, cambios de nivel y outliers aditivos:

\[ x_t^\top\beta=C_t+L_t+A_t, \]

donde \(C_t\) reúne Weekday y Leap.Year, \(L_t\) la suma de los LS y \(A_t\) la suma de los AO.

efectos_por_regresor <- sweep(
  as.matrix(matriz_x),
  MARGIN = 2,
  STATS = beta,
  FUN = "*"
)

colnames(efectos_por_regresor) <- make.names(colnames(efectos_por_regresor))

componentes_efecto_regresion <- efectos_por_regresor |>
  as.data.frame() |>
  as_tibble() |>
  mutate(date = pib_original$date, .before = 1) |>
  pivot_longer(
    cols = -date,
    names_to = "regresor",
    values_to = "value"
  ) |>
  mutate(
    componente = case_when(
      regresor %in% c("Leap.Year", "Weekday") ~ "C(t): calendario",
      startsWith(regresor, "LS") ~ "L(t): cambios de nivel",
      startsWith(regresor, "AO") ~ "A(t): outliers aditivos",
      .default = NA_character_
    )
  )

regresores_no_clasificados <- componentes_efecto_regresion |>
  filter(is.na(componente)) |>
  distinct(regresor) |>
  pull(regresor)

if (length(regresores_no_clasificados) > 0L) {
  stop(
    "Falta clasificar en C(t), L(t) o A(t): ",
    paste(regresores_no_clasificados, collapse = ", ")
  )
}

componentes_efecto_regresion <- componentes_efecto_regresion |>
  summarise(value = sum(value), .by = c(date, componente))

efecto_regresion_total <- componentes_efecto_regresion |>
  summarise(value = sum(value), .by = date)

if (!isTRUE(all.equal(
  efecto_regresion_total$value,
  efecto_x_beta_x13,
  tolerance = 1e-8
))) {
  stop("La suma C(t) + L(t) + A(t) no coincide con el efecto total.")
}

efecto_regresion_total <- efecto_regresion_total |>
  mutate(componente = "X beta: efecto total")

componentes_efecto_regresion <- bind_rows(
  efecto_regresion_total,
  componentes_efecto_regresion
) |>
  mutate(
    componente = factor(
      componente,
      levels = c(
        "X beta: efecto total",
        "C(t): calendario",
        "L(t): cambios de nivel",
        "A(t): outliers aditivos"
      )
    )
  )

componentes_efecto_regresion
# A tibble: 488 × 3
   date        value componente          
   <date>      <dbl> <fct>               
 1 1996-01-01 0.102  X beta: efecto total
 2 1996-04-01 0.0942 X beta: efecto total
 3 1996-07-01 0.0949 X beta: efecto total
 4 1996-10-01 0.0949 X beta: efecto total
 5 1997-01-01 0.0910 X beta: efecto total
 6 1997-04-01 0.0942 X beta: efecto total
 7 1997-07-01 0.0949 X beta: efecto total
 8 1997-10-01 0.0949 X beta: efecto total
 9 1998-01-01 0.0910 X beta: efecto total
10 1998-04-01 0.0942 X beta: efecto total
# ℹ 478 more rows
Código del gráfico
ggplot(
  componentes_efecto_regresion,
  aes(x = date, y = value)
) +
  geom_hline(yintercept = 0, colour = "grey70") +
  geom_line(linewidth = 0.6) +
  facet_wrap(vars(componente), ncol = 2) +
  labs(
    title = "Descomposición del efecto de regresión",
    subtitle = "En cada trimestre, X beta = C(t) + L(t) + A(t)",
    x = NULL,
    y = escala_modelo
  )

Los cuatro paneles comparten la escala vertical: el primero es exactamente la suma de los otros tres. Los escalones de \(L_t\) aparecen porque cada LS cambia de régimen en una fecha distinta; \(C_t\) oscila con el calendario y \(A_t\) solo se aparta de cero en el trimestre del outlier aditivo.

5.4 ¿Qué es la serie linealizada?

Al despejar \(z_t\) de la ecuación regARIMA obtenemos:

\[ z_t = y_t - x_t^\top\beta. \]

Se resta \(X\beta\) para neutralizar temporalmente efectos identificables que podrían confundirse con la estacionalidad: calendario, cambios de nivel y outliers. El resultado es una serie ajustada por esos efectos de regresión, no una observación nueva ni una estimación causal. Su función es evitar que SEATS los confunda con la estacionalidad que debe estimar.

Como \(y_t=\log(Y_t)\), la resta en logaritmos equivale a una división en niveles:

\[ z_t = \log(Y_t)-x_t^\top\beta = \log\left(\frac{Y_t}{\exp(x_t^\top\beta)}\right). \]

La columna Reg.Resids tiene un nombre potencialmente confuso: contiene los residuos de retirar la regresión, es decir, la serie linealizada. No contiene las innovaciones finales del ARIMA.

pib_linealizado <- efectos_regarima[, "Reg.Resids"]
pib_linealizado
          Qtr1      Qtr2      Qtr3      Qtr4
1996  9.814760  9.827540  9.787471  9.877262
1997  9.880191  9.893418  9.865851  9.962853
1998  9.953980  9.960551  9.904417  9.998454
1999  9.977616  9.976997  9.946675 10.050864
2000 10.023098 10.033778  9.997709 10.086306
2001 10.065012 10.075012 10.024969 10.106684
2002 10.079129 10.100561 10.064112 10.153388
2003 10.131090 10.150750 10.109005 10.191431
2004 10.166304 10.205644 10.180770 10.275731
2005 10.237232 10.263370 10.236429 10.331736
2006 10.292809 10.324784 10.293413 10.392443
2007 10.346608 10.378085 10.340922 10.437366
2008 10.397220 10.426939 10.373587 10.444827
2009 10.385229 10.398294 10.366387 10.457540
2010 10.404371 10.460447 10.438707 10.529399
2011 10.488404 10.521136 10.482754 10.585470
2012 10.536009 10.584080 10.549185 10.637029
2013 10.587391 10.617847 10.575809 10.664481
2014 10.611643 10.634750 10.588776 10.681455
2015 10.636323 10.654993 10.611186 10.699567
2016 10.654870 10.667557 10.630892 10.710294
2017 10.659999 10.674664 10.651851 10.740952
2018 10.704850 10.726733 10.675496 10.774361
2019 10.719552 10.741762 10.702005 10.794694
2020 10.753514 10.772808 10.738530 10.831521
2021 10.804449 10.832284 10.800422 10.902197
2022 10.822660 10.839152 10.801788 10.884313
2023 10.819938 10.835542 10.817729 10.901196
2024 10.845532 10.847212 10.839073 10.939099
2025 10.885246 10.883877 10.855976 10.954484
2026 10.882001 10.882012                    
preajuste_regarima_ancho <- tibble(
  date = pib_original$date,
  `Original usada por X-13` = as.numeric(original(modelo)),
  `Efectos de regresión` = as.numeric(efectos_regarima[, "Total.Reg"]),
  `Linealizada = Original - efectos de regresión` = as.numeric(
    pib_linealizado
  )
)

preajuste_regarima <- series_a_formato_largo(preajuste_regarima_ancho) |>
  mutate(
    serie = factor(
      serie,
      levels = c(
        "Original usada por X-13",
        "Linealizada = Original - efectos de regresión",
        "Efectos de regresión"
      )
    )
  )
preajuste_regarima
# A tibble: 366 × 3
   date       serie                                          value
   <date>     <fct>                                          <dbl>
 1 1996-01-01 Original usada por X-13                       9.92  
 2 1996-01-01 Efectos de regresión                          0.102 
 3 1996-01-01 Linealizada = Original - efectos de regresión 9.81  
 4 1996-04-01 Original usada por X-13                       9.92  
 5 1996-04-01 Efectos de regresión                          0.0942
 6 1996-04-01 Linealizada = Original - efectos de regresión 9.83  
 7 1996-07-01 Original usada por X-13                       9.88  
 8 1996-07-01 Efectos de regresión                          0.0949
 9 1996-07-01 Linealizada = Original - efectos de regresión 9.79  
10 1996-10-01 Original usada por X-13                       9.97  
# ℹ 356 more rows
Código del gráfico
ggplot(
  preajuste_regarima,
  aes(x = date, y = value)
) +
  geom_hline(
    data = filter(
      preajuste_regarima,
      serie == "Efectos de regresión"
    ) |>
      distinct(serie) |>
      mutate(yintercept = 0),
    aes(yintercept = yintercept),
    inherit.aes = FALSE,
    colour = "grey70"
  ) +
  geom_line(linewidth = 0.6) +
  facet_wrap(vars(serie), ncol = 1, scales = "free_y") +
  labs(title = "Preajuste regARIMA antes de SEATS",
    subtitle = "Escalas verticales independientes; el efecto de regresión aparece abajo",
    x = NULL, y = escala_modelo,
    caption = "La línea cero se muestra solo para los efectos de regresión")

Los paneles tienen la misma altura, pero usan escalas verticales independientes. Esto permite ver la forma temporal de cada serie; no corresponde comparar directamente la amplitud visual entre paneles. El efecto de regresión queda al final porque es el término que se resta de la serie original para construir la linealizada.

La linealización no desestacionaliza la serie. Por construcción, este ajuste elimina solamente los efectos modelados en \(X\beta\), como calendario y outliers, y debe conservar el componente estacional \(S_t\). SEATS lo estimará y retirará recién en el paso siguiente. Las dos líneas se separan especialmente donde existen cambios de nivel u outliers; su diferencia es exactamente \(X\beta\).

Código del gráfico
comparacion_original_linealizada_ancho <- tibble(
  date = pib_original$date,
  `log(PIB) observado` = as.numeric(original(modelo)),
  `Serie linealizada` = as.numeric(pib_linealizado)
)

comparacion_original_linealizada <- series_a_formato_largo(
  comparacion_original_linealizada_ancho
)

ggplot(
  comparacion_original_linealizada,
  aes(
    x = date,
    y = value,
    colour = serie,
    linetype = serie
  )
) +
  geom_line(linewidth = 0.7) +
  labs(title = "log(PIB) observado y serie linealizada",
    subtitle = "La linealización retira X beta, no la estacionalidad",
    x = NULL, y = escala_modelo, colour = NULL, linetype = NULL)

Nota¿Por qué ambas series coinciden al final y no al comienzo?

La explicación está en la codificación de los cambios de nivel. Para un LS en \(t_0\), X-13 utiliza:

\[ LS_t= \begin{cases} -1, & t<t_0,\\ 0, & t\ge t_0. \end{cases} \]

La serie linealizada se calcula como \(z_t=y_t-X_t\beta\). Si, por ejemplo, \(\beta_{LS}=-0.05\), entonces:

\[ z_t= \begin{cases} y_t-0.05, & t<t_0,\\ y_t, & t\ge t_0. \end{cases} \]

El coeficiente negativo indica que el nivel posterior es 0.05 log-puntos menor que el anterior. La linealización lleva el tramo anterior al nivel del régimen posterior; por eso aparece debajo de la observada antes del cambio y coincide con ella después.

Cuando existen varios LS, cada columna deja de contribuir al cruzar su propia fecha. Después del último cambio de nivel todas las columnas LS valen cero. Si los demás efectos de regresión también son pequeños, las dos líneas terminan casi superpuestas. Esto es una convención de referencia: la serie linealizada queda anclada al régimen más reciente; no significa que los outliers solo afecten al pasado.

Un AO funciona distinto: solo separa las líneas durante el trimestre puntual en que su columna vale 1.

5.5 ¿Para qué se usa el ARIMA después de linealizar?

La serie linealizada ya fue obtenida en el paso anterior:

\[ z_t=y_t-x_t^\top\beta. \]

Este paso no vuelve a construirla ni genera otra serie corregida. El ARIMA describe la dependencia temporal que todavía existe en \(z_t\): cuánto ayuda su propio pasado a explicar su valor actual. Al aplicar esa dinámica quedan las innovaciones \(\varepsilon_t\), es decir, la información nueva que el pasado del modelo no pudo anticipar.

El propósito es doble:

  1. comprobar que las innovaciones se parezcan a ruido blanco, como diagnóstico de que el ARIMA capturó adecuadamente la dependencia temporal;
  2. entregar a SEATS la estructura ARIMA que utilizará para estimar tendencia, estacionalidad e irregular.

Por tanto, la secuencia conceptual es \(y_t \rightarrow z_t \rightarrow \varepsilon_t\): primero se retira \(X\beta\), después se modela la dinámica de la serie resultante. La estimación de \(\beta\) y del ARIMA sí ocurre conjuntamente dentro de regARIMA; esta secuencia solo sirve para leer el resultado.

residuos_regarima <- residuals(modelo)

# La diferenciación ARIMA puede dejar sin innovación el primer trimestre. Los
# índices ts permiten alinear las 121 innovaciones con los 122 datos originales.
residuos_df <- tibble(
  date = pib_original$date,
  periodo = as.numeric(stats::time(pib_modelo_ts))
) |>
  left_join(
    tibble(
      periodo = as.numeric(stats::time(residuos_regarima)),
      value = as.numeric(residuos_regarima)
    ),
    by = "periodo"
  ) |>
  select(-all_of("periodo"))
residuos_df
# A tibble: 122 × 2
   date          value
   <date>        <dbl>
 1 1996-01-01 NA      
 2 1996-04-01  0.00157
 3 1996-07-01 -0.00166
 4 1996-10-01 -0.00325
 5 1997-01-01  0.0123 
 6 1997-04-01  0.00137
 7 1997-07-01  0.0115 
 8 1997-10-01  0.00530
 9 1998-01-01 -0.00458
10 1998-04-01 -0.00585
# ℹ 112 more rows

Aquí residuals(modelo) devuelve las innovaciones del modelo completo. El primer trimestre puede quedar como NA debido a la diferenciación ARIMA; no es una observación faltante del PIB.

Las innovaciones deberían fluctuar alrededor de cero y no mostrar autocorrelaciones persistentes.

Código del gráfico
ggplot(
  residuos_df,
  aes(x = date, y = value)
) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_line(linewidth = 0.6) +
  labs(title = "Innovaciones del modelo regARIMA",
    subtitle = "Deben fluctuar alrededor de cero sin patrones persistentes",
    x = NULL, y = "Innovación")

Código del gráfico
acf_innovaciones <- stats::acf(
  stats::na.omit(residuos_regarima),
  plot = FALSE
)

autoplot(acf_innovaciones) +
  labs(
    title = "Autocorrelación de las innovaciones regARIMA",
    subtitle = "Las barras deberían permanecer dentro de las bandas de confianza",
    x = "Rezago",
    y = "Autocorrelación"
  )

NotaInnovación regARIMA e irregular de SEATS no son lo mismo

residuals(modelo) devuelve \(\varepsilon_t\), el error de predicción del regARIMA. irregular(modelo) devuelve \(I_t\), una componente latente de la descomposición realizada posteriormente por SEATS.

6. La descomposición SEATS

SEATS sí trabaja sobre la serie linealizada. Para distinguir esta descomposición interna de los componentes que X-13 entrega al final, usaremos el superíndice lin:

\[ z_t = T_t^{lin} + S_t^{lin} + I_t^{lin}, \]

donde \(T_t^{lin}\) es la tendencia-ciclo, \(S_t^{lin}\) la estacionalidad e \(I_t^{lin}\) el irregular estimados a partir de la serie sin efectos de regresión. SEATS llama tendencia-ciclo al componente persistente no estacional: reúne tanto la trayectoria de largo plazo como movimientos económicos de mediano plazo, porque no busca separarlos en dos componentes distintos. En esta etapa, el ajuste estacional interno elimina \(S_t^{lin}\). El superíndice \(SA\) abrevia seasonally adjusted, es decir, «desestacionalizada»:

\[ z_t^{SA}=z_t-S_t^{lin}=T_t^{lin}+I_t^{lin}. \]

Este resultado interno todavía no es lo mismo que final(modelo). Como vimos al descomponer el efecto de regresión, en esta especificación \(x_t^\top\beta=C_t+L_t+A_t\). El componente de calendario es:

\[ C_t =\beta_{\text{Weekday}}\,\text{Weekday}_t +\beta_{\text{Leap.Year}}\,\text{Leap.Year}_t. \]

Weekday controla la composición de días hábiles y fines de semana de cada trimestre; Leap.Year controla el día adicional de febrero en los años bisiestos.

Después de que SEATS estima los componentes sobre \(z_t\), X-13 distribuye los outliers que había retirado: los LS se reincorporan a la tendencia-ciclo final y los AO al irregular final. El episodio del COVID-19 permite ver esta recomposición: la caída de \(L_t\) desplaza la tendencia-ciclo final hacia abajo respecto de la estimada sobre la serie linealizada, mientras que el pulso positivo de \(A_t\) en 2020.3 desplaza el irregular final hacia arriba:

\[ T_t^{final}=T_t^{lin}+L_t, \qquad I_t^{final}=I_t^{lin}+A_t. \]

El gráfico muestra esos aportes, no los componentes finales completos. Por eso \(A_t>0\) no implica necesariamente \(I_t^{final}>0\): indica que el AO eleva el irregular respecto de \(I_t^{lin}\). En este ajuste concreto, sin embargo, los aportes son suficientemente fuertes para conservar la dirección del movimiento: el AO mantiene el alza del irregular final en 2020.3 y los LS mantienen el desplazamiento descendente de la tendencia-ciclo final durante la crisis. El signo positivo del AO no convierte la crisis en un episodio positivo; señala que 2020.3 quedó por encima de la trayectoria condicionada a los cambios de nivel y al resto del modelo. El irregular de SEATS tampoco debe confundirse con la innovación \(\varepsilon_t\) del regARIMA.

El efecto de calendario \(C_t\) no se reincorpora porque precisamente se desea retirarlo de la serie desestacionalizada. Por tanto, no se suma nuevamente todo \(x_t^\top\beta\), sino solamente los efectos de outliers asignados a componentes no estacionales:

\[ \begin{aligned} y_t^{SA,\,final} &=z_t^{SA}+L_t+A_t \\ &=(T_t^{lin}+L_t)+(I_t^{lin}+A_t) \\ &=T_t^{final}+I_t^{final}. \end{aligned} \]

Por eso la desestacionalizada final conserva episodios como la crisis del COVID-19. Desestacionalizar significa retirar movimientos recurrentes de la época del año y efectos de calendario; no significa borrar shocks económicos no estacionales. Una serie corregida también por outliers sería otro producto, distinto de la serie desestacionalizada final utilizada aquí.

pib_desestacionalizado <- final(modelo)
pib_tendencia <- trend(modelo)
pib_estacional <- series(modelo, "seats.seasonal")
pib_irregular <- irregular(modelo)
Código del gráfico
comparacion_original_desestacionalizada_ancho <- tibble(
  date = pib_original$date,
  `PIB observado` = as.numeric(original(modelo)),
  `PIB desestacionalizado` = as.numeric(pib_desestacionalizado)
)

comparacion_original_desestacionalizada <- series_a_formato_largo(
  comparacion_original_desestacionalizada_ancho
)
comparacion_original_desestacionalizada
# A tibble: 244 × 3
   date       serie                  value
   <date>     <chr>                  <dbl>
 1 1996-01-01 PIB observado           9.92
 2 1996-01-01 PIB desestacionalizado  9.90
 3 1996-04-01 PIB observado           9.92
 4 1996-04-01 PIB desestacionalizado  9.92
 5 1996-07-01 PIB observado           9.88
 6 1996-07-01 PIB desestacionalizado  9.93
 7 1996-10-01 PIB observado           9.97
 8 1996-10-01 PIB desestacionalizado  9.94
 9 1997-01-01 PIB observado           9.97
10 1997-01-01 PIB desestacionalizado  9.97
# ℹ 234 more rows
Código del gráfico
ggplot(
  comparacion_original_desestacionalizada,
  aes(
    x = date,
    y = value,
    colour = serie,
    linetype = serie
  )
) +
  geom_line(linewidth = 0.7) +
  labs(title = "PIB observado y desestacionalizado",
    subtitle = "Se retira la estacionalidad; los AO y LS se conservan por defecto",
    x = NULL, y = escala_modelo, colour = NULL, linetype = NULL)

La diferencia regular entre ambas líneas corresponde principalmente a la estacionalidad y a los efectos de calendario retirados. Los episodios extraordinarios permanecen en la serie final mediante el mecanismo descrito al inicio de esta sección. Por eso la corrección de outliers es temporal: facilita la estimación de SEATS, pero no borra la crisis del resultado final.

componentes_ancho <- tibble(
  date = pib_original$date,
  `Serie observada` = as.numeric(original(modelo)),
  `Desestacionalizada = Tendencia + Irregular` = as.numeric(
    pib_desestacionalizado
  ),
  Tendencia = as.numeric(pib_tendencia),
  `Estacional (alrededor de 0)` = as.numeric(pib_estacional),
  `Irregular (alrededor de 0)` = as.numeric(pib_irregular)
)

componentes <- series_a_formato_largo(componentes_ancho)
componentes
# A tibble: 610 × 3
   date       serie                                          value
   <date>     <chr>                                          <dbl>
 1 1996-01-01 Serie observada                             9.92    
 2 1996-01-01 Desestacionalizada = Tendencia + Irregular  9.90    
 3 1996-01-01 Tendencia                                   9.90    
 4 1996-01-01 Estacional (alrededor de 0)                 0.00875 
 5 1996-01-01 Irregular (alrededor de 0)                 -0.000424
 6 1996-04-01 Serie observada                             9.92    
 7 1996-04-01 Desestacionalizada = Tendencia + Irregular  9.92    
 8 1996-04-01 Tendencia                                   9.92    
 9 1996-04-01 Estacional (alrededor de 0)                 0.00573 
10 1996-04-01 Irregular (alrededor de 0)                  0.000872
# ℹ 600 more rows

La identidad desestacionalizada \(=T_t+I_t\) corresponde a las componentes finales, que ya contienen los outliers reincorporados. La serie observada se presenta como referencia y no como una suma visual de las demás facetas: para reconstruirla también habría que incluir los efectos de calendario retirados.

Código del gráfico
ggplot(
  componentes,
  aes(x = date, y = value)
) +
  geom_line(linewidth = 0.6) +
  facet_wrap(vars(serie),
    ncol = 2, scales = "free_y") +
  labs(title = "Descomposición aditiva del PIB",
    subtitle = paste(escala_modelo, "· escalas verticales independientes"),
    x = NULL, y = NULL)

Nota¿Por qué el estacional y el irregular están alrededor de cero?

El PIB fue transformado a logaritmos y X-13 recibió transform.function = "none". Por eso SEATS trabaja de manera aditiva. En una descomposición multiplicativa en niveles, ambos factores fluctuarían alrededor de uno.

ImportanteEl resultado principal de X-13

El objetivo del ajuste estacional se alcanza con pib_desestacionalizado <- final(modelo). Esta es la serie que se usa para comparar trimestres consecutivos sin que el patrón estacional ni los efectos de calendario oculten los movimientos económicos de corto plazo:

\[ y_t^{SA,\,final}=T_t^{final}+I_t^{final}. \]

No es una tendencia suavizada: todavía contiene el irregular y los episodios extraordinarios reincorporados por X-13. La tendencia sirve para mirar la trayectoria más persistente; la desestacionalizada sirve para analizar el movimiento trimestral observado una vez retirados los efectos recurrentes.

7. Ya tenemos la serie desestacionalizada, ¿y ahora qué?

Al volver a niveles, denotamos por \(Y_t^{SA}\) la serie desestacionalizada final. Ahora sí tiene sentido comparar cada trimestre con el inmediatamente anterior:

\[ g_t =100\left(\frac{Y_t^{SA}}{Y_{t-1}^{SA}}-1\right). \]

Como el objeto del modelo permanece en logaritmos, el código evalúa este mismo cociente mediante la identidad \(Y_t^{SA}/Y_{t-1}^{SA}=\exp(y_t^{SA}-y_{t-1}^{SA})\).

variacion_trimestral <- tibble(
  date = pib_original$date,
  pib_observado = as.numeric(original(modelo)),
  pib_desestacionalizado = as.numeric(pib_desestacionalizado)
) |>
  mutate(
    pib_observado_anterior = lag(pib_observado),
    pib_desestacionalizado_anterior = lag(pib_desestacionalizado),
    variacion_observada_pct = if (usar_log) {
      100 * (
        exp(pib_observado - pib_observado_anterior) - 1
      )
    } else {
      100 * (
        pib_observado / pib_observado_anterior - 1
      )
    },
    variacion_desestacionalizada_pct = if (usar_log) {
      100 * (
        exp(
          pib_desestacionalizado - pib_desestacionalizado_anterior
        ) - 1
      )
    } else {
      100 * (
        pib_desestacionalizado / pib_desestacionalizado_anterior - 1
      )
    },
    efecto_ajuste_pp = variacion_observada_pct -
      variacion_desestacionalizada_pct
  )
variacion_trimestral
# A tibble: 122 × 8
   date       pib_observado pib_desestacionalizado pib_observado_anterior
   <date>             <dbl>                  <dbl>                  <dbl>
 1 1996-01-01          9.92                   9.90                  NA   
 2 1996-04-01          9.92                   9.92                   9.92
 3 1996-07-01          9.88                   9.93                   9.92
 4 1996-10-01          9.97                   9.94                   9.88
 5 1997-01-01          9.97                   9.97                   9.97
 6 1997-04-01          9.99                   9.98                   9.97
 7 1997-07-01          9.96                  10.0                    9.99
 8 1997-10-01         10.1                   10.0                    9.96
 9 1998-01-01         10.0                   10.0                   10.1 
10 1998-04-01         10.1                   10.0                   10.0 
# ℹ 112 more rows
# ℹ 4 more variables: pib_desestacionalizado_anterior <dbl>,
#   variacion_observada_pct <dbl>, variacion_desestacionalizada_pct <dbl>,
#   efecto_ajuste_pp <dbl>
Código del gráfico
variacion_trimestral_ancho <- variacion_trimestral |>
  transmute(
    date,
    `Sin ajuste estacional` = variacion_observada_pct,
    Desestacionalizada = variacion_desestacionalizada_pct
  )

variacion_trimestral_larga <- series_a_formato_largo(
  variacion_trimestral_ancho
)
variacion_trimestral_larga
# A tibble: 244 × 3
   date       serie                   value
   <date>     <chr>                   <dbl>
 1 1996-01-01 Sin ajuste estacional NA     
 2 1996-01-01 Desestacionalizada    NA     
 3 1996-04-01 Sin ajuste estacional  0.513 
 4 1996-04-01 Desestacionalizada     1.59  
 5 1996-07-01 Sin ajuste estacional -3.87  
 6 1996-07-01 Desestacionalizada     1.25  
 7 1996-10-01 Sin ajuste estacional  9.39  
 8 1996-10-01 Desestacionalizada     1.06  
 9 1997-01-01 Sin ajuste estacional -0.0911
10 1997-01-01 Desestacionalizada     2.74  
# ℹ 234 more rows
Código del gráfico
ggplot(
  variacion_trimestral_larga,
  aes(x = date, y = value,
    colour = serie, linetype = serie)
) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_line(linewidth = 0.65) +
  labs(title = "Variación trimestral con y sin ajuste estacional",
    subtitle = "La serie original mezcla el movimiento económico con efectos recurrentes",
    x = NULL, y = "Variación trimestral (%)", colour = NULL, linetype = NULL)

Código del gráfico
ggplot(
  variacion_trimestral,
  aes(x = date, y = efecto_ajuste_pp)
) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_linerange(
    aes(ymin = 0, ymax = efecto_ajuste_pp),
    colour = x13_navy,
    linewidth = 0.45,
    alpha = 0.75
  ) +
  geom_point(colour = x13_navy, size = 0.8) +
  labs(title = "Diferencia producida por el ajuste estacional",
    subtitle = "Cada segmento parte de cero; la diferencia no representa un error",
    x = NULL, y = "Diferencia (puntos porcentuales)")

resumen_ajuste <- variacion_trimestral |>
  summarise(
    volatilidad_sin_ajuste = sd(
      variacion_observada_pct,
      na.rm = TRUE
    ),
    volatilidad_desestacionalizada = sd(
      variacion_desestacionalizada_pct,
      na.rm = TRUE
    ),
    diferencia_absoluta_media_pp = mean(abs(efecto_ajuste_pp), na.rm = TRUE)
  )
resumen_ajuste
# A tibble: 1 × 3
  volatilidad_sin_ajuste volatilidad_desestacionalizada diferencia_absoluta_me…¹
                   <dbl>                          <dbl>                    <dbl>
1                   5.75                           1.83                     4.80
# ℹ abbreviated name: ¹​diferencia_absoluta_media_pp

El segundo gráfico no mide un error de pronóstico ni permite saber cuál valor es “verdadero”. Muestra cuánto cambia la lectura trimestral al retirar los efectos estacionales y de calendario. El resultado se evalúa por la reducción de oscilaciones recurrentes y por los diagnósticos del modelo, no simplemente porque disminuya la volatilidad. Por eso resumen_ajuste es descriptivo y no una medida de precisión.

La desestacionalización hace la comparación más homogénea y reproducible, no la convierte en una verdad libre de supuestos. El resultado se puede replicar si se mantienen la misma versión de los datos, la especificación y la versión de X-13. Además, las últimas observaciones pueden revisarse cuando se agregan nuevos datos y vuelven a estimarse el modelo y los factores estacionales.

7.1 Comparación con la serie oficial del Banco Central

La serie oficial ofrece un contraste útil, pero no funciona como una respuesta única contra la cual validar mecánicamente el ajuste estimado. El procedimiento desarrollado aquí aplica un solo regARIMA y SEATS al PIB agregado. Como este análisis no reconstruye las especificaciones, los insumos ni las políticas de revisión de la serie oficial, no se supone que ambas implementaciones sean idénticas.

La serie oficial utilizada es F032.PIB.FLU.R.CLP.EP18.Z.Z.1.T, correspondiente al PIB encadenado, desestacionalizado y con referencia 2018.

pib_desestacionalizado_bcch <- get_series(
  codigo_pib_desestacionalizado_bcch,
  from = min(pib_original$date),
  to = fecha_corte
) |>
  transmute(
    date = as.Date(date),
    nivel_bcch = as.numeric(value)
  ) |>
  filter(!is.na(nivel_bcch)) |>
  arrange(date)

comparacion_bcch <- tibble(
  date = pib_original$date,
  nivel_estimado = if (usar_log) {
    exp(as.numeric(pib_desestacionalizado))
  } else {
    as.numeric(pib_desestacionalizado)
  }
) |>
  inner_join(pib_desestacionalizado_bcch, by = "date") |>
  arrange(date) |>
  mutate(
    variacion_estimada_pct = 100 * (
      nivel_estimado / lag(nivel_estimado) - 1
    ),
    variacion_bcch_pct = 100 * (
      nivel_bcch / lag(nivel_bcch) - 1
    ),
    diferencia_nivel_pct = 100 * (nivel_estimado / nivel_bcch - 1),
    diferencia_variacion_pp = variacion_estimada_pct - variacion_bcch_pct
  )

comparacion_bcch
# A tibble: 122 × 7
   date       nivel_estimado nivel_bcch variacion_estimada_pct
   <date>              <dbl>      <dbl>                  <dbl>
 1 1996-01-01         19935.     20074.                 NA    
 2 1996-04-01         20252.     20273.                  1.59 
 3 1996-07-01         20504.     20475.                  1.25 
 4 1996-10-01         20722.     20769.                  1.06 
 5 1997-01-01         21290.     21240.                  2.74 
 6 1997-04-01         21640.     21682.                  1.64 
 7 1997-07-01         22178.     22172.                  2.49 
 8 1997-10-01         22545.     22644.                  1.65 
 9 1998-01-01         22934.     22826.                  1.73 
10 1998-04-01         23152.     23213.                  0.947
# ℹ 112 more rows
# ℹ 3 more variables: variacion_bcch_pct <dbl>, diferencia_nivel_pct <dbl>,
#   diferencia_variacion_pp <dbl>
Código del gráfico
comparacion_bcch_larga <- bind_rows(
  comparacion_bcch |>
    transmute(
      date,
      medida = "Nivel desestacionalizado",
      `Ajuste estimado con X-13` = nivel_estimado,
      `Serie oficial del Banco Central` = nivel_bcch
    ),
  comparacion_bcch |>
    transmute(
      date,
      medida = "Variación respecto del trimestre anterior (%)",
      `Ajuste estimado con X-13` = variacion_estimada_pct,
      `Serie oficial del Banco Central` = variacion_bcch_pct
    )
) |>
  pivot_longer(
    cols = -all_of(c("date", "medida")),
    names_to = "serie",
    values_to = "value"
  )

ggplot(
  comparacion_bcch_larga,
  aes(
    x = date,
    y = value,
    colour = serie,
    linetype = serie
  )
) +
  geom_hline(
    data = tibble(
      medida = "Variación respecto del trimestre anterior (%)",
      yintercept = 0
    ),
    aes(yintercept = yintercept),
    inherit.aes = FALSE,
    colour = "grey70"
  ) +
  geom_line(linewidth = 0.65) +
  facet_wrap(vars(medida), ncol = 1, scales = "free_y") +
  scale_colour_manual(values = c(
    "Ajuste estimado con X-13" = x13_navy,
    "Serie oficial del Banco Central" = "#0072B2"
  )) +
  labs(
    title = "Dos estimaciones del PIB desestacionalizado",
    subtitle = "La cercanía es informativa; las especificaciones no se suponen idénticas",
    x = NULL,
    y = NULL,
    colour = NULL,
    linetype = NULL
  )

Código del gráfico
ggplot(
  comparacion_bcch,
  aes(x = date, y = diferencia_nivel_pct)
) +
  geom_hline(yintercept = 0, colour = "grey70") +
  geom_line(colour = x13_navy, linewidth = 0.65) +
  labs(
    title = "Diferencia entre las series desestacionalizadas",
    subtitle = "Ajuste estimado con X-13 respecto de la serie oficial",
    x = NULL,
    y = "Diferencia (%)"
  )

Para resumir la distancia se calculan RMSE y MAPE. El RMSE conserva la unidad de las series y el MAPE expresa su distancia media en porcentaje. En este contexto no miden precisión predictiva: comparan dos estimaciones contemporáneas del PIB desestacionalizado. La serie oficial no se trata como un valor verdadero ni como un conjunto de prueba.

resumen_comparacion_bcch <- comparacion_bcch |>
  summarise(
    observaciones_nivel = sum(complete.cases(nivel_estimado, nivel_bcch)),
    observaciones_variacion = sum(complete.cases(
      variacion_estimada_pct,
      variacion_bcch_pct
    )),
    desde = min(date),
    hasta = max(date),
    rmse_nivel = sqrt(mean(
      (nivel_estimado - nivel_bcch)^2,
      na.rm = TRUE
    )),
    mape_nivel_pct = mean(abs(diferencia_nivel_pct), na.rm = TRUE),
    diferencia_absoluta_media_variacion_pp = mean(
      abs(diferencia_variacion_pp),
      na.rm = TRUE
    )
  )

resumen_comparacion_bcch
# A tibble: 1 × 7
  observaciones_nivel observaciones_variacion desde      hasta      rmse_nivel
                <int>                   <int> <date>     <date>          <dbl>
1                 122                     121 1996-01-01 2026-04-01       136.
# ℹ 2 more variables: mape_nivel_pct <dbl>,
#   diferencia_absoluta_media_variacion_pp <dbl>

La comparación utiliza 122 niveles trimestrales y 121 variaciones, desde 1996-01 hasta 2026-04. En ese intervalo, el RMSE de 136 está expresado en la misma unidad que el PIB: miles de millones de pesos encadenados. Por construcción, esta métrica da más peso a las diferencias grandes antes de tomar la raíz. El MAPE de 0,263% expresa la diferencia absoluta entre los niveles como porcentaje de la serie oficial (no corresponde a 26%, pues el resultado ya está expresado en porcentaje). Finalmente, la distancia entre las tasas de variación trimestral es de 0,383 puntos porcentuales en promedio.

Las diferencias pueden reflejar especificaciones, insumos, revisiones y decisiones de modelación distintas; no deben interpretarse automáticamente como errores de una de las dos series. La metodología oficial se describe en la Base de Datos Estadísticos del Banco Central.

8. Pronóstico como extensión opcional

El pronóstico normalmente no se realiza sobre el irregular. Si el irregular se comporta como ruido blanco, su esperanza futura es aproximadamente cero. El regARIMA pronostica la serie completa utilizando los regresores conocidos y la dinámica ARIMA estimada.

Nota¿Qué significa que el pronóstico sea opcional?

Publicar un pronóstico no es el objetivo del ajuste estacional: el producto principal ya se obtuvo con final(modelo). Sin embargo, X-13 puede usar pronósticos y retropredicciones internamente para extender los extremos de la serie y mejorar la estimación de los filtros. Aquí se muestran doce trimestres futuros solo como un uso adicional del regARIMA estimado.

8.1 El SARIMA seleccionado

orden_arima <- modelo$model$arima$model
theta_estacional <- stats::coef(modelo)[
  grepl("^MA-Seasonal", names(stats::coef(modelo)))
]

if (!identical(orden_arima, "(0 1 0)(0 1 1)")) {
  stop(
    "El SARIMA seleccionado cambió a ", orden_arima,
    "; actualice la explicación de esta sección."
  )
}

orden_arima
[1] "(0 1 0)(0 1 1)"
theta_estacional
MA-Seasonal-04 
     0.5886613 

Con los datos actuales, X-13 seleccionó un SARIMA \((0,1,0)(0,1,1)_4\). El subíndice \(4\) corresponde a la frecuencia trimestral. El modelo aplica una primera diferencia regular (\(d=1\)), una primera diferencia estacional (\(D=1\)) e incluye un término de media móvil estacional en el rezago 4 (\(Q=1\)); no incorpora términos AR ni MA regulares.

No es solamente “el trimestre anterior”. Primero se aplica la diferencia regular y luego la estacional:

\[ (1-B)(1-B^4)z_t =z_t-z_{t-1}-z_{t-4}+z_{t-5}, \]

donde \(Bz_t=z_{t-1}\). Según la convención de signos de X-13, el modelo completo es:

\[ \underbrace{(1-B)(1-B^4)z_t}_{\text{transformación de los datos}} = \underbrace{(1-\Theta B^4)\varepsilon_t}_{\text{modelo probabilístico}} =\varepsilon_t-\Theta\varepsilon_{t-4}. \]

Esta igualdad no es una identidad algebraica producida por las diferencias. El lado izquierdo es la transformación aplicada a los datos; el lado derecho es el supuesto probabilístico con que el SARIMA modela la serie transformada. El coeficiente \(\Theta\) se estima con los datos; no procede de la operación de diferenciar.

El gráfico siguiente aplica las dos diferencias a la serie linealizada. La primera reduce la evolución persistente del nivel; la segunda atenúa la persistencia que se repite cada cuatro trimestres.

arima_paso_a_paso <- tibble(
  date = pib_original$date,
  linealizada = as.numeric(pib_linealizado)
) |>
  mutate(
    diferencia_regular = linealizada - lag(linealizada),
    diferencia_regular_estacional = diferencia_regular -
      lag(diferencia_regular, 4L)
  )
arima_paso_a_paso
# A tibble: 122 × 4
   date       linealizada diferencia_regular diferencia_regular_estacional
   <date>           <dbl>              <dbl>                         <dbl>
 1 1996-01-01        9.81           NA                           NA       
 2 1996-04-01        9.83            0.0128                      NA       
 3 1996-07-01        9.79           -0.0401                      NA       
 4 1996-10-01        9.88            0.0898                      NA       
 5 1997-01-01        9.88            0.00293                     NA       
 6 1997-04-01        9.89            0.0132                       0.000448
 7 1997-07-01        9.87           -0.0276                       0.0125  
 8 1997-10-01        9.96            0.0970                       0.00721 
 9 1998-01-01        9.95           -0.00887                     -0.0118  
10 1998-04-01        9.96            0.00657                     -0.00666 
# ℹ 112 more rows
Código del gráfico
arima_paso_a_paso_ancho <- arima_paso_a_paso |>
  transmute(
    date,
    `"1. Serie linealizada:"~z[t]` = linealizada,
    `"2. Primera diferencia regular:"~(1-B)*z[t]` = diferencia_regular,
    `"3. Primera diferencia regular y estacional:"~(1-B)*(1-B^4)*z[t]` =
      diferencia_regular_estacional
  )

arima_paso_a_paso_largo <- series_a_formato_largo(arima_paso_a_paso_ancho)

ggplot(
  arima_paso_a_paso_largo,
  aes(x = date, y = value)
) +
  geom_line(linewidth = 0.6) +
  facet_wrap(vars(serie),
    ncol = 1, scales = "free_y", labeller = label_parsed) +
  labs(title = "Cómo las diferencias transforman la serie linealizada",
    subtitle = "Cada panel agrega una diferencia; las escalas verticales son independientes",
    x = NULL, y = NULL)

Lo esperado es que en el primer panel todavía se observen el nivel creciente y la dinámica estacional. La primera diferencia regular, mostrada en el segundo panel, elimina gran parte de la evolución persistente del nivel, pero puede dejar visible una estructura que se repite cada cuatro trimestres. Esa persistencia motiva la primera diferencia estacional del tercer panel, que compara el cambio trimestral con el observado cuatro trimestres antes:

\[ (1-B^4)(1-B)z_t =(z_t-z_{t-1})-(z_{t-4}-z_{t-5}). \]

Después de ambas diferencias, se espera que esa estructura estacional se atenúe y que la serie fluctúe de manera más estable alrededor de cero. El tercer panel todavía no contiene directamente las innovaciones \(\varepsilon_t\). Bajo el modelo seleccionado, la serie completa de ese panel sigue un MA estacional de orden 1:

\[ w_t=(1-B)(1-B^4)z_t =\varepsilon_t-\Theta\varepsilon_{t-4}. \]

Por tanto, no vemos una componente MA aislada, sino una serie diferenciada cuya dependencia restante en el rezago 4 es modelada por \(Q=1\). En otras palabras, la tercera faceta muestra una realización de un MA estacional de orden 1 con periodo 4. Visto como un polinomio MA ordinario, también puede describirse como un MA(4) restringido:

\[ w_t =\varepsilon_t +0\varepsilon_{t-1} +0\varepsilon_{t-2} +0\varepsilon_{t-3} -\Theta\varepsilon_{t-4}. \]

No es un MA(4) general porque los coeficientes de los rezagos 1, 2 y 3 están fijados en cero y solo se estima el coeficiente del rezago estacional 4. Esta relación suele apreciarse mejor en la autocorrelación que en el gráfico temporal. Una vez ajustado ese término, las innovaciones resultantes deberían parecer ruido blanco.

Al invertir las diferencias para pronosticar se obtiene, esquemáticamente:

\[ \widehat z_t =\widehat z_{t-1}+\widehat z_{t-4}-\widehat z_{t-5} -\Theta\widehat\varepsilon_{t-4}, \]

porque la esperanza de una innovación futura es cero. La recurrencia con los rezagos 4 y 5 hace que el pronóstico de la serie completa conserve un patrón trimestral marcado. La persistencia del patrón procede principalmente de \(D=1\): el cambio de cada trimestre queda relacionado con el cambio del mismo trimestre del año anterior. El término \(Q=1\) corrige además los primeros horizontes mediante la innovación ocurrida cuatro trimestres atrás; no significa que simplemente se copie el trimestre inmediatamente anterior.

X-13 seleccionó esa estructura porque encontró dependencia estacional en la serie linealizada. El patrón futuro procede del SARIMA y, en menor medida, de los regresores futuros de calendario. No se obtiene sumando nuevamente la componente estacional estimada por SEATS. Si una actualización de los datos lleva a X-13 a seleccionar otro orden, esta ecuación también debe actualizarse.

8.2 Pronóstico de la serie completa

Para un horizonte \(h\), la predicción puntual puede escribirse como:

\[ \widehat y_{T+h\mid T} =x_{T+h}^\top\widehat\beta +\widehat z_{T+h\mid T}, \]

donde \(x_{T+h}^\top\widehat\beta\) aporta los efectos de regresión conocidos para el periodo futuro y \(\widehat z_{T+h\mid T}\) es el pronóstico ARIMA de la serie linealizada. La innovación futura no se suma a la predicción puntual porque su esperanza condicional es cero; su incertidumbre sí se refleja en los intervalos de pronóstico.

El objeto matriz_x guardado anteriormente sirve solo para inspeccionar la linealización; no se entrega manualmente a update(). Esta llamada vuelve a ejecutar la especificación regARIMA almacenada y X-13 reconstruye la matriz de regresores, incluidas sus filas para el horizonte futuro. En esas filas, un AO histórico vale cero porque su efecto se concentra en un trimestre. Con la codificación centrada utilizada aquí, un LS también vale cero desde su fecha y, por tanto, en el horizonte futuro. El cambio de régimen sí permanece en el nivel y la historia de la serie linealizada que alimentan al ARIMA; no se suma como una contribución futura distinta de cero de esa columna de \(X\). Ambos tipos de intervención influyeron previamente en la estimación de la dinámica.

horizonte_pronostico <- 12L

modelo_pronostico <- stats::update(
  modelo,
  forecast.maxlead = horizonte_pronostico,
  forecast.save = "fct"
)

pronostico <- series(
  modelo_pronostico,
  "forecast.forecasts",
  reeval = FALSE
)
pronostico
        forecast  lowerci  upperci
2026 Q3 10.85929 10.83750 10.88108
2026 Q4 10.95442 10.92362 10.98523
2027 Q1 10.88780 10.85009 10.92551
2027 Q2 10.89468 10.85115 10.93820
2027 Q3 10.87196 10.81865 10.92526
2027 Q4 10.96709 10.90556 11.02862
2028 Q1 10.91132 10.84242 10.98023
2028 Q2 10.90734 10.83204 10.98264
2028 Q3 10.88237 10.79720 10.96755
2028 Q4 10.97750 10.88355 11.07146
2029 Q1 10.91538 10.81338 11.01738
2029 Q2 10.92001 10.81061 11.02940

Como la serie entregada a X-13 ya estaba en logaritmos y la transformación interna se desactivó, estos resultados permanecen en la escala logarítmica. La exponencial los devuelve a niveles; para interpretar la media condicional en niveles sería necesario considerar además la corrección por sesgo.

pronostico_df <- as_tibble(as.data.frame(pronostico)) |>
  slice_head(n = horizonte_pronostico) |>
  mutate(
    date = seq.Date(
      from = max(pib_original$date),
      by = "3 months",
      length.out = horizonte_pronostico + 1L
    )[-1L],
    .before = 1
  )
pronostico_df
# A tibble: 12 × 4
   date       forecast lowerci upperci
   <date>        <dbl>   <dbl>   <dbl>
 1 2026-07-01     10.9    10.8    10.9
 2 2026-10-01     11.0    10.9    11.0
 3 2027-01-01     10.9    10.9    10.9
 4 2027-04-01     10.9    10.9    10.9
 5 2027-07-01     10.9    10.8    10.9
 6 2027-10-01     11.0    10.9    11.0
 7 2028-01-01     10.9    10.8    11.0
 8 2028-04-01     10.9    10.8    11.0
 9 2028-07-01     10.9    10.8    11.0
10 2028-10-01     11.0    10.9    11.1
11 2029-01-01     10.9    10.8    11.0
12 2029-04-01     10.9    10.8    11.0
historial_pronostico <- tibble(
  date = pib_original$date,
  value = as.numeric(original(modelo))
) |>
  slice_tail(n = 60L)
historial_pronostico
# A tibble: 60 × 2
   date       value
   <date>     <dbl>
 1 2011-07-01  10.5
 2 2011-10-01  10.6
 3 2012-01-01  10.6
 4 2012-04-01  10.6
 5 2012-07-01  10.6
 6 2012-10-01  10.7
 7 2013-01-01  10.6
 8 2013-04-01  10.7
 9 2013-07-01  10.6
10 2013-10-01  10.7
# ℹ 50 more rows
trayectoria_pronostico <- bind_rows(
  historial_pronostico |>
    slice_tail(n = 1L) |>
    transmute(
      date,
      value,
      serie = "Pronóstico regARIMA"
    ),
  pronostico_df |>
    transmute(
      date,
      value = forecast,
      serie = "Pronóstico regARIMA"
  )
)
Código del gráfico
ggplot() +
  geom_ribbon(
    data = pronostico_df,
    aes(x = date, ymin = lowerci, ymax = upperci),
    fill = "#D55E00", alpha = 0.15
  ) +
  geom_line(
    data = historial_pronostico,
    aes(x = date, y = value,
      colour = "log(PIB) observado"),
    linewidth = 0.7
  ) +
  geom_line(
    data = trayectoria_pronostico,
    aes(x = date, y = value, colour = serie),
    linewidth = 0.7
  ) +
  scale_colour_manual(values = c(
    "log(PIB) observado" = x13_navy,
    "Pronóstico regARIMA" = "#D55E00"
  )) +
  labs(title = "log(PIB) observado y pronóstico regARIMA",
    subtitle = "Últimos quince años observados y doce trimestres pronosticados",
    x = NULL, y = escala_modelo, colour = NULL)

Se muestra el pronóstico junto a log(PIB) observado porque ambos representan la serie completa \(y_t\): X-13 vuelve a sumar \(X\beta\) al pronóstico ARIMA. La serie linealizada es una entrada intermedia y la desestacionalizada es otro producto; ninguna de ellas es el objetivo de forecast.forecasts.

ImportantePara pronosticar no se necesita la descomposición SEATS

Casi basta con llegar a la serie linealizada, pero todavía falta estimar su ARIMA. El recorrido mínimo es:

  1. estimar conjuntamente \(\beta\) y el ARIMA mediante regARIMA;
  2. construir \(z_t=y_t-X_t\beta\);
  3. pronosticar \(z_t\) con el ARIMA;
  4. sumar los efectos futuros \(X_{T+h}\beta\) para recuperar el pronóstico de \(y_{T+h}\).

SEATS no es necesario para ese pronóstico: su función es descomponer la serie y producir el ajuste estacional. Tampoco debe suponerse que el regARIMA elegido por X-13 será siempre el mejor modelo predictivo. Esa pregunta se responde mediante una evaluación fuera de muestra frente a alternativas como otros SARIMA, modelos estructurales, métodos multivariados o modelos de nowcasting.

Nota¿Qué aportan los outliers al pronóstico?

Las intervenciones evitan que el ARIMA tenga que explicar observaciones excepcionales como si fueran parte de su dinámica regular. Los LS dejan la serie linealizada referida al régimen más reciente y los AO aíslan pulsos puntuales. Con la codificación de este modelo, sus columnas históricas valen cero en el horizonte, pero ambas influyeron en la estimación del ARIMA. X-13 puede así evitar que los episodios pasados distorsionen el pronóstico; no puede anticipar una crisis futura que todavía no ha sido especificada.

Un regARIMA puede ampliarse con regresores externos \(w_t\), como actividad mundial, precio del cobre o tasas externas. Sin embargo, conviene separar los objetivos. Para ajustar estacionalmente, incluir una variable porque predice el PIB podría retirar parte del ciclo económico que se desea conservar. Para pronosticar, en cambio, puede construirse otro modelo dinámico sobre la serie desestacionalizada.

Antes de incorporar esas variables deben revisarse su transformación, estacionariedad o cointegración, rezagos, fecha real de publicación y disponibilidad futura. Su aporte se evalúa fuera de muestra, no mediante una correlación contemporánea aislada.

Resumen conceptual

La receta completa puede leerse así en escala logarítmica, sin olvidar que los parámetros de la regresión y del ARIMA se estiman conjuntamente. En esta especificación, \(x_t^\top\beta=C_t+L_t+A_t\) reúne calendario, cambios de nivel y outliers aditivos:

\[ y_t \xrightarrow{\;-x_t^\top\beta\;} z_t \xrightarrow{\;\text{SEATS}\;} T_t^{lin}+S_t^{lin}+I_t^{lin}. \]

SEATS retira \(S_t^{lin}\) para producir su ajuste interno. Después X-13 reincorpora los efectos no estacionales:

\[ z_t^{SA}=T_t^{lin}+I_t^{lin} \xrightarrow{\;+(L_t+A_t)\;} y_t^{SA,\,final}=T_t^{final}+I_t^{final}. \]

Así, el resultado principal final(modelo) deja fuera la estacionalidad y el calendario, pero conserva los shocks no estacionales. La contabilidad completa de la observada en esta especificación es:

\[ y_t=C_t+S_t^{lin}+T_t^{final}+I_t^{final}. \]

  1. Serie observada y transformación. X-13 puede escoger automáticamente entre trabajar en niveles o aplicar logaritmos mediante transform.function = "auto". Aquí se aplica log() antes de llamar a X-13 y luego se usa transform.function = "none" para hacer explícita la escala: \(y_t=\log(Y_t)\). En esta escala la descomposición se expresa mediante sumas y las componentes estacional e irregular fluctúan alrededor de cero, lo que facilita su visualización pedagógica.
  2. Preajuste regARIMA. X-13 representa la observada como \(y_t=x_t^\top\beta+z_t\), donde \(X\beta\) contiene calendario y outliers.
  3. Serie linealizada. Se retiran temporalmente esos efectos: \(z_t=y_t-x_t^\top\beta\). Esta serie todavía conserva la estacionalidad.
  4. Dinámica ARIMA. El ARIMA modela la dependencia temporal de \(z_t\) y deja innovaciones \(\varepsilon_t\) que deberían parecer ruido blanco.
  5. Descomposición SEATS. Sobre la linealizada se estima \(z_t=T_t^{lin}+S_t^{lin}+I_t^{lin}\).
  6. Serie desestacionalizada interna. Se elimina solamente la estacionalidad: \(z_t^{SA}=z_t-S_t^{lin}=T_t^{lin}+I_t^{lin}\).
  7. Resultado final de X-13. Los AO y LS seleccionados en esta especificación se reincorporan porque forman parte de la historia económica. No se reincorporan la estacionalidad ni los efectos de calendario, como Weekday, Leap.Year y feriados: esos permanecen retirados para permitir comparaciones entre trimestres consecutivos.

\[ y_t^{SA,\,final} =z_t^{SA}+L_t+A_t =T_t^{final}+I_t^{final}. \]

  1. Uso del resultado. Al volver a niveles, llamamos \(Y_t^{SA}\) a la serie desestacionalizada final. Su variación respecto del trimestre anterior es:

\[ g_t =100\left(\frac{Y_t^{SA}}{Y_{t-1}^{SA}}-1\right). \]

El pronóstico es una extensión posterior. Para generarlo se vuelve a la serie completa, no al irregular ni a la serie desestacionalizada:

\[ \widehat y_{T+h\mid T} =x_{T+h}^\top\widehat\beta +\widehat z_{T+h\mid T}. \]

Versión práctica sin el paso a paso

En una aplicación real no es necesario reconstruir \(X\beta\), las diferencias ARIMA ni cada identidad mostrada en este documento. Una vez limpia y ordenada la serie, las funciones esenciales son:

Función Uso práctico
stats::ts() Declara la frecuencia y el periodo inicial de la serie.
seas() Estima conjuntamente el regARIMA y ejecuta SEATS.
final() Extrae el resultado principal: la serie desestacionalizada final.
summary() Revisa regresores, outliers, orden ARIMA y coeficientes.
plot() Compara la serie original y la desestacionalizada, y etiqueta las intervenciones.
trend() y irregular() Extraen componentes para análisis adicional.
series() Recupera tablas específicas de X-13 y pronósticos guardados.

El flujo mínimo, dejando que X-13 seleccione automáticamente la transformación, es el siguiente:

pib_niveles_ts <- stats::ts(
  pib_original$value,
  start = stats::start(pib_modelo_ts),
  frequency = 4
)

modelo_real <- seas(
  pib_niveles_ts,
  transform.function = "auto",
  outlier.types = c("ao", "ls")
)

summary(modelo_real)

Call:
seas(x = pib_niveles_ts, transform.function = "auto", outlier.types = c("ao", 
    "ls"))

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
LS1998.4       -0.04913    0.01078  -4.558 5.16e-06 ***
LS2019.4       -0.04865    0.01083  -4.494 6.99e-06 ***
LS2020.2       -0.15378    0.01063 -14.468  < 2e-16 ***
LS2020.3        0.06073    0.01084   5.604 2.09e-08 ***
LS2020.4        0.04926    0.01083   4.548 5.41e-06 ***
LS2021.3        0.04017    0.01085   3.702 0.000214 ***
MA-Seasonal-04  0.62306    0.07417   8.400  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

SEATS adj.  ARIMA: (0 1 0)(0 1 1)  Obs.: 122  Transform: log
AICc:  1772, BIC:  1793  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 28.45   Shapiro (normality): 0.9776 *
pib_desestacionalizado_real <- final(modelo_real)
pib_desestacionalizado_real
         Qtr1     Qtr2     Qtr3     Qtr4
1996 20117.16 20285.40 20522.21 20718.89
1997 21257.24 21665.15 22194.30 22549.99
1998 22901.73 23171.30 23066.36 22221.49
1999 22348.77 22429.62 22893.26 23401.39
2000 23675.01 23723.87 24016.34 24198.31
2001 24530.22 24688.06 24643.99 24754.92
2002 24891.50 25306.52 25645.82 25926.08
2003 26279.32 26598.55 26787.68 26909.17
2004 27587.86 28104.49 28745.57 29222.00
2005 29360.18 29784.21 30381.71 30789.75
2006 31156.33 31662.04 32092.76 32687.98
2007 32913.37 33381.02 33652.43 34256.11
2008 34963.86 35034.54 34809.63 34493.95
2009 34246.19 34049.48 34525.73 34910.91
2010 34956.66 36224.02 37109.53 37496.01
2011 38015.65 38485.88 38815.98 39570.57
2012 40276.94 40986.46 41408.22 41749.45
2013 41914.50 42422.41 42642.98 42901.66
2014 42885.63 43187.41 43225.11 43642.97
2015 43881.25 44119.24 44211.21 44453.48
2016 45135.63 44721.47 45078.67 44827.31
2017 44967.68 45066.79 45945.14 46193.33
2018 47029.34 47469.53 47065.09 47838.47
2019 47639.23 48210.68 48424.20 46463.83
2020 47472.96 41008.17 43902.94 46455.06
2021 47738.81 48634.09 50951.43 51756.83
2022 50638.14 51087.47 50957.41 50647.92
2023 50601.91 51051.40 51620.00 51430.64
2024 52310.82 51786.88 52865.63 53444.93
2025 53819.79 53785.79 53785.39 54208.75
2026 53661.37 53708.54                  
resultado_real <- tibble(
  date = pib_original$date,
  pib_observado = as.numeric(original(modelo_real)),
  pib_desestacionalizado = as.numeric(pib_desestacionalizado_real)
) |>
  mutate(
    variacion_trimestral_pct = 100 * (
      pib_desestacionalizado / lag(pib_desestacionalizado) - 1
    )
  )
resultado_real
# A tibble: 122 × 4
   date       pib_observado pib_desestacionalizado variacion_trimestral_pct
   <date>             <dbl>                  <dbl>                    <dbl>
 1 1996-01-01        20265.                 20117.                   NA    
 2 1996-04-01        20369.                 20285.                    0.836
 3 1996-07-01        19581.                 20522.                    1.17 
 4 1996-10-01        21421.                 20719.                    0.958
 5 1997-01-01        21401.                 21257.                    2.60 
 6 1997-04-01        21756.                 21665.                    1.92 
 7 1997-07-01        21178.                 22194.                    2.44 
 8 1997-10-01        23335.                 22550.                    1.60 
 9 1998-01-01        23040.                 22902.                    1.56 
10 1998-04-01        23266.                 23171.                    1.18 
# ℹ 112 more rows

El método gráfico de seasonal resume directamente la serie original, la desestacionalizada y las etiquetas AO y LS seleccionadas:

local({
  parametros_graficos <- par(
    family = plot_font_family,
    font.main = 1,
    col.main = "#24364B",
    col.lab = "#24364B",
    col.axis = "#5F6873"
  )
  on.exit(par(parametros_graficos))

  plot(
    modelo_real,
    outliers = TRUE,
    main = "PIB original y desestacionalizado",
    xlab = "",
    ylab = "PIB (escala original)"
  )
})

Esta versión parte del PIB en niveles y usa transform.function = "auto". Por eso X-13 puede escoger una transformación y un modelo distintos de los del recorrido pedagógico, donde el logaritmo se aplicó explícitamente y se fijó transform.function = "none" para poder observar una descomposición aditiva.

En términos operativos, stats::ts() + seas() + final() constituyen el núcleo. El resto sirve para diagnosticar, interpretar o reutilizar resultados.

El pronóstico continúa siendo opcional:

modelo_real_pronostico <- stats::update(
  modelo_real,
  forecast.maxlead = 12L,
  forecast.save = "fct"
)

pronostico_real <- series(
  modelo_real_pronostico,
  "forecast.forecasts",
  reeval = FALSE
)
utils::head(pronostico_real, 12L)
        forecast  lowerci  upperci
2026 Q3 51970.13 50779.86 53188.30
2026 Q4 57127.61 55287.31 59029.16
2027 Q1 53624.74 51517.83 55817.83
2027 Q2 53933.38 51494.54 56487.74
2027 Q3 52676.11 49795.93 55722.88
2027 Q4 57903.65 54279.42 61769.87
2028 Q1 54353.20 50574.87 58413.80
2028 Q2 54666.03 50525.15 59146.29
2028 Q3 53391.68 48861.70 58341.64
2028 Q4 58690.23 53237.19 64701.82
2029 Q1 55091.55 49570.05 61228.08
2029 Q2 55408.63 49481.81 62045.35

Cómo citar

BibTeX
@online{kunst2026,
  author = {Kunst, Joshua},
  title = {Del PIB original al ajuste estacional con X-13},
  date = {2026-09-01},
  url = {https://jkunst.com/blog/posts/2026-09-01-desestacionalizar-pib-x13/},
  langid = {es}
}
Por favor, cita este trabajo como:
Kunst, Joshua. 2026. “Del PIB original al ajuste estacional con X-13.” September 1. https://jkunst.com/blog/posts/2026-09-01-desestacionalizar-pib-x13/.