Pairs trading y arbitraje estadístico: distancia, cointegración y Kalman.

Un diferencial que vuelve

El pairs trading es la forma más simple del arbitraje estadístico. Se buscan dos activos cuyo diferencial (la diferencia entre el logaritmo del precio de uno y una cantidad fija del logaritmo del otro) oscila alrededor de un nivel estable. Cuando se separa mucho, se compra el que se ha quedado atrás y se vende el que se ha adelantado; cuando vuelve, se cierra. La posición es casi neutral al mercado: si todo sube o baja a la vez, las dos patas se compensan.

Si $s_t$ es estacionario, $y$ y $x$ están cointegrados (lección 2.1.5): cada uno es un paseo aleatorio, pero una combinación lineal de los dos no lo es. Esa es la condición que hace funcionar el par. Que los rendimientos de los dos estén muy correlados no basta: dos paseos aleatorios con rendimientos correlados al 0,85 pueden separarse para siempre.

Los logaritmos no son un detalle. En logaritmos el diferencial es un porcentaje: una σ de 0,02 es un desajuste típico del 2 % entre las dos patas, cuesten 20 o 200 USD. En precios ($y_t - \alpha - \beta x_t$), la σ del diferencial en dólares crece con el nivel de los precios, y un z calculado con la σ de la formación deja de significar lo mismo cuando los precios se duplican.

Gatev, Goetzmann y Rouwenhorst estudiaron en datos reales de EE. UU. de 1962 a 1997 una regla que emparejaba acciones por la mínima distancia entre sus precios normalizados, y encontraron rendimientos anualizados en exceso de hasta el 12 % en carteras autofinanciadas de los mejores pares, que superaban una estimación conservadora de los costes en la mayor parte del periodo11Nota 1 · E. Gatev, W. N. Goetzmann y K. G. Rouwenhorst, «Pairs Trading: Performance of a Relative-Value Arbitrage Rule», Review of Financial Studies 19(3), 2006, 797-827. DOI 10.1093/rfs/hhj020. Versión de trabajo NBER w7032.. Esta lección aplica ese método y el de cointegración a datos donde la respuesta está escrita.

flowchart TB
  U[Universo] --> F[Formación 2014-17]
  F --> S[Selección de pares]
  S --> P[Parámetros fijos]
  P --> N[Negocia 2018-25]
  N --> R[Resultado por par]

Dieciséis activos, ciento veinte pares

import itertools
import warnings

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.stattools import coint

pares = pd.read_parquet("pares2.parquet")
w = pares.pivot(index="fecha", columns="activo", values="precio")
lw = np.log(w)                         # las relaciones son en log
verdad = pd.read_csv("pares2_verdad.csv").set_index("par")
FIN_FORMACION = "2017-12-31"
form, lform = w.loc[:FIN_FORMACION], lw.loc[:FIN_FORMACION]
print(w.shape[1], "activos,", len(form), "sesiones de formación,",
      len(w.loc["2018":]), "de negociación")
print("formulación:", verdad["formulacion"].unique().tolist())
verdad[["activo_y", "activo_x", "tipo", "beta",
        "vida_media_sesiones", "sigma_spread"]]

16 activos, 1007 sesiones de formación, 2011 de negociación
formulación: ['log']
    activo_y activo_x                      tipo  beta  vida_media_sesiones  sigma_spread
par
A       PZAY     PZAX               cointegrado   1.5                  5.0         0.020
B       PZBY     PZBX               cointegrado   0.8                 15.0         0.030
C       PZCY     PZCX               cointegrado   1.1                 40.0         0.045
D       PZDY     PZDX             beta_variable   NaN                 10.0         0.025
E       PZEY     PZEX  correlado_no_cointegrado   1.0                  NaN           NaN
F       PZFY     PZFX  correlado_no_cointegrado   0.9                  NaN           NaN
G       PZGY     PZGX  correlado_no_cointegrado   1.2                  NaN           NaN
H       PZHY     PZHX                   ruptura   0.9                 10.0         0.025

Dieciséis activos, 1.007 sesiones de formación (2014-2017) y 2.011 de negociación (2018-2025). La tabla de la verdad solo se mira para corregir: la selección y los parámetros saldrán de la formación como si no existiera. Todas las relaciones son en logaritmos (formulacion = log): ln y = α + β·ln x + u, con u estacionario de σ entre 0,020 y 0,045 en los pares verdaderos y un paseo aleatorio en E, F y G. Por eso los contrastes, las regresiones y el diferencial se calculan sobre lw, los logaritmos de los precios.

Selección por distancia

El método de la distancia normaliza los precios de cada activo a 1 el primer día de la formación y, para cada pareja, suma los cuadrados de la diferencia entre sus trayectorias. Los pares con menor suma son los que más se han parecido. No hay modelo ni contraste: es una medida de parecido.

norm = form / form.iloc[0]
todos = list(itertools.combinations(w.columns, 2))
verdaderos = {(r.activo_x, r.activo_y) for r in verdad.itertuples()}

