Estacionariedad: ADF, KPSS, diferenciación y diferenciación fraccional.

Qué significa estacionaria

Una serie es estacionaria (en sentido débil) si su media, su varianza y la covarianza entre dos observaciones separadas $k$ periodos no dependen de la fecha. Dicho de otra forma: un tramo de 2005 y uno de 2023 son dos muestras de la misma distribución. Es lo que da sentido a estimar una media, una volatilidad o una regresión con datos del pasado y usarla en el futuro.

Los precios no son estacionarios. Un índice que cotizaba a 1.000 en 2003 y a 5.000 en 2024 no tiene «una media» a la que volver; su varianza crece con el horizonte. El modelo más sencillo que lo recoge es el paseo aleatorio: el precio de hoy es el de ayer más un choque, $p_t = p_{t-1} + \varepsilon_t$. Su diferencia, la rentabilidad, sí es estacionaria. En la jerga, el precio es integrado de orden 1, $I(1)$, y la rentabilidad, $I(0)$.

flowchart LR
  P[Precio I1] --> L[Logaritmo]
  L --> D1[Diferencia d=1]
  L --> DF[Diferencia fracc.]
  D1 --> R[Rentabilidad I0]
  DF --> M[I0 con memoria]

Por qué importa en un sistema de trading:

• Regresiones espurias. Dos series integradas independientes dan correlaciones y t enormes por azar (lección 0.7.2). Una regresión de un precio sobre otro solo tiene sentido si están cointegrados (lección 2.1.5).
• Reversión a la media. Una regla que compra cuando el precio está «por debajo de su media» supone que hay una media a la que volver. En un precio $I(1)$ no la hay; en un spread estacionario, sí (lección 2.3.3).
• Variables de un modelo. Un modelo entrenado con el nivel de un precio aprende la relación de un periodo que no se repetirá (nivel 3).

Los datos de esta lección son sintéticos (README del nivel 2): un futuro de índice que es un paseo aleatorio con deriva, futuros de materias primas que vuelven muy despacio a un ancla (vidas medias de años) y pares cuyo spread verdadero es estacionario o no lo es, según el par. Se sabe la respuesta de antemano; eso permite ver dónde aciertan los contrastes y dónde no.

Dickey-Fuller: contrastar la raíz unitaria

La idea del contraste de Dickey-Fuller es estimar cuánto «vuelve» la serie. En la regresión

un $\gamma$ negativo significa que, cuando $y$ está alto, tiende a bajar: reversión. Si $\gamma = 0$, la serie es un paseo aleatorio (tiene raíz unitaria). El contraste ADF (augmented, por los retardos $\Delta y_{t-i}$, que absorben la autocorrelación de corto plazo) usa la nula $H_0: \gamma = 0$ frente a $\gamma < 0$.

La trampa es que, bajo la nula, $y_{t-1}$ no es estacionario y el estadístico t de $\gamma$ no sigue una t de Student. Sigue la distribución de Dickey-Fuller, más desplazada a la izquierda. Se ve simulando paseos aleatorios con semilla y calculando el estadístico en cada uno:

import numpy as np
import pandas as pd
from scipy import stats

rng = np.random.default_rng(212)
sims, n = 4000, 1000
y = rng.standard_normal((sims, n)).cumsum(axis=1)  # paseos
dy, ylag = np.diff(y, axis=1), y[:, :-1]

def t_gamma(dy, ylag):
    # MCO de dy = a + g·ylag, fila a fila
    xm = ylag - ylag.mean(1, keepdims=True)
    g = (xm * dy).sum(1) / (xm**2).sum(1)
    a = dy.mean(1) - g * ylag.mean(1)
    e = dy - a[:, None] - g[:, None] * ylag
    s2 = (e**2).sum(1) / (dy.shape[1] - 2)
    return g / np.sqrt(s2 / (xm**2).sum(1))

t = t_gamma(dy, ylag)
print(f"percentil 5 % de la t de Dickey-Fuller: "
      f"{np.percentile(t, 5):.2f}  (normal: "
      f"{stats.norm.ppf(0.05):.2f})")
