Introducción

La serie temporal que analizaremos es AirPassengers, la cual registra el número de pasajeros internacionales aéreos (en miles) entre 1949 y 1960. Los datos son mensuales y provienen de la librería datasets de R. Esta serie es útil para análisis de series temporales debido a su marcada tendencia al alza y su patrón estacional anual.

El gráfico muestra claramente una tendencia ascendente a lo largo del tiempo, indicando un crecimiento constante en el número de pasajeros aéreos. Se observa también una estacionalidad pronunciada, con picos que se repiten anualmente, generalmente en los meses de verano, y valles en los meses de invierno.

La amplitud de las fluctuaciones estacionales parece aumentar con el tiempo, lo que sugiere una heterocedasticidad, es decir, varianza no constante. La serie no es estacionaria debido a la presencia de tendencia y estacionalidad.

Transformaciones necesarias para lograr la estacionariedad

Para lograr la estacionariedad, aplicaremos primero una transformación logarítmica para estabilizar la varianza y luego una diferenciación para eliminar la tendencia y la estacionalidad.

# Transformación logarítmica
log_ap <- log(ap_ts)

# Gráfico de la serie logarítmica
autoplot(log_ap) +
  ggtitle("Logaritmo del Número de Pasajeros Aéreos") +
  xlab("Año") +
  ylab("Log(Número de Pasajeros)") +
  theme_minimal()

# Aplicar diferenciación estacional (diferencia de orden 1 con frecuencia 12)
diff_log_ap <- diff(log_ap, lag = 12)

# Gráfico de la serie diferenciada estacionalmente
autoplot(diff_log_ap) +
  ggtitle("Serie Diferenciada Estacionalmente (Log-Transformed)") +
  xlab("Año") +
  ylab("Diferencia Estacional") +
  theme_minimal()

# Aplicar diferenciación no estacional (diferencia de orden 1)
diff_log_ap_final <- diff(diff_log_ap, lag = 1)

# Gráfico de la serie doblemente diferenciada
autoplot(diff_log_ap_final) +
  ggtitle("Serie Doblemente Diferenciada (Log-Transformed)") +
  xlab("Año") +
  ylab("Diferencia Final") +
  theme_minimal()

# Prueba de Dickey-Fuller Aumentada para verificar estacionariedad
library(tseries)
adf.test(diff_log_ap_final)
## Warning in adf.test(diff_log_ap_final): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_log_ap_final
## Dickey-Fuller = -5.1993, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Interpretación:

El resultado del test de Dickey-Fuller aumentado (ADF) aplicado a la serie transformada mediante logaritmo y diferenciación () arrojó un estadístico de prueba de \(-5.1993\) con un rezago de orden 5 y un valor \(p < 0.01\). Estos resultados permiten rechazar la hipótesis nula de presencia de raíz unitaria y concluir que la serie diferenciada es estacionaria. La advertencia sobre el valor \(p\) (“”) indica que este es aún menor al umbral de impresión estándar, reforzando la evidencia a favor de la estacionariedad. En consecuencia, puede afirmarse que las transformaciones aplicadas han sido efectivas para estabilizar la media y la varianza de la serie.

Identificación del modelo

Analizamos las funciones de Autocorrelación (ACF) y Autocorrelación Parcial (PACF) de la serie estacionaria (diff_log_ap_final) para identificar los órdenes (p, q) y (P, Q) de los modelos SARIMA.

# Graficar ACF y PACF de la serie estacionaria
par(mfrow=c(1,2))
acf(diff_log_ap_final, lag.max = 48, main="ACF de la Serie Estacionaria")
pacf(diff_log_ap_final, lag.max = 48, main="PACF de la Serie Estacionaria")

par(mfrow=c(1,1)) # Resetear el layout

Interpretación de ACF y PACF:

La inspección conjunta de las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) de la serie estacionaria muestra que, tras la transformación aplicada, los valores de autocorrelación caen rápidamente dentro de las bandas de significancia, indicando que la dependencia serial se ha reducido de manera considerable. En la ACF se observa un decaimiento abrupto luego del rezago inicial, lo cual es consistente con la presencia de un componente de media móvil de bajo orden. Por su parte, la PACF exhibe algunos picos aislados dentro de los primeros rezagos, pero sin una estructura clara y persistente, lo que sugiere la posible inclusión de términos autorregresivos adicionales de bajo orden. En conjunto, este comportamiento respalda la hipótesis de estacionariedad previamente confirmada mediante pruebas estadísticas y aporta evidencia para la especificación de modelos ARIMA o SARIMA, donde un modelo con términos MA y posiblemente algún componente AR de orden reducido podría capturar de manera adecuada la dinámica residual de la serie.