def etiqueta(a, b):
    for r in verdad.itertuples():
        if {a, b} == {r.activo_x, r.activo_y}:
            return f"{r.Index} ({r.tipo})"
    return "cruzado"

dist = pd.Series({(a, b): ((norm[a] - norm[b]) ** 2).sum()
                  for a, b in todos}).sort_values()
top = dist.head(10).round(2).rename("distancia").to_frame()
top["par"] = [etiqueta(a, b) for a, b in top.index]
top

           distancia                           par
PZDX PZDY       0.70             D (beta_variable)
PZBX PZBY       1.31               B (cointegrado)
PZHX PZHY       1.89                   H (ruptura)
PZCX PZCY       3.24               C (cointegrado)
PZDX PZEX       3.83                       cruzado
PZDY PZEX       4.89                       cruzado
PZEX PZEY       6.75  E (correlado_no_cointegrado)
PZFX PZHX       6.93                       cruzado
     PZHY       7.57                       cruzado
PZBY PZHY       9.01                       cruzado

Los cuatro pares más cercanos son pares de verdad (D, B, H y C). El quinto y el sexto son cruzados (PZDX con PZEX y PZDY con PZEX: activos de pares distintos que comparten el factor de mercado y en estos cuatro años fueron juntos), y el séptimo es E, correlado y no cointegrado. Falta A, el par con la reversión más rápida del universo: su beta es 1,5, así que y se mueve un 50 % más que x y sus trayectorias normalizadas se separan aunque el diferencial en logaritmos vuelva cada pocos días. La distancia supone, sin decirlo, una cobertura de uno a uno. Mide parecido: no distingue «se parecen» de «están obligados a volver», y no ve un par cuya cobertura no sea uno a uno.

Selección por cointegración

El contraste de Engle-Granger estima por mínimos cuadrados $\ln y = \alpha + \beta \ln x + s$ y contrasta si el residuo $s$ tiene raíz unitaria con un ADF con valores críticos propios (lección 2.1.5). statsmodels.tsa.stattools.coint lo hace en una llamada. Se aplica a los logaritmos de los 120 pares del universo en la formación:

filas = []
for a, b in todos:
    _, p, _ = coint(lform[a], lform[b])
    filas.append({"y": a, "x": b, "p": p, "par": etiqueta(a, b)})
eg = pd.DataFrame(filas).sort_values("p")
print("p < 0,05:", (eg["p"] < 0.05).sum(), "de", len(eg),
      "  esperados por azar:", round(0.05 * len(eg), 1))
print("p < 0,05/120 (Bonferroni):", (eg["p"] < 0.05 / 120).sum())
eg.head(8).round(4)

p < 0,05: 18 de 120   esperados por azar: 6.0
p < 0,05/120 (Bonferroni): 4
        y     x       p                par
0    PZAX  PZAY  0.0000    A (cointegrado)
29   PZBX  PZBY  0.0000    B (cointegrado)
75   PZDX  PZDY  0.0000  D (beta_variable)
119  PZHX  PZHY  0.0003        H (ruptura)
109  PZFX  PZHY  0.0008            cruzado
101  PZEY  PZGX  0.0011            cruzado
49   PZBY  PZFY  0.0034            cruzado
53   PZBY  PZHY  0.0067            cruzado

Dieciocho pares pasan el 5 %, tres veces los 6 que se esperan por azar entre 120 contrastes. Los cuatro primeros son de verdad (A, B, D y H, con p de 0,0003 o menos), y la corrección de Bonferroni (p < 0,05/120 ≈ 0,0004) deja exactamente esos cuatro. Los otros catorce son falsos positivos, casi todos cruzados.

¿Por qué tantos? El contraste no está mal: el 5 % es una frecuencia media sobre muchas muestras, y en una sola los 120 contrastes no son independientes, porque comparten series. Si PZFX parece cointegrado con PZHY, también lo parecerá con PZHX, que va pegado a PZHY; en la tabla, PZBY aparece con PZFY y con PZHY. Los falsos positivos llegan en racimos, y su número varía de una muestra a otra mucho más que el de 120 contrastes independientes. Contarlos en las otras dos ventanas de cuatro años lo enseña:

reales = {frozenset((r.activo_x, r.activo_y))
          for r in verdad.itertuples() if not r.tipo.startswith("corr")}
for a, b in (("2014", "2017"), ("2018", "2021"), ("2022", "2025")):
    tramo = lw.loc[a:b]
    sig = [frozenset(p) for p in todos
           if coint(tramo[p[0]], tramo[p[1]])[1] < 0.05]
    falsos = sum(p not in reales for p in sig)
    print(f"{a}-{b}: {len(sig):2d} con p < 0,05, {falsos:2d} falsos")

2014-2017: 18 con p < 0,05, 14 falsos
2018-2021:  8 con p < 0,05,  4 falsos
2022-2025:  5 con p < 0,05,  2 falsos