print(f"paseos «estacionarios» con el umbral de la normal: "
      f"{np.mean(t < stats.norm.ppf(0.05)):.1%}")

percentil 5 % de la t de Dickey-Fuller: -2.88  (normal: -1.64)
paseos «estacionarios» con el umbral de la normal: 45.3%

El 5 % inferior del estadístico simulado está en −2,88 (el valor crítico tabulado para la regresión con constante y esta longitud es −2,86; la diferencia es el ruido de 4.000 simulaciones), y no en el −1,64 de la normal. Quien use el umbral de la t de Student declarará estacionario un 45,3 % de los paseos aleatorios: casi la mitad. statsmodels.tsa.stattools.adfuller ya usa los valores críticos correctos (de MacKinnon) y elige el número de retardos $p$ por criterio de información (autolag="AIC").

Dos nulas opuestas: ADF y KPSS

El ADF tiene un problema de diseño: la nula es la raíz unitaria, así que «no rechazar» puede significar que la serie es un paseo aleatorio o que la muestra es demasiado corta para verlo. El contraste KPSS (Kwiatkowski, Phillips, Schmidt y Shin, 1992) invierte los papeles: su nula es que la serie es estacionaria alrededor de una constante, y la rechaza cuando la suma acumulada de las desviaciones crece demasiado. Usar los dos juntos da más información que cualquiera por separado.

import warnings
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tools.sm_exceptions import InterpolationWarning

# kpss avisa cuando el estadístico cae fuera de su tabla de
# p-valores (los recorta a 0,01 o 0,10): el aviso es esperado
warnings.filterwarnings("ignore", category=InterpolationWarning)

def contrastes(x):
    x = pd.Series(x).dropna()
    a = adfuller(x, regression="c", autolag="AIC",
                 result_object=True)
    k = kpss(x, regression="c", nlags="auto", result_object=True)
    return {"n": len(x), "adf": a.statistic, "p_adf": a.pvalue,
            "kpss": k.statistic, "p_kpss": k.pvalue}

cont = pd.read_parquet("futuros2_continuos.parquet")
ida = np.log(cont[cont["mercado"] == "IDA"]
             .set_index("fecha")["cierre_ajustado"])
pd.DataFrame({"log precio IDA": contrastes(ida),
              "rentabilidad IDA": contrastes(ida.diff())}).T.round(3)

                       n     adf  p_adf    kpss  p_kpss
log precio IDA    6539.0   0.176  0.971  10.817    0.01
rentabilidad IDA  6538.0 -26.296  0.000   0.264    0.10

Con result_object=True, las dos funciones devuelven un objeto con nombres (statistic, pvalue, critical_values) en lugar de una tupla cuya longitud cambia según los argumentos; statsmodels avisa con un FutureWarning de que ese será el comportamiento por defecto, y pedirlo ya lo silencia. kpss devuelve p-valores interpolados en una tabla que va de 0,01 a 0,10; fuera de ese rango los recorta y emite un InterpolationWarning. El aviso no indica un error, sino que el p-valor real es menor que 0,01 o mayor que 0,10; se silencia solo esa categoría para que la salida sea limpia. Los valores críticos del KPSS con constante son 0,347 (10 %), 0,463 (5 %) y 0,739 (1 %).

El logaritmo del precio de IDA, un índice ficticio de renta variable construido como paseo aleatorio con deriva, da un ADF de 0,18 (p = 0,97): no se rechaza la raíz unitaria. El KPSS vale 10,82, muy por encima de 0,463: se rechaza la estacionariedad. Los dos dicen lo mismo con lenguajes opuestos. La rentabilidad diaria da lo contrario: ADF de −26,3 (p ≈ 0) y KPSS de 0,26, por debajo del umbral. Es el caso de libro.

Spreads con respuesta conocida

Los pares de pares2.parquet permiten ir más allá. Todas sus relaciones son en logaritmos: en los pares A, B y C, el spread $u = \ln y - \alpha - \beta \ln x$ es un proceso de Ornstein-Uhlenbeck estacionario con vida media de 5, 15 y 40 sesiones (un porcentaje, con desviación típica del 2 % al 4,5 %); en E, el mismo spread es un paseo aleatorio, aunque los dos activos tienen una correlación de rendimientos de 0,86. Aquí se usan los $\alpha$ y $\beta$ verdaderos de pares2_verdad.csv: estimarlos cambia los valores críticos (Engle-Granger, lección 2.1.5).

