Modelos ARMA y ARIMA: autocorrelación simple (ACF) y parcial (PACF).

Memoria de una serie

La lección 1.6.2 midió la autocorrelación de las rentabilidades y concluyó que la dirección casi no tiene memoria. La lección 2.1.2 separó las series estacionarias de las que no lo son. Esta lección da el paso siguiente: describir con un modelo cuánta memoria tiene una serie estacionaria y usarlo para prever.

Los modelos ARMA combinan dos ingredientes:

• AR(p), autorregresivo: el valor de hoy depende de los $p$ valores anteriores, $y_t = c + \phi_1 y_{t-1} + \dots + \phi_p y_{t-p} + \varepsilon_t$.
• MA(q), media móvil (de los errores, no una media móvil de precios): el valor de hoy depende de los $q$ choques anteriores, $y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q}$.

Un ARMA(p, q) tiene los dos. Un ARIMA(p, d, q) es un ARMA sobre la serie diferenciada $d$ veces: un ARIMA(1,1,1) del logaritmo del precio es un ARMA(1,1) de sus rentabilidades. El procedimiento clásico de Box y Jenkins tiene cuatro pasos que se repiten hasta que el diagnóstico es limpio:

flowchart TB
  I[Identificar] --> E[Estimar]
  E --> D[Diagnosticar]
  D -- residuos limpios --> P[Prever]
  D -- queda memoria --> I

Se usan dos series sintéticas del nivel 2 con respuesta conocida (README de los datos):

• el spread verdadero del par A, $u = \ln y + 2{,}1 - 1{,}5\,\ln x$ (en el código, multiplicado por 100 para leerlo en %), un Ornstein-Uhlenbeck discreto: un AR(1) con $\phi = 0{,}5^{1/5} = 0{,}871$ (vida media de 5 sesiones);
• las rentabilidades diarias del futuro IDA, un índice ficticio de renta variable sin autocorrelación propia salvo la que inducen la deriva por tramos del generador y un momentum de series temporales pequeño.

import warnings
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import acf, pacf

cont = pd.read_parquet("futuros2_continuos.parquet")
ida = np.log(cont[cont["mercado"] == "IDA"]
             .set_index("fecha")["cierre_ajustado"])
r = 100 * ida.diff().dropna()               # rentabilidad en %

pares = pd.read_parquet("pares2.parquet").pivot(
    index="fecha", columns="activo", values="precio")
# spread verdadero en log (ln y = -2,1 + 1,5·ln x + u), en %
sA = 100 * (np.log(pares["PZAY"]) + 2.1
            - 1.5 * np.log(pares["PZAX"]))

print(f"IDA: {len(r)} sesiones; spread A: {len(sA)} sesiones")
print(f"φ verdadero del spread A: {0.5 ** (1 / 5):.4f}")

IDA: 6538 sesiones; spread A: 3018 sesiones
φ verdadero del spread A: 0.8706

ACF y PACF: la huella de cada modelo

La ACF (función de autocorrelación) es la correlación de la serie con ella misma $k$ periodos antes. La PACF (autocorrelación parcial) es la correlación en el retardo $k$ descontando los retardos intermedios: el último coeficiente de una regresión de $y_t$ sobre $y_{t-1}, \dots, y_{t-k}$. Cada modelo deja una huella distinta en las dos:

La banda de ruido es $\pm 1{,}96/\sqrt{n}$: si la serie no tuviera memoria, un 95 % de las autocorrelaciones caería dentro. Leída como un árbol de decisión, la tabla queda así:

flowchart TB
  A{¿ACF dentro banda?} -- sí --> W[Sin memoria]
  A -- no --> B{¿PACF se corta?}
  B -- sí --> AR[AR de orden p]
  B -- no --> C{¿ACF se corta?}
  C -- sí --> MA[MA de orden q]
  C -- no --> AM[ARMA p, q]

En la práctica las huellas se mezclan con el ruido de muestreo y la rejilla de criterios de información (más abajo) decide entre los candidatos.

def huella(x, k=5):
    x = pd.Series(x).dropna()
    return pd.DataFrame({"acf": acf(x, nlags=k)[1:],
                         "pacf": pacf(x, nlags=k)[1:]},
                        index=range(1, k + 1))

banda = {"spread A": 1.96 / np.sqrt(len(sA)),
         "IDA": 1.96 / np.sqrt(len(r))}
