1 Introducción

Este documento presenta un análisis completo de la serie temporal del PIB trimestral utilizando modelos SARIMA (Seasonal Autoregressive Integrated Moving Average). El objetivo es identificar patrones estacionales, modelar la serie y generar pronósticos confiables.

1.1 Objetivos

  • Analizar la estacionariedad de la serie temporal del PIB
  • Identificar componentes estacionales y de tendencia
  • Seleccionar el mejor modelo SARIMA
  • Generar pronósticos para los próximos 8 trimestres
library(readxl)
library(tidyverse)
library(forecast)
library(tseries)
library(ggplot2)
library(urca)
library(gridExtra)
library(kableExtra)

2 Carga y Exploración de Datos

2.1 Carga de Datos

serie <- read_excel("macro.xlsx")

# Explorar la estructura de los datos
cat("Estructura de los datos:\n")
## Estructura de los datos:
str(serie)
## tibble [80 × 4] (S3: tbl_df/tbl/data.frame)
##  $ Periodo  : num [1:80] 2005 2005 2005 2005 2006 ...
##  $ Trimestre: num [1:80] 1 2 3 4 1 2 3 4 1 2 ...
##  $ Fecha    : chr [1:80] "2005-TI" "2005-TII" "2005-TIII" "2005-TIV" ...
##  $ PIB      : num [1:80] 121053 125547 129214 139039 128776 ...

2.2 Exploración Inicial

# Mostrar primeras y últimas observaciones
kable(head(serie, 10), caption = "Primeras 10 observaciones") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Primeras 10 observaciones
Periodo Trimestre Fecha PIB
2005 1 2005-TI 121053.1
2005 2 2005-TII 125546.7
2005 3 2005-TIII 129214.5
2005 4 2005-TIV 139038.8
2006 1 2006-TI 128775.7
2006 2 2006-TII 132230.6
2006 3 2006-TIII 139102.8
2006 4 2006-TIV 149325.8
2007 1 2007-TI 137419.1
2007 2 2007-TII 141303.0 
kable(tail(serie, 10), caption = "Últimas 10 observaciones") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Últimas 10 observaciones
Periodo Trimestre Fecha PIB
2022 3 2022-TIII 247207.2
2022 4 2022-TIV 257033.8
2023 1 2023-TI 235531.0
2023 2 2023-TII 240292.1
2023 3 2023-TIII 245602.4
2023 4 2023-TIV 258158.5
2024 1 2024-TI 236922.8
2024 2 2024-TII 244915.2
2024 3 2024-TIII 250697.4
2024 4 2024-TIV 264130.2 
# Resumen estadístico
summary(serie)
##     Periodo       Trimestre       Fecha                PIB        
##  Min.   :2005   Min.   :1.00   Length:80          Min.   :121053  
##  1st Qu.:2010   1st Qu.:1.75   Class :character   1st Qu.:159427  
##  Median :2014   Median :2.50   Mode  :character   Median :196212  
##  Mean   :2014   Mean   :2.50                      Mean   :191274  
##  3rd Qu.:2019   3rd Qu.:3.25                      3rd Qu.:216727  
##  Max.   :2024   Max.   :4.00                      Max.   :264130

Este análisis estadístico descriptivo de la serie temporal del PIB trimestral revela características fundamentales de los datos que serán cruciales para el modelado SARIMA posterior. La serie comprende 80 observaciones trimestrales que abarcan desde 2005 hasta 2024, proporcionando una ventana temporal de 20 años que es suficientemente amplia para capturar patrones estacionales y tendencias de largo plazo en la actividad económica.

El PIB presenta una distribución con un valor mínimo de 121,053 millones y un máximo de 264,130 millones, evidenciando un crecimiento económico sustancial durante el período analizado. La media de 191,274 millones es ligeramente inferior a la mediana de 196,212 millones, sugiriendo una distribución con una ligera asimetría negativa, lo que indica que los valores más bajos del PIB en los primeros años del período ejercen influencia en el promedio general. Esta característica es típica en series económicas de países en crecimiento, donde los valores iniciales son considerablemente menores que los más recientes.

La diferencia entre el primer cuartil (159,427 millones) y el tercer cuartil (216,727 millones) muestra una dispersión considerable de 57,300 millones, reflejando la variabilidad natural del crecimiento económico a lo largo de dos décadas. Esta variabilidad será un elemento clave a considerar en el modelado SARIMA, ya que puede indicar la presencia de heteroscedasticidad que podría requerir transformaciones logarítmicas para estabilizar la varianza. El rango completo de los datos (143,077 millones) confirma la presencia de una tendencia creciente pronunciada que deberá ser tratada mediante diferenciación en el proceso de modelado para lograr estacionariedad.

2.3 Creación de la Serie Temporal

# Crear la serie temporal con frecuencia trimestral
ts_pib <- ts(serie$PIB, 
             start = c(min(serie$Periodo), 1), 
             frequency = 4)

cat("Información de la serie temporal:\n")
## Información de la serie temporal:
cat("Inicio:", paste(start(ts_pib), collapse = "-"), "\n")
## Inicio: 2005-1
cat("Final:", paste(end(ts_pib), collapse = "-"), "\n")
## Final: 2024-4
cat("Frecuencia:", frequency(ts_pib), "(trimestral)\n")
## Frecuencia: 4 (trimestral)
cat("Longitud:", length(ts_pib), "observaciones\n")
## Longitud: 80 observaciones

La configuración exitosa de la serie temporal confirma la estructura de datos apropiada para el análisis SARIMA. El período de 20 años (2005-2024) con 80 observaciones trimestrales proporciona una muestra robusta que permite identificar patrones estacionales de 4 trimestres y estimar modelos con componentes estacionales confiables. La frecuencia trimestral es óptima para capturar ciclos económicos y patrones estacionales típicos en variables macroeconómicas como el PIB, donde factores como efectos calendario, patrones de consumo estacional y ciclos de producción se manifiestan claramente cada cuatro trimestres

2.4 Visualización Inicial

# Visualización inicial
par(mfrow = c(2, 2))

# 1. Serie temporal completa
plot(ts_pib, main = "PIB Trimestral", 
     xlab = "Año", ylab = "PIB (millones)", 
     col = "blue", lwd = 2, type = "o")
grid()

# 2. Serie en logaritmos
log_pib <- log(ts_pib)
plot(log_pib, main = "Log(PIB) Trimestral", 
     xlab = "Año", ylab = "Log(PIB)", 
     col = "red", lwd = 2, type = "o")
grid()

# 3. Tasa de crecimiento trimestral
growth_rate <- diff(log_pib) * 100
plot(growth_rate, main = "Tasa de Crecimiento Trimestral (%)", 
     xlab = "Año", ylab = "Crecimiento %", 
     col = "green", lwd = 2, type = "o")
abline(h = 0, col = "gray", lty = 2)
grid()

# 4. Boxplot por trimestre
boxplot(PIB ~ Trimestre, data = serie, 
        main = "PIB por Trimestre",
        xlab = "Trimestre", ylab = "PIB",
        col = c("lightblue", "lightgreen", "lightyellow", "lightcoral"))

par(mfrow = c(1, 1))

La visualización múltiple revela patrones cruciales en la serie del PIB trimestral. El gráfico superior izquierdo muestra una clara tendencia creciente con fluctuaciones estacionales regulares y una notable caída en 2020, evidenciando el impacto de la pandemia. La transformación logarítmica (superior derecha) linealiza la tendencia y estabiliza la varianza, confirmando que esta transformación será apropiada para el modelado SARIMA.

La tasa de crecimiento trimestral (inferior izquierda) exhibe un patrón estacional marcado con alternancia regular entre valores positivos y negativos, y una volatilidad extrema durante 2020. Este patrón sugiere la necesidad de componentes estacionales en el modelo. El boxplot por trimestre (inferior derecha) confirma la estacionalidad, mostrando que el cuarto trimestre sistemáticamente presenta valores más altos, mientras que el primer trimestre tiende a ser más bajo

2.5 Análisis Estacional Preliminar