pares = pd.read_parquet("pares2.parquet").pivot(
    index="fecha", columns="activo", values="precio")
verdad = pd.read_csv("pares2_verdad.csv").set_index("par")

def spread(par):
    # todas las relaciones son en logaritmos (formulacion = "log")
    v = verdad.loc[par]
    y, x = pares[v["activo_y"]], pares[v["activo_x"]]
    return np.log(y) - v["alfa"] - v["beta"] * np.log(x)

tabla = {}
for par in ["A", "B", "C", "E"]:
    vm = verdad.loc[par, "vida_media_sesiones"]
    s = spread(par)
    tabla[f"{par}, vida media {vm:g}"] = contrastes(s)
    tabla[f"{par}, último año"] = contrastes(s.iloc[-252:])
pd.DataFrame(tabla).T.round(3)

                        n     adf  p_adf   kpss  p_kpss
A, vida media 5    3018.0 -14.929  0.000  0.108   0.100
A, último año       252.0  -4.251  0.001  0.287   0.100
B, vida media 15   3018.0  -8.975  0.000  0.385   0.084
B, último año       252.0  -2.825  0.055  0.133   0.100
C, vida media 40   3018.0  -5.562  0.000  0.826   0.010
C, último año       252.0  -1.643  0.461  0.333   0.100
E, vida media nan  3018.0  -2.012  0.281  4.775   0.010
E, último año       252.0  -2.371  0.150  0.586   0.024

Con los doce años completos (3.018 sesiones), el ADF acierta en los cuatro: rechaza la raíz unitaria en A, B y C (estadísticos de −14,9, −9,0 y −5,6) y no la rechaza en E (p = 0,28). El KPSS acierta en A y en E, pero rechaza la estacionariedad de C (0,826 > 0,739, p ≤ 0,01), que sí es estacionario. Con un solo año de datos, la imagen cambia: el ADF ya no rechaza ni en B (p = 0,055) ni en C (p = 0,46). Con una vida media de 40 sesiones, un año contiene unos seis ciclos de reversión, y eso no basta para distinguirlos de un paseo.

Las cuatro combinaciones

Las dos filas de arriba son la respuesta clara. Las dos de abajo dicen que el contraste no puede responder con estos datos, y son más frecuentes de lo que parece con series cortas o persistentes, que son justo las que interesan para la reversión a la media.

flowchart TB
  A{¿ADF rechaza?} -- sí --> K1{¿KPSS rechaza?}
  A -- no --> K2{¿KPSS rechaza?}
  K1 -- no --> E[Estacionaria]
  K1 -- sí --> P[Muy persistente]
  K2 -- sí --> U[Raíz unitaria]
  K2 -- no --> M[Falta muestra]

h = spread("H")
ruptura = verdad.loc["H", "fecha_ruptura"]
pd.DataFrame({"completa": contrastes(h),
              "antes de la ruptura": contrastes(h[:ruptura].iloc[:-1]),
              }).T.round(3)

                          n    adf  p_adf   kpss  p_kpss
completa             3018.0 -1.448  0.559  6.643   0.010
antes de la ruptura  1875.0 -5.702  0.000  0.407   0.074

Diferenciar: lo que se gana y lo que se pierde

Diferenciar una vez ($d = 1$) convierte el precio en rentabilidad y resuelve la estacionariedad. El coste es la memoria: la rentabilidad de hoy no dice nada del nivel en que está el precio. Se mide con la correlación entre la serie transformada y la original:

r = ida.diff().dropna()
print(f"correlación de la rentabilidad con el log precio: "
      f"{np.corrcoef(r, ida.loc[r.index])[0, 1]:.3f}")

