--- title: "2.1.5 Cointegration: Engle-Granger, Johansen and half-life" description: "Correlation is not cointegration: Engle-Granger and Johansen on synthetic pairs with a known answer, Ornstein-Uhlenbeck half-life and a break." url: "https://pythonparatrading.com/en/python-trading-course/level-2/cointegration-engle-granger-johansen-half-life/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-2/cointegration-engle-granger-johansen-half-life.png" --- # 2.1.5 · Cointegration: Engle-Granger, Johansen and the half-life of mean reversion > Two prices can move together day by day and drift apart without limit, or move with little correlation and always come back to the same distance. This lesson separates the two with eight synthetic pairs whose answer is known: which ones are cointegrated, with what half-life and when the relationship breaks. - Course: [Level 2: Systematic strategies and rigorous backtesting](https://pythonparatrading.com/en/python-trading-course/level-2/) - Module: 2.1 · Statistics and econometrics for trading - Duration: 70 min - Prerequisites: [0.7.2 · Probability and distributions, correlation versus causation and simple regression](https://pythonparatrading.com/en/python-trading-course/level-0/probability-distributions-and-regression/), [1.6.1 · Returns, volatility, rolling correlations and drawdown](https://pythonparatrading.com/en/python-trading-course/level-1/returns-volatility-correlation-and-drawdown/) ## By the end of this lesson - Tell correlation of returns from cointegration of prices, and explain why the first does not imply the second. - Apply the two-step Engle-Granger test and the Johansen test with `statsmodels`, and read their statistics and critical values. - Estimate the half-life of a spread's mean reversion as an Ornstein-Uhlenbeck process and measure the bias of the estimator. - Recognize the usual mistakes: levels versus logarithms, screens of many pairs, static βs that pass and short stretches after a break. ## Correlation is not cointegration A drunk leaving a bar walks aimlessly: each step is independent of the previous one and, over time, he can end up anywhere. If he has a dog off the leash, the dog wanders too. But if the dog is on a long leash, the two wander aimlessly and, even so, never get further apart than the leash. Murray (1994) used this image to explain **cointegration**: two series that are random walks on their own but whose distance is stationary. In trading, the two walkers are two prices. The question that matters for a pair (lesson 2.3.3) is not whether they go up and down together, but whether the distance between them returns to a level. They are different questions: - **Correlation of returns** measures whether the two prices move in the same direction on the same day. It is a short-term property. Two random walks with a common factor have a high correlation and can drift apart without limit. - **Cointegration** says that there is a $\beta$ such that the spread $\ln y_t - \beta \ln x_t$ is stationary: it returns to its mean at a finite speed. It is a long-term property, and it does not require daily returns to be highly correlated. Formally: two series integrated of order 1 (their differences are stationary, lesson 2.1.2) are cointegrated when there is a vector $(1, -\beta)$ whose combination is stationary. With prices one usually works with their logarithms, because what returns is a ratio and the spread is expressed as a percentage: $$ \ln y_t = \alpha + \beta \ln x_t + u_t, \qquad u_t \sim I(0). $$ ```mermaid flowchart TB P["I(1) prices"] --> R{Correlated} P --> S{"I(0) spread"} R -- yes --> C1[Short-term hedge] R -- no --> N1[Little relationship] S -- yes --> C2[Mean-reverting pair] S -- no --> N2[They drift apart] ``` The two branches are independent: a pair can be in any of the four combinations. The figure at the top shows the difference with data: above, the spread of pair B (cointegrated) crosses zero constantly; below, the residual of control E (more correlated than B, not cointegrated) spends whole years on one side of zero. ## The eight pairs of the course The data for this lesson are eight synthetic pairs (`pares2.parquet`), NYSE sessions from 2014 to 2025, with the **known answer** in `pares2_verdad.csv`. All the relationships were planted between logarithms. Three pairs are cointegrated with different half-lives (A, B and C); D has a β that changes over time (lesson 2.1.6); E, F and G are controls: their returns are highly correlated, but the relationship between their logarithms carries a random walk, so they are not cointegrated; H is cointegrated until June 15, 2021 and stops being so from that day on. They are made-up data with the effect planted: they serve to check whether a method recovers what is known to be there, not as evidence that a real pair works. ```python import numpy as np import pandas as pd pares = pd.read_parquet("pares2.parquet") verdad = pd.read_csv("pares2_verdad.csv", index_col="par") precios = pares.pivot(index="fecha", columns="activo", values="precio") lp = np.log(precios) # the relationships are between logarithms print(precios.shape, precios.index[0].date(), precios.index[-1].date()) verdad[["tipo", "formulacion", "beta", "vida_media_sesiones", "sigma_spread"]] ``` ```output (3018, 16) 2014-01-02 2025-12-31 tipo formulacion beta vida_media_sesiones sigma_spread par A cointegrado log 1.5 5.0 0.020 B cointegrado log 0.8 15.0 0.030 C cointegrado log 1.1 40.0 0.045 D beta_variable log NaN 10.0 0.025 E correlado_no_cointegrado log 1.0 NaN NaN F correlado_no_cointegrado log 0.9 NaN NaN G correlado_no_cointegrado log 1.2 NaN NaN H ruptura log 0.9 10.0 0.025 ``` Each pair has an asset `x` and an asset `y` (`PZAX` and `PZAY` for A). The `formulacion` column says the relationship is between logarithms (`log`), and `sigma_spread` is the standard deviation of the log spread: 0.02 is a typical 2% distance between the two prices. The first thing is to measure what an analyst in a hurry usually measures, correlation, next to the Engle-Granger test explained in the following section: ```python from statsmodels.tsa.stattools import coint filas = [] for par, v in verdad.iterrows(): y, x = lp[v.activo_y], lp[v.activo_x] t, p, _ = coint(y, x) filas.append({"par": par, "tipo": v.tipo, "corr_rend": y.diff().corr(x.diff()), "corr_logprecios": y.corr(x), "t_eg": t, "p_eg": p}) tabla = pd.DataFrame(filas).set_index("par") tabla.round(3) ``` ```output tipo corr_rend corr_logprecios t_eg p_eg par A cointegrado 0.912 1.000 -14.971 0.000 B cointegrado 0.792 0.990 -8.979 0.000 C cointegrado 0.875 0.992 -5.586 0.000 D beta_variable 0.871 0.994 -3.472 0.035 E correlado_no_cointegrado 0.859 0.978 -2.599 0.237 F correlado_no_cointegrado 0.856 0.976 -2.036 0.510 G correlado_no_cointegrado 0.856 0.931 -2.371 0.339 H ruptura 0.635 0.941 -2.647 0.219 ``` *Fig. — Daily return correlation against the Engle-Granger statistic for the eight pairs: correlation does not rank cointegration.* Scatter of eight pairs. Horizontal axis, return correlation from 0.6 to 0.95; vertical, Engle-Granger t statistic from −16 to 0, with the 5% critical value at −3.34. Below it, cointegrated: A (0.91; −15.0), B, highlighted (0.79; −9.0), C (0.88; −5.6) and, narrowly, D (0.87; −3.5). Above it, controls E, F and G, with 0.86 correlation (−2.6, −2.0 and −2.4), and H (0.64; −2.6), whose break also sinks its correlation. The three controls have a correlation of returns of 0.86, more than pair B (0.79), and the correlation of their logarithms is above 0.93. If the criterion were correlation, E would be a better candidate than B. The test says the opposite: B rejects the hypothesis of “no cointegration” comfortably (t = −8.98, far below the 5% critical value of −3.34) and E, F and G do not reject it (p between 0.24 and 0.51). The correlation of prices is an even worse guide: two trending series have correlations of levels close to 1 even if they have nothing to do with each other (spurious regression). D and H are the trick cases. H does not reject in the full sample (p = 0.22): the relationship existed and broke, and the jump of the break drags down even its correlation of returns (0.64). D does reject (p = 0.035), and it is a misleading pass: its true β moves and the full-sample regression estimates an average one; the residual of that fixed β reverts, but much more slowly than the true spread (it shows in the half-life, below). Both are dealt with in this lesson and in 2.1.6. ## Engle-Granger in two steps The Engle and Granger (1987) test is the most direct way to check the definition: 1. Estimate $\ln y_t = \alpha + \beta \ln x_t + u_t$ by least squares. 2. Test whether the residual $\hat u_t$ has a unit root with an augmented Dickey-Fuller regression, $\Delta \hat u_t = \gamma \hat u_{t-1} + \sum_{i=1}^{k} c_i \Delta \hat u_{t-i} + e_t$. The null hypothesis is $\gamma = 0$: the residual is a random walk and there is **no** cointegration. ```mermaid flowchart TB A["ln y, ln x"] --> B[OLS: y on x] B --> C[Residual u] C --> D[ADF of u] D --> E{t below critical} E -- yes --> F[Cointegrated] E -- no --> G[No evidence] ``` The detail that gets forgotten is the critical value. In the first step, OLS chooses the β that minimizes the variance of the residual, that is, the one that makes it look most stationary. That is why the ADF of the residual cannot be compared with the ADF tables for an observed series: more negative critical values are needed, MacKinnon's for cointegration residuals, which are the ones `coint` uses. Doing it by hand shows it: ```python import statsmodels.api as sm from statsmodels.tsa.stattools import adfuller def dos_pasos(y, x): mco = sm.OLS(y, sm.add_constant(x)).fit() adf = adfuller(mco.resid, result_object=True) eg = coint(y, x) return {"alfa": mco.params.iloc[0], "beta": mco.params.iloc[1], "t": adf.statistic, "p_adf": adf.pvalue, "crit5_adf": adf.critical_values["5%"], "p_coint": eg[1], "crit5_coint": eg[2][1]} pd.DataFrame({par: dos_pasos(lp[f"PZ{par}Y"], lp[f"PZ{par}X"]) for par in ["B", "E"]}).round(3) ``` ```output B E alfa 0.683 0.913 beta 0.805 0.821 t -8.978 -2.598 p_adf 0.000 0.093 crit5_adf -2.862 -2.862 p_coint 0.000 0.237 crit5_coint -3.338 -3.338 ``` In B, OLS recovers the planted relationship ($\alpha$ = 0.7 and $\beta$ = 0.8; it gets 0.68 and 0.805) and the residual is clearly stationary with any table. In E the statistic is the same in both columns (−2.60), but the plain ADF gives p = 0.09, a rejection at 10%, and `coint`, p = 0.24: the 5% critical value goes from −2.86 to −3.34. In a screen of many pairs, using the wrong table multiplies the false positives. Engle-Granger has two known weaknesses. The first is that it is not symmetric: regressing `y` on `x` or `x` on `y` gives different statistics and p-values, and in a borderline case each direction can lead to a different conclusion. The second is that with more than two series it only finds one relationship, even if there are several. Johansen solves both. ## Johansen: how many relationships there are among several series Johansen (1991) starts from a vector error correction model (VECM) for a vector $Y_t$ of $n$ series, here the logarithms of the prices: $$ \Delta Y_t = \Pi Y_{t-1} + \sum_{i=1}^{k} \Gamma_i \Delta Y_{t-i} + \varepsilon_t . $$ If the series are not cointegrated, $\Pi = 0$ and the model is a VAR in differences. If there are $r$ cointegrating relationships, $\Pi = \alpha \beta'$ has rank $r$: the $r$ columns of $\beta$ are the stationary combinations and $\alpha$ says how strongly each series returns to equilibrium. The trace test asks in sequence “rank ≤ 0?”, “rank ≤ 1?”… and the estimated rank is the first one that is not rejected. With the four series of A and B, the known answer is $r = 2$: ```python import warnings from statsmodels.tsa.vector_ar.vecm import coint_johansen # With numpy 2.5, the eigenvectors come out as complex numbers with a # zero imaginary part and statsmodels warns (ComplexWarning) when # casting them to real: the warning does not change any number warnings.filterwarnings("ignore", category=np.exceptions.ComplexWarning) def johansen(datos, k_ar_diff=1): j = coint_johansen(datos, det_order=0, k_ar_diff=k_ar_diff) n = datos.shape[1] tabla = pd.DataFrame( {"traza": j.lr1, "crit95_traza": j.cvt[:, 1], "max_autov": j.lr2, "crit95_max": j.cvm[:, 1]}, index=[f"rank <= {r}" for r in range(n)]) rango = int((j.lr1 > j.cvt[:, 1]).cumprod().sum()) return tabla, rango, np.real(j.evec) cuatro = lp[["PZAY", "PZAX", "PZBY", "PZBX"]] tabla_j, rango, vectores = johansen(cuatro) print("Estimated rank:", rango) tabla_j.round(1) ``` ```output Estimated rank: 2 traza crit95_traza max_autov crit95_max rank <= 0 296.2 47.9 209.6 27.6 rank <= 1 86.6 29.8 80.0 21.1 rank <= 2 6.6 15.5 5.6 14.3 rank <= 3 0.9 3.8 0.9 3.8 ``` The trace rejects “rank ≤ 0” (296.2 versus 47.9) and “rank ≤ 1” (86.6 versus 29.8), and does not reject “rank ≤ 2” (6.6 versus 15.5): two relationships, the planted ones. The cointegrating vectors are not unique (any combination of two stationary vectors is stationary too), but once normalized they are recognizable: ```python v1 = vectores[:, 0] / vectores[0, 0] # normalized by PZAY v2 = vectores[:, 1] / vectores[2, 1] # normalized by PZBY pd.DataFrame({"vector_1": v1, "vector_2": v2}, index=cuatro.columns).round(3) ``` ```output vector_1 vector_2 activo PZAY 1.000 -0.040 PZAX -1.501 0.065 PZBY 0.070 1.000 PZBX -0.059 -0.812 ``` The first vector is almost exactly the relationship of A, $\ln y - 1.5 \ln x$ (−1.501), with small weights on B; the second, almost that of B, $\ln y - 0.8 \ln x$ (−0.812), with small weights on A. Here they come out clean because the two relationships have very different reversion strengths, but that is not guaranteed: any combination of two stationary vectors is stationary too, and what Johansen identifies uniquely is the number of relationships and the space they span, not a particular basis. If a series that has nothing to do with a cointegrated pair is added to it, Johansen should find a single relationship and give zero weight to the unrelated series. And with series that are not cointegrated with each other it should give rank 0, except by chance (at 5%, in one case in twenty). This last check depends on the formulation: ```python from itertools import combinations tabla_3, rango_3, vec_3 = johansen(lp[["PZAY", "PZAX", "PZEX"]]) print("A + PZEX: rank", rango_3, "vector", (vec_3[:, 0] / vec_3[0, 0]).round(3)) controles = ["PZEY", "PZEX", "PZFY", "PZFX", "PZGY", "PZGX"] trios = list(combinations(controles, 3)) falsos = {"levels": 0, "logs": 0} for trio in trios: datos = precios[list(trio)] falsos["levels"] += johansen(datos)[1] > 0 falsos["logs"] += johansen(lp[list(trio)])[1] > 0 print(len(trios), "control trios; with rank > 0:", falsos) crec = precios[controles].iloc[-1] / precios[controles].iloc[0] print("Final / initial price:", crec.round(1).to_dict()) ``` ```output A + PZEX: rank 1 vector [ 1. -1.51 0.008] 20 control trios; with rank > 0: {'levels': 4, 'logs': 0} Final / initial price: {'PZEY': 2.4, 'PZEX': 3.1, 'PZFY': 3.3, 'PZFX': 3.4, 'PZGY': 5.1, 'PZGX': 3.2} ``` With the unrelated series, rank 1 and an almost zero weight (0.008) for `PZEX`: correct. With the controls, in logs none of the 20 trios gives a rank above zero, as expected (chance would give one); in levels, 4. The reason is in the last line: the controls multiply by between 2.4 and 5.1 in twelve years. The test assumes errors with constant variance, and in levels the daily change of a price is proportional to its level: a price that has multiplied by 5 has changes five times larger at the end than at the start, and the distribution of the trace is no longer that of the tables. In logs the daily change is a percentage and its size does not depend on the level. ## Half-life: the spread as an Ornstein-Uhlenbeck process Knowing that a spread reverts is not enough to trade it: you need to know how fast. The standard model is the **Ornstein-Uhlenbeck** process, a random walk with a force that pulls toward the mean: $$ dS_t = \theta (\mu - S_t)\,dt + \sigma\, dW_t . $$ $\theta$ is the speed of reversion: the expected distance to the mean decays like $e^{-\theta t}$. Sampled once per session, it is an AR(1), $S_t = c + \varphi S_{t-1} + \varepsilon_t$ with $\varphi = e^{-\theta}$, and the **half-life** (the time in which the expected distance to the mean is halved) is $$ h = \frac{\ln 2}{\theta} = -\frac{\ln 2}{\ln \varphi} . $$ It is estimated with the regression of the change in the spread on its previous level, $\Delta S_t = c + b\, S_{t-1} + e_t$, with $\varphi = 1 + b$. In the cointegrated pairs it is applied to the Engle-Granger residual: ```python def vida_media(spread): """Half-life in sessions of the AR(1) dS_t = c + b * S_(t-1).""" s = pd.Series(spread).dropna() X = sm.add_constant(s.shift(1).iloc[1:]) b = sm.OLS(s.diff().iloc[1:], X).fit().params.iloc[1] return -np.log(2) / np.log(1 + b) def residuo_eg(par): y, x = lp[f"PZ{par}Y"], lp[f"PZ{par}X"] return sm.OLS(y, sm.add_constant(x)).fit().resid vidas = pd.DataFrame({ "estimada": {p: vida_media(residuo_eg(p)) for p in "ABCDEFG"}, "verdadera": verdad["vida_media_sesiones"]}).dropna( subset=["estimada"]) vidas.round(1) ``` ```output estimada verdadera A 4.7 5.0 B 13.0 15.0 C 35.4 40.0 D 38.4 10.0 E 162.1 NaN F 259.4 NaN G 270.0 NaN ``` In A, B and C the estimate is close to the truth but always below it: 4.7 versus 5, 13.0 versus 15 and 35.4 versus 40. In E, F and G the function returns half-lives of 162 to 270 sessions, even though there are none: an AR(1) always gives a $\varphi$ below 1 in a finite sample, so the half-life of a random walk comes out finite. **A half-life only makes sense after a cointegration test that supports it.** And D shows the opposite: the test passed, but the residual of its fixed β has a half-life of 38.4 sessions when that of the true spread is 10. What reverts in D is something other than what was estimated: with a β that moves, the residual of an average β mixes the spread with the slow drift of β (lesson 2.1.6). Is it bad luck or bias? It can be answered by simulation: generate many Ornstein-Uhlenbeck processes with the true half-life and the length of this data, and estimate each one the same way. ```python from scipy.signal import lfilter rng = np.random.default_rng(2028) n, caminos = len(lp), 400 filas = {} for h in [5, 15, 40]: phi = 0.5 ** (1 / h) e = rng.standard_normal((caminos, n)) * np.sqrt(1 - phi**2) s = lfilter([1.0], [1.0, -phi], e, axis=1) ds, s0 = np.diff(s, axis=1), s[:, :-1] s0c = s0 - s0.mean(axis=1, keepdims=True) b = (ds * s0c).sum(axis=1) / (s0c**2).sum(axis=1) est = -np.log(2) / np.log(1 + b) filas[h] = np.percentile(est, [5, 50, 95]) pd.DataFrame(filas, index=["p5", "mediana", "p95"]).T.round(1) ``` ```output p5 mediana p95 5 4.4 4.9 5.6 15 11.8 14.7 17.8 40 27.3 37.2 51.2 ``` *Fig. — Half-life of the cointegrated pairs: estimated in the data (diamond), true value (tick) and 90% band of 400 simulations of the same length.* Three rows on a log scale from 3 to 70 sessions. A: true 5, estimated 4.7, band from 4.4 to 5.6. B: true 15, estimated 13.0, band from 11.8 to 17.8. C, highlighted: true 40, estimated 35.4, band from 27.3 to 51.2. All three estimates fall below the truth and inside their band, which widens the slower the reversion. The median of the simulations is below the truth, and more so the slower the reversion: it is the finite-sample bias of the autoregressive coefficient, which pushes $\hat\varphi$ down. The three estimates from the data fall inside the 90% band, so they are what one would expect. The band is what matters most: with twelve years of data, the half-life of 40 sessions could have been estimated at 27 or at 51. A pair with a half-life of 40 sessions has barely 75 half-lives in the sample, few independent observations of the reversion. The half-life is used later for two decisions of the system in lesson 2.3.3: how long to hold a position before calling it failed and which window to use to estimate the mean and deviation of the spread. A badly estimated half-life spoils both. ## When cointegration breaks: pair H Pair H is cointegrated (half-life 10, spread with a standard deviation of 0.025 in logs) until June 15, 2021. That day the spread jumps five standard deviations and from then on it is a random walk. It is the most dangerous real-world situation: a relationship that worked for years, that any test confirms with past data, and that ceases to exist. ```python y, x = lp["PZHY"], lp["PZHX"] antes = lp.index < "2021-06-15" for nombre, m in [("full sample", slice(None)), ("before the break", antes), ("after the break", ~antes)]: t, p, _ = coint(y[m], x[m]) print(f"{nombre:22} sessions {len(y[m]):5} t {t:6.2f} " f"p {p:.4f}") ``` ```output full sample sessions 3018 t -2.65 p 0.2191 before the break sessions 1875 t -5.91 p 0.0000 after the break sessions 1143 t -3.66 p 0.0205 ``` *Fig. — Spread of pair H with the β estimated before the break, 2014–2025: cointegrated until 15 June 2021, a random walk since then.* Line from 2014 to 2025. Until June 2021 the spread oscillates around zero within a few percentage points (Engle-Granger, p < 0.0001). The 15 June 2021 break, highlighted, makes it jump, and from there it wanders without returning, between 12 and 47%; even so, Engle-Granger on that stretch gives p = 0.021, a false positive. Before the break the cointegration is unambiguous (t = −5.91; p < 0.0001). The full sample no longer detects it (p = 0.22). And the later stretch, which is a random walk by construction, gives t = −3.66 and p = 0.021: a **false positive**. An analyst who looked only at the last four years would conclude that the relationship is still alive. A 5% test rejects one time in twenty when there is nothing; if you look at a single stretch, chance can strike. One way to see the whole problem is to test in rolling one-year windows, as a system that re-estimates the relationship every month would: ```python def ventanas(y, x, n=252, paso=21): return np.array([coint(y.iloc[i:i + n], x.iloc[i:i + n], maxlag=10)[1] for i in range(0, len(y) - n + 1, paso)]) casos = {"A (half-life 5)": "A", "C (half-life 40)": "C", "E (control)": "E", "F (control)": "F", "G (control)": "G"} filas = {} for nombre, par in casos.items(): p = ventanas(lp[f"PZ{par}Y"], lp[f"PZ{par}X"]) filas[nombre] = (len(p), (p < 0.05).mean()) p = ventanas(y[~antes], x[~antes]) filas["H after the break"] = (len(p), (p < 0.05).mean()) pd.DataFrame(filas, index=["ventanas", "rechazo_5pct"]).T.round(3) ``` ```output ventanas rechazo_5pct A (half-life 5) 132.0 0.977 C (half-life 40) 132.0 0.091 E (control) 132.0 0.038 F (control) 132.0 0.098 G (control) 132.0 0.061 H after the break 43.0 0.070 ``` Three readings. In A, with a half-life of 5 sessions, one year of data is enough and the test rejects in 98% of the windows. In C, truly cointegrated but with a half-life of 40, one year contains barely six half-lives and the test rejects in 9%, almost as if there were nothing: it is a problem of **power**, not of the absence of a relationship. In the controls the rejection rate ranges from 3.8% to 9.8%, around the 5% expected from chance (the windows overlap, so those rates are very noisy), and in H after the break it is 7% (3 of 43 windows). The short window does not tell a slow pair from a non-existent one well; the long one mixes stretches with and without a relationship. There is no right window: you have to know what is lost with each one. Detecting the break itself is the subject of lesson 2.1.7. ## Screening pairs: selection bias In practice nobody tests one pair: hundreds are screened and the ones that pass are kept. With the sixteen assets in the data there are 120 possible pairings, and only those of A, B and C (and H for part of the sample) have a relationship. The 112 pairings of assets from different pairs are not cointegrated by construction: ```python filas = [] for a, b in combinations(precios.columns, 2): filas.append({ "mismo_par": a[:3] == b[:3], "p_logs": coint(lp[a], lp[b], maxlag=10)[1], "p_niveles": coint(precios[a], precios[b], maxlag=10)[1]}) criba = pd.DataFrame(filas) ajenas = criba[~criba["mismo_par"]] print("Pairings from different pairs:", len(ajenas)) (ajenas[["p_logs", "p_niveles"]] < 0.05).sum() ``` ```output Pairings from different pairs: 112 p_logs 4 p_niveles 2 dtype: int64 ``` In logs, 4 of the 112 unrelated pairings pass the test at 5%, and in levels, 2: both numbers are close to the 5.6 that chance would give, and there is no lesson about the formulation between them (the pairings share a market factor and are not independent tests, so the number of rejections varies a lot from one sample to another). What matters is that they exist: any of those four pairings, presented on its own, would look like a finding. The problem of picking the best of many tests has a name and corrections (lesson 2.1.8); the simplest defense is to require an economic reason before looking at the p-value and to validate outside the sample in which the pair was chosen (lesson 2.5.2). ## Common mistakes with cointegration - **Estimating β on the full sample and trading on it.** The β and the mean of the spread in a backtest have to be estimated with data prior to each trade; with the full sample, the spread reverts “by construction” (lesson 2.5.1). - **Screening without correcting.** Of 112 unrelated pairings, 4 pass in logs. With hundreds of candidates, the “pairs” that pass are largely chance. - **Confusing correlation and cointegration.** The controls E, F and G have more correlation than B and do not revert. - **Trusting a pass with a fixed β.** D passes Engle-Granger with an average β, but the residual that reverts has a half-life four times longer than the real spread. - **Choosing the wrong formulation.** In levels, Johansen gives false positives and Engle-Granger misses pair A. - **Ignoring power.** A slow pair looks non-cointegrated in short windows; a fast pair looks cointegrated in almost any window. - **Assuming it lasts.** The relationship of H worked for seven years and broke in a day. Monitoring the spread and the stability of β (lessons 2.1.6 and 2.1.7) is part of the system, not an extra. - **Unadjusted prices.** An unadjusted split or dividend puts a jump in one of the two prices and breaks a cointegration that exists (lesson 1.4.3). ## Exercises 1. Repeat the Engle-Granger test for pairs A, B, C and for the controls E, F and G with prices in levels instead of in logarithms. Does any conclusion change? **Solution.** ```python for par in "ABCEFG": t, p, _ = coint(precios[f"PZ{par}Y"], precios[f"PZ{par}X"]) print(par, f"t {t:6.2f} p {p:.3f}") ``` ```output A t -3.30 p 0.055 B t -8.56 p 0.000 C t -5.44 p 0.000 E t -2.32 p 0.366 F t -1.85 p 0.607 G t -2.04 p 0.508 ``` The controls still do not reject, B and C still reject, but A stops doing so: in levels its statistic goes from −14.97 to −3.30 and the p-value to 0.055. A is the pair whose `y` has grown the most (it multiplies by 5.5 in twelve years), and a relationship between logarithms with β = 1.5 is not linear in levels: the money distance between the two prices grows with the level and the residual of the regression in levels looks like a random walk. It is the other side of the formulation error: not only false positives, but also real relationships that go unseen. 2. Estimate the half-life of pair B with the first three years (2014-2016) and with the last three (2023-2025), using in each stretch its own β. How much does it vary? Does it fall inside the simulated band for a half-life of 15? **Solution.** ```python yb, xb = lp["PZBY"], lp["PZBX"] for desde, hasta in [("2014", "2016"), ("2023", "2025")]: ys, xs = yb.loc[desde:hasta], xb.loc[desde:hasta] u = sm.OLS(ys, sm.add_constant(xs)).fit().resid print(desde, hasta, f"sessions {len(u)}", f"half-life {vida_media(u):.1f}") ``` ```output 2014 2016 sessions 756 half-life 9.7 2023 2025 sessions 752 half-life 13.5 ``` With three years, some 750 sessions, the half-life comes out at 9.7 and 13.5 sessions, both below 15 and the first outside the simulated band for 3,018 sessions (from 11.8 to 17.8): with a quarter of the data the band widens and the downward bias grows. In addition, each stretch estimates its own β, and a β estimated on a small sample absorbs part of the spread and makes it look faster. A half-life is an estimate with a large error: use it as an order of magnitude. 3. A colleague screens 300 pairings of stocks from the same sector in levels, finds 24 with p < 0.05 in Engle-Granger and proposes trading all of them. With what you have seen in this lesson, what would you tell them? **Solution.** That with 300 tests at 5%, chance alone already produces about 15 rejections (in the course data, 4 of 112 unrelated pairings), and that in levels the test does not even measure what it should: with prices that grow multiplicatively, it gives false positives and misses real relationships. The 24 may be almost all false. Before trading: fix the formulation (logarithms), require an economic reason for each pair, correct for multiple testing (lesson 2.1.8), estimate the half-life and check that it is compatible with the system's horizon, and validate outside the sample in which the screen was run, monitoring the stability of the relationship as with pair H. ## Self-assessment 1. The daily returns of pair E have a correlation of 0.86 and those of pair B, 0.79. What can be concluded about their cointegration? - A. Nothing: correlation of returns measures the joint move of a single day, not whether the distance between the prices returns to a level - B. That E is more cointegrated than B - C. That neither is cointegrated, because the correlation does not reach 0.9 - D. That both are cointegrated, because the correlation is above 0.8 *Answer A: Nothing: correlation of returns measures the joint move of a single day, not whether the distance between the prices returns to a level.* In the course data B is cointegrated (Engle-Granger gives t = −8.98) and E is not (t = −2.60, p = 0.24), even though E has more correlation. Correlation of returns and cointegration measure different things: the first, the short term; the second, a long-run equilibrium. 2. Why can't Engle-Granger use the usual ADF critical values in the second step? - A. Because the residual does not have zero mean - B. Because the ADF only works with returns - C. Because the β of the first step is chosen by OLS so that the residual has the smallest possible variance: the residual looks more stationary than it is, and the critical values have to be more negative - D. Because the residual is always a random walk *Answer C: Because the β of the first step is chosen by OLS so that the residual has the smallest possible variance: the residual looks more stationary than it is, and the critical values have to be more negative.* With the same statistic of −2.60, the plain ADF on the residual of pair E gives p = 0.09 and `coint`, with MacKinnon's critical values for residuals, p = 0.24. The 5% critical value goes from −2.86 to −3.34. 3. Johansen on four series (the two assets of A and the two of B) gives rank 2. What does it mean? - A. That two of the four series are stationary - B. That there are two independent linear combinations of the four series that are stationary: the relationships of A and of B - C. That the series have to be differenced twice - D. That the test has failed *Answer B: That there are two independent linear combinations of the four series that are stationary: the relationships of A and of B.* The rank of the VECM's Π matrix is the number of cointegrating relationships. With four non-stationary log prices and two independent pairs, the known answer is 2, and that is what comes out. 4. The half-life estimator from the regression of ΔS on S is biased. In which direction, and when does it matter most? - A. Upward, and more the faster the reversion - B. It is unbiased if the sample has more than a thousand sessions - C. Downward only if the spread has fat tails - D. Downward (the spread seems to revert faster), and more the longer the half-life relative to the sample *Answer D: Downward (the spread seems to revert faster), and more the longer the half-life relative to the sample.* It is the finite-sample bias of the autoregressive coefficient: the estimated φ ends up below the true one. With 3,018 sessions, pair C's half-life of 40 comes out at 35.4, and the simulation gives a median of 37.2 with 90% of cases between 27.3 and 51.2. 5. After the 2021 break, the final stretch of pair H gives p = 0.021 in Engle-Granger. What is it? - A. Proof that the cointegration is back - B. A data error - C. A false positive: since the break the spread is a random walk by construction, and a stretch of a few years can pass the test by chance - D. The sign that β has changed and has to be re-estimated *Answer C: A false positive: since the break the spread is a random walk by construction, and a stretch of a few years can pass the test by chance.* The known answer is that since 2021-06-15 the spread of H is a random walk. A 5% test rejects by chance one time in twenty; in one-year windows after the break, it rejects in 3 of 43, close to what chance would produce. ## Key points - Correlation of returns and cointegration measure different things: the controls E, F and G have a correlation of 0.86 and are not cointegrated; B has 0.79 and is. - Engle-Granger: OLS of one log price on the other and an ADF of the residual with its own critical values (MacKinnon), stricter than the ADF ones. - Johansen tests the number of cointegrating relationships among several series: rank 2 in the four series of A and B, rank 1 with an unrelated series. - The half-life of the spread, ln 2 / θ, comes from an AR(1); the estimator underestimates it, more so the slower the reversion. - Mistakes come from the formulation (in levels, Johansen gives a rank above 0 in 4 of 20 control trios; in logs, in none), from screening many pairs, from static βs that pass with a slow spread (D) and from short stretches after a break (H). ## Resources ### Research - [Stochastic Processes - Random Walks](https://pythonparatrading.com/en/stochastic-processes/) · Research-0005 · 11 min ### Documentation - [statsmodels: coint (Engle-Granger)](https://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.coint.html) · statsmodels.org - [statsmodels: coint_johansen](https://www.statsmodels.org/stable/generated/statsmodels.tsa.vector_ar.vecm.coint_johansen.html) · statsmodels.org ### Books - [Engle, R. F. and Granger, C. W. J. (1987): Co-Integration and Error Correction: Representation, Estimation, and Testing. Econometrica 55(2), 251-276](https://doi.org/10.2307/1913236) · doi.org - [Johansen, S. (1991): Estimation and Hypothesis Testing of Cointegration Vectors in Gaussian Vector Autoregressive Models. Econometrica 59(6), 1551-1580](https://doi.org/10.2307/2938278) · doi.org - [Murray, M. P. (1994): A Drunk and Her Dog: An Illustration of Cointegration and Error Correction. The American Statistician 48(1), 37-39](https://doi.org/10.1080/00031305.1994.10476017) · doi.org ### Data - [Level 2 data: eight synthetic pairs, 2014-2025 (Parquet) and their truth table](https://pythonparatrading.com/curso/datos/nivel-2/pares2.parquet) · Download · pares2.parquet Previous lesson: [2.1.4 · Volatility models: ARCH, GARCH, EGARCH and GJR with the arch library](https://pythonparatrading.com/en/python-trading-course/level-2/volatility-models-arch-garch-egarch-gjr/) Next lesson: [2.1.6 · Kalman filter for dynamic hedge ratios](https://pythonparatrading.com/en/python-trading-course/level-2/kalman-filter-dynamic-hedge-ratio/)