# Estadísticas por trimestre
medias_trim <- tapply(serie$PIB, serie$Trimestre, mean, na.rm = TRUE)
sd_trim <- tapply(serie$PIB, serie$Trimestre, sd, na.rm = TRUE)
cv_trim <- sd_trim / medias_trim * 100

stats_preliminar <- data.frame(
  Trimestre = 1:4,
  Media = round(medias_trim, 2),
  Desv_Std = round(sd_trim, 2),
  CV_Porcentaje = round(cv_trim, 2)
)

kable(stats_preliminar, caption = "Estadísticas por Trimestre") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Estadísticas por Trimestre
Trimestre Media Desv_Std CV_Porcentaje
1 181777.6 35778.87 19.68
2 185537.6 36151.79 19.48
3 192882.6 37171.96 19.27
4 204897.9 38726.31 18.90

Observación: El trimestre 4 muestra consistentemente los valores más altos del PIB, indicando un fuerte patrón estacional.

Las estadísticas por trimestre confirman cuantitativamente la presencia de un patrón estacional robusto y consistente en el PIB. El trimestre 4 presenta la media más alta (204,897.9 millones), superando al trimestre 1 por aproximadamente 23,120 millones, lo que representa una diferencia del 12.7% que es económicamente significativa. Esta progresión ascendente desde el primer trimestre (181,777.6) hasta el cuarto trimestre evidencia un patrón estacional aditivo donde cada trimestre sucesivo tiende a mostrar mayor actividad económica.

Los coeficientes de variación son notablemente similares entre trimestres (oscilando entre 18.90% y 19.68%), indicando que la variabilidad relativa es estable a lo largo del año. Esta homogeneidad en la dispersión relativa sugiere que la estacionalidad es predominantemente aditiva rather than multiplicativa, lo que simplifica la especificación del modelo SARIMA. La consistencia de estos patrones a través de 20 años de datos proporciona confianza en que la estacionalidad es estructural y predecible, justificando la inclusión de componentes estacionales en el modelo final.

3 Descomposición y Análisis de Componentes

3.1 Descomposición de la Serie

# Descomposición aditiva
decomp_add <- decompose(ts_pib, type = "additive")
plot(decomp_add)
title("Descomposición Aditiva del PIB", line = -1, outer = TRUE)

# Descomposición multiplicativa
decomp_mult <- decompose(ts_pib, type = "multiplicative")
plot(decomp_mult)
title("Descomposición Multiplicativa del PIB", line = -1, outer = TRUE)

La descomposición de series temporales es una técnica metodológica fundamental que permite separar los componentes no observables de tendencia, estacionalidad y ruido aleatorio para evaluar cuál modelo (aditivo o multiplicativo) describe mejor la estructura de los datos. Esta separación es crucial para la especificación correcta del modelo SARIMA, ya que determina si se requieren transformaciones logarítmicas y cómo interpretar los componentes estacionales.

La descomposición aditiva muestra una tendencia claramente creciente y lineal hasta 2020, con una caída pronunciada durante la pandemia seguida de recuperación. El componente estacional exhibe un patrón regular y estable con amplitud constante a lo largo del tiempo, confirmando que la estacionalidad es aditiva. El componente aleatorio presenta variabilidad relativamente constante, excepto por el shock de 2020, sugiriendo homocedasticidad general.

La descomposición multiplicativa revela patrones similares, pero el componente estacional muestra amplitud proporcional al nivel de la serie, mientras que el componente aleatorio presenta mayor estabilidad en términos relativos. La comparación entre ambas descomposiciones indica que el modelo aditivo es más apropiado para esta serie, ya que la amplitud estacional permanece relativamente constante en términos absolutos. Esta conclusión metodológica sugiere que no será necesario aplicar transformaciones logarítmicas únicamente por razones de estacionalidad multiplicativa, aunque podrían considerarse para estabilizar la varianza de la tendencia.

3.2 Análisis Estacional Detallado

# Gráficos estacionales
par(mfrow = c(2, 2))

# Patrón estacional por trimestre
monthplot(ts_pib, main = "Patrón Estacional por Trimestre", 
          xlab = "Trimestre", ylab = "PIB")

# Gráfico de subseries estacionales
seasonplot(ts_pib, main = "PIB por Trimestre y Año", 
           col = rainbow(20), year.labels = TRUE, 
           xlab = "Trimestre", ylab = "PIB")

# Distribución por trimestre
boxplot(PIB ~ Trimestre, data = serie, 
        main = "Distribución del PIB por Trimestre",
        xlab = "Trimestre", ylab = "PIB",
        col = c("lightblue", "lightgreen", "lightyellow", "lightcoral"))

# Estadísticas detalladas por trimestre
stats_trimestre <- serie %>%
  group_by(Trimestre) %>%
  summarise(
    Media = mean(PIB, na.rm = TRUE),
    Mediana = median(PIB, na.rm = TRUE),
    Desv_Std = sd(PIB, na.rm = TRUE),
    Min = min(PIB, na.rm = TRUE),
    Max = max(PIB, na.rm = TRUE),
    CV = (Desv_Std/Media)*100,
    .groups = 'drop'
  )

plot(stats_trimestre$Trimestre, stats_trimestre$Media, 
     type = "b", pch = 19, col = "red", lwd = 2,
     main = "Media del PIB por Trimestre",
     xlab = "Trimestre", ylab = "PIB Promedio")
grid()

par(mfrow = c(1, 1))

El análisis estacional detallado profundiza en la caracterización de los patrones cíclicos mediante múltiples perspectivas visuales que son esenciales para validar la especificación de componentes estacionales en el modelo SARIMA. El patrón estacional por trimestre (superior izquierdo) confirma la consistencia del ciclo anual, mostrando una progresión ascendente sistemática desde Q1 hasta Q4, lo que valida la necesidad de incluir términos estacionales autorregresivos y de media móvil en la especificación del modelo.

El gráfico de subseries estacionales (superior derecho) revela que cada trimestre mantiene su posición relativa a lo largo de todos los años, con Q4 consistentemente superior y Q1 inferior, confirmando la estabilidad temporal del patrón estacional. Esta estabilidad es crucial para la validez de los modelos SARIMA, que asumen patrones estacionales constantes. La distribución por trimestre (inferior izquierdo) muestra que Q4 presenta mayor variabilidad y valores atípicos superiores, mientras que Q1 es más compacto, información relevante para evaluar la homocedasticidad estacional.

La progresión lineal de medias trimestrales (inferior derecho) cuantifica la magnitud del efecto estacional, mostrando un incremento aproximadamente lineal de 23,000 millones entre Q1 y Q4. Esta linearidad sugiere que un modelo SARIMA con componentes estacionales de primer orden será suficiente para capturar la estructura estacional, sin necesidad de órdenes estacionales más complejos que podrían introducir sobreajuste en el modelo final.

3.3 Análisis de Tendencia

# Modelo de tendencia lineal
tiempo <- 1:length(ts_pib)
modelo_tendencia <- lm(as.vector(ts_pib) ~ tiempo)

# Gráfico con tendencia
plot(ts_pib, main = "PIB con Tendencia Lineal", 
     xlab = "Año", ylab = "PIB", col = "blue", lwd = 2)
lines(ts(fitted(modelo_tendencia), start = c(2005, 1), frequency = 4), 
      col = "red", lwd = 2, lty = 2)
legend("topleft", legend = c("PIB Original", "Tendencia Lineal"), 
       col = c("blue", "red"), lwd = 2, lty = c(1, 2))
grid()

# Resumen del modelo de tendencia
summary(modelo_tendencia)
## 
## Call:
## lm(formula = as.vector(ts_pib) ~ tiempo)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -44945  -5281   -349   5723  18134 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 128589.55    2276.43   56.49   <2e-16 ***
## tiempo        1547.76      48.83   31.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10090 on 78 degrees of freedom
## Multiple R-squared:  0.928,  Adjusted R-squared:  0.927 
## F-statistic:  1005 on 1 and 78 DF,  p-value: < 2.2e-16

