Cointegration: The Johansen Method

The previous article left a loose end: the deepest limitation of Engle-Granger is the cointegration rank. A system of $K$ variables can carry up to $K - 1$ cointegrating relations, and Engle-Granger tests exactly one at a time — by picking a variable for the left-hand side of the regression, you are effectively betting that the system contains exactly one stationary combination. With two variables that bet is safe. With more, there are more ways to lose it: the test loses power, and the “which variable goes on the left” normalization stops being harmless, because different choices land on different directions inside the cointegrating space.

Johansen (1988, 1991, 1995) turns the problem around. Instead of guessing one relation and testing its residual, write the whole system as a single model and let maximum likelihood answer three questions at once: how many cointegrating relations there are (the rank test), what they are (the vector estimates), and who does the adjusting (the adjustment matrix). The price is a heavier derivation — reduced rank regression and a family of non-standard asymptotic distributions — but what you buy is formal inference on the long-run parameters, which Engle-Granger cannot give you. That is why Johansen is the default for multivariate empirical work.

This article derives the Johansen approach from the VECM, derives the trace and max eigenvalue tests, and then deals with a question that is often skipped but decides whether your results mean anything: identification of the cointegrating vectors. The textbook references are Hamilton (1994, Ch. 20), Lütkepohl (2005, Ch. 6-7), and Tsay (2010, Ch. 8).

From VAR to VECM: Everything Rides on $\Pi$

Start with a $K$-variable VAR($p$) in levels:

where $\mathbf{D}_t$ holds deterministic terms (constant, trend). Subtract $\mathbf{y}_{t-1}$ from both sides and rearrange the lags:

where $\Pi = \sum_{j=1}^p \mathbf{A}_j - \mathbf{I}$ and $\Gamma_j = -\sum_{i=j+1}^p \mathbf{A}_i$. This is the VECM (Vector Error Correction Model). Note that it is just a re-parameterization of the VAR — no new information — but it turns cointegration from an implicit property into explicit structure. The left side, $\Delta \mathbf{y}_t$, is stationary, and so are the lagged differences on the right, which forces $\Pi \mathbf{y}_{t-1}$ to be stationary too. A matrix that maps a non-stationary vector into something stationary can only carry so much information, and how much is exactly its rank.

Granger’s representation theorem (Engle & Granger, 1987) makes this precise: if $\mathbf{y}_t$ is I(1) and cointegrated with rank $r$, then $\Pi = \boldsymbol{\alpha} \boldsymbol{\beta}'$, where $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are $K \times r$ matrices of full column rank. The columns of $\boldsymbol{\beta}$ are the cointegrating vectors; the columns of $\boldsymbol{\alpha}$ are the adjustment speeds. Substituting back:

The product $\boldsymbol{\beta}' \mathbf{y}_{t-1}$ is an $r$-vector of stationary cointegrating relations — the multivariate version of the scalar error-correction term $\hat{u}_{t-1}$ from the previous article. The matrix $\boldsymbol{\alpha}$ says how much each variable corrects for a deviation in each relation.

The rank of $\Pi$ tells you which kind of system you are holding:

So the first job in a Johansen analysis is not estimating cointegrating vectors — it is answering “what is $r$”. That is what the trace and max eigenvalue tests do. But to test the rank you first need the likelihood under a given rank constraint, so estimation comes first.

Maximum Likelihood: Reduced Rank Regression

Johansen’s central move: treat $\Pi$ as the unknown matrix and maximize the Gaussian likelihood subject to “rank at most $r$”. That sounds like a nasty constrained optimization; it actually has a closed-form solution, via concentrated likelihood and reduced rank regression.

First concentrate out the short-run dynamics. Look at the VECM again:

This is a multivariate regression of $\Delta \mathbf{y}_t$ on $\mathbf{y}_{t-1}$, the lagged differences, and $\mathbf{D}_t$. Wash out the nuisance terms in two steps: regress $\Delta \mathbf{y}_t$ on the lagged differences and $\mathbf{D}_t$, keeping the residuals $\mathbf{R}_{0t}$; then regress $\mathbf{y}_{t-1}$ on the same regressors, keeping the residuals $\mathbf{R}_{1t}$. The problem is now compressed into the relationship between $\mathbf{R}_{0t}$ and $\mathbf{R}_{1t}$. Define the moment matrices $\mathbf{S}_{01} = T^{-1} \sum \mathbf{R}_{0t} \mathbf{R}_{1t}'$ and $\mathbf{S}_{11} = T^{-1} \sum \mathbf{R}_{1t} \mathbf{R}_{1t}'$, and analogously $\mathbf{S}_{00}$ and $\mathbf{S}_{10}$.

Without a rank constraint, the ML estimate is just OLS:

with rank $K$. To impose rank $r$, solve the generalized eigenvalue problem:

The resulting eigenvalues $\hat{\lambda}_1 > \hat{\lambda}_2 > \cdots > \hat{\lambda}_K$ are the squared canonical correlations between $\mathbf{R}_{0t}$ and $\mathbf{R}_{1t}$ — each $\hat{\lambda}_i$ measures how much explanatory power some linear combination of $\mathbf{y}_{t-1}$ has for $\Delta \mathbf{y}_t$. The cointegrating relations are precisely the directions where that explanatory power does not vanish, so the rank-$r$ ML estimate takes the first $r$ eigenvectors:

The adjustment speeds follow:

One warning: the eigenvalue decomposition hands you $\hat{\boldsymbol{\beta}}$ under a normalization that is mathematically convenient and economically arbitrary (typically $\hat{\boldsymbol{\beta}}' \mathbf{S}_{11} \hat{\boldsymbol{\beta}} = \mathbf{I}$). To make the coefficients economically readable, you have to re-normalize — that is the identification problem, treated in its own section below.

The Rank Tests: Trace and Max Eigenvalue

With the estimation problem solved, rank testing comes almost for free: the maximized likelihood under any rank constraint can be written in terms of the eigenvalues, and comparing them is a likelihood ratio test.

The trace test runs a sequential test of the null that the rank is at most $r$:

with statistic:

The reading is intuitive: under the null, the eigenvalues past position $r$ should all be zero; the statistic sums up these “should-be-zero” eigenvalues, and you reject when the sum is too large. The sequential procedure walks from the bottom up — test $H_0(0)$, and a rejection means at least one cointegrating relation; test $H_0(1)$, and a rejection means at least two; keep going until the first failure to reject, and the rank is the $r$ of that first non-rejected hypothesis. The asymptotic distribution is non-standard (a functional of Brownian motion), tabulated by Johansen (1988) and Osterwald-Lenum (1992); the critical values depend on $K$ and on the deterministic specification.

The max eigenvalue test takes one step at a time:

testing $H_0$: rank $= r$ against $H_1$: rank $= r + 1$. Its asymptotic distribution is also non-standard, and it is not the same as the trace distribution. The two tests usually agree but can disagree: the trace test sums all remaining eigenvalues and has more power against alternatives with several small eigenvalues; the max eigenvalue test looks only at the next one and has more power against a single dominant eigenvalue. When they disagree there is no agreed arbiter; in practice most applied work treats the trace test as the primary verdict and the max eigenvalue as corroboration.

That is the formal procedure. Now the honest part: sequential testing carries a built-in multiple-testing problem, and the non-standard statistics over-reject in small samples (the workflow section returns to this). So do not fix the rank from the test table alone. Plot the eigenvalues — a clean cliff between large and small ones is a strong hint about the rank; five eigenvalues at 0.9, 0.85, 0.1, 0.05, 0.02 point at rank 2. Estimate with $r = 1, 2, 3$ and check which set of cointegrating vectors makes economic sense: vectors with huge coefficients and no economic reading usually mean you over-ranked. And when theory makes a sharp prediction, use it — UIP plus PPP in an open economy predicts two cointegrating relations, so test $H_0(2)$ directly. Some authors select the rank with AIC or BIC across VECMs with different $r$, but this is far less standard than it is for VAR order selection.

Identification: The Likelihood Only Pins Down a Rotation

Now the question that decides whether your results can be read at all. The ML estimate $\hat{\boldsymbol{\beta}}$ is identified only up to a rotation: for any $r \times r$ non-singular matrix $\mathbf{H}$, $\hat{\boldsymbol{\beta}} \mathbf{H}$ spans the same cointegrating space and produces exactly the same likelihood. The likelihood recognizes the space, not the particular vectors. Without restrictions, the columns of $\hat{\boldsymbol{\beta}}$ carry no economic meaning — every coefficient you report could be rotated away.

When $r = 1$ the problem barely exists: the cointegrating vector is unique up to scale, so pick a variable, normalize its coefficient to 1, done — this is exactly the Engle-Granger normalization. The trouble starts at $r > 1$. The standard just-identifying scheme normalizes the leading $r \times r$ block of $\boldsymbol{\beta}$ to the identity: each relation is normalized on one variable, and that variable is excluded from the other relations. Count the degrees of freedom: $K \times r$ parameters minus $r^2$ restrictions leaves $r(K - r)$ free parameters — just identified, nothing to spare. But note what this buys you: mathematical uniqueness only. Which variable anchors which relation is still your choice, and if that choice has no economic rationale, the estimates are unique and still unreadable.

If you want more structure, there are several routes: long-run restrictions in the Blanchard-Quah style (certain shocks have zero long-run effect on certain variables — requires a fully structural model), sign restrictions as in SVARs, or prior information from theory, earlier empirical work, or institutional knowledge written directly into constraints. In practice, the vast majority of applied work uses the normalization scheme and hopes theory is informative enough to interpret the result — fine at $r = 1$; at $r > 1$, either add structure or accept that the reported vectors will not read.

Once the identification scheme is set, theory can be put on trial. What theory typically predicts is not just “cointegration exists” but the exact shape of the vector. The classics from finance and macro: money demand says $\log M_t - \log P_t - \log Y_t$ is stationary ($M$ money supply, $P$ price level, $Y$ income), cointegrating vector $(1, -1, -1)$; PPP says $\log E_t - \log P_t + \log P_t^*$ is stationary ($E$ nominal exchange rate, $P$ and $P^*$ domestic and foreign price levels), vector $(1, -1, 1)$; the term structure says the spread $R_t - r_t$ is stationary, vector $(1, -1)$ — the financial example in the previous article is its bivariate special case; the present-value model says $\log P_t - \log D_t$ is stationary ($P$ stock price, $D$ dividend), vector $(1, -1)$.

Every one of these predictions can be put to a vote. A restricted model has fewer free parameters, so its maximized log-likelihood is lower or equal, and the likelihood ratio statistic:

where $q$ is the number of restrictions; reject the restrictions when $LR$ is too large. The procedure: the unrestricted Johansen estimate tells you what $\boldsymbol{\beta}$ the data prefers, theory tells you what $\boldsymbol{\beta}$ should look like, and the LR test rules on whether the two are compatible. Say you believe the cointegrating vector is $(1, -1, 0, 0)$ — the first two variables enter with opposite signs and unit coefficients, the others do not enter — that is $q$ testable restrictions. The asymptotic distribution is a plain $\chi^2$, which is one of the places the Johansen framework earns its keep over Engle-Granger: hypothesis tests on long-run parameters are routine here.

Weak Exogeneity: Who Gets to Ignore the Equilibrium

The previous article introduced weak exogeneity from a single-equation angle; here it acquires matrix form. A variable $y_{it}$ is weakly exogenous for the cointegrating parameters if the $i$-th row of $\boldsymbol{\alpha}$ is zero: $y_{it}$ does not respond to any equilibrium deviation, and the equilibrium is restored entirely by the other variables. The test is direct — a $t$-test (or Wald test) on $\alpha_i = 0$ inside the VECM. If you cannot reject, $y_{it}$ can be treated as a conditioning variable and the cointegrating parameters estimated in a partial system that drops its dynamic equation.

What is the point? Efficiency. Johansen (1992) supplies the formal framework: when the weak exogeneity restriction is true, the partial-system estimates are more efficient — in small samples, estimating fewer equations genuinely shrinks the standard errors. When the restriction is false, the long-run parameter estimates typically remain consistent, but the efficiency gain is lost and inference based on the partial system is no longer reliable. So this is a “test first, then use” tool, not a default.

A Code Snippet: Johansen in Python

import pandas as pd
from statsmodels.tsa.vector_ar.vecm import coint_johansen


result = coint_johansen(y, det_order=0, k_ar_diff=1)  # det_order=-1: no deterministic terms, =0: constant, =1: linear trend

print("Eigenvalues:", result.eig)
print("Trace statistics:", result.lr1)  # for r=0,1,..,K-1
print("Trace critical values (90%, 95%, 99%):\n", result.cvt)
print("Max eigenvalue stats:", result.lr2)
print("Max eigenvalue critical values:\n", result.cvm)

# Cointegrating vectors (eigenvectors)
print("Cointegrating vectors (columns):\n", result.evec)

# Plot for trace test
import matplotlib.pyplot as plt
plt.plot(result.lr1, label='trace stat')
plt.plot(result.cvt[:, 1], label='95% crit')
plt.xlabel('r'); plt.legend(); plt.show()

Mind the det_order encoding: $-1$ means no deterministic terms, $0$ a constant, $1$ a linear trend. This argument is not decoration — get the deterministic specification wrong and the critical values are wrong, which is the next section’s business.

A Workflow That Holds Up

Johansen results are sensitive to the specification, and the sensitivity concentrates in three places: the lag order, the deterministic terms, and the sample length. Rather than doubting the output afterwards, put the checks in front. A workable order:

1. Confirm integration orders. Run ADF on each variable. The method assumes the variables are I(1) and does not test this premise itself.
2. Choose the lag order. Select $p$ for the VAR in levels via information criteria (AIC, BIC); the corresponding VECM carries $p - 1$ lagged differences.
3. Choose the deterministic specification. Plot the data: mean-zero variables get no deterministic terms; non-zero means put a constant inside the cointegrating relation; linear trends call for constant plus trend, with the trend coefficient tested.
4. Run Johansen: eigenvalues, trace and max eigenvalue statistics, against critical values.
5. Fix the rank. Sequential tests, the eigenvalue cliff, and the economic sanity of the vectors — cross-examine all three.
6. Identify and test restrictions. Normalize $\boldsymbol{\beta}$, test theoretical restrictions with LR.
7. Estimate the full VECM and check whether the signs and sizes in $\boldsymbol{\alpha}$ make sense.
8. Diagnose. Residual autocorrelation (Portmanteau), normality, parameter stability.

Step 2 is not a formality: the cointegration rank can change with $p$; too few lags leave autocorrelated residuals, too many burn degrees of freedom in small samples. Same for step 3 — the constant and trend directly change the test’s distribution, and the bias runs both ways: including them when they should be absent, or excluding them when they should be present. A common middle case is a constant inside the cointegrating relation with no trend: $\Delta \mathbf{y}_t = \boldsymbol{\alpha}(\boldsymbol{\beta}' \mathbf{y}_{t-1} + \mu) + \sum \Gamma_j \Delta \mathbf{y}_{t-j} + \boldsymbol{\varepsilon}_t$, which is “Case 2” of Johansen (1995).

At step 5, remember the small-sample problem: the trace and max eigenvalue tests are known to over-reject in small samples, which inflates the rank. With short samples, use the Bartlett correction or re-size the critical values by simulation, and lean harder on the eigenvalue cliff and vector sanity as corroborating evidence.

A few situations break the whole procedure and are worth ruling out before starting. Structural breaks are the most treacherous: a one-time level shift is “stationary” on either side of the break, and the tests will mistake it for a cointegrating relation — the “the Fed changes its policy framework” risk flagged in the previous article’s term-structure example is exactly this; if you suspect breaks, use break-aware tests like Gregory-Hansen (1996). If some variables are I(2) — stationary only after differencing twice — the I(1) framework described here does not apply at all, and you need a dedicated I(2) cointegration framework. Monthly or quarterly data with strong seasonality should be deseasonalized first, or given seasonal dummies. And a single outlier can dominate the cointegrating vector estimate — the eigenvalue problem is built on second moments — so pre-process outliers or use robust estimation.

Side by Side with Engle-Granger

The previous article ended with the boundary of what Engle-Granger can do; now the two methods can be lined up in full. With two variables, a single cointegrating relation, and a strong theoretical prior about the normalization, Engle-Granger is entirely adequate — one OLS plus one ADF, at nearly zero cost, and it works fine as a quick screen in exploratory analysis. Three or more variables, suspected multiple cointegrating relations, formal inference on the cointegrating parameters, or structural analysis where identification matters — those are precisely the things Engle-Granger structurally cannot do, and that is Johansen territory. For multivariate systems in finance and macro (term structure, money demand, exchange-rate determination), Johansen is the default.

With two variables and one cointegrating vector, the two methods often give similar results — no surprise, they are estimating the same object. And Johansen is not a free upgrade: the ML estimates of the cointegrating vectors are biased in finite samples too (a smaller bias than Engle-Granger’s, but not zero), and the identification problem that Engle-Granger quietly buries inside its normalization choice must be answered head-on here. Putting the problem on the table is not the same as solving it for you.

References

• Engle, R. F., & Granger, C. W. J. (1987). Co-integration and error correction: Representation, estimation, and testing. Econometrica, 55(2), 251-276.
• Gregory, A. W., & Hansen, B. E. (1996). Residual-based tests for cointegration in models with regime shifts. Journal of Econometrics, 70(1), 99-126.
• Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press. Ch. 20.
• Johansen, S. (1988). Statistical analysis of cointegration vectors. Journal of Economic Dynamics and Control, 12(2-3), 231-254.
• Johansen, S. (1991). Estimation and hypothesis testing of cointegration vectors in Gaussian vector autoregressive models. Econometrica, 59(6), 1551-1580.
• Johansen, S. (1995). Likelihood-Based Inference in Cointegrated Vector Autoregressive Models. Oxford University Press.
• Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer. Ch. 6-7.
• Osterwald-Lenum, M. (1992). A note with quantiles of the asymptotic distribution of the likelihood ratio cointegration rank test statistics. Oxford Bulletin of Economics and Statistics, 54(3), 461-472.
• Park, J. Y., & Hansen, B. E. (1990). Inference for I(1) cointegrated regression models with trend. Working paper, University of Rochester.
• Tsay, R. S. (2010). Analysis of Financial Time Series (3rd ed.). Wiley. Ch. 8.