print({k: round(v, 3) for k, v in banda.items()})
pd.concat({"spread A": huella(sA), "rent. IDA": huella(r)},
          axis=1).round(3)

{'spread A': np.float64(0.036), 'IDA': np.float64(0.024)}
  spread A        rent. IDA
       acf   pacf       acf   pacf
1    0.863  0.864     0.020  0.020
2    0.742 -0.013    -0.015 -0.016
3    0.640  0.007     0.017  0.018
4    0.549 -0.007     0.060  0.059
5    0.469 -0.011     0.025  0.024

La figura de la cabecera dibuja los 12 primeros retardos. El spread A tiene la huella de libro de un AR(1): la ACF baja como una potencia de 0,86 (0,74 es casi 0,86²; 0,64, casi 0,86³) y la PACF vale 0,86 en el retardo 1 y nada después. Las rentabilidades de IDA tienen autocorrelaciones minúsculas, de −0,015 a 0,06; la del retardo 4 (0,060) y la del 5 (0,025) se salen de la banda (±0,024), pero con 6.538 observaciones la banda es estrecha y cualquier dependencia diminuta la supera.

Un MA se ve en un sitio familiar: las rentabilidades solapadas de la lección 2.1.1. La suma de 5 rentabilidades diarias calculada cada día es un MA(4): dos sumas separadas $k$ días comparten $5-k$ días, así que su ACF teórica es $(5-k)/5$ hasta el retardo 4 y cero después.

r5 = r.rolling(5).sum().dropna()
ma = huella(r5, k=7)
ma["acf_teorica"] = [max(5 - k, 0) / 5 for k in ma.index]
ma.round(3)

     acf   pacf  acf_teorica
1  0.814  0.815          0.8
2  0.628 -0.105          0.6
3  0.454 -0.079          0.4
4  0.278 -0.125          0.2
5  0.084 -0.193          0.0
6  0.072  0.402          0.0
7  0.078 -0.021          0.0

La ACF sigue la teoría (0,81 frente a 0,80; 0,63 frente a 0,60…) y cae casi a cero tras el retardo 4; lo que queda en 5, 6 y 7 (de 0,07 a 0,08) es la pequeña autocorrelación propia de IDA, sobre todo la del retardo 4, que entra en varias sumas a la vez. La PACF, en cambio, no se corta: cambia de signo y vuelve a subir en el retardo 6, la huella de un MA. Este es el motivo técnico de que un objetivo solapado infle los t: es una serie MA por construcción.

Estimar: el AR(1) del spread

statsmodels.tsa.arima.model.ARIMA estima por máxima verosimilitud (con un modelo de espacio de estados). El orden se da como (p, d, q) y trend="c" añade la constante. Se pasan arrays de NumPy: con un índice de fechas sin frecuencia regular (las sesiones de la NYSE no la tienen), statsmodels avisa de que no puede prever fechas futuras.

from statsmodels.tsa.arima.model import ARIMA

res_a = ARIMA(sA.to_numpy(), order=(1, 0, 0), trend="c").fit()
tabla = pd.DataFrame({"estimado": res_a.params, "error": res_a.bse},
                     index=["media", "phi", "sigma2"])
tabla["verdad"] = [0, 0.5 ** (1 / 5), np.nan]
tabla.round(4)

        estimado   error  verdad
media     0.2329  0.1297  0.0000
phi       0.8655  0.0090  0.8706
sigma2    0.9475  0.0245     NaN

El $\phi$ estimado, 0,865 con un error típico de 0,009, está a menos de un error típico del verdadero (0,871). La media del spread (0,23 % con error 0,13) no se distingue de cero al 5 %, como debe ser con los $\alpha$ y $\beta$ verdaderos. $\sigma^2$ (0,95, en %²) es la varianza del choque diario; la del spread es mayor, $\sigma^2/(1-\phi^2)$, porque los choques se acumulan.

Elegir el orden: AIC y BIC

Más parámetros siempre ajustan mejor dentro de la muestra. Los criterios de información penalizan cada parámetro:

con $L$ la verosimilitud, $k$ el número de parámetros y $n$ las observaciones. Se elige el menor. Con $n$ grande, la penalización del BIC ($\ln n \approx 8$ con 3.000 datos) es cuatro veces la del AIC, y por eso elige modelos más sencillos.

from statsmodels.tools.sm_exceptions import (ConvergenceWarning,
                                             EstimationWarning)