El análisis de tendencia mediante regresión lineal es metodológicamente fundamental para cuantificar el componente de largo plazo y determinar el orden de integración necesario en el modelo SARIMA. El modelo de tendencia lineal muestra un ajuste excepcional con R² = 0.928, indicando que el 92.8% de la variabilidad del PIB se explica por una tendencia temporal lineal. El coeficiente de tendencia de 1,547.76 millones por trimestre es altamente significativo (t = 31.70, p < 2e-16), confirmando un crecimiento promedio de aproximadamente 6,191 millones anuales.

La visualización superpuesta revela que la tendencia lineal captura efectivamente el comportamiento de largo plazo, pero los residuos muestran patrones sistemáticos que incluyen fluctuaciones estacionales y desviaciones cíclicas, particularmente la caída pronunciada en 2020. Esta evidencia confirma que la serie es no estacionaria debido a la presencia de una tendencia determinística fuerte, justificando la necesidad de diferenciación (d = 1) en la especificación SARIMA.

El error estándar residual de 10,090 millones y la distribución de residuos (con rango de -44,945 a 18,134) indican que, aunque la tendencia lineal es dominante, existe variabilidad sustancial no explicada que incluye componentes estacionales y choques aleatorios. Esta variabilidad residual será capturada por los componentes AR, MA y estacionales del modelo SARIMA, mientras que la diferenciación eliminará la tendencia determinística para lograr estacionariedad.

4 Pruebas de Estacionariedad

4.1 Preparación de Series

# Crear diferentes transformaciones de la serie
ts_pib <- ts(serie$PIB, start = c(2005, 1), frequency = 4)
log_pib <- log(ts_pib)
diff_pib <- diff(ts_pib)
diff_log_pib <- diff(log_pib)
diff_seasonal_pib <- diff(ts_pib, lag = 4)

4.2 Pruebas Augmented Dickey-Fuller (ADF)

cat("=== PRUEBAS AUGMENTED DICKEY-FULLER ===\n\n")
## === PRUEBAS AUGMENTED DICKEY-FULLER ===
# Serie original
adf_original <- ur.df(ts_pib, type = "trend", lags = 8, selectlags = "AIC")
cat("ADF - Serie Original:\n")
## ADF - Serie Original:
summary(adf_original)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -37837  -1582    714   1919  16499 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.219e+04  1.582e+04   3.931 0.000210 ***
## z.lag.1     -4.722e-01  1.233e-01  -3.830 0.000295 ***
## tt           7.114e+02  1.891e+02   3.762 0.000368 ***
## z.diff.lag1  1.045e-01  1.376e-01   0.760 0.450209    
## z.diff.lag2  8.460e-04  1.186e-01   0.007 0.994331    
## z.diff.lag3 -1.103e-01  1.026e-01  -1.075 0.286365    
## z.diff.lag4  6.782e-01  9.422e-02   7.198 8.34e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5968 on 64 degrees of freedom
## Multiple R-squared:  0.7947, Adjusted R-squared:  0.7755 
## F-statistic:  41.3 on 6 and 64 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -3.8296 6.5917 7.3331 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2  6.50  4.88  4.16
## phi3  8.73  6.49  5.47
# Serie en logaritmos
adf_log <- ur.df(log_pib, type = "trend", lags = 8, selectlags = "AIC")
cat("\nADF - Serie en Logaritmos:\n")
## 
## ADF - Serie en Logaritmos:
summary(adf_log)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.190007 -0.011062  0.003692  0.010362  0.098400 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.8849933  1.2770561   3.042  0.00340 ** 
## z.lag.1     -0.3276500  0.1081090  -3.031  0.00352 ** 
## tt           0.0025097  0.0008763   2.864  0.00565 ** 
## z.diff.lag1 -0.0415648  0.1306707  -0.318  0.75145    
## z.diff.lag2 -0.1131532  0.1151149  -0.983  0.32933    
## z.diff.lag3 -0.1836080  0.1028858  -1.785  0.07907 .  
## z.diff.lag4  0.6076844  0.0964267   6.302 3.06e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03092 on 64 degrees of freedom
## Multiple R-squared:  0.7764, Adjusted R-squared:  0.7555 
## F-statistic: 37.05 on 6 and 64 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -3.0307 4.8929 4.7114 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2  6.50  4.88  4.16
## phi3  8.73  6.49  5.47
# Primera diferencia
adf_diff <- ur.df(diff_pib, type = "drift", lags = 8, selectlags = "AIC")
cat("\nADF - Primera Diferencia:\n")
## 
## ADF - Primera Diferencia:
summary(adf_diff)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -37805  -1580     92   1798  18788 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3516.8747  1057.1262   3.327  0.00148 ** 
## z.lag.1       -2.2521     0.4735  -4.756 1.22e-05 ***
## z.diff.lag1    1.0026     0.4274   2.346  0.02220 *  
## z.diff.lag2    0.7090     0.3836   1.848  0.06931 .  
## z.diff.lag3    0.3841     0.3461   1.110  0.27139    
## z.diff.lag4    0.7666     0.2582   2.970  0.00424 ** 
## z.diff.lag5    0.4650     0.1938   2.399  0.01944 *  
## z.diff.lag6    0.2165     0.1260   1.718  0.09081 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6273 on 62 degrees of freedom
## Multiple R-squared:  0.9121, Adjusted R-squared:  0.9022 
## F-statistic: 91.91 on 7 and 62 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -4.7561 11.3108 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1  6.70  4.71  3.86
# Diferencia de logaritmos
adf_diff_log <- ur.df(diff_log_pib, type = "drift", lags = 8, selectlags = "AIC")
cat("\nADF - Diferencia de Logaritmos:\n")
## 
## ADF - Diferencia de Logaritmos:
summary(adf_diff_log)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.192954 -0.007835  0.000044  0.009597  0.083641 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.018171   0.005545   3.277  0.00172 ** 
## z.lag.1     -2.207311   0.475903  -4.638 1.86e-05 ***
## z.diff.lag1  0.938277   0.432170   2.171  0.03376 *  
## z.diff.lag2  0.630426   0.388291   1.624  0.10954    
## z.diff.lag3  0.325286   0.348956   0.932  0.35487    
## z.diff.lag4  0.698396   0.261158   2.674  0.00956 ** 
## z.diff.lag5  0.435067   0.192719   2.258  0.02751 *  
## z.diff.lag6  0.214952   0.122729   1.751  0.08482 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03181 on 62 degrees of freedom
## Multiple R-squared:  0.9091, Adjusted R-squared:  0.8988 
## F-statistic: 88.58 on 7 and 62 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -4.6381 10.7662 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1  6.70  4.71  3.86

Las pruebas Augmented Dickey-Fuller constituyen el método estándar para evaluar formalmente la presencia de raíces unitarias y determinar el orden de integración necesario para la especificación SARIMA. Los resultados proporcionan evidencia concluyente sobre la naturaleza no estacionaria de la serie original y la efectividad de la diferenciación para lograr estacionariedad.

Serie Original: El estadístico ADF de -3.8296 no supera el valor crítico de -4.04 al 1%, indicando que no se puede rechazar la hipótesis nula de raíz unitaria. Esto confirma que la serie en niveles es no estacionaria, como se esperaba por la presencia de tendencia creciente observada en el análisis gráfico previo.

Serie en Logaritmos: El estadístico de -3.0307 es aún menos negativo que el de la serie original, confirmando que la transformación logarítmica por sí sola no elimina la no estacionariedad. Aunque estabiliza la varianza, la tendencia determinística persiste, requiriendo diferenciación adicional.

Primera Diferencia: El estadístico ADF de -4.7561 supera ampliamente el valor crítico de -3.51 al 1%, rechazando contundentemente la hipótesis nula de raíz unitaria. Esto confirma que la primera diferencia es estacionaria, estableciendo que d = 1 en la especificación SARIMA.

Diferencia de Logaritmos: Con un estadístico de -4.6381, también rechaza la hipótesis nula al 1%, confirmando estacionariedad. Esta transformación combina estabilización de varianza y eliminación de tendencia, siendo la más apropiada para el modelado SARIMA posterior. Los resultados validan que la especificación requerirá d = 1 y que la transformación logarítmica será beneficiosa para la estabilidad del modelo.