Catorce falsos en 2014-2017, cuatro en 2018-2021 y dos en 2022-2025: una media de unos siete por ventana, cerca de los ≈ 5,8 que se esperan entre los 115 pares que no son verdaderos, pero con una ventana que da más del doble. La cifra «esperados por azar» es una media, no un techo. El par C, que está cointegrado, no aparece entre los significativos de la formación: con una vida media de 40 sesiones, cuatro años no bastan para que el contraste distinga su diferencial de un paseo aleatorio. Un contraste con poca potencia también falla hacia el otro lado.

Los ocho pares designados, uno por uno, con la correlación de sus rendimientos, el contraste en la formación y en toda la muestra, y la σ y la vida media estimadas del diferencial. La vida media sale de una regresión AR(1) del diferencial: $\Delta s_t = \lambda s_{t-1} + e_t$, y $h = -\ln 2 / \ln(1 + \lambda)$ (lección 2.1.5).

def formacion(y, x):
    m = sm.OLS(lform[y], sm.add_constant(lform[x])).fit()
    a, b = m.params.iloc[0], m.params.iloc[1]
    s = lform[y] - a - b * lform[x]
    ar = sm.OLS(s.diff().dropna(), s.shift().dropna()).fit()
    lam = ar.params.iloc[0]
    vida = -np.log(2) / np.log(1 + lam) if lam < 0 else np.inf
    return a, b, s.mean(), s.std(), vida

resumen = []
for par, r in verdad.iterrows():
    y, x = r["activo_y"], r["activo_x"]
    a, b, mu, sd, vida = formacion(y, x)
    resumen.append({
        "par": par,
        "corr_rend": lform[[y, x]].diff().corr().iloc[0, 1],
        "p_formacion": coint(lform[y], lform[x])[1],
        "p_total": coint(lw[y], lw[x])[1],
        "beta_est": b, "beta_real": r["beta"], "sd_est": sd,
        "vida_est": vida, "vida_real": r["vida_media_sesiones"]})
pd.DataFrame(resumen).set_index("par").round(3)

     corr_rend  p_formacion  p_total  beta_est  beta_real  sd_est  vida_est  vida_real
par
A        0.900        0.000    0.000     1.503        1.5   0.019     4.423        5.0
B        0.724        0.000    0.000     0.981        0.8   0.025     9.331       15.0
C        0.850        0.187    0.000     1.092        1.1   0.045    43.827       40.0
D        0.830        0.000    0.035     1.204        NaN   0.024     9.749       10.0
E        0.825        0.207    0.237     1.237        1.0   0.057    60.003        NaN
F        0.821        0.314    0.510     0.757        0.9   0.056    69.790        NaN
G        0.821        0.308    0.339     2.117        1.2   0.105    58.862        NaN
H        0.677        0.000    0.219     0.867        0.9   0.027    12.339       10.0

Lo que dice la tabla, par a par:

• La correlación no separa nada. Todos los pares, verdaderos y falsos, tienen rendimientos correlados entre 0,68 y 0,90, y dos de los verdaderos (H, 0,68, y B, 0,72) están entre los menos correlados.
• A, B, D y H pasan el contraste en la formación con p ≈ 0. A y H tienen betas casi exactas (1,503 frente a 1,5; 0,867 frente a 0,9) y vidas medias cercanas (4,4 y 12,3 sesiones, frente a 5 y 10). B se estima peor: 0,981 frente a 0,8, y una vida media de 9 sesiones frente a 15. La σ estimada del diferencial (sd_est) queda cerca de la verdadera en todos (0,019, 0,025 y 0,027, frente a 0,020, 0,030 y 0,025).
• C no pasa en la formación (p = 0,19) y sí en toda la muestra (p ≈ 0): la vida media estimada, 44 sesiones, es buena, pero con pocos ciclos el contraste no tiene potencia.
• D, con beta variable, pasa en la formación con una vida media de 9,7 sesiones, casi la verdadera: en esos cuatro años su beta se movió poco y una beta fija de 1,204 sirve. En toda la muestra el contraste ya duda (p = 0,035): la beta ha seguido moviéndose y el diferencial con una beta fija revierte mucho más despacio que el verdadero (el README mide una vida media de 38 sesiones con toda la muestra).
• E, F y G no están cointegrados en ninguno de los dos periodos, como debe ser (p de 0,21 a 0,51). Sus «vidas medias» de 59 a 70 sesiones son lo que da una regresión AR(1) sobre un paseo aleatorio en una muestra finita: un número que parece razonable y no significa nada. La beta de G, 2,117 frente a 1,2, tampoco.
• H pasa en la formación (p ≈ 0) y no en la muestra total (p = 0,22): la ruptura de 2021 está en la parte que no se podía ver al elegirlo.

Reglas de entrada y salida

Con los parámetros de la formación fijos, la regla clásica abre cuando el diferencial se separa más de dos desviaciones típicas y cierra cuando vuelve a su media. Dos salidas más protegen de que la relación haya dejado de existir: un stop por nivel (|z| > 4) y un stop por tiempo (más de 60 sesiones abierta, unas cuatro veces la vida media de un par típico). Tras un stop, la regla no vuelve a entrar hasta que el diferencial regrese a |z| < 0,5: sin esa condición, volvería a abrir al día siguiente en el mismo sitio.