# sobrediferenciar un spread que ya era estacionario
ds = spread("A").diff().dropna()
print(f"autocorrelación en el retardo 1 del spread A: "
      f"{spread('A').autocorr(1):.3f}; de su diferencia: "
      f"{ds.autocorr(1):.3f}")

correlación de la rentabilidad con el log precio: 0.025
autocorrelación en el retardo 1 del spread A: 0.865; de su diferencia: -0.055

La rentabilidad tiene una correlación de 0,025 con el nivel del precio: toda la información sobre dónde está el mercado ha desaparecido. Y diferenciar lo que ya era estacionario estropea la serie: el spread A tiene una autocorrelación de 0,865 en el retardo 1 (la teórica del Ornstein-Uhlenbeck es $0{,}5^{1/5} = 0{,}871$); su diferencia tiene −0,055, cerca del valor de un AR(1) sobrediferenciado ($-(1-\phi)/2 = -0{,}065$): la persistencia positiva que delataba la reversión ha desaparecido y queda una serie casi sin memoria. La regla: diferenciar lo justo, no por reflejo.

Diferenciación fraccional

¿Hay algo entre $d = 0$ (el precio, con toda la memoria y no estacionario) y $d = 1$ (la rentabilidad, estacionaria y sin memoria)? Sí: el operador de diferencia $(1-B)^d$, con $B$ el operador de retardo ($B\,y_t = y_{t-1}$), admite cualquier $d$ real. Su desarrollo binomial da una suma ponderada infinita del presente y el pasado,

Con $d = 1$ los pesos son $1, -1, 0, 0, \ldots$: la diferencia de siempre. Con $0 < d < 1$ todos los pesos pasados son negativos y decaen despacio, como una potencia: la serie recuerda el nivel durante mucho tiempo.

def pesos(d, umbral=1e-4, kmax=100_000):
    w = [1.0]
    for k in range(1, kmax):
        wk = -w[-1] * (d - k + 1) / k
        if abs(wk) < umbral:
            break
        w.append(wk)
    return np.array(w)

for d in [0.1, 0.3, 0.5, 1.0]:
    w = pesos(d)
    print(f"d={d:.1f}: {len(w):4d} pesos; primeros "
          f"{np.round(w[:5], 3).tolist()}")

d=0.1:  503 pesos; primeros [1.0, -0.1, -0.045, -0.028, -0.021]
d=0.3:  388 pesos; primeros [1.0, -0.3, -0.105, -0.06, -0.04]
d=0.5:  200 pesos; primeros [1.0, -0.5, -0.125, -0.062, -0.039]
d=1.0:    2 pesos; primeros [1.0, -1.0]

La versión de ventana fija (fixed-width window fractional differencing, López de Prado, 2018, capítulo 5) corta la suma donde los pesos caen por debajo de un umbral (aquí $10^{-4}$) y aplica la misma ventana a todas las fechas. Así todos los valores tienen el mismo peso de historia y ninguno usa datos futuros: el valor de la fecha $t$ solo combina $y_t, y_{t-1}, \ldots$

def ffd(x, d, umbral=1e-4):
    w = pesos(d, umbral)
    # convolución «valid»: el valor en t usa y_t … y_(t-L+1)
    v = np.convolve(x.to_numpy(), w, mode="valid")
    return pd.Series(v, index=x.index[len(w) - 1:])

# comprobación: d = 1 es la rentabilidad logarítmica
print(np.allclose(ffd(ida, 1.0), ida.diff().dropna()))
x3 = ffd(ida, 0.3)
print(f"d = 0,3: empieza el {x3.index[0]:%Y-%m-%d}, "
      f"{len(ida) - len(x3)} sesiones después que el precio")

True
d = 0,3: empieza el 2001-07-17, 387 sesiones después que el precio

Ahora se barre $d$ y, para cada valor, se calcula el ADF de la serie transformada y su correlación con el logaritmo del precio. Siguiendo a López de Prado, el ADF lleva un solo retardo (maxlag=1, sin selección automática), para que todas las $d$ se comparen con la misma regresión:

filas = []
for d in np.round(np.arange(0, 1.01, 0.05), 2):
    x = ffd(ida, d)
    a = adfuller(x, maxlag=1, regression="c", autolag=None,
                 result_object=True)
    crit = a.critical_values
    filas.append({"d": d, "adf": a.statistic, "p": a.pvalue,
                  "corr": np.corrcoef(x, ida.loc[x.index])[0, 1]})
barrido = pd.DataFrame(filas).set_index("d")
d_min = barrido.index[barrido["adf"] < crit["5%"]].min()
print(f"valor crítico al 5 %: {crit['5%']:.2f}; d mínima: {d_min}")
barrido.loc[[0.0, 0.1, 0.2, 0.25, 0.3, 0.5, 1.0]].round(3)

valor crítico al 5 %: -2.86; d mínima: 0.25
         adf      p   corr
d
0.00   0.393  0.981  1.000
0.10  -0.877  0.795  0.998
0.20  -2.166  0.219  0.987
0.25  -3.152  0.023  0.977
0.30  -4.173  0.001  0.963
0.50 -12.092  0.000  0.845
1.00 -57.513  0.000  0.025

La figura de la cabecera dibuja el barrido completo. El ADF baja de forma continua al subir $d$ y cruza el valor crítico (−2,86) en $d = 0{,}25$ (ADF de −3,15; con $d = 0{,}20$ aún vale −2,17); con esa $d$, la serie es estacionaria según el ADF y conserva una correlación de 0,977 con el logaritmo del precio. La rentabilidad ($d = 1$) es muchísimo más estacionaria de lo necesario (ADF de −57,5) y ha perdido toda la memoria (0,025). Entre las dos hay un abanico de series que conservan información del nivel sin los problemas de un regresor no estacionario.

Lo que esto cambia en un sistema

Tres advertencias. La primera, que «no rechazar» no es «aceptar»: un ADF que no rechaza con un año de datos no prueba que el spread sea un paseo aleatorio. La segunda, que los contrastes de esta lección miden la estacionariedad del pasado: el par H era estacionario hasta que dejó de serlo. La tercera, que barrer $d$, ventanas o pares hasta que uno pase el ADF es una búsqueda múltiple como cualquier otra (lección 2.1.8).

Ejercicios

1. Repite el barrido de $d$ con los futuros de oro (MEA) y de crudo ligero (ENA). ¿Qué $d$ mínima sale en cada uno y qué correlación conserva?

def d_minima(serie, crit=-2.86):
    for d in np.round(np.arange(0, 1.01, 0.05), 2):
        x = ffd(serie, d)
        adf = adfuller(x, maxlag=1, regression="c", autolag=None,
                       result_object=True).statistic
        if adf < crit:
            return d, adf, np.corrcoef(x, serie.loc[x.index])[0, 1]

for m in ["IDA", "MEA", "ENA"]:
    s = np.log(cont[cont["mercado"] == m]
               .set_index("fecha")["cierre_ajustado"])
    d, adf, c = d_minima(s)
    print(f"{m}: d mínima {d:.2f}, ADF {adf:.2f}, correlación {c:.3f}")

IDA: d mínima 0.25, ADF -3.15, correlación 0.977
MEA: d mínima 0.15, ADF -3.05, correlación 0.981
ENA: d mínima 0.10, ADF -3.40, correlación 0.988

2. El ADF tiene más potencia cuanto más larga es la muestra. Para el spread B (vida media de 15 sesiones), calcula el p-valor del ADF con los últimos 126, 252, 504 y 1.008 días. ¿Desde qué longitud rechaza la raíz unitaria al 5 %?

sb = spread("B")
for n in [126, 252, 504, 1008]:
    p = adfuller(sb.iloc[-n:], regression="c", autolag="AIC",
                 result_object=True).pvalue
    print(f"{n:5d} sesiones: p = {p:.4f}")

  126 sesiones: p = 0.2449
  252 sesiones: p = 0.0548
  504 sesiones: p = 0.0090
 1008 sesiones: p = 0.0000

3. Un compañero propone una regla de reversión: comprar IDA cuando su precio está un 10 % por debajo de su media de 200 sesiones. ¿Qué dicen los contrastes de esta lección sobre esa regla, y qué contrastarías en su lugar?