def rejilla(x, pmax=2, qmax=2):
    filas = []
    for p in range(pmax + 1):
        for q in range(qmax + 1):
            with warnings.catch_warnings():
                # algunos MA no convergen del todo o arrancan de
                # valores no válidos; no cambia el orden elegido
                warnings.simplefilter("ignore", ConvergenceWarning)
                warnings.simplefilter("ignore", EstimationWarning)
                f = ARIMA(x, order=(p, 0, q), trend="c").fit()
            filas.append({"p": p, "q": q, "aic": f.aic, "bic": f.bic})
    t = pd.DataFrame(filas)
    return t.assign(d_aic=t.aic - t.aic.min(),
                    d_bic=t.bic - t.bic.min())

ra = rejilla(sA.to_numpy())
rr = rejilla(r.to_numpy())
cols = ["p", "q", "d_aic", "d_bic"]
print("mínimo del AIC en el spread (p, q):",
      ra.loc[ra.aic.idxmin(), ["p", "q"]].tolist())
print("spread A, por BIC:")
print(ra.sort_values("bic")[cols].head(3).round(1).to_string(
    index=False))
print("rentabilidades IDA, por BIC (y el ARMA(1,1)):")
orden = rr.sort_values("bic")
orden = pd.concat([orden.head(4),
                   orden[(orden.p == 1) & (orden.q == 1)]])
print(orden[cols].round(1).to_string(index=False))

mínimo del AIC en el spread (p, q): [1, 0]
spread A, por BIC:
 p  q  d_aic  d_bic
 1  0    0.0    0.0
 1  1    1.5    7.5
 2  0    1.5    7.5
rentabilidades IDA, por BIC (y el ARMA(1,1)):
 p  q  d_aic  d_bic
 0  0   26.7    0.0
 2  2    0.0    0.4
 1  2   12.2    5.8
 0  1   26.1    6.2
 1  1   30.0   16.9

En el spread, los dos criterios eligen el AR(1), la respuesta correcta, pero con márgenes distintos: el BIC le da 7,5 puntos de ventaja sobre el siguiente y el AIC solo 1,5, porque penaliza menos el parámetro de más. En las rentabilidades de IDA los dos criterios discrepan: el BIC elige el modelo sin memoria (0,0), y el AIC, un ARMA(2,2) con 26,7 puntos de ventaja sobre él. Con 6.538 observaciones, cuatro parámetros más cuestan 8 puntos de AIC y 35 de BIC; unas dependencias minúsculas pagan lo primero y no lo segundo. El ARMA(1,1), el punto de partida habitual, queda por detrás del (0,0) en los dos (30,0 y 16,9 puntos). Los coeficientes explican por qué:

res_r = ARIMA(r.to_numpy(), order=(1, 0, 1), trend="c").fit()
pd.Series(res_r.params, index=["media", "phi", "theta",
                               "sigma2"]).round(4)

media     0.0251
phi      -0.9201
theta     0.9240
sigma2    1.6084
dtype: float64

$\phi = -0{,}920$ y $\theta = 0{,}924$: casi se anulan. En $(1 - \phi B)\,y_t = (1 + \theta B)\,\varepsilon_t$, con $\phi \approx -\theta$ los dos factores se cancelan y queda $y_t \approx \varepsilon_t$: ruido blanco con parámetros de más. Es la forma típica en que un ARMA ajustado a rentabilidades dice «aquí no hay nada». Las dependencias minúsculas que sí hay (la deriva por tramos del generador, con años alcistas y crisis de fechas fijas, y el momentum plantado, que es lento) no tienen la forma de un ARMA de orden bajo.

Diagnosticar: Ljung-Box

Si el modelo es bueno, sus residuos no deben tener memoria. El contraste de Ljung-Box resume las $m$ primeras autocorrelaciones $\hat\rho_k$ de una serie en un estadístico,

que, sin autocorrelación, se reparte como una $\chi^2$ con $m$ grados de libertad (menos los parámetros ARMA estimados, si se aplica a residuos).

from statsmodels.stats.diagnostic import acorr_ljungbox

def lb(x, gl=0):
    t = acorr_ljungbox(x, lags=[10], model_df=gl)
    return t.iloc[0].round(4).to_dict()

pd.DataFrame({
    "spread A": lb(sA),
    "residuos AR(1)": lb(res_a.resid, gl=1),
    "rent. IDA": lb(r),
    "residuos ARMA(1,1)": lb(res_r.resid, gl=2)}).T

                      lb_stat  lb_pvalue
