División de la serie

Aantes de iniciar dividimos la serie de tiempo o DataFrame en conjuntos de entrenamiento (train) y prueba (test), usando un porcentaje definido por \(train\_size\). El conjunto de entrenamiento contiene los primeros datos (orden temporal) y el de prueba el resto, evitando mezclar información futura. Es adecuada para modelos de series de tiempo, donde no debe aplicarse un muestreo aleatorio.

def train_test_split(data, train_size=0.90):
  size = int(len(data)*train_size)
  # split data
  X_train = data.iloc[:size]
  X_test = data.iloc[size:]
  return X_train, X_test

# split
X_train, X_test = train_test_split(xt, train_size=0.85)

print('Train : ', len(X_train))

## Train :  107

print('Test : ', len(X_test))

## Test :  20

print('Total : ', len(xt))

## Total :  127

La imagen abajo muestra un crecimiento sostenido de los homicidios entre 2015 y 2018, seguido de un periodo de alta volatilidad y estancamiento entre 2019 y 2022. Al final de la serie se observa una tendencia a la baja hacia el final del periodo.

plt.figure(figsize=(12,5))
plt.plot(X_train, label='Train')
plt.plot(X_test, label = 'Test')
plt.title('Homicidios Registrados en México')
plt.xlabel('Fecha')
plt.ylabel('Total')
plt.legend()
#mplcyberpunk.add_glow_effects()
plt.show()

El gráfico abajo evidencia una estacionalidad mensual marcada: cada inicio de año los valores parten bajos y aumentan de forma sistemática dentro de cada mes. Las líneas rojas representan el promedio mensual, mostrando que marzo, mayo y diciembre concentran niveles más altos de homicidios, mientras que febrero y junio tienden a ser más bajos. En conjunto, se observa un patrón recurrente estable a lo largo de los años, lo que justifica el uso de modelos con componente estacional.

fig = month_plot(X_train.Total);
fig.set_size_inches((12, 5))
plt.show()

El gráfico de autocorrelación (ACF) muestra valores positivos y estadísticamente significativos en múltiples rezagos, lo que indica una fuerte dependencia temporal en la serie de homicidios. La disminución gradual de la autocorrelación sugiere la presencia de tendencia y persistencia, lo que puede suponer que la serie no es estacionaria.

fig = plot_acf(X_train.Total)
fig.set_size_inches((12, 4))
plt.show()

La descomposición muestra claramente que la serie de homicidios está compuesta por tendencia, estacionalidad y un componente irregular bien definidos. La tendencia crece de forma sostenida hasta alrededor de 2018–2019 y posteriormente presenta una ligera estabilización y descenso, lo que podría sugierir un cambio estructural en el nivel del fenómeno. La estacionalidad es marcada y recurrente a lo largo de los años, con variaciones mensuales sistemáticas, lo que confirma un patrón anual consistente.

dec = seasonal_decompose(xt, model="additive")

series = [dec.observed,dec.trend, dec.seasonal, dec.resid]
titulos= ['Serie observada', 'Componente de tendencia', 'Componente estacional', 'Residuo']
fig, axes = plt.subplots(2,2, sharex=False, figsize=(16, 8))
k=0
for i in range(2):
  for j in range(2):
    serie = series[k]
    axes[i, j].plot(serie)
    axes[i, j].set_title(titulos[k])
    axes[i, j].set_xlabel('Fecha')
    axes[i, j].set_ylabel('Total')
    k +=1
plt.tight_layout()
plt.show()

Modelo 1

El modelo \(SARIMAX(0,1,1)(0,1,1)12\) se ajusta adecuadamente a la serie Total, capturando tanto la dinámica de corto plazo como la estacionalidad anual. Los coeficientes MA(1) y MA estacional (12) son estadísticamente significativos, lo que confirma la presencia de dependencia temporal y estacional. El diagnóstico de residuos (Ljung-Box, Jarque-Bera y heterocedasticidad) indican que estos se comportan como ruido blanco, sin autocorrelación ni problemas graves de normalidad o varianza. En conjunto, el modelo es estadísticamente consistente y adecuado para pronóstico.

