Pairs trading and statistical arbitrage: distance, cointegration and Kalman.

A spread that comes back

Pairs trading is the simplest form of statistical arbitrage. You look for two assets whose spread (the difference between the logarithm of one's price and a fixed amount of the logarithm of the other's) oscillates around a stable level. When it widens a lot, you buy the one that has fallen behind and sell the one that has moved ahead; when it comes back, you close. The position is almost market neutral: if everything rises or falls at once, the two legs offset each other.

If $s_t$ is stationary, $y$ and $x$ are cointegrated (lesson 2.1.5): each one is a random walk, but a linear combination of the two is not. That is the condition that makes the pair work. Having highly correlated returns is not enough: two random walks with returns correlated at 0.85 can drift apart forever.

The logarithms are not a detail. In logarithms the spread is a percentage: a σ of 0.02 is a typical 2% mismatch between the two legs, whether they cost 20 or 200 USD. In prices ($y_t - \alpha - \beta x_t$), the σ of the spread in dollars grows with the price level, and a z computed with the formation σ stops meaning the same thing when prices double.

Gatev, Goetzmann and Rouwenhorst studied, on real US data from 1962 to 1997, a rule that paired stocks by the minimum distance between their normalized prices, and found annualized excess returns of up to 12% on self-financing portfolios of the best pairs, which exceeded a conservative estimate of costs over most of the period11Note 1 · E. Gatev, W. N. Goetzmann and 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. NBER working paper version w7032.. This lesson applies that method and the cointegration one to data where the answer is written down.

flowchart TB
  U[Universe] --> F[Formation 2014-17]
  F --> S[Pair selection]
  S --> P[Fixed parameters]
  P --> N[Trading 2018-25]
  N --> R[Result per pair]

Sixteen assets, a hundred and twenty pairs

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)                         # the relationships are in logs
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], "assets,", len(form), "formation sessions,",
      len(w.loc["2018":]), "trading sessions")
print("formulation:", verdad["formulacion"].unique().tolist())
verdad[["activo_y", "activo_x", "tipo", "beta",
        "vida_media_sesiones", "sigma_spread"]]

16 assets, 1007 formation sessions, 2011 trading sessions
formulation: ['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

Sixteen assets, 1,007 formation sessions (2014-2017) and 2,011 trading sessions (2018-2025). The truth table is only looked at to grade the work: the selection and the parameters will come from the formation period as if it did not exist. All relationships are in logarithms (formulacion = log): ln y = α + β·ln x + u, with u stationary with σ between 0.020 and 0.045 in the true pairs and a random walk in E, F and G. That is why the tests, the regressions and the spread are computed on lw, the logarithms of the prices.

Selection by distance

The distance method normalizes each asset's prices to 1 on the first day of the formation period and, for each pair, sums the squares of the difference between their paths. The pairs with the smallest sum are the ones that have been most alike. There is no model and no test: it is a measure of similarity.

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 "cross"

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                         cross
PZDY PZEX       4.89                         cross
PZEX PZEY       6.75  E (correlado_no_cointegrado)
PZFX PZHX       6.93                         cross
     PZHY       7.57                         cross
PZBY PZHY       9.01                         cross

The four closest pairs are true pairs (D, B, H and C). The fifth and sixth are cross pairs (PZDX with PZEX and PZDY with PZEX: assets from different pairs that share the market factor and moved together over these four years), and the seventh is E, correlated and not cointegrated. A is missing, the pair with the fastest reversion in the universe: its beta is 1.5, so y moves 50% more than x and their normalized paths drift apart even though the spread in logarithms comes back every few days. Distance assumes, without saying so, a one-to-one hedge. It measures similarity: it does not distinguish “they look alike” from “they are forced to come back”, and it does not see a pair whose hedge is not one to one.

Selection by cointegration

The Engle-Granger test estimates $\ln y = \alpha + \beta \ln x + s$ by least squares and tests whether the residual $s$ has a unit root with an ADF test with its own critical values (lesson 2.1.5). statsmodels.tsa.stattools.coint does it in one call. It is applied to the logarithms of the universe's 120 pairs in the formation period:

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(), "of", len(eg),
      "  expected by chance:", 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 of 120   expected by chance: 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              cross
101  PZEY  PZGX  0.0011              cross
49   PZBY  PZFY  0.0034              cross
53   PZBY  PZHY  0.0067              cross

Eighteen pairs pass at 5%, three times the 6 expected by chance among 120 tests. The first four are true (A, B, D and H, with p of 0.0003 or less), and the Bonferroni correction (p < 0.05/120 ≈ 0.0004) leaves exactly those four. The other fourteen are false positives, almost all of them cross pairs.

Why so many? The test is not wrong: 5% is an average frequency over many samples, and in a single one the 120 tests are not independent, because they share series. If PZFX looks cointegrated with PZHY, it will also look so with PZHX, which tracks PZHY closely; in the table, PZBY appears with PZFY and with PZHY. False positives arrive in clusters, and their number varies from one sample to another much more than that of 120 independent tests. Counting them in the other two four-year windows shows it:

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} with p < 0.05, {falsos:2d} false")