spread A            8347.1553     0.0000
residuos AR(1)         6.8739     0.6502
rent. IDA             50.5221     0.0000
residuos ARMA(1,1)    49.6918     0.0000

El spread A tiene una memoria enorme (Q = 8.347); los residuos del AR(1) ya no (p = 0,65): el modelo recoge toda la estructura. Las rentabilidades de IDA rechazan la ausencia de autocorrelación (Q = 50,5, p < 0,0001), y los residuos del ARMA(1,1) casi igual (Q = 49,7): el modelo no ha quitado nada, como anunciaban sus coeficientes. La dependencia que detecta Ljung-Box está sobre todo en el retardo 4 (0,060), fuera del alcance de un ARMA(1,1). Se podría subir el orden hasta limpiar los residuos, pero la pregunta que importa es otra: si algo de esto sirve para prever.

Prever fuera de la muestra

Se estima cada modelo con una parte de la historia y se prevé el resto sin reestimar: rentabilidades de IDA con 2000-2019 y previsión de 2020-2025; spread A con 2014-2021 y previsión de 2022-2025. res.apply(nuevos_datos) aplica los parámetros ya estimados a datos nuevos y devuelve en fittedvalues la previsión a un paso de cada fecha, hecha con la información del día anterior. La referencia es la previsión ingenua: cero para la rentabilidad (o la media de entrenamiento) y el último valor para el spread (paseo aleatorio).

import statsmodels.api as sm

tr, te = r[:"2019"].to_numpy(), r["2020":].to_numpy()
m = ARIMA(tr, order=(1, 0, 1), trend="c").fit()
f1 = m.apply(te).fittedvalues            # E[r_t | hasta t-1]

def ecm(real, prev):
    return np.mean((real - prev) ** 2)

print(f"ECM ARMA / ECM cero: {ecm(te, f1) / ecm(te, 0):.4f}")
print(f"ECM ARMA / ECM media: {ecm(te, f1) / ecm(te, tr.mean()):.4f}")
print(f"aciertos de signo: {np.mean(np.sign(f1) == np.sign(te)):.1%}"
      f" (días al alza: {np.mean(te > 0):.1%})")

# Diebold-Mariano: ¿la diferencia de pérdidas es distinta de cero?
dif = (te - f1) ** 2 - te ** 2
dm = sm.OLS(dif, np.ones(len(dif))).fit(
    cov_type="HAC", cov_kwds={"maxlags": 5})
print(f"Diebold-Mariano: t = {dm.tvalues[0]:.2f}, "
      f"p = {dm.pvalues[0]:.3f}")

ECM ARMA / ECM cero: 1.0027
ECM ARMA / ECM media: 1.0030
aciertos de signo: 53.4% (días al alza: 53.2%)
Diebold-Mariano: t = 1.11, p = 0.265

El ARMA(1,1) aumenta el error cuadrático medio un 0,3 % frente a prever cero. Acierta el signo el 53,4 % de los días, pero los días al alza son el 53,2 %: una regla «siempre comprado» acierta lo mismo. El contraste de Diebold y Mariano (1995) toma la diferencia de pérdidas de cada día y contrasta si su media es cero, con errores de Newey-West (lección 2.1.1): t = 1,11, p = 0,27. La diferencia no se distingue del azar, y va en contra del modelo. Este es el resultado habitual de un ARMA sobre rentabilidades: alguna autocorrelación significativa dentro de la muestra, nada aprovechable fuera.

En el spread pasa lo contrario. La previsión a $h$ pasos de un AR(1) tiene forma cerrada, $\hat u_{t+h} = \mu + \phi^h (u_t - \mu)$, y se compara con el paseo aleatorio ($\hat u_{t+h} = u_t$) en varios horizontes; para las rentabilidades, se compara la previsión de la suma de las $h$ siguientes con cero:

ta, ea = sA[:"2021"].to_numpy(), sA["2022":].to_numpy()
mu, phi, _ = ARIMA(ta, order=(1, 0, 0), trend="c").fit().params
todo, n0 = np.r_[ta, ea], len(ta)
c, phi_r, _, _ = m.params