model = SARIMAX(X_train, order=(0,1,1), seasonal_order=(0,1, 1, 12)).fit()
print(model.summary())

##                                      SARIMAX Results                                      
## ==========================================================================================
## Dep. Variable:                              Total   No. Observations:                  107
## Model:             SARIMAX(0, 1, 1)x(0, 1, 1, 12)   Log Likelihood                -606.648
## Date:                          lun., 12 ene. 2026   AIC                           1219.297
## Time:                                    22:20:02   BIC                           1226.927
## Sample:                                01-01-2015   HQIC                          1222.379
##                                      - 11-01-2023                                         
## Covariance Type:                              opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ma.L1         -0.5944      0.092     -6.489      0.000      -0.774      -0.415
## ma.S.L12      -0.7678      0.129     -5.958      0.000      -1.020      -0.515
## sigma2      2.096e+04   2779.778      7.540      0.000    1.55e+04    2.64e+04
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.06   Jarque-Bera (JB):                 2.99
## Prob(Q):                              0.80   Prob(JB):                         0.22
## Heteroskedasticity (H):               1.06   Skew:                            -0.41
## Prob(H) (two-sided):                  0.86   Kurtosis:                         3.27
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

Modelo 2

Probamos un segundo modelo \(SARIMAX(1,1,1)(1,1,1)12\), al igual que el modelo anterior los residuos no presentan autocorrelaciones ni problemas de normalidad o varianza gravez. Por otro lado, si observamos algunas métricas como el AIC, y BIC, podemos notar que el modelo anterior es más simple, es decir, posee menos parámetros, lo cual nos inclina hacia la selección del primer modelo

model2 = SARIMAX(X_train, order=(1,1,1), seasonal_order=(1,1, 1, 12)).fit()
print(model2.summary())

##                                      SARIMAX Results                                      
## ==========================================================================================
## Dep. Variable:                              Total   No. Observations:                  107
## Model:             SARIMAX(1, 1, 1)x(1, 1, 1, 12)   Log Likelihood                -606.616
## Date:                          lun., 12 ene. 2026   AIC                           1223.232
## Time:                                    22:20:03   BIC                           1235.949
## Sample:                                01-01-2015   HQIC                          1228.369
##                                      - 11-01-2023                                         
## Covariance Type:                              opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ar.L1          0.0365      0.203      0.180      0.857      -0.361       0.434
## ma.L1         -0.6131      0.172     -3.555      0.000      -0.951      -0.275
## ar.S.L12       0.0036      0.170      0.021      0.983      -0.329       0.336
## ma.S.L12      -0.7665      0.202     -3.797      0.000      -1.162      -0.371
## sigma2      2.098e+04   3081.211      6.807      0.000    1.49e+04     2.7e+04
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.16   Jarque-Bera (JB):                 2.84
## Prob(Q):                              0.69   Prob(JB):                         0.24
## Heteroskedasticity (H):               1.07   Skew:                            -0.41
## Prob(H) (two-sided):                  0.85   Kurtosis:                         3.22
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

Prueba de autocorrelación Ljung-Box Test

A veces es importante probar si los valores de una serie temporal, digamos \(u_t\), son independientes. Podemos decir que existe un patrón de autocorrelación si las series están seriamente correlacionadas. Una prueba muy utilizada para verificar si una serie de valores está correlacionada es el estadístico de Ljung Box, que se basa en la autocorrelación como conjunto. Supongamos que nos dan \(m\) valores autocorrelacionados. Probamos la siguiente hipótesis nula conjunta sobre las correlaciones entre los valores

\[ \begin{equation} H_0 : \rho_1 = ... = \rho_k = 0 \end{equation} \] el estadístico \[ \begin{equation} Q = n(n+2) \sum_{k=1}^{m} \frac{\hat{\rho}_k^2}{(n-k)} \end{equation} \] donde \(m\) es el número de rezagos (lags), and \(\hat{\rho}_k^2\) representa \(\hat{Cov} (u_t, u_{t+k})\) for \(k=1,...,m\).