4.3 Pruebas Phillips-Perron

cat("=== PRUEBAS PHILLIPS-PERRON ===\n\n")
## === PRUEBAS PHILLIPS-PERRON ===
# Serie original
pp_original <- ur.pp(ts_pib, type = "Z-tau", model = "trend")
cat("Phillips-Perron - Serie Original:\n")
## Phillips-Perron - Serie Original:
summary(pp_original)
## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept and trend 
## 
## 
## Call:
## lm(formula = y ~ y.l1 + trend)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -41350  -6003    783   5462  18290 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.424e+05  2.112e+04   6.743 2.67e-09 ***
## y.l1        2.583e-01  1.113e-01   2.322   0.0229 *  
## trend       1.142e+03  1.777e+02   6.428 1.03e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9819 on 76 degrees of freedom
## Multiple R-squared:  0.9303, Adjusted R-squared:  0.9285 
## F-statistic: 507.3 on 2 and 76 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic, type: Z-tau  is: -6.5132 
## 
##            aux. Z statistics
## Z-tau-mu              8.2373
## Z-tau-beta            6.1877
## 
## Critical values for Z statistics: 
##                      1pct      5pct     10pct
## critical values -4.077136 -3.466583 -3.159722
# Primera diferencia
pp_diff <- ur.pp(diff_pib, type = "Z-tau", model = "constant")
cat("\nPhillips-Perron - Primera Diferencia:\n")
## 
## Phillips-Perron - Primera Diferencia:
summary(pp_diff)
## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -38300  -8701   3047   8811  27390 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 2216.3534  1362.4827   1.627    0.108  
## y.l1          -0.2645     0.1113  -2.377    0.020 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11920 on 76 degrees of freedom
## Multiple R-squared:  0.06918,    Adjusted R-squared:  0.05694 
## F-statistic: 5.649 on 1 and 76 DF,  p-value: 0.01998
## 
## 
## Value of test-statistic, type: Z-tau  is: -17.5223 
## 
##          aux. Z statistics
## Z-tau-mu            2.5014
## 
## Critical values for Z statistics: 
##                      1pct      5pct     10pct
## critical values -3.515218 -2.898577 -2.586272
# Diferencia de logaritmos
pp_diff_log <- ur.pp(diff_log_pib, type = "Z-tau", model = "constant")
cat("\nPhillips-Perron - Diferencia de Logaritmos:\n")
## 
## Phillips-Perron - Diferencia de Logaritmos:
summary(pp_diff_log)
## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19489 -0.04248  0.01718  0.04530  0.13060 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.01209    0.00692   1.747   0.0847 .
## y.l1        -0.27349    0.11055  -2.474   0.0156 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06043 on 76 degrees of freedom
## Multiple R-squared:  0.07452,    Adjusted R-squared:  0.06235 
## F-statistic:  6.12 on 1 and 76 DF,  p-value: 0.0156
## 
## 
## Value of test-statistic, type: Z-tau  is: -18.0655 
## 
##          aux. Z statistics
## Z-tau-mu             2.737
## 
## Critical values for Z statistics: 
##                      1pct      5pct     10pct
## critical values -3.515218 -2.898577 -2.586272

Las pruebas Phillips-Perron complementan el análisis de estacionariedad mediante una metodología robusta ante heterocedasticidad y autocorrelación, utilizando correcciones no paramétricas que las hacen menos sensibles a la especificación del número de rezagos. Estos resultados refuerzan las conclusiones obtenidas con las pruebas ADF y proporcionan mayor confianza en la especificación del orden de integración.

Serie Original: El estadístico Z-tau de -6.5132 supera el valor crítico de -4.077 al 1%, lo que contradice el resultado de la prueba ADF. Esta discrepancia metodológica es interesante y sugiere que la no estacionariedad detectada por ADF podría estar influenciada por la estructura de autocorrelación específica. Sin embargo, la presencia visual clara de tendencia y el consenso general en la literatura económica sugieren mantener cautela y proceder con diferenciación.

Primera Diferencia: El estadístico Z-tau de -17.5223 es extremadamente significativo, superando ampliamente el valor crítico de -3.515 al 1%. Este resultado contundente confirma que la primera diferencia es fuertemente estacionaria, validando la especificación d = 1 para el modelo SARIMA.

Diferencia de Logaritmos: Con un estadístico Z-tau de -18.0655, también rechaza la hipótesis nula al 1% con mayor contundencia que la primera diferencia simple. Esta evidencia refuerza la decisión de utilizar la transformación logarítmica y la primera diferencia para el modelado, ya que proporciona la serie más claramente estacionaria según ambas metodologías de prueba.

La concordancia entre las pruebas ADF y Phillips-Perron para las series diferenciadas proporciona evidencia robusta de que d = 1 es la especificación apropiada, mientras que la ligera superioridad de la diferencia logarítmica sugiere que esta transformación será óptima para el modelado SARIMA, combinando estabilización de varianza y eliminación de tendencia.

4.4 Resumen de Estacionariedad

Conclusiones de las pruebas:

  • Serie Original: No estacionaria (requiere diferenciación)
  • Primera Diferencia: Estacionaria (d = 1)
  • Diferencia de Logaritmos: Estacionaria (recomendada para modelado)
  • Diferenciación Estacional: No necesaria (D = 0)

5 Análisis de Funciones de Autocorrelación

5.1 Gráficos ACF y PACF

# Trabajar con la diferencia de logaritmos
diff_log_pib <- diff(log_pib)

# Gráficos de ACF y PACF
par(mfrow = c(2, 3))

# ACF y PACF de la serie original
acf(ts_pib, lag.max = 20, main = "ACF - Serie Original")
pacf(ts_pib, lag.max = 20, main = "PACF - Serie Original")

# ACF y PACF de la serie diferenciada
acf(diff_log_pib, lag.max = 20, main = "ACF - Diff(Log(PIB))")
pacf(diff_log_pib, lag.max = 20, main = "PACF - Diff(Log(PIB))")

# ACF y PACF con más lags para detectar estacionalidad
acf(diff_log_pib, lag.max = 40, main = "ACF - 40 lags")
pacf(diff_log_pib, lag.max = 40, main = "PACF - 40 lags")

par(mfrow = c(1, 1))

El análisis de las funciones de autocorrelación constituye el núcleo metodológico para la identificación de modelos ARIMA, aplicando los principios de Box-Jenkins para determinar los órdenes apropiados de los componentes autorregresivos y de media móvil. La comparación entre la serie original y la serie diferenciada revela claramente el impacto de la transformación en la estructura de correlación temporal.

Serie Original: La ACF muestra el patrón típico de no estacionariedad con decaimiento lento y persistente, confirmando la presencia de tendencia. La PACF presenta valores significativos en los primeros lags, pero la interpretación es problemática debido a la no estacionariedad subyacente, validando la necesidad de diferenciación antes del análisis de identificación.

Serie Diferenciada (Diff(Log(PIB))): La ACF muestra un corte abrupto después del lag 1, sugiriendo un componente MA(1), mientras que presenta valores significativos en lags estacionales (4, 8). La PACF exhibe significancia en el lag 1 y decaimiento posterior, indicando un posible componente AR(1). Crucialmente, la PACF muestra significancia en el lag 4, sugiriendo un componente SAR(1).

Análisis Extendido (40 lags): La extensión del análisis confirma la presencia de correlaciones estacionales significativas en múltiplos de 4 (lags 4, 8, 12), validando la necesidad de componentes estacionales. El patrón sugiere que un modelo SARIMA(1,1,1)(1,0,1)4 podría ser apropiado, combinando componentes regulares AR(1) y MA(1) con componentes estacionales SAR(1) y SMA(1). Esta identificación preliminar guiará la especificación de modelos candidatos en la fase de estimación.

5.2 Análisis Detallado de Correlaciones

# Calcular ACF y PACF
acf_values <- acf(diff_log_pib, lag.max = 20, plot = FALSE)
pacf_values <- pacf(diff_log_pib, lag.max = 20, plot = FALSE)

# Límite de significancia
n <- length(diff_log_pib)
sig_level <- 1.96/sqrt(n)