filas = {}
for h in [1, 2, 5, 10, 20]:
    t0 = np.arange(n0, len(todo) - h + 1)  # se prevé desde t0 - 1
    ultimo, real = todo[t0 - 1], todo[t0 + h - 1]
    ar = mu + phi ** h * (ultimo - mu)
    # rentabilidades: suma prevista de las h siguientes frente a 0
    idx = np.arange(len(te) - h + 1)
    prev = np.array([sum(c + phi_r ** j * (f1[i] - c)
                         for j in range(h)) for i in idx])
    suma = np.array([te[i:i + h].sum() for i in idx])
    filas[h] = {"spread: AR/paseo": ecm(real, ar) / ecm(real, ultimo),
                "rent.: ARMA/cero": ecm(suma, prev) / ecm(suma, 0)}
pd.DataFrame(filas).T.rename_axis("h").round(3)

    spread: AR/paseo  rent.: ARMA/cero
h
1              0.924             1.003
2              0.857             1.002
5              0.710             1.003
10             0.587             1.001
20             0.535             0.999

A un día, el AR(1) del spread solo reduce el error un 8 % frente al paseo aleatorio: con $\phi = 0{,}87$, el último valor ya es casi la mejor previsión. A 10 sesiones, la reducción es del 41 %, porque el spread ha recorrido buena parte del camino de vuelta ($\phi^{10} \approx 0{,}25$) y el paseo aleatorio no lo sabe. Con las rentabilidades, el cociente se queda entre 0,999 y 1,003 en todos los horizontes: ninguna ganancia. La diferencia entre las dos líneas es la diferencia entre una estructura que se puede operar (lección 2.3.3) y una que no.

ARIMA sobre precios

Con $d = 1$, ARIMA diferencia por dentro. Un ARIMA(1,1,1) del logaritmo del precio es el ARMA(1,1) de sus rentabilidades, sin constante salvo que se pida una tendencia:

lp = 100 * ida.to_numpy()
a111 = ARIMA(lp, order=(1, 1, 1), trend="t").fit()
pd.DataFrame({"ARIMA(1,1,1) log precio": a111.params[:3],
              "ARMA(1,1) rentabilidad": res_r.params[:3]},
             index=["deriva", "phi", "theta"]).round(4)

        ARIMA(1,1,1) log precio  ARMA(1,1) rentabilidad
deriva                   0.0251                  0.0251
phi                     -0.9201                 -0.9201
theta                    0.9240                  0.9240

Los coeficientes coinciden casi exactamente (con trend="t" la deriva del precio es la media de la rentabilidad). La lección 2.1.2 explica cuándo diferenciar; ARIMA solo lo hace por ti. Para el trading, la parte útil del ARIMA está casi siempre en el ARMA de algo estacionario: un spread, una base, un diferencial de tipos.

Ejercicios

1. Ajusta un AR(1) al spread verdadero del par B ($\ln y - 0{,}7 - 0{,}8\,\ln x$, en %) y calcula su vida media, $\ln 0{,}5 / \ln \phi$. Compárala con la verdadera, 15 sesiones.

sB = 100 * (np.log(pares["PZBY"]) - 0.7
            - 0.8 * np.log(pares["PZBX"]))
fb = ARIMA(sB.to_numpy(), order=(1, 0, 0), trend="c").fit()
phi_b = fb.params[1]
print(f"phi = {phi_b:.4f} (verdad {0.5 ** (1 / 15):.4f}); "
      f"vida media = {np.log(0.5) / np.log(phi_b):.1f} sesiones")

phi = 0.9479 (verdad 0.9548); vida media = 13.0 sesiones

2. Simula con semilla 2131 un AR(2) de 3.000 observaciones con $\phi_1 = 0{,}5$ y $\phi_2 = 0{,}3$. Comprueba que la PACF se corta tras el retardo 2 y que el BIC elige el orden correcto en la rejilla.

rng = np.random.default_rng(2131)
e = rng.standard_normal(3000)
x = np.zeros(3000)
for t in range(2, 3000):
    x[t] = 0.5 * x[t - 1] + 0.3 * x[t - 2] + e[t]
print("PACF:", np.round(pacf(x, nlags=4)[1:], 3))
g = rejilla(x)
print(g.sort_values("bic")[["p", "q", "d_bic"]].head(3)
      .round(1).to_string(index=False))

PACF: [0.677 0.306 0.02  0.001]
 p  q  d_bic
 2  0    0.0
 2  1    6.9
 1  2    7.6

3. Un compañero ajusta un ARMA(3,3) a las rentabilidades diarias de una acción, obtiene el menor AIC de su rejilla y un acierto de signo del 54 % dentro de la muestra. ¿Qué tres comprobaciones le pedirías antes de creérselo?