Referencia: Pankratz(\(1983\)).

Prueba de normalidad Jarque-Bera

La prueba de normalidad Jarque-Bera es una prueba estadística que se emplea para ver si una determinada serie de tiempo se distribuye de forma normal.

El estadístico de prueba Jarque-Bera se define como \[ \begin{equation*} JB =[(n-k+1) / 6] * [S^2 + (0,25*(C-3)2)] \end{equation*} \]

Bajo la hipótesis nula de normalidad, la prueba Jarque-Bera(JB) \(\sim \chi^2(2)\)

donde \(n\) denota el número de observaciones en la muestra, \(k\) denota el número de regresores (\(k=1\) si no se utiliza en una regresión), \(S\) denota la asimetría de la muestra y \(C\) denota la curtosis de la muestra.

La hipótesis nula es \[ \begin{equation} H_0 : \text{La serie sigue una distribución normal} \end{equation} \]

Prueba de heterocedasticidad Goldfeld-Quandt

La heterocedasticidad en la situación en la que la variabilidad de una serie es desigual para un rango de valores distinto.

La prueba Goldfeld-Quandt se utiliza para probar la presencia de heteroscedasticidad en los datos dados.

La hipótesis nula es \[ \begin{equation} H_0 : \text{No hay presencia de heterocedasticidad} \end{equation} \]

Observación:

Una prueba de hipótesis se rechaza mediante: - El estadístico de prueba. - El p-valor. Se rechaza si p-valor < \(\alpha\)

En el gr+afico abajo, se observa que los residuos estandarizados oscilan alrededor de cero sin patrón visible, lo que indica una buena especificación del modelo. El histograma muestra una distribución cercana a la normal, con ligeras desviaciones en las colas.El Q-Q plot confirma una normalidad razonable, con pequeños desajustes en los extremos. Por otra parte, el correlograma no presenta autocorrelaciones significativas, sugiriendo residuos tipo ruido blanco y un ajuste adecuado del modelo.

from statsmodels.tsa.statespace.sarimax import SARIMAXResults

model.plot_diagnostics(figsize=(14,6))
plt.tight_layout()
plt.show()

Confirmación

Por mera curiosidad vamos a corroborar el modelo mediante auto_arima

import pmdarima as pm

model = pm.auto_arima(X_train, start_p=1, start_q=1,
                      test='adf',       #adftest para encontrar el valor óptimo de 'd'
                      max_p=5, max_q=5, # máx p y q
                      m=12,              # frecuencia de la serie
                      d=1,              # diferenciación de la parte regular (no estacional)
                      seasonal=True,    #
                      start_P=0, start_Q=0,  # valores de inicio de P y Q de la parte estacional
                      D=1,                   # Diferenciación de la parte estacional
                      suppress_warnings=True,
                      stepwise=True)

print(model.summary())

##                                      SARIMAX Results                                      
## ==========================================================================================
## Dep. Variable:                                  y   No. Observations:                  107
## Model:             SARIMAX(0, 1, 1)x(0, 1, 1, 12)   Log Likelihood                -606.648
## Date:                          lun., 12 ene. 2026   AIC                           1219.297
## Time:                                    22:20:14   BIC                           1226.927
## Sample:                                01-01-2015   HQIC                          1222.379
##                                      - 11-01-2023                                         
## Covariance Type:                              opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ma.L1         -0.5944      0.092     -6.489      0.000      -0.774      -0.415
## ma.S.L12      -0.7678      0.129     -5.958      0.000      -1.020      -0.515
## sigma2      2.096e+04   2779.778      7.540      0.000    1.55e+04    2.64e+04
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.06   Jarque-Bera (JB):                 2.99
## Prob(Q):                              0.80   Prob(JB):                         0.22
## Heteroskedasticity (H):               1.06   Skew:                            -0.41
## Prob(H) (two-sided):                  0.86   Kurtosis:                         3.27
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).