2014-2017: 18 with p < 0.05, 14 false
2018-2021:  8 with p < 0.05,  4 false
2022-2025:  5 with p < 0.05,  2 false

Fourteen false in 2014-2017, four in 2018-2021 and two in 2022-2025: an average of about seven per window, close to the ≈ 5.8 expected among the 115 pairs that are not true, but with one window giving more than twice that. The “expected by chance” figure is an average, not a ceiling. Pair C, which is cointegrated, does not appear among the significant ones in the formation period: with a half-life of 40 sessions, four years are not enough for the test to distinguish its spread from a random walk. A test with little power also fails in the other direction.

The eight designated pairs, one by one, with the correlation of their returns, the test in the formation period and over the whole sample, and the estimated σ and half-life of the spread. The half-life comes from an AR(1) regression of the spread: $\Delta s_t = \lambda s_{t-1} + e_t$, and $h = -\ln 2 / \ln(1 + \lambda)$ (lesson 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

What the table says, pair by pair:

• Correlation separates nothing. All the pairs, true and false, have returns correlated between 0.68 and 0.90, and two of the true ones (H, 0.68, and B, 0.72) are among the least correlated.
• A, B, D and H pass the test in the formation period with p ≈ 0. A and H have almost exact betas (1.503 against 1.5; 0.867 against 0.9) and close half-lives (4.4 and 12.3 sessions, against 5 and 10). B is estimated worse: 0.981 against 0.8, and a half-life of 9 sessions against 15. The estimated σ of the spread (sd_est) stays close to the true one in all of them (0.019, 0.025 and 0.027, against 0.020, 0.030 and 0.025).
• C does not pass in the formation period (p = 0.19) but does over the whole sample (p ≈ 0): the estimated half-life, 44 sessions, is good, but with few cycles the test has no power.
• D, with a time-varying beta, passes in the formation period with a half-life of 9.7 sessions, almost the true one: in those four years its beta moved little and a fixed beta of 1.204 works. Over the whole sample the test is already in doubt (p = 0.035): the beta has kept moving and the spread with a fixed beta reverts much more slowly than the true one (the README measures a half-life of 38 sessions with the whole sample).
• E, F and G are not cointegrated in either period, as it should be (p from 0.21 to 0.51). Their “half-lives” of 59 to 70 sessions are what an AR(1) regression gives on a random walk in a finite sample: a number that looks reasonable and means nothing. G's beta, 2.117 against 1.2, does not either.
• H passes in the formation period (p ≈ 0) and not over the whole sample (p = 0.22): the 2021 break is in the part that could not be seen when it was chosen.

Entry and exit rules

With the formation parameters fixed, the classic rule opens when the spread moves more than two standard deviations away and closes when it returns to its mean. Two more exits protect against the relationship having stopped existing: a level stop (|z| > 4) and a time stop (more than 60 sessions open, about four times the half-life of a typical pair). After a stop, the rule does not re-enter until the spread returns to |z| < 0.5: without that condition, it would reopen the next day in the same place.

stateDiagram-v2
  [*] --> FLAT
  FLAT --> LONG: z below -2
  FLAT --> SHORT: z above 2
  LONG --> FLAT: z crosses 0
  SHORT --> FLAT: z crosses 0
  LONG --> LOCKOUT: stop
  SHORT --> LOCKOUT: stop
  LOCKOUT --> FLAT: abs z below 0.5

“Long the spread” means buying $V$ USD of $y$ and selling $\beta V$ USD of $x$: with the spread in logarithms, that portfolio earns $V \cdot \Delta s$ (for small moves), because $\Delta \ln y$ and $\Delta \ln x$ are returns. The size is set so that one standard deviation of the spread is worth 1,000 USD: $V = 1{,}000 / \sigma_s$, about 50,000 USD of notional in $y$ with a σ of 0.02. The shares are computed with the prices at the close at which the position is opened and are not touched until it closes. Costs are per share and side (0.02 USD: a commission of 0.005, a half spread of 0.01 and slippage of 0.005, for stocks with a one-cent tick) plus 0.5% a year of borrowing on the notional sold short. As in the whole level, the signal from the close of t is executed at the close of 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:                       # crossed the mean
                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)               # at the next close
    Y, X = w[y].loc[z.index], w[x].loc[z.index]
    usd = pos * 1000 / sd                        # notional of y
    nuevo = pos.diff().fillna(pos) != 0      # shares fixed at opening
    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: "long", -1.0: "short", 0.0: "close"})
      .head(8).to_string())
