Filtro de Kalman para hedge ratios dinámicos.

Un hedge ratio que no se está quieto

En la lección 2.1.5 el par D aprobaba el contraste de Engle-Granger con un β fijo, pero el residuo de ese β revertía cuatro veces más despacio que el spread verdadero: su relación es $\ln y_t = a_t + \beta_t \ln x_t + u_t$ con un $\beta_t$ que cambia cada día. Una regresión con un β fijo para doce años describe una relación media que nunca existió. Es una situación real frecuente: la sensibilidad de una acción a su sector, la de un ETF a su índice o la de un bono a otro de distinto plazo cambian con la estructura de la empresa, con los tipos o con la composición del índice.

El hedge ratio importa por dos motivos. Con la relación en logaritmos, β es una elasticidad: para cubrir una posición de 1 USD en y hay que vender unos $\beta$ USD de x, y un β equivocado deja riesgo sin cubrir. Para operar un par (lección 2.3.3), el spread $\ln y - \beta \ln x$ es la señal: un β equivocado mete en el spread movimientos que no son reversión.

Los datos de esta lección son el par D de pares2.parquet, sintético, con la serie del β verdadero de cada sesión en pares2_beta_real.csv: un paseo aleatorio reflejado con paso diario de desviación típica 0,006, y un spread de Ornstein-Uhlenbeck con vida media 10 y desviación típica 0,025 (un 2,5 %). La relación gira alrededor del primer precio de x, $x_0$: $\ln y_t = \alpha + \ln x_0 + \beta_t \ln(x_t / x_0) + u_t$, así que en la forma habitual la constante $a_t = \alpha + (1 - \beta_t) \ln x_0$ también se mueve con β. Todo se trabaja en logaritmos. Saber la verdad permite medir el error de cada método, algo que con datos reales nunca se puede hacer.

import numpy as np
import pandas as pd
import statsmodels.api as sm

pares = pd.read_parquet("pares2.parquet")
precios = pares.pivot(index="fecha", columns="activo",
                      values="precio")
y, x = np.log(precios["PZDY"]), np.log(precios["PZDX"])
beta = pd.read_csv("pares2_beta_real.csv", index_col="fecha",
                   parse_dates=True)["beta_real"]
print(f"beta verdadero: de {beta.min():.3f} a {beta.max():.3f}, "
      f"paso diario {beta.diff().std():.4f}")

fijo = sm.OLS(y, sm.add_constant(x)).fit().params
error = np.sqrt(((fijo.iloc[1] - beta) ** 2).mean())
print(f"MCO de toda la muestra: alfa {fijo.iloc[0]:.2f}, "
      f"beta {fijo.iloc[1]:.3f}, error cuadrático medio {error:.3f}")

beta verdadero: de 0.976 a 1.489, paso diario 0.0061
MCO de toda la muestra: alfa -1.03, beta 1.337, error cuadrático medio 0.166

El β verdadero recorre de 0,98 a 1,49. La regresión de toda la muestra da un β de 1,337 y una constante de −1,03: la regresión ajusta la relación media de los precios, no estima la de ningún día concreto. Su error cuadrático medio frente al β verdadero es de 0,166. Todo lo que sigue se mide con esa vara: la raíz del error cuadrático medio (RMSE) del β estimado frente al verdadero.

La regresión móvil y su ventana

La respuesta habitual es una regresión en una ventana móvil: con los últimos $n$ días se estima un β, y al día siguiente se desplaza la ventana. Hay dos formas de hacerla, en precios o en variaciones de precio, y se comportan de forma muy distinta:

def rmse(est):
    """Error del beta disponible el día anterior frente al verdadero."""
    e = (est.shift(1) - beta).iloc[500:]
    return np.sqrt((e ** 2).mean())

dy, dx = y.diff(), x.diff()
filas = {}
for n in [10, 20, 40, 60, 120, 250, 500]:
    en_precios = y.rolling(n).cov(x) / x.rolling(n).var()
    en_variac = ((dy * dx).rolling(n).sum()
                 / (dx ** 2).rolling(n).sum())
    filas[n] = {"precios": rmse(en_precios),
                "variaciones": rmse(en_variac)}
