--- title: "2.3.3 Pairs trading and statistical arbitrage" description: "Pair selection by distance and cointegration with controls, z-score rules, Kalman hedge ratio, the broken pair and why stops are needed." url: "https://pythonparatrading.com/en/python-trading-course/level-2/pairs-trading-statistical-arbitrage/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-2/pairs-trading-statistical-arbitrage.png" --- # 2.3.3 · Pairs trading and statistical arbitrage: distance, cointegration and Kalman > Two assets that move together and a spread that returns to its mean: the oldest idea in statistical arbitrage. This lesson selects pairs among 16 synthetic assets where we know which are cointegrated, which are only correlated and which one breaks, and checks which method finds them, which rules trade them and what happens when the relationship stops existing. - Course: [Level 2: Systematic strategies and rigorous backtesting](https://pythonparatrading.com/en/python-trading-course/level-2/) - Module: 2.3 · Strategy families - Duration: 75 min - Prerequisites: [1.6.1 · Returns, volatility, rolling correlations and drawdown](https://pythonparatrading.com/en/python-trading-course/level-1/returns-volatility-correlation-and-drawdown/), [1.1.5 · Costs, settlement, leverage and margin](https://pythonparatrading.com/en/python-trading-course/level-1/trading-costs-settlement-leverage-and-margin/), [1.4.4 · Universes without survivorship bias: delistings, mergers and symbol changes](https://pythonparatrading.com/en/python-trading-course/level-1/survivorship-bias-and-point-in-time-data/) ## By the end of this lesson - Select pairs with the distance method and with the Engle-Granger test, and count the false positives of each in a universe with a known answer. - Estimate the hedge ratio and the half-life of the spread in a formation period and trade it outside it with z-score rules, per-share costs, borrowing and execution at the next close. - Explain why a broken pair requires level and time stops, and measure how much they avoid. - Track a time-varying hedge ratio with a Kalman filter and understand the trade-off between following the beta and eating the spread. - Read the result of a pairs portfolio separating the true pairs from the false ones. ## 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. $$ s_t = \ln y_t - \alpha - \beta \ln x_t, \qquad z_t = \frac{s_t - \mu_s}{\sigma_s} $$ 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 period[^1]. This lesson applies that method and the cointegration one to data where the answer is written down. ```mermaid 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 ```python 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"]] ``` ```output 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. ```python 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 ``` ```output 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: ```python 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) ``` ```output 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: ```python 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") ``` ```output 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). ```python 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) ``` ```output 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. ```mermaid 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`. ```python 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") ``` ```output 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 ```python 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 ``` ```output 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 ``` *Fig. — Net result of each pair from 2018 to 2025, without stops and with level (|z| > 4) and time (60 sessions) stops.* Horizontal bars per pair, no stops on top and with stops below, in USD thousands per USD 1,000 of spread σ. The three cointegrated pairs win with and without stops: A 57.2 and 57.2; B 27.0 and 7.3; C 13.3 and 8.1. D, the varying-beta pair, 32.7 and 23.1. The three false pairs lose without stops (E −11.7, F −4.4, G −11.4) and end between −1.5 and +2.9 with stops. H, the broken pair (highlighted), goes from 3.8 to 6.5. 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. ```python 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()) ``` ```output 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 ``` *Fig. — Broken pair H: spread z-score (top) and cumulative result with and without stops (bottom), 2018-2025.* Top, the weekly z-score of the H log spread moves between −2.3 and +3 until 15 June 2021, jumps above the stop at 4 and drifts away to 18. Bottom, two lines of cumulative result in USD thousands: without stops it rises to about 15 by mid-2021, sinks with the break to −7, recovers some with the noise of the random walk and ends at +3.8; with stops (highlighted) it reaches about 12, loses the jump of the break, is closed the next session and stays flat at +6.5. 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): $$ \ln y_t = \alpha + \ln x_0 + \beta_t \, (\ln x_t - \ln x_0) + u_t $$ 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}$. ```python 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) ``` ```output 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 ``` *Fig. — Pair D hedge ratio: the generator's true beta, a centred Kalman filter's and the fixed beta of a formation-period regression, 2015-2025.* Three weekly lines from 2015 to 2025. The true beta moves between 1.02 and 1.49 and rises mostly from 2022. The centred Kalman filter with step variance 10 to the −6 (highlighted) follows it with some lag and noise; its mean error is 0.043. The fixed formation beta, dashed, is 1.20 throughout: right in 2014-2017 and too low once the true beta rises. 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. ```python 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 ``` ```output 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: ```python 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()})") ``` ```output 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? **Solution.** ```python 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 ``` ```output 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 ``` In A, lowering the threshold to 1.5 doubles the trades (57 against 28) and almost doubles the net (98,600 against 57,200); raising it to 2.5 leaves it at 10 trades and 26,700. In a truly cointegrated pair with a well-estimated beta, entering earlier collects more reversion cycles, and the cost per share is small relative to the gain of each cycle. In B the threshold barely matters (6,500, 7,300 and 7,500, with 20, 14 and 12 trades): with its poorly estimated beta the spread drifts and the result is decided by the stops, not the entry. With higher costs per trade, or on a false pair, the conclusion would be the opposite: each extra entry pays costs without an edge. Choosing the threshold with these results would be, once again, fitting to the trading sample. 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? **Solution.** ```python 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()) ``` ```output 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 ``` The rolling test shows cointegration in every window up to the one on June 8, 2021 and loses it in the next one, on July 7, three weeks after the break, and does not recover it in the following year. As monitoring it works, but it is not an immediate alarm: the 250-session window needs the broken stretch to weigh enough. And in the other direction, a rolling test on a false pair will show cointegration in some window by chance: at 20 sessions per step, thirteen tests are run a year. 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? **Solution.** How many pairs were tested to arrive at 200 (with 500 stocks there are almost 125,000 possible pairs, and at 5% about 6,000 would come out “cointegrated” by chance) and what multiple-testing correction was applied. Whether the selection was made with the same data with which the result is measured or with an earlier formation period. With what universe (the index's current members have survivorship bias). Whether there is an economic reason for each pair (same business, same commodity, two share classes of the same company) or just a statistical coincidence. How long each relationship lasts out of sample and what stops and per-pair limits were used. And the real costs: bid-ask spread, commissions and, above all, borrowing and availability of the short leg. [^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. ## Self-assessment 1. In the formation period, the Engle-Granger test at 5% on the 120 possible pairs of 16 assets gives 18 significant pairs. Only four are true pairs. What explains the other fourteen? - A. That the test is coded wrong - B. That those pairs are cointegrated but the generator does not say so - C. That with 120 tests at 5% about 6 false positives are expected even with no relationship at all, and since the tests share series they arrive in clusters: in a given sample they can be twice as many or more - D. That the correlation of their returns is high *Answer C: That with 120 tests at 5% about 6 false positives are expected even with no relationship at all, and since the tests share series they arrive in clusters: in a given sample they can be twice as many or more.* A 5% level means that, without a relationship, one test in twenty comes out significant on average. If PZFX looks cointegrated with PZHY, it will also look so with PZHX, which tracks PZHY closely: false positives are not independent. In the 2018-2021 and 2022-2025 windows there are 4 and 2. With the Bonferroni correction (0.05/120) only the four true ones survive (lesson 2.1.8). 2. The distance method leaves out of its top ten pair A, which is cointegrated and has the shortest half-life in the universe. Why? - A. Because its beta is 1.5: distance compares normalized prices one to one, and with a beta other than 1 the two paths drift apart even though the spread in logarithms reverts - B. Because it is not cointegrated in the formation period - C. Because its half-life is too short to measure - D. Because of a rounding error *Answer A: Because its beta is 1.5: distance compares normalized prices one to one, and with a beta other than 1 the two paths drift apart even though the spread in logarithms reverts.* Distance is a measure of similarity with an implicit one-to-one hedge, not of reversion. In exchange, it puts two cross pairs and E, correlated and without cointegration, in its top seven. Gatev, Goetzmann and Rouwenhorst use it with short trading periods (six months), partly because the relationship of a distance pair need not last. 3. Pair H breaks on June 15, 2021: the spread jumps from under 2 to 6.5 σ in one session and from then on is a random walk. The rule was short. What did the level stop (|z| > 4) do? - A. It avoided the loss from the jump - B. Nothing: with stops it loses the same - C. It closed the position before the break - D. It did not avoid the loss from the jump, which went past the stop in a single session, but it closed the position in the next session, and the lockout prevented re-entering while the spread drifted away to 18.4 σ *Answer D: It did not avoid the loss from the jump, which went past the stop in a single session, but it closed the position in the next session, and the lockout prevented re-entering while the spread drifted away to 18.4 σ.* With stops, 2021 closes at −3,200 USD and there are no more trades afterwards; without stops, −15,900 in 2021 and a short position open until the end of 2025 on a relationship that no longer exists. A stop does not protect from the jump, but from what comes after. 4. In pair D, a Kalman filter with a fixed constant and the generator's step variance (q = 3.6e-5) leaves a residual spread with almost no memory. Why? - A. Because the filter is coded wrong - B. Because in that model the true constant moves with the beta: the filter compensates for the error by moving the beta and explains every deviation of the spread as a change of beta; it eats the reversion you wanted to trade - C. Because the spread of pair D does not revert - D. Because Kalman only works for stationary series *Answer B: Because in that model the true constant moves with the beta: the filter compensates for the error by moving the beta and explains every deviation of the spread as a change of beta; it eats the reversion you wanted to trade.* With a fixed constant, the residual's autocorrelation drops to 0.33 (half-life of 0.6 sessions) and the rule earns 12,800 USD in 78 trades. With the regressor centered on the first price of x, the constant is truly fixed, the residual keeps its memory (0.86) and, with q = 1e-6, the rule earns 39,400 in 22 trades. 5. Why is the spread of these pairs defined in logarithms and not in prices? - A. Because logarithms always give more profit - B. Because statsmodels does not accept prices - C. Because in logarithms the spread is a percentage: its size does not change with the price level, whereas in prices its σ in dollars grows when prices rise and a z with the formation σ stops meaning the same thing - D. Because that way there are no costs *Answer C: Because in logarithms the spread is a percentage: its size does not change with the price level, whereas in prices its σ in dollars grows when prices rise and a z with the formation σ stops meaning the same thing.* The generator writes every relationship in logarithms (`ln y = α + β·ln x + u`). In prices, with volatility proportional to the price, cointegration tests reject more than they should (the README measures it with Johansen), and the beta in dollars changes with the level of both prices. ## Key points - A pair is a bet on a stationary spread in logarithms; the correlation of returns does not guarantee it: the three correlated pairs without cointegration have correlations of 0.82 to 0.83, more than two of the true ones (B and H). - With 120 possible pairs, the Engle-Granger test at 5% gives 14 false positives in the formation period, more than twice what is expected: they arrive in clusters. Bonferroni leaves only the four true ones. Distance assumes a beta of 1 and leaves out pair A. - Estimate in formation and trade outside it: hedge ratio, mean, σ and half-life from 2014-2017, trading from 2018 to 2025, with execution at the next close. - The true pairs win and the false ones lose without stops (−11,700, −4,400 and −11,400 USD per 1,000 USD of σ). Stops cost the true pairs almost what they save on the false ones: they do not improve the mean, they trim the tail. - The broken pair H teaches the tail risk of the family: a relationship that works for years can disappear in a day, and the stop does not avoid the jump, but what comes after. - Kalman tracks a time-varying beta if the model is well specified; with a fixed constant and a very agile beta it eats the spread. The step variance is a trading decision, not just a fitting one. ## Resources ### Research - [Stochastic Processes - Random Walks](https://pythonparatrading.com/en/stochastic-processes/) · Research-0005 · 11 min ### Documentation - [Gatev, Goetzmann and Rouwenhorst: “Pairs Trading: Performance of a Relative Value Arbitrage Rule” (NBER w7032; Review of Financial Studies 19(3), 2006)](https://www.nber.org/papers/w7032) · nber.org - [Engle and Granger (1987): “Co-Integration and Error Correction: Representation, Estimation, and Testing”, Econometrica 55(2)](https://doi.org/10.2307/1913236) · doi.org - [statsmodels: Engle-Granger cointegration test (coint)](https://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.coint.html) · statsmodels.org ### Books - [Ernest P. Chan: Algorithmic Trading: Winning Strategies and Their Rationale (Wiley, 2013), chapter 3 (Kalman filter as dynamic regression)](https://www.wiley.com/en-us/Algorithmic+Trading%3A+Winning+Strategies+and+Their+Rationale-p-9781118460146) · wiley.com ### Data - [Level 2 data: eight synthetic pairs, 2014-2025 (Parquet)](https://pythonparatrading.com/curso/datos/nivel-2/pares2.parquet) · Download · pares2.parquet - [Level 2 data: the truth for each pair (CSV)](https://pythonparatrading.com/curso/datos/nivel-2/pares2_verdad.csv) · Download · pares2_verdad.csv Previous lesson: [2.3.2 · Mean reversion: RSI and Bollinger, intraday reversal and gaps](https://pythonparatrading.com/en/python-trading-course/level-2/mean-reversion-and-gaps/) Next lesson: [2.3.4 · Carry in currencies, futures and bonds; seasonality and calendar effects](https://pythonparatrading.com/en/python-trading-course/level-2/carry-seasonality-calendar-effects/)