ACF: La ACF muestra decaimiento lento en los lags no estacionales y picos significativos en los lags estacionales (12, 24, 36…). Esto sugiere la necesidad de un componente MA estacional (Q). PACF: La PACF muestra un decaimiento más rápido en los lags no estacionales, pero también picos significativos en los lags estacionales. Esto sugiere la necesidad de un componente AR estacional (P).

Basándonos en los gráficos ACF y PACF de la serie estacionaria, y considerando la estacionalidad anual (periodo 12), proponemos dos modelos SARIMA.

Modelo 1 (SARIMA (1,1,1)(1,1,1)): Este modelo es un punto de partida común. El (1,1,1) no estacional indica un componente MA(1) y una diferenciación de orden 1. El (1,1,1)[12] estacional indica un componente MA(1) estacional y una diferenciación estacional de orden 1.

Modelo 2 (SARIMA (2,1,1)(0,1,1)): Este modelo incluye términos AR y MA tanto no estacionales como estacionales, además de las diferenciaciones. Estimación de parámetros Estimaremos los parámetros de ambos modelos utilizando la función Arima de la librería forecast.

# Ajustar Modelo 1: SARIMA(1,1,1)(1,1,1)[12]
modelo1 <- Arima(ap_ts, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))
summary(modelo1)
## Series: ap_ts 
## ARIMA(1,1,1)(1,1,1)[12] 
## 
## Coefficients:
##           ar1      ma1     sar1    sma1
##       -0.1386  -0.2028  -0.9228  0.8329
## s.e.   0.5865   0.6128   0.2387  0.3519
## 
## sigma^2 = 135:  log likelihood = -506.15
## AIC=1022.3   AICc=1022.78   BIC=1036.68
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.3153117 10.91024 8.003304 0.04379749 2.848886 0.2498666
##                      ACF1
## Training set 0.0004692493
# Ajustar Modelo 2: SARIMA(2,1,1)(0,1,1)[12]
modelo2 <- Arima(ap_ts, order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
summary(modelo2)
## Series: ap_ts 
## ARIMA(2,1,1)(0,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1     sma1
##       0.6031  0.1998  -0.9773  -0.1065
## s.e.  0.0898  0.0898   0.0294   0.0864
## 
## sigma^2 = 131.7:  log likelihood = -504.18
## AIC=1018.36   AICc=1018.84   BIC=1032.74
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 1.361857 10.77916 7.756231 0.4160733 2.744797 0.2421529 -0.0036558

Comparación de modelos ARIMA

Característica ARIMA(1,1,1)(1,1,1)[12] ARIMA(2,1,1)(0,1,1)[12]
Coeficientes estimados ar1 = -0.1386 (s.e. = 0.5865)
ma1 = -0.2028 (s.e. = 0.6128)
sar1 = -0.9228 (s.e. = 0.2387)
sma1 = 0.8329 (s.e. = 0.3519)		 ar1 = 0.6031 (s.e. = 0.0898)
ar2 = 0.1998 (s.e. = 0.0898)
ma1 = -0.9773 (s.e. = 0.0294)
sma1 = -0.1065 (s.e. = 0.0864)
Varianza residual (σ²) 135.0 131.7
Log-verosimilitud -506.15 -504.18
Criterios de información AIC = 1022.30
AICc = 1022.78
BIC = 1036.68		 AIC = 1018.36
AICc = 1018.84
BIC = 1032.74
Errores de entrenamiento ME = 0.3153
RMSE = 10.9102
MAE = 8.0033
MPE = 0.0438
MAPE = 2.8489%
MASE = 0.2499
ACF1 = 0.0005		 ME = 1.3619
RMSE = 10.7792
MAE = 7.7562
MPE = 0.4161
MAPE = 2.7448%
MASE = 0.2422
ACF1 = -0.0037

Interpretación

  1. Ajuste de parámetros

    • En el modelo ARIMA(1,1,1)(1,1,1)[12], los coeficientes no son estadísticamente fuertes, ya que presentan errores estándar relativamente altos, lo que indica cierta inestabilidad en la estimación.
    • En el modelo ARIMA(2,1,1)(0,1,1)[12], los parámetros tienen errores estándar bajos, especialmente el MA(1), lo que sugiere mayor precisión y significancia estadística en la estimación.

  2. Bondad del ajuste

    • Ambos modelos presentan varianzas residuales similares, aunque el modelo ARIMA(2,1,1)(0,1,1)[12] tiene una varianza más baja (131.7 vs. 135.0), lo que refleja un ajuste ligeramente mejor.
    • Los criterios de información (AIC, AICc, BIC) favorecen al ARIMA(2,1,1)(0,1,1)[12], ya que presenta valores consistentemente menores, lo que indica mejor parsimonia y ajuste.

  3. Errores de predicción en el entrenamiento

    • El ARIMA(2,1,1)(0,1,1)[12] presenta menores valores de RMSE (10.77 vs. 10.91) y MAE (7.76 vs. 8.00), además de un MAPE más bajo (2.74% vs. 2.85%), lo cual indica un desempeño predictivo superior en términos de precisión relativa.
    • El MASE de ambos modelos es bajo (<0.25), lo que evidencia un ajuste adecuado, pero nuevamente el segundo modelo muestra un desempeño ligeramente mejor.

  4. Autocorrelación residual (ACF1)

    • En ambos modelos, la autocorrelación del primer rezago residual es cercana a cero, lo que sugiere que los errores son aproximadamente ruido blanco.
    • El ARIMA(2,1,1)(0,1,1)[12] presenta un valor aún más cercano a cero (-0.0037), reforzando la calidad de su ajuste.

El modelo ARIMA(2,1,1)(0,1,1)[12] es el que mejor representa la serie ap_ts, dado que presenta:

  • Parámetros más significativos estadísticamente.
  • Mejores indicadores de ajuste (AIC, BIC, varianza residual).
  • Menores errores de predicción.
# Verificación de residuos del Modelo 1
checkresiduals(modelo1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(1,1,1)[12]
## Q* = 38.83, df = 20, p-value = 0.007
## 
## Model df: 4.   Total lags used: 24
# Verificación de residuos del Modelo 2
checkresiduals(modelo2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(0,1,1)[12]
## Q* = 38.937, df = 20, p-value = 0.006789
## 
## Model df: 4.   Total lags used: 24

Interpretación de la verificación de residuos:

Modelo Estadístico Q* Grados de libertad (df) p-valor Interpretación
ARIMA(1,1,1)(1,1,1)[12] 38.83 20 0.007 El p-valor < 0.05 indica que se rechaza la hipótesis nula de independencia de los residuos. Esto sugiere que aún queda autocorrelación significativa en los residuos, por lo que el modelo no capta completamente la estructura de la serie.
ARIMA(2,1,1)(0,1,1)[12] 38.94 20 0.0068 De manera similar, el p-valor < 0.05 lleva a rechazar la hipótesis nula. Esto implica que también persiste autocorrelación en los residuos, por lo que el ajuste no es del todo adecuado.

Interpretación

  1. Hipótesis del test

    • H0: Los residuos son ruido blanco (no autocorrelación).
    • H1: Los residuos presentan autocorrelación.

  2. Resultados

    • Ambos modelos presentan un p-valor menor a 0.05, lo que significa que los residuos no son ruido blanco.
    • En consecuencia, ninguno de los dos modelos captura completamente la dinámica de la serie.

  3. Implicaciones

    • Aunque el modelo ARIMA(2,1,1)(0,1,1)[12] tiene mejores métricas de ajuste (AIC, BIC, errores de predicción), desde el punto de vista diagnóstico, todavía presenta problemas de autocorrelación en los residuos.

    • Esto sugiere que sería recomendable:

      • Ajustar modelos alternativos (por ejemplo, revisar órdenes SARIMA más complejos o modelos multiplicativos distintos).
      • Considerar transformaciones adicionales (Box-Cox, diferenciación estacional alternativa).
      • Explorar enfoques complementarios (modelos ARIMAX con variables externas, o modelos basados en Machine Learning si la serie es más compleja).

Ambos modelos presentan buen ajuste en términos de AIC/BIC y precisión, pero el test de Ljung-Box evidencia autocorrelación en los residuos, lo que limita su validez como modelos definitivos. El ARIMA(2,1,1)(0,1,1)[12] sigue siendo preferible, pero debería refinarse o compararse con modelos alternativos.

# Resumen de métricas de los modelos
print("Métricas del Modelo 1: SARIMA(1,1,1)(1,1,1)[12]:")
## [1] "Métricas del Modelo 1: SARIMA(1,1,1)(1,1,1)[12]:"
print(modelo1)
## Series: ap_ts 
## ARIMA(1,1,1)(1,1,1)[12] 
## 
## Coefficients:
##           ar1      ma1     sar1    sma1
##       -0.1386  -0.2028  -0.9228  0.8329
## s.e.   0.5865   0.6128   0.2387  0.3519
## 
## sigma^2 = 135:  log likelihood = -506.15
## AIC=1022.3   AICc=1022.78   BIC=1036.68
print("Métricas del Modelo 2: SARIMA(2,1,1)(0,1,1)[12]")
## [1] "Métricas del Modelo 2: SARIMA(2,1,1)(0,1,1)[12]"
print(modelo2)
## Series: ap_ts 
## ARIMA(2,1,1)(0,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1     sma1
##       0.6031  0.1998  -0.9773  -0.1065
## s.e.  0.0898  0.0898   0.0294   0.0864
## 
## sigma^2 = 131.7:  log likelihood = -504.18
## AIC=1018.36   AICc=1018.84   BIC=1032.74
# Comparación explícita de AIC y BIC
print(paste("AIC Modelo 1:", round(AIC(modelo1), 3)))
## [1] "AIC Modelo 1: 1022.3"
print(paste("BIC Modelo 1:", round(BIC(modelo1), 3)))
## [1] "BIC Modelo 1: 1036.676"
print(paste("AIC Modelo 2:", round(AIC(modelo2), 3)))
## [1] "AIC Modelo 2: 1018.36"
print(paste("BIC Modelo 2:", round(BIC(modelo2), 3)))
## [1] "BIC Modelo 2: 1032.736"

Comparación de Modelos SARIMA

Característica Modelo 1: SARIMA(1,1,1)(1,1,1)[12] Modelo 2: SARIMA(2,1,1)(0,1,1)[12]
Coeficientes estimados ar1 = -0.1386 (s.e. 0.5865)
ma1 = -0.2028 (s.e. 0.6128)
sar1 = -0.9228 (s.e. 0.2387)
sma1 = 0.8329 (s.e. 0.3519)		 ar1 = 0.6031 (s.e. 0.0898)
ar2 = 0.1998 (s.e. 0.0898)
ma1 = -0.9773 (s.e. 0.0294)
sma1 = -0.1065 (s.e. 0.0864)
Varianza residual (σ²) 135.0 131.7
Log-verosimilitud -506.15 -504.18
Criterios de información AIC = 1022.30
AICc = 1022.78
BIC = 1036.68		 AIC = 1018.36
AICc = 1018.84
BIC = 1032.74

Interpretación

  1. Precisión de los coeficientes

    • En el Modelo 1, los errores estándar de ar1 y ma1 son elevados, lo que sugiere baja significancia estadística de esos parámetros.
    • En el Modelo 2, todos los coeficientes muestran errores estándar pequeños, lo que indica mayor robustez en las estimaciones.

  2. Ajuste del modelo

    • La varianza residual del Modelo 2 es menor (131.7 vs. 135.0), lo que sugiere un mejor ajuste a los datos.
    • La log-verosimilitud es más alta en el Modelo 2 (-504.18 > -506.15), lo cual también respalda su superioridad.

  3. Criterios de información

    • Tanto el AIC como el BIC son menores en el Modelo 2 (AIC = 1018.36; BIC = 1032.74), lo que indica mayor parsimonia y capacidad de generalización.

Conclusión

El SARIMA(2,1,1)(0,1,1)[12] (Modelo 2) es estadística y computacionalmente superior al SARIMA(1,1,1)(1,1,1)[12] (Modelo 1).

  • Presenta parámetros más significativos, menor varianza residual y mejores valores de AIC y BIC.
  • Sin embargo, al evaluar los tests de Ljung-Box, ambos modelos aún presentan autocorrelación en los residuos, por lo que, aunque el Modelo 2 es preferible, podría ser necesario explorar ajustes adicionales para garantizar un modelo más adecuado.

Serie de entrenamiento y prueba.

# Dividir la serie en entrenamiento y prueba (aprox. 80% train, 20% test)
n_obs <- length(ap_ts)
train_end_year <- 1949 + floor((n_obs * 0.8) / 12) - 1
train_end_month <- ceiling((n_obs * 0.8) %% 12)
if(train_end_month == 0) train_end_month <- 12 # Ajuste para el último mes del año

train_data <- window(ap_ts, end = c(train_end_year, train_end_month))
test_data <- window(ap_ts, start = c(train_end_year, train_end_month + 1))
if(length(test_data) == 0) { # Caso si la división resulta en un set de test vacío
  test_data <- window(ap_ts, start = c(train_end_year + 1, 1))
}

# Ajustar el modelo seleccionado (SARIMA(2,1,1)(0,1,1)) en los datos de entrenamiento
# Reajustamos por si acaso 'modelo1' no era el seleccionado en la lógica anterior
modelo_final <- Arima(train_data, order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))

# Realizar pronósticos para el periodo de prueba
pronostico_out_of_sample <- forecast(modelo_final, h = length(test_data))

# Comparar pronósticos con los datos de prueba
p_oos <- autoplot(test_data) +
  autolayer(pronostico_out_of_sample, series="Pronóstico Out-of-Sample") +
  ggtitle("Pronóstico Out-of-Sample vs Datos Reales") +
  xlab("Año") +
  ylab("Número de Pasajeros (miles)") +
  theme_minimal()
print(p_oos)

# Calcular métricas de evaluación (ECM y MAPE)
eval_out <- accuracy(pronostico_out_of_sample, test_data)
print("--- Métricas de Evaluación Out-of-Sample ---")
## [1] "--- Métricas de Evaluación Out-of-Sample ---"
print(eval_out)
##                       ME     RMSE       MAE         MPE      MAPE      MASE
## Training set   0.7410824  8.84234  6.440203   0.2029285  2.819947 0.2130524
## Test set     -65.9353373 71.54077 65.935337 -16.3724560 16.372456 2.1812481
##                      ACF1 Theil's U
## Training set -0.008271525        NA
## Test set      0.627318041  1.594536
# Pronóstico in-sample (para comparar la bondad de ajuste en el set de entrenamiento)
pronostico_in_sample <- forecast(modelo_final, h = length(train_data))
p_is <- autoplot(train_data) +
  autolayer(pronostico_in_sample, series="Pronóstico In-Sample") +
  ggtitle("Pronóstico In-Sample vs Datos Reales (Entrenamiento)") +
  xlab("Año") +
  ylab("Número de Pasajeros (miles)") +
  theme_minimal()
print(p_is)

Métricas de Evaluación Out-of-Sample

Conjunto ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 0.741 8.842 6.440 0.203% 2.82% 0.213 -0.0083 NA
Test set -65.935 71.541 65.935 -16.372% 16.37% 2.181 0.6273 1.595

Interpretación

1. Rendimiento en entrenamiento (in-sample)

  • Precisión elevada:

    • RMSE = 8.84 y MAE = 6.44 → errores pequeños en relación al nivel de la serie.
    • MAPE = 2.82% → error porcentual muy bajo, excelente capacidad predictiva in-sample.
    • MASE = 0.213 (<1) → el modelo es mucho mejor que un pronóstico ingenuo en entrenamiento.

  • Diagnóstico de residuos:

    • ACF1 = -0.008 → residuos prácticamente no autocorrelacionados, lo cual es positivo.

2. Rendimiento en prueba (out-of-sample)

  • Deterioro importante:

    • RMSE = 71.54 y MAE = 65.94 → errores muy elevados en el set de prueba.
    • MAPE = 16.37% → se multiplica por casi 6 respecto al entrenamiento (2.82%), evidencia de sobreajuste.
    • MASE = 2.18 (>1) → en test el modelo es peor que un pronóstico ingenuo, lo que limita su utilidad práctica.

  • Residuos fuera de muestra:

    • ACF1 = 0.627 → fuerte autocorrelación, lo que indica que el modelo no logra capturar patrones futuros relevantes.

  • Benchmarking:

    • Theil’s U = 1.595 (>1) → el modelo es inferior a un random walk / naive forecast en el test set.

Conclusión

  • El modelo presenta un ajuste sobresaliente dentro de la muestra, pero sufre un grave deterioro fuera de muestra.
  • Los indicadores (MAPE, MASE, Theil’s U) muestran que el modelo no generaliza y es superado por pronósticos ingenuos en validación.
  • Este comportamiento revela sobreajuste (overfitting) o un cambio estructural en la serie en el periodo de prueba.