pd.DataFrame(filas).T.rename_axis("ventana").round(4)

         precios  variaciones
ventana
10        0.5309       0.2841
20        0.4458       0.1845
40        0.3591       0.1336
60        0.3189       0.1060
120       0.2379       0.0753
250       0.1500       0.0704
500       0.1495       0.0900

El error se mide a partir de la sesión 500, para que todas las ventanas tengan datos, y siempre con el β del día anterior (shift(1)): la cobertura de hoy se decide con lo que se sabía ayer.

En precios (sus logaritmos, con constante), la regresión falla: el error va de 0,53 con 10 sesiones a 0,15 con 500, nunca mejor que un β fijo. El motivo es la constante. En una ventana de 120 sesiones el logaritmo del precio de PZDX se mueve con una desviación típica de un 5 % alrededor de un nivel de 4, y en ese tramo corto muchas combinaciones de α y β ajustan casi igual de bien: un β más alto con un α más bajo. Las dos estimaciones se compensan y cada una por separado es ruido:

b120 = y.rolling(120).cov(x) / x.rolling(120).var()
a120 = y.rolling(120).mean() - b120 * x.rolling(120).mean()
print(f"correlación entre alfa y beta móviles: {a120.corr(b120):.2f}")
print(f"beta móvil de {b120.min():.2f} a {b120.max():.2f}; "
      f"alfa de {a120.min():.0f} a {a120.max():.0f}")

correlación entre alfa y beta móviles: -1.00
beta móvil de 0.45 a 1.96; alfa de -4 a 3

En variaciones (la regresión de $\Delta y$ sobre $\Delta x$, sin constante) la constante desaparece y el error baja mucho, con la forma de U típica de toda estimación móvil. Con ventanas cortas el β sigue bien los cambios pero se estima con pocos datos: varianza. Con ventanas largas se estima con muchos datos pero mezcla βs de hace meses: sesgo. El mínimo, 0,070 con 250 sesiones (0,075 con 120), es el mejor equilibrio para este par; otro par con otro ritmo de cambio tendría otro.

El defecto de fondo de la ventana es que trata igual el dato de hace 249 días que el de ayer, y al de hace 251 no le da ningún peso. El filtro de Kalman sustituye esa ventana rectangular por una actualización con pesos que salen de un modelo.

El modelo de espacio de estados

Un modelo de espacio de estados separa lo que no se ve (el estado, aquí $\beta_t$) de lo que se observa (las variaciones de los logaritmos de los precios, es decir, los rendimientos). Para el modelo en variaciones:

La primera ecuación dice que el β de hoy es el de ayer más un cambio aleatorio; la segunda, que la variación de y es β veces la de x más un ruido. Las dos varianzas son las únicas decisiones: $q$, cuánto puede moverse β de un día a otro; $r$, cuánto ruido tiene cada observación. En estos datos la verdad es $q = 0{,}006^2 = 3{,}6 \times 10^{-5}$; el ruido de observación no es exactamente blanco (incluye la variación del spread y el efecto del cambio de β sobre el precio), así que $r$ es una aproximación: el modelo es útil aunque no sea exacto. Como orden de magnitud, $r$ es la varianza diaria de lo que la cobertura no explica: con un spread de desviación 0,025 y vida media 10, del orden de $10^{-4}$, que es el valor que se toma abajo.

El filtro de Kalman calcula, para cada día, la mejor estimación de $\beta_t$ con los datos hasta ese día y su varianza $P_t$. Lo hace en dos fases:

flowchart TB
  A["Ayer: beta, P"] --> B[Predicción]
  B --> C{"Llega el precio"}
  C --> D[Innovación e]
  D --> E[Ganancia K]
  E --> F["Hoy: beta, P"]
  F --> B