print(f"net of pair B, 2018-2025: {neto_b.sum():,.0f} USD")

fecha
2018-05-02     long
2018-09-18    close
2018-10-09     long
2019-03-07    close
2019-03-28     long
2019-04-16    close
2019-05-03    short
2019-07-30    close
net of pair B, 2018-2025: 26,962 USD

Pair B has few, long trades: the first lasts from May to September 2018 and the second, from October 2018 to March 2019, something odd for a spread with a half-life of 15 sessions. The cause is the beta error: 0.981 instead of 0.8. With the estimated beta, the spread is $u_t - 0.181 \ln x_t$ plus a constant, and when $x$ moves, that term drifts and the spread stays for months on one side of the formation mean: a 25% rise in x moves $\ln x$ by 0.22 and the spread by 0.04, more than one and a half σ. The pair earns 26,962 USD net in eight years, without stops. In real data the problem is the same: a small error in the hedge ratio turns into a trend in the spread when prices move a lot.

Results per pair: true, false and 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

The figures are dollars per 1,000 USD of spread σ, from 2018 to 2025, net of costs:

• The true pairs win: A, 57,200 with 28 trades; B, 27,000; C, 13,300 with only 6 (its 40-session half-life gives few cycles). D, the one with the time-varying beta, also wins 32,700 with the fixed formation beta: it is a true pair, and that beta serves it for years, albeit with a drawdown of 15,200.
• The false ones lose without stops: E, −11,700; F, −4,400; G, −11,400, with drawdowns of up to 22,500. Each one makes a single trade in eight years: it opens and stays open while the spread drifts away. It is exactly the behavior of a random walk traded as if it reverted.
• Stops trim the tail, they do not improve the mean. With |z| > 4 and 60 sessions as the limit, the three false ones end between −1,500 and +2,900, with drawdowns of 1,600 to 3,000: they save 27,500 USD among the three. In the true ones, stops cost almost the same, 34,500 across B, C and D (A does not change): they cut trades that would have come back. B is the extreme case, from 27,000 to 7,300, because with its poorly estimated beta the trades last months and the time stop closes them before they come back.
• H goes from 3,800 to 6,500: the next section explains it.

The figure brings together both results. Added up, the stops barely change the total for this sample; what they change is the worst drawdown of a pair, which falls from 22,500 to 5,900. A stop on pairs is insurance with a premium: it is paid on the pairs that work and protects from the ones that do not, which in real data do not come labeled.

The broken pair

Pair H meets everything one can ask of a pair in the formation period: cointegrated with p ≈ 0, short half-life and a stable hedge ratio. During trading it works for three years. On June 15, 2021 its spread jumps five standard deviations and from that day on it is a random walk: the relationship no longer exists. It is what happens in real data when one of the two companies changes business, is acquired or suffers a problem of its own.

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("maximum z after the break:",
      round(z_h.loc["2021-06-15":].max(), 1))