stateDiagram-v2
  [*] --> PLANO
  PLANO --> LARGO: z menor de -2
  PLANO --> CORTO: z mayor de 2
  LARGO --> PLANO: z cruza 0
  CORTO --> PLANO: z cruza 0
  LARGO --> BLOQUEO: stop
  CORTO --> BLOQUEO: stop
  BLOQUEO --> PLANO: abs z menor de 0,5

«Largo el diferencial» es comprar $V$ USD de $y$ y vender $\beta V$ USD de $x$: con el diferencial en logaritmos, esa cartera gana $V \cdot \Delta s$ (para movimientos pequeños), porque $\Delta \ln y$ y $\Delta \ln x$ son rendimientos. El tamaño se fija para que una desviación típica del diferencial valga 1.000 USD: $V = 1.000 / \sigma_s$, unos 50.000 USD de nominal en $y$ con una σ de 0,02. Las acciones se calculan con los precios del cierre en que se abre y no se tocan hasta cerrar. Los costes son por acción y lado (0,02 USD: comisión de 0,005, media horquilla de 0,01 y deslizamiento de 0,005, para acciones de un céntimo de tick) más un 0,5 % anual de préstamo sobre el nominal vendido. Como en todo el nivel, la señal del cierre de t se ejecuta al cierre de t+1.

COSTE_ACCION = 0.02
PRESTAMO = 0.005

def maquina(z, entrada=2.0, stop=None, tiempo=None):
    est, cur, dias, armado = np.zeros(len(z)), 0, 0, True
    for i, v in enumerate(z):
        if not armado and abs(v) < 0.5:
            armado = True
        if cur == 0 and armado:
            cur = -1 if v > entrada else (1 if v < -entrada else 0)
            dias = 0
        elif cur != 0:
            dias += 1
            if cur * v > 0:                       # cruzó la media
                cur = 0
            elif (stop and abs(v) > stop) or (tiempo and dias > tiempo):
                cur, armado = 0, False
        est[i] = cur
    return est

def opera(y, x, z, beta, sd, **reglas):
    est = pd.Series(maquina(z.to_numpy(), **reglas), index=z.index)
    pos = est.shift(1).fillna(0.0)               # al cierre siguiente
    Y, X = w[y].loc[z.index], w[x].loc[z.index]
    usd = pos * 1000 / sd                        # nominal de y
    nuevo = pos.diff().fillna(pos) != 0      # acciones fijas al abrir
    acc_y = (usd / Y).where(nuevo).ffill().where(pos != 0, 0.0)
    acc_x = (-usd * beta / X).where(nuevo).ffill().where(pos != 0, 0.0)
    bruto = acc_y.shift(1) * Y.diff() + acc_x.shift(1) * X.diff()
    giro = (acc_y.diff().abs() + acc_x.diff().abs()) * COSTE_ACCION
    corto = (acc_y.clip(upper=0).abs() * Y
             + acc_x.clip(upper=0).abs() * X) * PRESTAMO / 252
    neto = (bruto - giro - corto).fillna(0.0)
    return est, neto, bruto.fillna(0.0)

def par(nombre, **reglas):
    r = verdad.loc[nombre]
    y, x = r["activo_y"], r["activo_x"]
    a, b, mu, sd, _ = formacion(y, x)
    z = ((lw[y] - a - b * lw[x] - mu) / sd).loc["2018":]
    est, neto, bruto = opera(y, x, z, b, sd, **reglas)
    return z, est, neto, bruto

z_b, est_b, neto_b, _ = par("B")
cambios = est_b[est_b.diff().fillna(est_b) != 0]
print(cambios.map({1.0: "largo", -1.0: "corto", 0.0: "cierra"})
      .head(8).to_string())
print(f"neto del par B, 2018-2025: {neto_b.sum():,.0f} USD")

fecha
2018-05-02     largo
2018-09-18    cierra
2018-10-09     largo
2019-03-07    cierra
2019-03-28     largo
2019-04-16    cierra
2019-05-03     corto
2019-07-30    cierra
neto del par B, 2018-2025: 26,962 USD

El par B tiene pocas operaciones y largas: la primera dura de mayo a septiembre de 2018 y la segunda, de octubre de 2018 a marzo de 2019, algo raro para un diferencial con vida media de 15 sesiones. La causa es el error de la beta: 0,981 en vez de 0,8. Con la beta estimada, el diferencial es $u_t - 0{,}181 \ln x_t$ más una constante, y cuando $x$ se mueve, ese término deriva y el diferencial se queda meses a un lado de la media de la formación: una subida de x del 25 % mueve $\ln x$ en 0,22 y el diferencial en 0,04, más de una σ y media. El par gana 26.962 USD netos en ocho años, sin stops. En datos reales el problema es el mismo: un error pequeño en el hedge ratio se convierte en una tendencia del diferencial cuando los precios se mueven mucho.

Resultados por par: verdaderos, falsos y stops