• Predicción. Sin datos nuevos, la mejor estimación de $\beta_t$ es la de ayer, pero con más incertidumbre: $\beta_{t|t-1} = \beta_{t-1|t-1}$ y $P_{t|t-1} = P_{t-1|t-1} + q$.
• Actualización. Llega $\Delta y_t$. La innovación es la sorpresa, $e_t = \Delta y_t - \Delta x_t\,\beta_{t|t-1}$, con varianza $S_t = \Delta x_t^2 P_{t|t-1} + r$. La ganancia $K_t = P_{t|t-1}\,\Delta x_t / S_t$ dice qué parte de la sorpresa se atribuye a un cambio de β, y $\beta_{t|t} = \beta_{t|t-1} + K_t e_t$, $P_{t|t} = (1 - K_t \Delta x_t) P_{t|t-1}$.

La ganancia lo resume todo: si el estado es muy incierto ($P$ grande) o la observación muy informativa ($r$ pequeño, $\Delta x_t$ grande), el β se mueve mucho con la sorpresa; si no, poco. Con una estimación inicial muy incierta, lo que decide el comportamiento es el cociente $q/r$.

El filtro desde cero

El filtro cabe en un bucle. Esta versión admite un estado de cualquier dimensión, porque más adelante se usará con dos ($\alpha$ y $\beta$):

def kalman(z, H, q, r, p0=1e4):
    """Filtro de Kalman con estado de paseo aleatorio.

    z: observaciones (n,); H: fila de diseño de cada día (n, k);
    q: varianzas del estado (k,); r: varianza de la observación.
    Devuelve el estado filtrado (n, k), la innovación y su varianza.
    """
    n, k = H.shape
    est, P, Q = np.zeros(k), np.eye(k) * p0, np.diag(q)
    filtrado = np.full((n, k), np.nan)
    innov, var_innov = np.full(n, np.nan), np.full(n, np.nan)
    for t in range(n):
        P = P + Q                       # predicción
        if not np.isfinite(z[t]):
            continue
        h = H[t]
        S = h @ P @ h + r
        K = P @ h / S                   # ganancia
        e = z[t] - h @ est              # innovación
        est = est + K * e               # actualización
        P = P - np.outer(K, h) @ P
        filtrado[t], innov[t], var_innov[t] = est, e, S
    return filtrado, innov, var_innov

dz, h = dy.to_numpy(), dx.to_numpy()[:, None]
r = 1e-4
filas = {}
for q in [1e-7, 1e-6, 1e-5, 3.6e-5, 1e-4, 3e-4, 1e-3, 1e-2]:
    f, _, _ = kalman(dz, np.nan_to_num(h), [q], r)
    filas[f"{q:g}"] = rmse(pd.Series(f[:, 0], index=y.index))
pd.Series(filas, name="rmse").rename_axis("q").round(4)

q
1e-07      0.1164
1e-06      0.0984
1e-05      0.0699
3.6e-05    0.0614
0.0001     0.0641
0.0003     0.0782
0.001      0.1057
0.01       0.1912
Name: rmse, dtype: float64

El error frente a q tiene la misma forma de U que el de la ventana móvil, y por la misma razón. Con q muy pequeño ($10^{-7}$) el filtro casi no deja moverse a β y va siempre detrás de la verdad (0,116): sesgo. Con q muy grande ($10^{-2}$) cada sorpresa del precio mueve β y la estimación es ruido (0,191): varianza. El mínimo, 0,0614, cae en el q verdadero, $3{,}6 \times 10^{-5}$, y mejora la mejor ventana móvil (0,0704); con $q = 10^{-4}$ el error apenas sube (0,0641): el fondo de la U es plano. El filtro es óptimo cuando el modelo y sus varianzas son los verdaderos, y aquí, con un r razonable, lo son casi.

Sesgo y varianza, a la vista

La tabla resume el error, pero no enseña cómo se equivoca cada filtro. Un año del β estimado con un q unas treinta veces menor y otro unas treinta veces mayor que el del mínimo lo muestra:

tramo = slice("2023-01-01", "2023-12-31")
ver = pd.DataFrame({"verdad": beta[tramo]})
for q in [1e-6, 1e-3]:
    f, _, _ = kalman(dz, np.nan_to_num(h), [q], r)
    est = pd.Series(f[:, 0], index=y.index).shift(1)[tramo]
    ver[f"q={q:g}"] = est