print("position without stops at the end:", e_h0.iloc[-1],
      "  with 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
maximum z after the break: 18.4
position without stops at the end: -1.0   with 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

Without stops, H earns 12,500 USD from 2018 to 2020 and loses 15,900 in 2021. The rule was short before the break (the spread was hovering around +2), and on June 15 the z goes from 1.96 to 6.49 in a single session. From there it keeps drifting away, up to 18.4 standard deviations, with the position open until the end of 2025. What the version without stops earns or loses from 2022 to 2025 (+6,000, +2,100, +2,300 and −3,200) is noise: a fixed position on a random walk.

With stops, the level stop triggers at the close of June 15 (the con_stops column goes to 0 that day) and the position is closed in the next session. It did not avoid the loss from the jump, which went past the stop at 4 in a single session: 2021 closes at −3,200. What it avoided is everything that came after: the lockout prevents re-entering, because the spread does not return to |z| < 0.5, and from 2022 to 2025 there is nothing. The figure shows the two paths: the spread before and after the break, and the cumulative result with and without stops.

Time-varying hedge ratio: the Kalman filter

Pair D has a beta that moves like a reflected random walk between 0.6 and 1.6, with a daily step of 0.006, and the relationship pivots around the first price of x, $x_0$ (data README):

In the formation period, a fixed beta passes the test, but over the whole sample it is already in doubt (table of the eight pairs). A Kalman filter treats the beta as a hidden state that changes a little every day and estimates it recursively, with data up to the previous day (lesson 2.1.6). There are two ways to write the observation, and they are not equivalent:

• Centered: $\ln y_t = c + \beta_t (\ln x_t - \ln x_0) + e_t$. Here the constant $c = \alpha + \ln x_0$ is truly fixed, because the generator pivots the relationship exactly at $x_0$.
• Fixed alpha, the usual form: $\ln y_t = \alpha + \beta_t \ln x_t + e_t$ with $\alpha$ fixed. When the beta changes, the true constant, $\alpha + (1 - \beta_t)\ln x_0$, also changes, and the filter that assumes it fixed is misspecified; lesson 2.1.6 measures it in detail.

In both, the beta follows a random walk, $\beta_t = \beta_{t-1} + \eta_t$ with $\eta_t \sim N(0, q)$, and the variance of the observation noise is that of the formation spread ($0.024^2$). The step variance $q$ is the key decision; the generator's is $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])                 # the constant does not move
    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]                       # before seeing 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 = {"centered": np.column_stack([uno, lx - lx[0]]),
           "fixed alpha": 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     centered  1.0e-07       0.047    0.955    14.965
1     centered  1.0e-06       0.043    0.930     9.604
2     centered  1.0e-05       0.047    0.894     6.186
3     centered  3.6e-05       0.055    0.862     4.671
4  fixed alpha  1.0e-07       0.103    0.891     6.027
5  fixed alpha  1.0e-06       0.109    0.788     2.906
6  fixed alpha  1.0e-05       0.092    0.536     1.110
7  fixed alpha  3.6e-05       0.086    0.328     0.621

The table shows two things. The first is the specification: with the centered model, the mean beta error stays between 0.043 and 0.055 on values between 1.0 and 1.5; with the fixed constant, between 0.086 and 0.109, twice as much, because the filter compensates for the wrong constant by moving the beta. The second is the trade-off between following the beta and preserving the spread. In the centered model, as q rises the residual's autocorrelation falls little (from 0.955 to 0.862) and its half-life, from 15 to 4.7 sessions, brackets the true one of 10. In the fixed-alpha model it falls from 0.891 to 0.328: with the generator's q the residual has a half-life of 0.6 sessions. The filter attributes every deviation of the spread to a change of beta and eats the reversion you wanted to trade.

One detail: in the centered model the best beta error (0.043) is not at the generator's q, but at one a hundred times smaller. The filter assumes that the observation noise is independent from one day to the next, and the true spread is a process with memory (half-life of 10 sessions): a large q confuses those runs with changes of beta. The figure compares the true beta with that of the centered filter with $q = 10^{-6}$ and with the fixed formation beta.

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])  # the truth
variantes = {"fixed beta (OLS)": (lyD - a_f - b_f * lxD,
                                  pd.Series(b_f, index=w.index)),
             "true beta": (u_real, beta_real)}