cat("Límite de significancia (95%):", round(sig_level, 4), "\n")
## Límite de significancia (95%): 0.2205
cat("Número de observaciones:", n, "\n\n")
## Número de observaciones: 79
# Tabla de ACF
acf_table <- data.frame(
  Lag = 1:12,
  ACF = round(acf_values$acf[2:13], 4),
  Significativo = abs(acf_values$acf[2:13]) > sig_level,
  Estacional = c(rep("No", 3), "Sí", rep("No", 3), "Sí", rep("No", 3), "Sí")
)

kable(acf_table, caption = "Función de Autocorrelación (ACF)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Función de Autocorrelación (ACF)
Lag ACF Significativo Estacional
1 -0.2718 TRUE No
2 -0.3196 TRUE No
3 -0.2661 TRUE No
4 0.8096 TRUE Sí
5 -0.2590 TRUE No
6 -0.2977 TRUE No
7 -0.2409 TRUE No
8 0.7213 TRUE Sí
9 -0.2355 TRUE No
10 -0.2482 TRUE No
11 -0.2255 TRUE No
12 0.6896 TRUE Sí 
# Tabla de PACF
pacf_table <- data.frame(
  Lag = 1:12,
  PACF = round(pacf_values$acf[1:12], 4),
  Significativo = abs(pacf_values$acf[1:12]) > sig_level,
  Estacional = c(rep("No", 3), "Sí", rep("No", 3), "Sí", rep("No", 3), "Sí")
)

kable(pacf_table, caption = "Función de Autocorrelación Parcial (PACF)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Función de Autocorrelación Parcial (PACF)
Lag PACF Significativo Estacional
1 -0.2718 TRUE No
2 -0.4249 TRUE No
3 -0.6658 TRUE No
4 0.5494 TRUE Sí
5 -0.1366 FALSE No
6 -0.1393 FALSE No
7 -0.1514 FALSE No
8 0.0210 FALSE Sí
9 -0.1022 FALSE No
10 -0.0288 FALSE No
11 -0.1173 FALSE No
12 0.1036 FALSE Sí

La cuantificación precisa de las correlaciones mediante tablas numéricas permite aplicar sistemáticamente las reglas de identificación de Box-Jenkins, utilizando el límite de significancia de 0.2205 basado en el tamaño muestral de 79 observaciones. Esta metodología objetiva es fundamental para determinar los órdenes específicos de los componentes ARIMA y validar las interpretaciones visuales previas.

Función de Autocorrelación (ACF): Los resultados revelan un patrón mixto con correlaciones negativas significativas en los lags 1, 2, 3, 5, 6, 7, 9, 10 y 11, alternando con correlaciones positivas significativas en los lags estacionales 4 y 8. Este patrón de alternancia sugiere sobrediferenciación o la presencia de componentes de media móvil complejos. Las correlaciones estacionales positivas y significativas (0.8096 en lag 4, 0.7213 en lag 8) confirman fuertemente la necesidad de componentes estacionales SMA.

Función de Autocorrelación Parcial (PACF): La PACF muestra significancia en los lags 1, 2, 3 y 4, con el lag 4 siendo estacional. Después del lag 4, todas las correlaciones parciales se vuelven no significativas, sugiriendo un corte abrupto que es característico de procesos autorregresivos. El valor significativo en el lag 4 (0.5494) indica claramente un componente SAR(1), mientras que la significancia en los primeros tres lags sugiere la necesidad de componentes AR de orden superior o una combinación AR-MA.

La interpretación conjunta sugiere un modelo SARIMA con componentes regulares mixtos (posiblemente ARMA(2,1) o ARMA(3,1)) y componentes estacionales que incluyen tanto SAR(1) como SMA(1). El patrón complejo en los lags no estacionales requiere evaluación cuidadosa de múltiples especificaciones candidatas para determinar la combinación óptima que capture toda la estructura de correlación sin sobreajustar el modelo.

5.3 Análisis de Lags Estacionales

# Lags estacionales (4, 8, 12, 16)
seasonal_lags <- c(4, 8, 12, 16)

seasonal_acf <- data.frame(
  Lag = seasonal_lags,
  ACF = round(acf_values$acf[seasonal_lags + 1], 4),
  Significativo = abs(acf_values$acf[seasonal_lags + 1]) > sig_level
)

seasonal_pacf <- data.frame(
  Lag = seasonal_lags,
  PACF = round(pacf_values$acf[seasonal_lags], 4),
  Significativo = abs(pacf_values$acf[seasonal_lags]) > sig_level
)

kable(seasonal_acf, caption = "ACF en Lags Estacionales") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
ACF en Lags Estacionales
Lag ACF Significativo
4 0.8096 TRUE
8 0.7213 TRUE
12 0.6896 TRUE
16 0.6434 TRUE 
kable(seasonal_pacf, caption = "PACF en Lags Estacionales") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
PACF en Lags Estacionales
Lag PACF Significativo
4 0.5494 TRUE
8 0.0210 FALSE
12 0.1036 FALSE
16 0.0643 FALSE

El análisis específico de los lags estacionales es metodológicamente crucial para identificar los órdenes estacionales (P, Q) del modelo SARIMA de manera independiente de los componentes regulares. Esta separación permite aplicar las reglas de Box-Jenkins específicamente a la estructura estacional, que opera en múltiplos de la frecuencia temporal (cada 4 trimestres).

ACF Estacional: Los resultados muestran correlaciones positivas y altamente significativas en todos los lags estacionales analizados, con un patrón de decaimiento gradual desde 0.8096 (lag 4) hasta 0.6434 (lag 16). Este decaimiento lento y persistente es característico de procesos estacionales autorregresivos, sugiriendo que los valores del PIB en un trimestre específico están fuertemente correlacionados con los valores del mismo trimestre en años anteriores.

PACF Estacional: La función de autocorrelación parcial estacional muestra un patrón claramente diferenciado, con significancia únicamente en el lag 4 (0.5494) y valores no significativos en los lags 8, 12 y 16. Este corte abrupto después del primer lag estacional es la evidencia definitiva de un proceso SAR(1), donde la dependencia estacional se captura completamente con un solo término autorregresivo estacional.

La interpretación conjunta de ambas funciones confirma inequívocamente la presencia de un componente SAR(1) en la especificación del modelo. El decaimiento gradual en la ACF estacional combinado con el corte abrupto en la PACF estacional después del lag 4 es el patrón textbook para un proceso autorregresivo estacional de primer orden. Esta evidencia sugiere que la especificación estacional del modelo SARIMA debe incluir P=1, mientras que la necesidad de un componente SMA(1) deberá evaluarse en conjunto con los componentes regulares durante la fase de estimación y comparación de modelos.

Interpretación:

  • ACF: Significativo en lags estacionales (4, 8, 12, 16) con decaimiento lento
  • PACF: Significativo principalmente en lag 4, sugiere SAR(1)
  • Componente regular: PACF corta después del lag 3, sugiere AR(3) o AR(2)

6 Estimación y Selección de Modelos SARIMA

6.1 Modelos Candidatos

# Definir modelos candidatos
models_candidates <- list(
  "SARIMA(1,1,0)(0,0,0)4" = list(order = c(1,1,0), seasonal = list(order = c(0,0,0), period = 4)),
  "SARIMA(2,1,0)(0,0,0)4" = list(order = c(2,1,0), seasonal = list(order = c(0,0,0), period = 4)),
  "SARIMA(3,1,0)(0,0,0)4" = list(order = c(3,1,0), seasonal = list(order = c(0,0,0), period = 4)),
  "SARIMA(0,1,1)(0,0,0)4" = list(order = c(0,1,1), seasonal = list(order = c(0,0,0), period = 4)),
  "SARIMA(0,1,2)(0,0,0)4" = list(order = c(0,1,2), seasonal = list(order = c(0,0,0), period = 4)),
  "SARIMA(1,1,0)(1,0,0)4" = list(order = c(1,1,0), seasonal = list(order = c(1,0,0), period = 4)),
  "SARIMA(2,1,0)(1,0,0)4" = list(order = c(2,1,0), seasonal = list(order = c(1,0,0), period = 4)),
  "SARIMA(3,1,0)(1,0,0)4" = list(order = c(3,1,0), seasonal = list(order = c(1,0,0), period = 4)),
  "SARIMA(0,1,1)(0,0,1)4" = list(order = c(0,1,1), seasonal = list(order = c(0,0,1), period = 4)),
  "SARIMA(1,1,1)(1,0,0)4" = list(order = c(1,1,1), seasonal = list(order = c(1,0,0), period = 4)),
  "SARIMA(1,1,1)(1,0,1)4" = list(order = c(1,1,1), seasonal = list(order = c(1,0,1), period = 4)),
  "SARIMA(2,1,1)(1,0,1)4" = list(order = c(2,1,1), seasonal = list(order = c(1,0,1), period = 4)),
  "SARIMA(3,1,1)(1,0,1)4" = list(order = c(3,1,1), seasonal = list(order = c(1,0,1), period = 4))
)

cat("Modelos candidatos a evaluar:\n")
## Modelos candidatos a evaluar:
for(i in 1:length(models_candidates)) {
  cat(names(models_candidates)[i], "\n")
}
## SARIMA(1,1,0)(0,0,0)4 
## SARIMA(2,1,0)(0,0,0)4 
## SARIMA(3,1,0)(0,0,0)4 
## SARIMA(0,1,1)(0,0,0)4 
## SARIMA(0,1,2)(0,0,0)4 
## SARIMA(1,1,0)(1,0,0)4 
## SARIMA(2,1,0)(1,0,0)4 
## SARIMA(3,1,0)(1,0,0)4 
## SARIMA(0,1,1)(0,0,1)4 
## SARIMA(1,1,1)(1,0,0)4 
## SARIMA(1,1,1)(1,0,1)4 
## SARIMA(2,1,1)(1,0,1)4 
## SARIMA(3,1,1)(1,0,1)4

6.2 Función de Evaluación

# Función para ajustar y evaluar modelos
fit_and_evaluate <- function(model_spec, data, model_name) {
  tryCatch({
    # Ajustar modelo
    fit <- Arima(data, order = model_spec$order, 
                 seasonal = model_spec$seasonal,
                 method = "ML")
    
    # Calcular métricas
    aic <- fit$aic
    bic <- fit$bic
    aicc <- fit$aicc
    loglik <- fit$loglik
    
    # Calcular RMSE y MAE
    residuals <- residuals(fit)
    rmse <- sqrt(mean(residuals^2, na.rm = TRUE))
    mae <- mean(abs(residuals), na.rm = TRUE)
    
    # Test de Ljung-Box para autocorrelación
    lb_test <- Box.test(residuals, lag = 12, type = "Ljung-Box")
    lb_pvalue <- lb_test$p.value
    
    # Devolver resultados
    return(data.frame(
      Model = model_name,
      AIC = aic,
      BIC = bic,
      AICc = aicc,
      LogLik = loglik,
      RMSE = rmse,
      MAE = mae,
      LB_pvalue = lb_pvalue,
      Residual_OK = lb_pvalue > 0.05
    ))
  }, error = function(e) {
    cat("Error en modelo", model_name, ":", e$message, "\n")
    return(NULL)
  })
}

6.3 Ajuste de Modelos

# Ajustar todos los modelos
results_list <- list()

cat("Ajustando modelos...\n")
## Ajustando modelos...
for(i in 1:length(models_candidates)) {
  model_name <- names(models_candidates)[i]
  cat("Ajustando modelo:", model_name, "\n")
  
  result <- fit_and_evaluate(models_candidates[[i]], log_pib, model_name)
  
  if(!is.null(result)) {
    results_list[[i]] <- result
  }
}
## Ajustando modelo: SARIMA(1,1,0)(0,0,0)4 
## Ajustando modelo: SARIMA(2,1,0)(0,0,0)4 
## Ajustando modelo: SARIMA(3,1,0)(0,0,0)4 
## Ajustando modelo: SARIMA(0,1,1)(0,0,0)4 
## Ajustando modelo: SARIMA(0,1,2)(0,0,0)4 
## Ajustando modelo: SARIMA(1,1,0)(1,0,0)4 
## Ajustando modelo: SARIMA(2,1,0)(1,0,0)4 
## Ajustando modelo: SARIMA(3,1,0)(1,0,0)4 
## Ajustando modelo: SARIMA(0,1,1)(0,0,1)4 
## Ajustando modelo: SARIMA(1,1,1)(1,0,0)4 
## Ajustando modelo: SARIMA(1,1,1)(1,0,1)4 
## Ajustando modelo: SARIMA(2,1,1)(1,0,1)4 
## Ajustando modelo: SARIMA(3,1,1)(1,0,1)4
# Combinar resultados
results_df <- do.call(rbind, results_list)

# Ordenar por AIC
results_sorted_aic <- results_df[order(results_df$AIC), ]

6.4 Comparación de Modelos

kable(results_sorted_aic, digits = 4, caption = "Modelos Ordenados por AIC") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  row_spec(1, bold = TRUE, color = "white", background = "#D7261E")
Modelos Ordenados por AIC
Model AIC BIC AICc LogLik RMSE MAE LB_pvalue Residual_OK
11 SARIMA(1,1,1)(1,0,1)4 -327.2886 -315.4413 -326.4667 168.6443 0.0262 0.0136 0.8959 TRUE
12 SARIMA(2,1,1)(1,0,1)4 -325.4414 -311.2247 -324.2748 168.7207 0.0262 0.0137 0.9008 TRUE
13 SARIMA(3,1,1)(1,0,1)4 -324.2095 -307.6234 -322.6321 169.1048 0.0257 0.0136 0.9095 TRUE
10 SARIMA(1,1,1)(1,0,0)4 -308.8400 -299.3622 -308.2995 158.4200 0.0314 0.0173 0.4825 TRUE
6 SARIMA(1,1,0)(1,0,0)4 -303.1566 -296.0482 -302.8366 154.5783 0.0330 0.0194 0.6333 TRUE
7 SARIMA(2,1,0)(1,0,0)4 -301.5425 -292.0647 -301.0020 154.7713 0.0329 0.0195 0.5841 TRUE
8 SARIMA(3,1,0)(1,0,0)4 -299.8658 -288.0186 -299.0439 154.9329 0.0329 0.0194 0.5539 TRUE
9 SARIMA(0,1,1)(0,0,1)4 -261.2808 -254.1724 -260.9608 133.6404 0.0438 0.0342 0.0000 FALSE
3 SARIMA(3,1,0)(0,0,0)4 -253.0563 -243.5786 -252.5158 130.5282 0.0456 0.0369 0.0000 FALSE
2 SARIMA(2,1,0)(0,0,0)4 -224.6466 -217.5383 -224.3266 115.3233 0.0557 0.0453 0.0000 FALSE
4 SARIMA(0,1,1)(0,0,0)4 -224.1049 -219.3660 -223.9470 114.0525 0.0567 0.0456 0.0000 FALSE
5 SARIMA(0,1,2)(0,0,0)4 -223.5196 -216.4113 -223.1996 114.7598 0.0562 0.0480 0.0000 FALSE
1 SARIMA(1,1,0)(0,0,0)4 -214.7390 -210.0001 -214.5811 109.3695 0.0602 0.0510 0.0000 FALSE

La comparación sistemática de modelos candidatos mediante múltiples criterios de información proporciona evidencia objetiva para la selección del modelo óptimo, balanceando ajuste estadístico con parsimonia. Los resultados revelan una jerarquía clara de desempeño que valida las especificaciones sugeridas por el análisis de autocorrelación previo.

Modelo Óptimo: El SARIMA(1,1,1)(1,0,1)4 emerge como el mejor modelo según AIC (-327.2886), combinando componentes regulares AR(1) y MA(1) con componentes estacionales SAR(1) y SMA(1). Este modelo también presenta el mejor BIC (-315.4413) y AICc (-326.4667), confirmando que la complejidad adicional se justifica estadísticamente. El p-valor de Ljung-Box de 0.8959 indica ausencia de autocorrelación residual, validando que el modelo captura adecuadamente la estructura temporal.

Modelos Alternativos Competitivos: Los modelos SARIMA(2,1,1)(1,0,1)4 y SARIMA(3,1,1)(1,0,1)4 muestran AIC ligeramente superiores (-325.4414 y -324.2095), pero sus BIC más altos (-311.2247 y -307.6234) indican que la complejidad adicional no se justifica según el criterio de parsimonia. Todos mantienen residuos no autocorrelacionados, confirmando especificaciones válidas.

Modelos Sin Componentes Estacionales: Los modelos sin términos estacionales (filas inferiores) presentan AIC sustancialmente peores, confirmando la importancia crítica de los componentes estacionales identificados en el análisis previo. La diferencia de más de 70 puntos en AIC entre el mejor modelo con componentes estacionales versus el mejor sin ellos demuestra que la estacionalidad es fundamental para capturar la estructura de los datos.

La convergencia de criterios hacia el SARIMA(1,1,1)(1,0,1)4, combinada con residuos bien comportados, proporciona confianza en que esta especificación captura óptimamente la dinámica temporal del PIB trimestral sin sobreajustar.

6.5 Mejor Modelo

# Seleccionar el mejor modelo
best_model_name <- results_sorted_aic$Model[1]
best_model_spec <- models_candidates[[best_model_name]]

cat("Mejor modelo según AIC:", best_model_name, "\n\n")
## Mejor modelo según AIC: SARIMA(1,1,1)(1,0,1)4
# Ajustar el mejor modelo
best_model <- Arima(log_pib, 
                    order = best_model_spec$order, 
                    seasonal = best_model_spec$seasonal,
                    method = "ML")

print(summary(best_model))
## Series: log_pib 
## ARIMA(1,1,1)(1,0,1)[4] 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1
##       0.7191  -0.8838  0.9992  -0.9130
## s.e.  0.1745   0.1229  0.0031   0.1632
## 
## sigma^2 = 0.0007307:  log likelihood = 168.64
## AIC=-327.29   AICc=-326.47   BIC=-315.44
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE      MAPE     MASE
## Training set -0.001873423 0.02617374 0.01358272 -0.01517105 0.1118646 0.314605
##                    ACF1
## Training set 0.04035596

El modelo SARIMA(1,1,1)(1,0,1)4 seleccionado representa la especificación óptima que captura tanto la dinámica de corto plazo como los patrones estacionales del PIB trimestral. Los coeficientes estimados proporcionan información económica valiosa sobre la persistencia y los mecanismos de corrección de errores en la serie temporal.

Componentes Regulares: El coeficiente AR(1) de 0.7191 indica una persistencia moderada-alta en las desviaciones de corto plazo, sugiriendo que los choques al PIB tienen efectos que se disipan gradualmente a lo largo de aproximadamente 3-4 trimestres. El coeficiente MA(1) de -0.8838 representa un mecanismo de corrección de errores robusto, donde errores de pronóstico previos se incorporan fuertemente en las predicciones actuales con signo opuesto, estabilizando la serie.

Componentes Estacionales: El coeficiente SAR(1) de 0.9992 es extraordinariamente alto, indicando una persistencia estacional casi perfecta donde el PIB de cada trimestre está casi completamente determinado por el valor del mismo trimestre del año anterior. Esta persistencia extrema refleja la naturaleza estructural de los patrones estacionales en la economía. El coeficiente SMA(1) de -0.9130 proporciona un mecanismo de corrección estacional fuerte, ajustando por errores estacionales previos.

Calidad del Ajuste: La varianza residual (sigma²) de 0.0007307 es muy baja en escala logarítmica, indicando un ajuste excelente. Las métricas de error muestran un MAPE de 0.11%, confirmando precisión excepcional. El ACF1 residual de 0.04 está muy cerca de cero, validando que no hay autocorrelación residual significativa. La combinación de criterios de información favorables (AIC=-327.29, BIC=-315.44) con residuos bien comportados confirma que este modelo logra el balance óptimo entre ajuste y parsimonia para la serie del PIB trimestral.

6.6 Diagnóstico de Residuos

# Diagnóstico de residuos del mejor modelo
par(mfrow = c(2, 2))

# Gráfico de residuos
plot(residuals(best_model), main = paste("Residuos -", best_model_name),
     ylab = "Residuos", xlab = "Tiempo", type = "l")
abline(h = 0, col = "red", lty = 2)

# ACF de residuos
acf(residuals(best_model), main = paste("ACF Residuos -", best_model_name))

# PACF de residuos
pacf(residuals(best_model), main = paste("PACF Residuos -", best_model_name))

# QQ-plot
qqnorm(residuals(best_model), main = paste("Q-Q Plot -", best_model_name))
qqline(residuals(best_model), col = "red")

par(mfrow = c(1, 1))

El diagnóstico de residuos constituye la validación fundamental del modelo SARIMA, verificando que los supuestos estadísticos se cumplan y que el modelo haya capturado adecuadamente toda la estructura temporal de los datos. Los cuatro gráficos proporcionan evidencia complementaria sobre la calidad del ajuste y la validez de las inferencias estadísticas.

Gráfico Temporal de Residuos: Los residuos muestran un comportamiento generalmente aleatorio alrededor de cero, con varianza relativamente constante a lo largo del tiempo. El spike pronunciado en 2020 corresponde al impacto extraordinario de la pandemia, representando un choque exógeno que no puede ser capturado por la estructura ARIMA. Fuera de este evento excepcional, los residuos no muestran patrones sistemáticos, tendencias o cambios evidentes en la varianza, validando los supuestos de homocedasticidad y ausencia de estructura temporal residual.

ACF y PACF de Residuos: Ambas funciones muestran que prácticamente todas las correlaciones se mantienen dentro de las bandas de confianza, confirmando que los residuos se comportan como ruido blanco. Esta evidencia valida que el modelo SARIMA ha capturado exitosamente toda la estructura de autocorrelación temporal y estacional presente en los datos originales, cumpliendo el supuesto fundamental de independencia residual.

Q-Q Plot: El gráfico cuantil-cuantil revela desviaciones de la normalidad en las colas de la distribución, particularmente en los cuantiles extremos donde los residuos se alejan de la línea teórica. Esta no normalidad es común en series económicas y no invalida el modelo para propósitos de pronóstico, aunque sugiere cautela en la interpretación de intervalos de confianza muy estrechos. El comportamiento aproximadamente lineal en la región central confirma que la mayoría de los residuos siguen una distribución cercana a la normal.

6.7 Pruebas de Diagnóstico

cat("=== PRUEBAS DE DIAGNÓSTICO DEL MEJOR MODELO ===\n\n")
## === PRUEBAS DE DIAGNÓSTICO DEL MEJOR MODELO ===
# Test de Ljung-Box
lb_test <- Box.test(residuals(best_model), lag = 12, type = "Ljung-Box")
cat("Test de Ljung-Box (autocorrelación):\n")
## Test de Ljung-Box (autocorrelación):
cat("  Estadístico:", round(lb_test$statistic, 4), "\n")
##   Estadístico: 6.3763
cat("  p-valor:", round(lb_test$p.value, 4), "\n")
##   p-valor: 0.8959
cat("  Resultado:", ifelse(lb_test$p.value > 0.05, 
                           "No hay autocorrelación (Bueno)", 
                           "Hay autocorrelación (Malo)"), "\n\n")
##   Resultado: No hay autocorrelación (Bueno)
# Test de normalidad
if(require(tseries, quietly = TRUE)) {
  jb_test <- jarque.bera.test(residuals(best_model))
  cat("Test de Jarque-Bera (normalidad):\n")
  cat("  Estadístico:", round(jb_test$statistic, 4), "\n")
  cat("  p-valor:", round(jb_test$p.value, 4), "\n")
  cat("  Resultado:", ifelse(jb_test$p.value > 0.05, 
                             "Residuos normales (Bueno)", 
                             "Residuos no normales (Precaución)"), "\n\n")
}
## Test de Jarque-Bera (normalidad):
##   Estadístico: 3456.016 
##   p-valor: 0 
##   Resultado: Residuos no normales (Precaución)

Test de Ljung-Box: El estadístico de 6.3763 con p-valor de 0.8959 proporciona evidencia contundente de que no se puede rechazar la hipótesis nula de ausencia de autocorrelación en los residuos. Este resultado excelente confirma que el modelo SARIMA(1,1,1)(1,0,1)4 ha capturado exitosamente toda la estructura temporal presente en los datos, cumpliendo el supuesto fundamental de que los residuos constituyen ruido blanco. La validación de este supuesto es esencial para la validez de los pronósticos y garantiza que no existe información temporal adicional que el modelo haya omitido.

Test de Jarque-Bera: El estadístico extremadamente alto de 3456.016 con p-valor prácticamente cero rechaza categóricamente la hipótesis nula de normalidad de los residuos. Esta desviación significativa de la normalidad es común en series económicas y refleja la presencia de eventos extremos (como la crisis de 2020) y asimetría en la distribución de los choques económicos. Aunque la no normalidad no invalida el modelo para propósitos de pronóstico puntual, sugiere precaución en la interpretación de intervalos de confianza, especialmente en los extremos de la distribución.

7 Pronósticos y Validación Final

7.1 Modelo Final

# Modelo final
final_model <- Arima(log_pib, 
                     order = c(1,1,1), 
                     seasonal = list(order = c(1,0,1), period = 4),
                     method = "ML")

cat("MODELO FINAL: SARIMA(1,1,1)(1,0,1)4\n\n")
## MODELO FINAL: SARIMA(1,1,1)(1,0,1)4
print(final_model)
## Series: log_pib 
## ARIMA(1,1,1)(1,0,1)[4] 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1
##       0.7191  -0.8838  0.9992  -0.9130
## s.e.  0.1745   0.1229  0.0031   0.1632
## 
## sigma^2 = 0.0007307:  log likelihood = 168.64
## AIC=-327.29   AICc=-326.47   BIC=-315.44

7.2 Pronósticos

# Generar pronósticos para 8 trimestres (2 años)
forecast_periods <- 8
forecasts_log <- forecast(final_model, h = forecast_periods, level = c(80, 95))

# Convertir a escala original
forecasts_original <- forecasts_log
forecasts_original$mean <- exp(forecasts_log$mean)
forecasts_original$lower <- exp(forecasts_log$lower)
forecasts_original$upper <- exp(forecasts_log$upper)
forecasts_original$x <- exp(forecasts_log$x)

# Tabla de pronósticos
forecast_table <- data.frame(
  Periodo = c("2025-TI", "2025-TII", "2025-TIII", "2025-TIV",
              "2026-TI", "2026-TII", "2026-TIII", "2026-TIV"),
  Pronostico = round(as.numeric(forecasts_original$mean), 0),
  Limite_Inferior_80 = round(as.numeric(forecasts_original$lower[,1]), 0),
  Limite_Superior_80 = round(as.numeric(forecasts_original$upper[,1]), 0),
  Limite_Inferior_95 = round(as.numeric(forecasts_original$lower[,2]), 0),
  Limite_Superior_95 = round(as.numeric(forecasts_original$upper[,2]), 0)
)

kable(forecast_table, caption = "Pronósticos del PIB (millones)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Pronósticos del PIB (millones)
Periodo Pronostico Limite_Inferior_80 Limite_Superior_80 Limite_Inferior_95 Limite_Superior_95
2025-TI 243237 234931 251837 230649 256513
2025-TII 248197 237209 259693 231591 265992
2025-TIII 258752 245718 272478 239085 280037
2025-TIV 275029 260006 290920 252389 299701
2026-TI 252662 237778 268477 230258 277246
2026-TII 257348 241340 274418 233273 283908
2026-TIII 267939 250547 286538 241802 296901
2026-TIV 284516 265388 305024 255788 316471

7.3 Gráficos de Pronósticos

par(mfrow = c(2, 1))

# Gráfico en escala logarítmica
plot(forecasts_log, main = "Pronósticos del Log(PIB) - SARIMA(1,1,1)(1,0,1)4",
     xlab = "Año", ylab = "Log(PIB)", col = "blue")
grid()

# Gráfico en escala original
plot(forecasts_original, main = "Pronósticos del PIB - SARIMA(1,1,1)(1,0,1)4",
     xlab = "Año", ylab = "PIB (millones)", col = "blue")
grid()

par(mfrow = c(1, 1))

7.4 Análisis de Crecimiento Proyectado

# Calcular tasas de crecimiento
pib_actual_2024 <- serie$PIB[nrow(serie)]
pronosticos_pib <- as.numeric(forecasts_original$mean)

# Crecimiento trimestral
crecimiento_trimestral <- c(
  (pronosticos_pib[1] / pib_actual_2024 - 1) * 100,
  diff(pronosticos_pib) / pronosticos_pib[-length(pronosticos_pib)] * 100
)

# Crecimiento interanual
crecimiento_interanual <- c(
  (pronosticos_pib[1:4] / serie$PIB[(nrow(serie)-3):nrow(serie)] - 1) * 100,
  (pronosticos_pib[5:8] / pronosticos_pib[1:4] - 1) * 100
)

growth_table <- data.frame(
  Periodo = c("2025-TI", "2025-TII", "2025-TIII", "2025-TIV",
              "2026-TI", "2026-TII", "2026-TIII", "2026-TIV"),
  PIB_Pronosticado = round(pronosticos_pib, 0),
  Crecimiento_Trimestral = round(crecimiento_trimestral, 2),
  Crecimiento_Interanual = round(crecimiento_interanual, 2)
)

kable(growth_table, caption = "Tasas de Crecimiento Proyectadas (%)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Tasas de Crecimiento Proyectadas (%)
Periodo PIB_Pronosticado Crecimiento_Trimestral Crecimiento_Interanual
2025-TI 243237 -7.91 2.67
2025-TII 248197 2.04 1.34
2025-TIII 258752 4.25 3.21
2025-TIV 275029 6.29 4.13
2026-TI 252662 -8.13 3.87
2026-TII 257348 1.85 3.69
2026-TIII 267939 4.12 3.55
2026-TIV 284516 6.19 3.45

8 Conclusiones y Recomendaciones

8.1 Resumen del Modelo

El análisis SARIMA del PIB trimestral ha identificado el modelo SARIMA(1,1,1)(1,0,1)4 como el más adecuado para la serie temporal.

8.1.1 Características del Modelo:

  • AR(1): Captura la dependencia a corto plazo
  • I(1): Una diferenciación regular para eliminar la tendencia
  • MA(1): Corrección de errores de pronóstico
  • SAR(1): Dependencia estacional trimestral muy fuerte (coef ≈ 0.999)
  • SMA(1): Corrección de errores estacionales

8.2 Hallazgos Principales

  1. Estacionalidad Fuerte: El PIB muestra un patrón estacional muy marcado, con el trimestre 4 consistentemente más alto.
  2. Tendencia Creciente: La serie presenta una tendencia de crecimiento sostenido a largo plazo.
  3. Persistencia Estacional: El coeficiente SAR(1) ≈ 0.999 indica que los patrones estacionales son muy persistentes.
  4. Impactos Externos: La serie captura eventos como la crisis de 2020 (COVID-19).

8.3 Pronósticos Clave

  • PIB 2025-TI: 2.43237^{5} millones
  • PIB 2026-TIV: 2.84516^{5} millones
  • Crecimiento promedio proyectado: 3.24% interanual

8.4 Recomendaciones

  1. Uso del Modelo: Apropiado para pronósticos a corto-mediano plazo (1-2 años)
  2. Actualización: Actualizar trimestralmente con nuevos datos
  3. Monitoreo: Vigilar cambios estructurales en la economía
  4. Intervalos de Confianza: Utilizar los rangos de incertidumbre proporcionados
  5. Consideraciones Adicionales: Para análisis más sofisticados, considerar variables exógenas o modelos de cambio de régimen

8.5 Referencias

  • Box, G. E. P., Jenkins, G. M., & Reinsel, G. C. (2015). Time Series Analysis: Forecasting and Control. 5th Edition.
  • Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice. 3rd Edition.