cambio = ver.diff().abs().mean()
print("Variación media diaria del beta:", cambio.round(4).to_dict())
ver.iloc[[0, 60, 120, 180, 240]].round(3)

Variación media diaria del beta: {'verdad': 0.0047, 'q=1e-06': 0.0006, 'q=0.001': 0.0209}
            verdad  q=1e-06  q=0.001
fecha
2023-01-03   1.342    1.189    1.394
2023-03-30   1.356    1.189    1.284
2023-06-27   1.463    1.200    1.261
2023-09-21   1.435    1.223    1.473
2023-12-15   1.400    1.236    1.438

En 2023 el β verdadero pasa de 1,34 a 1,46 en junio y acaba el año en 1,40. Con $q = 10^{-6}$ el filtro va muy por detrás y por debajo (1,19 en enero, 1,24 en diciembre) y casi no se mueve de un día a otro: 0,0006 de variación media diaria, frente a 0,0047 de la verdad. Con $q = 10^{-3}$ sigue la verdad en media, pero se mueve 0,0209 al día, más de cuatro veces más que ella, y puede estar en 1,26 cuando la verdad está en 1,46. Un hedge ratio que salta cada día obliga a rehacer la cobertura cada día y paga comisiones y horquilla en cada ajuste (lección 2.4.7). El q óptimo para el error no es necesariamente el óptimo para la cuenta de resultados: un poco más de sesgo puede salir más barato.

Máxima verosimilitud con statsmodels

En datos reales no hay β verdadero para medir el error, así que q y r no se pueden elegir mirando el RMSE. Lo que sí hay es la verosimilitud: el filtro produce cada día una innovación con su varianza $S_t$, y si el modelo y sus varianzas son buenos, las innovaciones se comportan como normales independientes de esas varianzas. Maximizar esa verosimilitud estima q y r con los propios datos.

statsmodels tiene un marco de espacio de estados con filtro, suavizador y máxima verosimilitud. Para un β variable basta una clase pequeña: la «matriz de diseño» cambia cada día (es $\Delta x_t$) y las dos varianzas son los parámetros.

class BetaVariable(sm.tsa.statespace.MLEModel):
    """dy_t = beta_t * dx_t + e_t; beta_t paseo aleatorio."""

    def __init__(self, dy, dx):
        super().__init__(dy, k_states=1, initialization="diffuse")
        self["design"] = np.asarray(dx, float).reshape(1, 1, -1)
        self["transition"] = np.eye(1)
        self["selection"] = np.eye(1)

    @property
    def param_names(self):
        return ["r", "q"]

    @property
    def start_params(self):
        return [1e-4, 1e-5]

    def transform_params(self, u):      # varianzas positivas
        return u ** 2

    def untransform_params(self, c):
        return c ** 0.5

    def update(self, params, **kwargs):
        params = super().update(params, **kwargs)
        self["obs_cov", 0, 0] = params[0]
        self["state_cov", 0, 0] = params[1]

# arrays sin fechas: el índice de sesiones no tiene frecuencia fija
# y statsmodels avisaría (ValueWarning) sin que cambie nada
modelo = BetaVariable(dy.iloc[1:].to_numpy(), dx.iloc[1:].to_numpy())
ajuste = modelo.fit(disp=False)
r_mv, q_mv = ajuste.params
print(f"r = {r_mv:.2e}  q = {q_mv:.2e}  "
      f"(q verdadero {0.006 ** 2:.1e})")

r = 8.92e-05  q = 5.12e-05  (q verdadero 3.6e-05)

La máxima verosimilitud da $q = 5{,}1 \times 10^{-5}$, cerca del verdadero, y $r = 8{,}9 \times 10^{-5}$, cerca del $10^{-4}$ que se había tomado como orden de magnitud. Con datos reales sería el punto de partida razonable, sin haber mirado la verdad.

Dos comprobaciones más: que el filtro propio y el de statsmodels dan lo mismo con los mismos parámetros, y qué ocurre con el suavizado:

f_mv, _, _ = kalman(dz, np.nan_to_num(h), [q_mv], r_mv)
propio = pd.Series(f_mv[:, 0], index=y.index)
sm_filt = pd.Series(ajuste.filtered_state[0], index=y.index[1:])
sm_suav = pd.Series(ajuste.smoothed_state[0], index=y.index[1:])

dif = (propio.iloc[500:] - sm_filt.iloc[499:]).abs().max()
print(f"diferencia máxima propio - statsmodels: {dif:.1e}")
print(f"RMSE filtrado (día anterior): {rmse(sm_filt):.4f}")
e_suav = (sm_suav - beta).iloc[500:]
print(f"RMSE suavizado (mismo día):   "
      f"{np.sqrt((e_suav ** 2).mean()):.4f}")

diferencia máxima propio - statsmodels: 1.3e-08
RMSE filtrado (día anterior): 0.0614
RMSE suavizado (mismo día):   0.0382

El filtro propio y el de statsmodels coinciden hasta la octava cifra decimal (la diferencia viene del arranque: statsmodels usa una inicialización difusa exacta y el propio, una varianza inicial de 10.000), y el β filtrado con las varianzas de máxima verosimilitud tiene un error de 0,061, como el mejor de la tabla, sin haber mirado la verdad. El suavizado baja a 0,038, y es una trampa: el suavizador de Kalman (Rauch-Tung-Striebel) recorre la muestra hacia atrás desde el final, así que el β suavizado de cada día usa todos los precios posteriores. Sirve para estudiar cómo fue la relación; en un backtest es sesgo de anticipación, el mismo de las probabilidades suavizadas de un HMM (lección 2.1.7).

En precios: un α que también se mueve y un riesgo para el spread

El modelo en variaciones tira la información de los niveles. Si se supone que los logaritmos de los precios están ligados por $\ln y_t = \alpha + \beta_t \ln x_t + u_t$, el estado puede ser el vector $(\alpha, \beta_t)$ y la observación, el logaritmo del precio:

Es el planteamiento más extendido en pares. La función kalman sirve tal cual con un diseño de dos columnas y, como ruido de observación, la varianza del spread ($0{,}025^2$). Con α fijo ($q_\alpha = 0$), que es lo que se suele hacer, y varios $q_\beta$, se mide el error del β y algo que importa para operar el spread: la autocorrelación de la innovación, que es el spread que ve el filtro.

Hp = np.column_stack([np.ones(len(x)), x.to_numpy()])
filas = {}
for qb in [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]:
    f, e, S = kalman(y.to_numpy(), Hp, [0.0, qb], 0.025 ** 2)
    inn = pd.Series(e, index=y.index).iloc[500:]
    filas[f"{qb:g}"] = {
        "rmse_beta": rmse(pd.Series(f[:, 1], index=y.index)),
                 "alfa_final": f[-1, 0],
                 "ac1_innovacion": inn.autocorr()}
# a_t = alfa + (1 - beta_t) * ln x0, con alfa = 0,3 (pares2_verdad.csv)
print("a verdadero el último día:",
      round(0.3 + (1 - beta.iloc[-1]) * x.iloc[0], 3))
pd.DataFrame(filas).T.rename_axis("q_beta").round(3)

a verdadero el último día: -1.063
        rmse_beta  alfa_final  ac1_innovacion
q_beta
1e-07       0.131      -0.458           0.891
1e-06       0.136      -0.439           0.791
1e-05       0.120      -0.475           0.545
0.0001      0.112      -0.498           0.182
0.001       0.110      -0.503           0.009
0.01        0.110      -0.503          -0.018

En estos datos el modelo en precios estima el β peor que el de variaciones: el error no baja de 0,110, frente al 0,061 del filtro en variaciones. La razón es de especificación. Aquí la constante verdadera $a_t$ se mueve con β (el último día vale −1,063) y el modelo la supone fija; el filtro la estima en −0,44 a −0,50 y compensa el error con el β, igual que la regresión móvil. Un modelo en niveles solo es más preciso si su forma es la correcta.