def caida(neto):
    acum = neto.cumsum()
    return (acum - acum.cummax()).min()

filas = {}
for nombre in verdad.index:
    _, e0, n0, b0 = par(nombre)
    _, e1, n1, _ = par(nombre, stop=4.0, tiempo=60)
    filas[nombre] = {
        "tipo": verdad.loc[nombre, "tipo"],
        "operaciones": int((e0.diff().fillna(e0) != 0).sum() // 2),
        "bruto": b0.sum(), "neto": n0.sum(), "caida": caida(n0),
        "neto_stops": n1.sum(), "caida_stops": caida(n1)}
res = pd.DataFrame(filas).T
for c in ["bruto", "neto", "caida", "neto_stops", "caida_stops"]:
    res[c] = res[c].astype(float).round(-2).astype(int)
res

                       tipo operaciones  bruto   neto  caida  neto_stops  caida_stops
A               cointegrado          28  59900  57200  -2900       57200        -2900
B               cointegrado          12  29100  27000 -11700        7300        -5000
C               cointegrado           6  14000  13300  -2100        8100        -2100
D             beta_variable          15  34700  32700 -15200       23100        -5400
E  correlado_no_cointegrado           1 -10000 -11700 -22500       -1500        -2600
F  correlado_no_cointegrado           1  -3800  -4400 -10100        2900        -1600
G  correlado_no_cointegrado           1 -10200 -11400 -20300       -1400        -3000
H                   ruptura           6   5300   3800 -22300        6500        -5900

Las cifras son dólares por cada 1.000 USD de σ del diferencial, de 2018 a 2025, netos de costes:

• Los pares verdaderos ganan: A, 57.200 con 28 operaciones; B, 27.000; C, 13.300 con solo 6 (su vida media de 40 sesiones da pocos ciclos). D, el de la beta variable, también gana 32.700 con la beta fija de la formación: es un par de verdad, y esa beta le sirve durante años, aunque con una caída de 15.200.
• Los falsos pierden sin stops: E, −11.700; F, −4.400; G, −11.400, con caídas de hasta 22.500. Cada uno hace una sola operación en ocho años: abre y se queda abierta mientras el diferencial se aleja. Es exactamente el comportamiento de un paseo aleatorio operado como si revirtiera.
• Los stops recortan la cola, no mejoran la media. Con |z| > 4 y 60 sesiones como límite, los tres falsos quedan entre −1.500 y +2.900, con caídas de 1.600 a 3.000: ahorran 27.500 USD entre los tres. En los verdaderos, los stops cuestan casi lo mismo, 34.500 entre B, C y D (A no cambia): cortan operaciones que habrían vuelto. B es el caso extremo, de 27.000 a 7.300, porque con su beta mal estimada las operaciones duran meses y el stop por tiempo las cierra antes de que vuelvan.
• H pasa de 3.800 a 6.500: la sección siguiente lo explica.

La figura reúne los dos resultados. Sumados, los stops casi no cambian el total de esta muestra; lo que cambian es la peor caída de un par, que baja de 22.500 a 5.900. Un stop en pares es un seguro con prima: se paga en los pares que funcionan y protege de los que no, que en datos reales no vienen etiquetados.

El par roto

El par H cumple todo lo que se le puede pedir a un par en la formación: cointegrado con p ≈ 0, vida media corta y hedge ratio estable. Durante la negociación funciona tres años. El 15 de junio de 2021 su diferencial salta cinco desviaciones típicas y desde ese día es un paseo aleatorio: la relación ya no existe. Es lo que pasa en datos reales cuando una de las dos empresas cambia de negocio, es adquirida o sufre un problema propio.

z_h, e_h0, n_h0, _ = par("H")
_, e_h1, n_h1, _ = par("H", stop=4.0, tiempo=60)
anual = pd.DataFrame({"sin_stops": n_h0, "con_stops": n_h1})
print(anual.groupby(anual.index.year).sum().round(-2).astype(int).T
      .to_string())
print("z máximo tras la ruptura:",
      round(z_h.loc["2021-06-15":].max(), 1))
print("posición sin stops al final:", e_h0.iloc[-1],
      "  con stops:", e_h1.iloc[-1])
print(pd.DataFrame({"z": z_h, "sin_stops": e_h0, "con_stops": e_h1})
      .loc["2021-06-11":"2021-06-17"].round(2).to_string())

fecha      2018  2019  2020   2021  2022  2023  2024  2025
sin_stops  3100  1000  8400 -15900  6000  2100  2300 -3200
con_stops  3700  1000  4900  -3200     0     0     0     0
z máximo tras la ruptura: 18.4
posición sin stops al final: -1.0   con stops: 0.0
               z  sin_stops  con_stops
fecha
2021-06-11  2.02       -1.0       -1.0
2021-06-14  1.96       -1.0       -1.0
2021-06-15  6.49       -1.0        0.0
2021-06-16  6.42       -1.0        0.0
2021-06-17  6.20       -1.0        0.0

Sin stops, H gana 12.500 USD de 2018 a 2020 y pierde 15.900 en 2021. La regla estaba corta antes de la ruptura (el diferencial rondaba +2), y el 15 de junio el z pasa de 1,96 a 6,49 en una sola sesión. Desde ahí sigue alejándose, hasta 18,4 desviaciones típicas, con la posición abierta hasta el final de 2025. Lo que la versión sin stops gana o pierde de 2022 a 2025 (+6.000, +2.100, +2.300 y −3.200) es ruido: una posición fija sobre un paseo aleatorio.

Con stops, el de nivel salta en el cierre del 15 de junio (la columna con_stops pasa a 0 ese día) y la posición se cierra en la sesión siguiente. No evitó la pérdida del salto, que pasó por encima del stop de 4 en una sola sesión: 2021 cierra con −3.200. Lo que evitó es todo lo que vino después: el bloqueo impide volver a entrar, porque el diferencial no regresa a |z| < 0,5, y de 2022 a 2025 no hay nada. La figura enseña los dos caminos: el diferencial antes y después de la ruptura, y el resultado acumulado con y sin stops.

Hedge ratio variable: el filtro de Kalman

El par D tiene una beta que se mueve como un paseo aleatorio reflejado entre 0,6 y 1,6, con un paso diario de 0,006, y la relación gira alrededor del primer precio de x, $x_0$ (README de los datos):

En la formación, una beta fija pasa el contraste, pero en toda la muestra ya duda (tabla de los ocho pares). Un filtro de Kalman trata la beta como un estado oculto que cambia un poco cada día y la estima de forma recursiva, con los datos hasta la víspera (lección 2.1.6). Hay dos maneras de escribir la observación, y no son equivalentes:

• Centrado: $\ln y_t = c + \beta_t (\ln x_t - \ln x_0) + e_t$. Aquí la constante $c = \alpha + \ln x_0$ es fija de verdad, porque el generador hace girar la relación justo en $x_0$.
• Alfa fijo, la forma habitual: $\ln y_t = \alpha + \beta_t \ln x_t + e_t$ con $\alpha$ fija. Cuando la beta cambia, la constante verdadera, $\alpha + (1 - \beta_t)\ln x_0$, también cambia, y el filtro que la supone fija está mal especificado; la lección 2.1.6 lo mide con detalle.

En los dos, la beta sigue un paseo aleatorio, $\beta_t = \beta_{t-1} + \eta_t$ con $\eta_t \sim N(0, q)$, y la varianza del ruido de observación es la del diferencial de la formación ($0{,}024^2$). La varianza $q$ del paso es la decisión clave; la del generador es $0{,}006^2 = 3{,}6 \times 10^{-5}$.

beta_real = pd.read_csv("pares2_beta_real.csv", parse_dates=["fecha"]
                        ).set_index("fecha")["beta_real"]
_, _, _, sd_D, _ = formacion("PZDY", "PZDX")

def kalman(y, X, q, r=sd_D ** 2):
    th, P = np.zeros(2), np.eye(2) * 100.0
    Q = np.diag([0.0, q])                 # la constante no se mueve
    n = len(y)
    beta, error = np.zeros(n), np.zeros(n)
    for t in range(n):
        P = P + Q
        H = X[t]
        S = H @ P @ H + r
        beta[t] = th[1]                       # antes de ver y_t
        error[t] = y[t] - H @ th
        K = P @ H / S
        th = th + K * error[t]
        P = P - np.outer(K, H @ P)
    return beta, error

yD, lx = lw["PZDY"].to_numpy(), lw["PZDX"].to_numpy()
uno = np.ones(len(lx))
modelos = {"centrado": np.column_stack([uno, lx - lx[0]]),
           "alfa fijo": np.column_stack([uno, lx])}
filas = []
for nombre, X in modelos.items():
    for q in (1e-7, 1e-6, 1e-5, 0.006 ** 2):
        beta_k, err = kalman(yD, X, q)
        rho = pd.Series(err, index=w.index).loc["2015":].autocorr()
        e_b = (pd.Series(beta_k, w.index) - beta_real).loc["2015":]
        filas.append({"modelo": nombre, "q": f"{q:.1e}",
                      "error_beta": e_b.abs().mean(), "ac1_dif": rho,
                      "vida_dif": -np.log(2) / np.log(rho)})
pd.DataFrame(filas).round(3)

      modelo        q  error_beta  ac1_dif  vida_dif
0   centrado  1.0e-07       0.047    0.955    14.965
1   centrado  1.0e-06       0.043    0.930     9.604
2   centrado  1.0e-05       0.047    0.894     6.186
3   centrado  3.6e-05       0.055    0.862     4.671
4  alfa fijo  1.0e-07       0.103    0.891     6.027
5  alfa fijo  1.0e-06       0.109    0.788     2.906
6  alfa fijo  1.0e-05       0.092    0.536     1.110
7  alfa fijo  3.6e-05       0.086    0.328     0.621

La tabla enseña dos cosas. La primera es la especificación: con el modelo centrado, el error medio de la beta queda entre 0,043 y 0,055 sobre valores entre 1,0 y 1,5; con la constante fija, entre 0,086 y 0,109, el doble, porque el filtro compensa la constante equivocada moviendo la beta. La segunda es el compromiso entre seguir la beta y conservar el diferencial. En el modelo centrado, al subir q la autocorrelación del residuo baja poco (de 0,955 a 0,862) y su vida media, de 15 a 4,7 sesiones, rodea la verdadera de 10. En el de alfa fijo cae de 0,891 a 0,328: con la q del generador el residuo tiene una vida media de 0,6 sesiones. El filtro atribuye cada desviación del diferencial a un cambio de beta y se come la reversión que se quería operar.

Un detalle: en el modelo centrado el mejor error de la beta (0,043) no está en la q del generador, sino en una cien veces menor. El filtro supone que el ruido de observación es independiente de un día a otro, y el diferencial verdadero es un proceso con memoria (vida media de 10 sesiones): una q grande confunde esas rachas con cambios de beta. La figura compara la beta real con la del filtro centrado de $q = 10^{-6}$ y con la beta fija de la formación.

def opera_dif(dif, beta_t):
    sd = dif.rolling(250, min_periods=60).std().shift(1)
    z = (dif / sd).loc["2018":]
    est, neto, _ = opera("PZDY", "PZDX", z, beta_t.loc[z.index],
                         sd.loc[z.index], stop=4.0, tiempo=60)
    n_op = int((est.diff().fillna(est) != 0).sum() // 2)
    return n_op, round(neto.sum(), -2)

a_f, b_f, *_ = formacion("PZDY", "PZDX")
lyD, lxD = lw["PZDY"], lw["PZDX"]
u_real = lyD - 0.3 - lx[0] - beta_real * (lxD - lx[0])  # la verdad
variantes = {"beta fija (OLS)": (lyD - a_f - b_f * lxD,
                                 pd.Series(b_f, index=w.index)),
             "beta real": (u_real, beta_real)}
for nombre, q in (("centrado", 1e-6), ("centrado", 0.006 ** 2),
                  ("alfa fijo", 0.006 ** 2)):
    beta_k, err = kalman(yD, modelos[nombre], q)
    variantes[f"{nombre} q={q:.1e}"] = (pd.Series(err, w.index),
                                        pd.Series(beta_k, w.index))
pd.DataFrame({k: opera_dif(*v) for k, v in variantes.items()},
             index=["operaciones", "neto_usd"]).T

                     operaciones  neto_usd
beta fija (OLS)             11.0   18400.0
beta real                   23.0   56100.0
centrado q=1.0e-06          22.0   39400.0
centrado q=3.6e-05          36.0   29700.0
alfa fijo q=3.6e-05         78.0   12800.0

Con la beta real, el diferencial es el proceso de Ornstein-Uhlenbeck del generador y la regla gana 56.100 USD con 23 operaciones: es el techo, inalcanzable porque usa la verdad. Con la beta fija de la formación, 18.400 en 11 operaciones: funciona mientras la beta no se aleja mucho de la de 2014-2017. El Kalman centrado con $q = 10^{-6}$ gana 39.400 con 22 operaciones, y con la q del generador, 29.700 con 36: sigue la beta con más ruido y opera más. El de alfa fijo con esa misma q gana 12.800 con 78 operaciones, menos que la beta fija: opera un diferencial sin memoria y paga costes en cada vuelta. Entre el techo y el mejor filtro quedan 16.700 USD: lo que cuesta no conocer la beta, en el mejor de los casos.

Una cartera de pares elegidos en la formación

En la práctica nadie conoce pares2_verdad.csv: se eligen los pares con los datos de la formación y se operan todos. La cartera honesta es la de los dieciocho pares que pasan Engle-Granger al 5 % en la formación, falsos positivos incluidos, con stops:

def par_libre(y, x, **reglas):
    a, b, mu, sd, _ = formacion(y, x)
    z = ((lw[y] - a - b * lw[x] - mu) / sd).loc["2018":]
    return opera(y, x, z, b, sd, **reglas)[1]

elegidos = eg[eg["p"] < 0.05]
cartera = pd.DataFrame({f"{r.y}-{r.x}": par_libre(r.y, r.x, stop=4.0,
                                                   tiempo=60)
                        for r in elegidos.itertuples()})
total = (cartera.sum().round(-2).astype(int)
         .rename("neto_usd").to_frame())
total["par"] = elegidos["par"].to_numpy()
print(total.to_string())
falso = total["par"].str.contains("cruzado|correlado").to_numpy()
print(f"cartera: {cartera.sum().sum():,.0f} USD;  verdaderos: "
      f"{cartera.loc[:, ~falso].sum().sum():,.0f};  falsos: "
      f"{cartera.loc[:, falso].sum().sum():,.0f} ({falso.sum()})")

           neto_usd                           par
PZAX-PZAY     62700               A (cointegrado)
PZBX-PZBY     10800               B (cointegrado)
PZDX-PZDY     14500             D (beta_variable)
PZHX-PZHY      8000                   H (ruptura)
PZFX-PZHY     10300                       cruzado
PZEY-PZGX       400                       cruzado
PZBY-PZFY     -1600                       cruzado
PZBY-PZHY      4900                       cruzado
PZBY-PZHX       600                       cruzado
PZBX-PZHY     -8800                       cruzado
PZEX-PZEY      5000  E (correlado_no_cointegrado)
PZBY-PZFX       200                       cruzado
PZFX-PZHX      2600                       cruzado
PZEX-PZGX      1100                       cruzado
PZBX-PZFY     -2700                       cruzado
PZEY-PZGY     -1300                       cruzado
PZCX-PZGY     -2000                       cruzado
PZDY-PZGX     -2000                       cruzado
cartera: 102,610 USD;  verdaderos: 96,047;  falsos: 6,564 (14)

La cartera gana 102.610 USD, y casi todo viene de los cuatro pares verdaderos (96.047; A solo aporta 62.700). Los catorce falsos suman +6.564: con stops, un par falso es sobre todo ruido con costes, unas veces positivo (+10.300 en PZFX-PZHY) y otras negativo (−8.800 en PZBX-PZHY). En este universo los falsos positivos no hunden la cartera porque los stops los mantienen pequeños, pero tampoco aportan nada: son catorce de dieciocho posiciones, con su capital, su préstamo y su riesgo de cola, para un resultado que es cero más ruido. Con más candidatos y sin stops serían la mayoría de la cartera y la mayoría de las pérdidas.

A y B ganan aquí más que en la tabla de los ocho pares (62.700 frente a 57.200 y 10.800 frente a 7.300), y D menos (14.500 frente a 23.100), porque el contraste los eligió con el orden (x, y) invertido: la beta de una regresión de x sobre y no es la inversa exacta de la de y sobre x, y elegir el orden también es una decisión del investigador que mueve el resultado.

Lo que esto cambia en tus backtests

• Formación y negociación separadas. La beta, la media, la σ y la selección salen de un periodo y se operan en otro. Seleccionar los pares con toda la muestra es elegir los que ya se sabe que funcionaron.
• Cuenta los pares que probaste. Con cien o mil candidatos, un contraste al 5 % da decenas de pares «cointegrados» por azar; Bonferroni o FDR antes de operar (lección 2.1.8).
• Correlación no es cointegración. Los pares E, F y G tienen rendimientos correlados al 0,82-0,83 y ningún ancla.
• Logaritmos, no precios. El diferencial en logaritmos es un porcentaje y su σ no cambia con el nivel del precio; el tamaño se fija en nominal, no en acciones.
• Stops por nivel y por tiempo, y bloqueo tras el stop. El riesgo de un par es que la relación deje de existir, y eso no avisa.
• El coste de la pata corta. Préstamo sobre el nominal vendido y disponibilidad de los títulos (lección 2.4.7): en un par de nominal fijo por σ, el préstamo es una fracción constante del tamaño, y pesa más cuanto más dura cada operación.

Ejercicios

1. Repite la cartera de los ocho pares con stops cambiando el umbral de entrada a 1,5 y a 2,5. ¿Qué pasa con el número de operaciones y el neto de los pares A y B?

filas = {}
for entrada in (1.5, 2.0, 2.5):
    for nombre in ("A", "B"):
        _, e, n, _ = par(nombre, entrada=entrada, stop=4.0, tiempo=60)
        filas[(entrada, nombre)] = {
            "operaciones": int((e.diff().fillna(e) != 0).sum() // 2),
            "neto": round(n.sum(), -2)}
pd.DataFrame(filas).T

       operaciones     neto
1.5 A         57.0  98600.0
    B         20.0   6500.0
2.0 A         28.0  57200.0
    B         14.0   7300.0
2.5 A         10.0  26700.0
    B         12.0   7500.0

2. Vigila el par H con un contraste de Engle-Granger móvil: cada 20 sesiones, el p-valor de las 250 sesiones anteriores. ¿Cuándo detecta la ruptura y cuánto tarda?

fechas = w.loc["2020-01-01":].index[::20]
p_movil = pd.Series({d: coint(lw["PZHY"].loc[:d].iloc[-250:],
                              lw["PZHX"].loc[:d].iloc[-250:])[1]
                     for d in fechas})
print((p_movil.loc["2021-01-01":"2022-06-30"] < 0.05)
      .map({True: "coint.", False: "no"}).to_string())

2021-01-13    coint.
2021-02-11    coint.
2021-03-12    coint.
2021-04-12    coint.
2021-05-10    coint.
2021-06-08    coint.
2021-07-07        no
2021-08-04        no
2021-09-01        no
2021-09-30        no
2021-10-28        no
2021-11-26        no
2021-12-27        no
2022-01-25        no
2022-02-23        no
2022-03-23        no
2022-04-21        no
2022-05-19        no
2022-06-17        no

3. Un proveedor te ofrece una lista de «200 pares cointegrados del S&P 500» seleccionados con datos de los últimos diez años. ¿Qué preguntas harías antes de operarlos?