for nombre, q in (("centered", 1e-6), ("centered", 0.006 ** 2),
                  ("fixed alpha", 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
fixed beta (OLS)              11.0   18400.0
true beta                     23.0   56100.0
centered q=1.0e-06            22.0   39400.0
centered q=3.6e-05            36.0   29700.0
fixed alpha q=3.6e-05         78.0   12800.0

With the true beta, the spread is the generator's Ornstein-Uhlenbeck process and the rule earns 56,100 USD with 23 trades: it is the ceiling, unattainable because it uses the truth. With the fixed formation beta, 18,400 in 11 trades: it works as long as the beta does not drift far from that of 2014-2017. The centered Kalman with $q = 10^{-6}$ earns 39,400 with 22 trades, and with the generator's q, 29,700 with 36: it follows the beta with more noise and trades more. The fixed-alpha one with that same q earns 12,800 with 78 trades, less than the fixed beta: it trades a spread with no memory and pays costs on every round trip. Between the ceiling and the best filter there are 16,700 USD: what it costs not to know the beta, in the best case.

A portfolio of pairs chosen in the formation period

In practice nobody knows pares2_verdad.csv: pairs are chosen with the formation data and all of them are traded. The honest portfolio is that of the eighteen pairs that pass Engle-Granger at 5% in the formation period, false positives included, with 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("cross|correlado").to_numpy()
print(f"portfolio: {cartera.sum().sum():,.0f} USD;  true: "
      f"{cartera.loc[:, ~falso].sum().sum():,.0f};  false: "
      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                         cross
PZEY-PZGX       400                         cross
PZBY-PZFY     -1600                         cross
PZBY-PZHY      4900                         cross
PZBY-PZHX       600                         cross
PZBX-PZHY     -8800                         cross
PZEX-PZEY      5000  E (correlado_no_cointegrado)
PZBY-PZFX       200                         cross
PZFX-PZHX      2600                         cross
PZEX-PZGX      1100                         cross
PZBX-PZFY     -2700                         cross
PZEY-PZGY     -1300                         cross
PZCX-PZGY     -2000                         cross
PZDY-PZGX     -2000                         cross
portfolio: 102,610 USD;  true: 96,047;  false: 6,564 (14)

The portfolio earns 102,610 USD, and almost all of it comes from the four true pairs (96,047; A alone contributes 62,700). The fourteen false ones add up to +6,564: with stops, a false pair is mostly noise with costs, sometimes positive (+10,300 on PZFX-PZHY) and sometimes negative (−8,800 on PZBX-PZHY). In this universe the false positives do not sink the portfolio because the stops keep them small, but they do not contribute anything either: they are fourteen of eighteen positions, with their capital, their borrowing and their tail risk, for a result that is zero plus noise. With more candidates and without stops they would be most of the portfolio and most of the losses.

A and B earn more here than in the table of the eight pairs (62,700 against 57,200 and 10,800 against 7,300), and D less (14,500 against 23,100), because the test chose them with the (x, y) order reversed: the beta of a regression of x on y is not the exact inverse of that of y on x, and choosing the order is also a researcher's decision that moves the result.

What this changes in your backtests

• Separate formation and trading. The beta, the mean, the σ and the selection come from one period and are traded in another. Selecting the pairs with the whole sample is choosing the ones already known to have worked.
• Count the pairs you tested. With a hundred or a thousand candidates, a test at 5% gives dozens of “cointegrated” pairs by chance; Bonferroni or FDR before trading (lesson 2.1.8).
• Correlation is not cointegration. Pairs E, F and G have returns correlated at 0.82-0.83 and no anchor.
• Logarithms, not prices. The spread in logarithms is a percentage and its σ does not change with the price level; the size is set in notional, not in shares.
• Level and time stops, and a lockout after the stop. The risk of a pair is that the relationship stops existing, and that gives no warning.
• The cost of the short leg. Borrowing on the notional sold short and availability of the shares (lesson 2.4.7): in a pair with a fixed notional per σ, borrowing is a constant fraction of the size, and it weighs more the longer each trade lasts.

Exercises

1. Repeat the portfolio of the eight pairs with stops, changing the entry threshold to 1.5 and to 2.5. What happens to the number of trades and to the net of pairs A and 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. Monitor pair H with a rolling Engle-Granger test: every 20 sessions, the p-value of the previous 250 sessions. When does it detect the break and how long does it take?

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. A vendor offers you a list of “200 cointegrated pairs from the S&P 500” selected with data from the last ten years. What questions would you ask before trading them?