La última columna enseña el otro riesgo, que existe aunque el modelo esté bien planteado. Con q pequeño la innovación tiene autocorrelación 0,89: conserva la memoria del spread, que es lo que revierte. Con $q_\beta = 10^{-4}$ ya baja a 0,18, y con más, a cero o algo menos: el filtro explica cada desviación del spread como un cambio de β, y el residuo deja de revertir. El filtro se come el spread. Para cubrir puede valer; para operar la reversión, la señal ha desaparecido.

Por eso en un par no se elige q solo por el error del β (que además en datos reales no se conoce): se elige mirando también si el spread resultante conserva una vida media útil (lección 2.1.5), comprobando que la forma del modelo es razonable (aquí, dejar que α se mueva, o usar el modelo en variaciones, que no lo necesita) y validando fuera de la muestra.

flowchart TB
  Q{"q del estado"} -- pequeño --> L[beta lento: sesgo]
  Q -- grande --> R[beta ruidoso]
  R --> C[Más costes]
  R --> S[Spread absorbido]
  L --> P[Cobertura tardía]

Errores típicos con el filtro de Kalman

• Usar el estado actualizado con el precio de hoy para operar hoy. El β utilizable es el de ayer. En el modelo en precios, la innovación $e_t$ ya está calculada con el estado de ayer: es la versión sin anticipación del spread.
• Usar el suavizado. Siempre es más preciso y siempre es anticipación.
• Elegir q y r por el resultado del backtest. Es sobreajuste con otro nombre. Máxima verosimilitud o validación fuera de la muestra (lección 2.5.2).
• Olvidar el arranque. Con una varianza inicial grande, las primeras decenas de sesiones son ruido. Descártalas, como aquí las 500 primeras al medir el error.
• Confiar en el modelo. El ruido de observación de estos datos no es blanco y el filtro funciona igual, pero sus varianzas $S_t$ no son exactas: un z-score $e_t / \sqrt{S_t}$ como señal (lección 2.3.3) hereda ese error.
• Precios en escalas distintas. q depende de las unidades: si x se multiplica por 100, el β se divide por 100 y su q por 10.000. Reescalar los precios cambia el q adecuado.

Ejercicios

1. Repite la tabla del error frente a q del modelo en variaciones con $r = 2{,}5 \times 10^{-5}$ y con $r = 4 \times 10^{-4}$. ¿Se desplaza el mínimo? ¿En qué proporción?

for r_alt in [2.5e-5, 4e-4]:
    errores = {}
    for q in [2.25e-6, 9e-6, 3.6e-5, 1.44e-4, 5.76e-4]:
        f, _, _ = kalman(dz, np.nan_to_num(h), [q], r_alt)
        errores[q] = rmse(pd.Series(f[:, 0], index=y.index))
    mejor = min(errores, key=errores.get)
    print(f"r = {r_alt}: mejor q = {mejor:g}, "
          f"rmse {errores[mejor]:.4f}")

r = 2.5e-05: mejor q = 9e-06, rmse 0.0614
r = 0.0004: mejor q = 0.000144, rmse 0.0614

2. Con el modelo en variaciones y las varianzas de máxima verosimilitud, calcula el z-score de la innovación, $e_t / \sqrt{S_t}$. Si el modelo fuera correcto, ¿qué media, desviación típica y autocorrelación debería tener? ¿Las tiene?

_, e_mv, S_mv = kalman(dz, np.nan_to_num(h), [q_mv], r_mv)
z = pd.Series(e_mv / np.sqrt(S_mv), index=y.index).iloc[500:]
print(f"media {z.mean():.3f}  desviación {z.std():.3f}  "
      f"autocorrelación {z.autocorr():.3f}")
print(f"curtosis en exceso {z.kurt():.2f}  autocorrelación "
      f"de z al cuadrado {(z ** 2).autocorr():.3f}")

media 0.004  desviación 1.010  autocorrelación -0.025
curtosis en exceso 0.10  autocorrelación de z al cuadrado 0.003

3. Un compañero propone elegir el $\delta$ del filtro de su par probando veinte valores y quedándose con el que da mejor Sharpe en el backtest de 2015 a 2024. ¿Qué problemas tiene y qué harías en su lugar?