Cointegration

This page states the mathematics of cointegration that the repository relies on: the definitions of integrated and cointegrated series, the exact representations, the Gaussian rank inference of the vector error-correction model, the exact theory of a mean-reverting relation with a known cointegrating vector, the filtering theory for time-varying activity, the high-dimensional spectral laws, and the optimal stopping of a mean-reverting spread. Results are stated as results, with the regime in which they hold and the reductions that validate them.

The open questions of the research programme that uses this material, none of which has a proof in the literature as stated, are in the temporal-cointegration plan, and the precise statements that must be proven to close them are on Cointegration: Proof Obligations. This page contains only established mathematics and does not restate them.

Throughout, \(p_t \in \mathbb{R}^N\) is the vector of log-prices, \(N\) is the dimension, \(T\) the number of observations, and \(\gamma = N/T\) the aspect ratio. The Ornstein-Uhlenbeck process used below is the same object defined on Stochastic Processes; the cointegrating residual is an instance of it, with the symbols of that page.

The model assumptions are Gaussian innovations, parameters that are constant within a regime, a first-order finite-state activity chain, and log-prices that are (I(1)). The technical assumptions are finite innovation variance, an irreducible and ergodic chain wherever the filter is used, and finite moments wherever a spectral limit is invoked. Assumptions specific to a single result are stated with that result. The symbol (\Delta) denotes the sampling interval, whereas (\Delta x_t = x_t - x_{t-1}) denotes the first difference; the two are distinguished by the presence of an operand.

Integrated series and the cointegrating representation

A scalar process \(x_t\) is \(I(0)\) if it has a causal linear representation \(x_t = \sum_{j \ge 0} \psi_j \varepsilon_{t-j}\) with \(\varepsilon_t\) independent, mean zero and finite variance, and \(\sum_{j \ge 0} \lvert \psi_j \rvert < \infty\); equivalently, its partial sums converge under the functional limit

\[ T^{-1/2} \sum_{t=1}^{\lfloor Tr \rfloor} x_t \;\Rightarrow\; \sigma W(r), \qquad 0 \le r \le 1, \]

where \(W\) is standard Brownian motion and \(\Rightarrow\) denotes weak convergence in \(D[0,1]\). A process is \(I(1)\) if \(\Delta x_t := x_t - x_{t-1}\) is \(I(0)\) and \(x_t\) itself is not \(I(0)\); then \(T^{-1/2} x_{\lfloor Tr \rfloor} \Rightarrow \sigma W(r)\), and \(\sigma = \psi(1) \sigma_\varepsilon\) is the long-run standard deviation.

The vector \(p_t = (p_{1t}, \dots, p_{Nt})^\top\) is cointegrated of rank \(r\) if every component is \(I(1)\) and there exist \(r\) linearly independent vectors collected in an \(N \times r\) matrix \(\beta\) such that \(\beta^\top p_t\) is \(I(0)\). The columns of \(\beta\) are the cointegrating vectors, \(\mathrm{span}(\beta)\) is the cointegrating space, and \(r\) is the rank. If \(r = 0\) the system has no cointegrating relation; if \(r = N\) every component is \(I(0)\).

The representation that all of the following uses is the exact coordinate identity. Choose any complement \(\beta_\perp \in \mathbb{R}^{N \times (N-r)}\) with \(\beta^\top \beta_\perp = 0\). Then \(\beta\) and \(\beta_\perp\) together span \(\mathbb{R}^N\), and every state decomposes into a stationary part and a trend part,

\[ p_t = \beta (\beta^\top \beta)^{-1} z_t + \beta_\perp (\beta_\perp^\top \beta_\perp)^{-1} f_t, \qquad z_t = \beta^\top p_t, \qquad f_t = \beta_\perp^\top p_t . \]

Substituting \(p_t\) into \(\beta^\top p_t\) gives \(\beta^\top p_t = z_t\) exactly, because \(\beta^\top \beta_\perp = 0\) kills the trend term. The identity is therefore exact finite-dimensional linear algebra and involves no approximation: choosing the dynamics of the stationary coordinates \(z_t\) and the trend coordinates \(f_t\) determines \(p_t\) with exact cointegration by construction.

The converse direction is the Granger representation theorem. If \(p_t\) is \(I(1)\) with cointegration rank \(r\), then it has a vector error-correction representation

\[ \Delta p_t = \alpha \beta^\top p_{t-1} + \sum_{i=1}^{k-1} \Gamma_i \Delta p_{t-i} + \varepsilon_t, \qquad \Pi := \alpha \beta^\top, \]

with \(\alpha, \beta \in \mathbb{R}^{N \times r}\) of full column rank and \(\alpha_\perp^\top \Gamma \beta_\perp\) invertible, where \(\Gamma = I_N - \sum_{i=1}^{k-1} \Gamma_i\) and \(\alpha_\perp, \beta_\perp\) are complements (Engle and Granger, 1987). The matrix \(\Pi\) has rank exactly \(r\): it is the only channel through which the levels \(p_{t-1}\) enter, and its rank is the rank of the cointegrating space.

Gaussian inference and the rank test

For inference assume the innovations are Gaussian and the dynamics are contiguous over the sample, so no parameter changes within \([0,T]\). The model is

\[ \Delta p_t = \Pi p_{t-1} + \sum_{i=1}^{k-1} \Gamma_i \Delta p_{t-i} + \varepsilon_t, \qquad \varepsilon_t \sim \mathcal{N}(0, \Omega), \]

with fixed \(N\) and \(T \to \infty\). The Gaussian likelihood is exact, and maximizing it over \(\Pi = \alpha \beta^\top\) of rank \(r\) is a reduced-rank regression (Johansen, 1988, 1991). Reduce \(\Delta p_t\) and \(p_{t-1}\) on the lagged differences \((1, \Delta p_{t-1}, \dots, \Delta p_{t-k+1})\) by least squares to obtain residuals \(R_{0t}\) and \(R_{1t}\), and set

\[ S_{ij} = \frac{1}{T} \sum_{t=1}^{T} R_{it} R_{jt}^\top, \qquad i, j \in \lbrace 0, 1 \rbrace . \]

The profile likelihood depends on \(\beta\) only through the determinant \(\lvert \beta^\top S_{11} \beta \rvert / \lvert \beta^\top \beta \rvert\), so the maximizer solves the generalized eigenvalue problem

\[ \det\!\left( \lambda S_{11} - S_{10} S_{00}^{-1} S_{01} \right) = 0 . \]

Order the eigenvalues \(\hat\lambda_1 \ge \dots \ge \hat\lambda_N\). The maximum likelihood estimate of \(\beta\) is the matrix of the \(r\) eigenvectors belonging to the \(r\) largest eigenvalues, normalized by \(\hat\beta^\top S_{11} \hat\beta = I\), and \(\hat\alpha = S_{01}\hat\beta\). This estimator is an explicit algebraic functional of the data; it is closed form in the sense of a finite eigenproblem, not an iterative search.

The likelihood-ratio statistic for the null hypothesis of rank at most \(r\) is

\[ \mathrm{LR}(r) = -T \sum_{i=r+1}^{N} \ln\!\left( 1 - \hat\lambda_i \right) . \]

Under the null its limit is not the chi-square law, because the regressor \(p_{t-1}\) is \(I(1)\) and the estimation error of \(\beta\) does not vanish fast enough. It is the trace of a functional of an \((N-r)\)-dimensional standard Brownian motion,

\[ \mathrm{LR}(r) \;\Rightarrow\; \operatorname{tr}\!\lbrace \left( \int_0^1 W \, dW^\top \right) \left( \int_0^1 W W^\top \, du \right)^{-1} \left( \int_0^1 dW \, W^\top \right) \rbrace, \]

whose percentiles are tabulated (Johansen, 1991). Two features matter for the present problem and are used repeatedly below: the distribution depends on \(N\), and it is derived under fixed \(N\) with \(T \to \infty\), so it does not apply when \(N\) and \(T\) grow together.

Exact theory of a single mean-reverting relation

If \(\beta\) is known, the relation is observed directly and the problem becomes univariate. Model the residual \(j\) as an Ornstein-Uhlenbeck process with speed \(\theta_j > 0\), long-run mean \(\mu_j\), and volatility \(\sigma_j > 0\),

\[ d z_{j,t} = \theta_j (\mu_j - z_{j,t})\, dt + \sigma_j\, dW_{j,t}, \]

the same process as on Stochastic Processes. Its stationary law and the relaxation time constant are

\[ z_j(\infty) \sim \mathcal{N}\!\left( \mu_j, \frac{\sigma_j^2}{2\theta_j} \right), \qquad h_j = \frac{\ln 2}{\theta_j}, \]

where \(h_j\) is the half-life: the time for the expected deviation from \(\mu_j\) to halve. The stationary standard deviation is \(\varsigma_j = \sigma_j / \sqrt{2\theta_j}\), and it is the natural scale of a deviation and of any trading threshold.

The exact discrete law at sampling interval \(\Delta\) follows from the integrating factor. Solving the linear SDE over one step,

\[ z_{j,t+\Delta} = \mu_j + (z_{j,t} - \mu_j) e^{-\theta_j \Delta}

• \sigma_j \int_0^{\Delta} e^{-\theta_j (\Delta - s)} \, dW_{j,s}, \]

and the Ito isometry evaluates the variance of the stochastic integral as \(\sigma_j^2 (1 - e^{-2\theta_j \Delta}) / (2\theta_j)\). Hence the sampled process is exactly an AR(1),

\[ z_{j,t+\Delta} = \mu_j (1 - \phi_j) + \phi_j z_{j,t} + \epsilon_{j,t}, \qquad \phi_j = e^{-\theta_j \Delta}, \qquad \operatorname{Var}(\epsilon_{j,t}) = \frac{\sigma_j^2}{2\theta_j} \left( 1 - \phi_j^2 \right), \]

with \(\epsilon_{j,t}\) independent and Gaussian. Everything here is exact, not an Euler approximation: the Euler scheme would replace \(\phi_j\) by \(1 - \theta_j \Delta\) and is only first order.

Two consequences of the exact law are used as design criteria. First, the resolvability condition: as \(\theta_j \Delta \to 0\) the autoregressive root satisfies \(\phi_j \to 1\), so the sampled residual converges to a random walk and cannot be distinguished from one by any test of fixed size. A non-vanishing \(\theta_j \Delta\) is therefore necessary for the relation to be identifiable at the sampling interval. Second, the observability condition: the expected deviation from the mean after a lifespan \(D\) is \(e^{-\theta_j D} = 2^{-D/h_j}\) of its initial value, so observing the relaxation to within a fraction \(\rho\) of the initial deviation requires \(D / h_j \ge \log_2(1/\rho)\). For example, \(\rho = 0.1\) requires \(D/h_j \gtrsim 3.3\), and \(\rho = 0.03\) requires \(D/h_j \gtrsim 5\).

The Gaussian estimator of \(\theta_j\) is a deterministic function of the AR(1) estimate. Conditional on the first observation, the least-squares estimator of \(\phi_j\) is the demeaned autocorrelation

\[ \hat\phi_j = \frac{\sum_{t} (z_{j,t} - \bar z_j)(z_{j,t+\Delta} - \bar z_j)}{\sum_{t} (z_{j,t} - \bar z_j)^2}, \qquad \hat\theta_j = -\frac{\ln \hat\phi_j}{\Delta}, \]

and \(\hat\phi_j\) is downward biased in finite samples, of order \(1/T\) (Kendall, 1954; Stambaugh, 1999). The bias inflates the estimated speed and understates the half-life; it is a property of the estimator, reproduced by the model, and not a numerical defect.

Time-varying activity

Cointegration that holds only on part of the sample is modelled by switching the adjustment channel, not the cointegrating space. Let \(s_{j,t} \in \lbrace 0, 1 \rbrace\) be the activity state of relation \(j\), with

\[ d z_{j,t} = s_{j,t}\, \theta_j (\mu_j - z_{j,t})\, dt + \sigma_j\, dW_{j,t}. \]

When \(s_{j,t} = 1\) the residual mean-reverts; when \(s_{j,t} = 0\) it is a driftless random walk. This is exactly the \(\alpha_j \to 0\) limit of the error-correction representation on the inactive intervals, so the local rank of the cointegrating space equals the number of active relations. The switching capability defines the class; estimating it uses the filter below.

Take the state to be a first-order Markov chain on a finite set with transition matrix \(P\), \(P_{ij} = \Pr(s_t = j \mid s_{t-1} = i)\), and let \(\eta_{t,j} = f(y_t \mid s_t = j)\) be the density of the observation under state \(j\). The Hamilton (1989) filter computes the exact state distribution by the forward recursion

\[ \xi_{t \mid t} = \frac{\xi_{t \mid t-1} \odot \eta_t}{\mathbf{1}^\top (\xi_{t \mid t-1} \odot \eta_t)}, \qquad \xi_{t+1 \mid t} = P^\top \xi_{t \mid t}, \]

where \(\odot\) is elementwise multiplication and \(\mathbf{1}\) is the vector of ones. The log-likelihood is exact and equals the sum of the normalizing constants,

\[ \ell(\psi) = \sum_{t=1}^{T} \ln \mathbf{1}^\top \!\left( \xi_{t \mid t-1} \odot \eta_t \right), \]

and the exact regime posterior is the backward-smoothed distribution \(\xi_{t \mid T}\). The objective and the posterior are exact; the maximizer over the parameters \(\psi\) is not available in closed form, so estimation is a finite-dimensional numerical optimization of an exactly computable criterion. This is the reason the switching model is classified below as exact but not closed form.

If the regime path is known, for instance from a deterministic schedule, the model is conditionally Gaussian with a known regressor partition, and the exact conditional maximum likelihood estimator is block least squares: ordinary least squares applied separately to each regime. No filter is needed, and this is the reference estimator the switching filter must reproduce.

For a single change at an unknown time, the break fraction estimator converges at rate \(T\) to an argmax of a two-sided Brownian motion with a parabolic drift (the Chernoff distribution), and the supremum of the likelihood-ratio statistic over the change date converges to a known functional of Brownian motion (Andrews, 1993; Bai, 1997). For a cointegrated system, the tests that separate a change in the adjustment \(\alpha\), in the cointegrating space \(\beta\), in the rank, and in the short-run dynamics are exactly constructed with tabulated limits (Hansen, 2003). These are the tools that distinguish the three meanings of temporary cointegration.

High-dimensional spectral theory

When \(N\) is of order \(10^3\) and \(\gamma = N/T\) is of order one, the eigenvalues of a sample covariance are governed by the Wishart ensemble. If \(X\) has \(N \times T\) independent entries and \(S = T^{-1} X X^\top\), the eigenvalues of \(S\) follow the Marchenko-Pastur law

\[ \rho_{\mathrm{MP}}(d\lambda) = \frac{\sqrt{(\lambda_+ - \lambda)(\lambda - \lambda_-)}}{2\pi \gamma \lambda} \quad \text{on } [\lambda_-, \lambda_+], \qquad \lambda_\pm = (1 \pm \sqrt{\gamma})^2 , \]

concentrated on \([\lambda_-, \lambda_+]\), with an atom at the origin when \(\gamma > 1\) (Marchenko and Pastur, 1967). The largest eigenvalue fluctuates around the upper edge \(\lambda_+\) on the scale \(T^{-2/3}\) and follows the Tracy-Widom law: \(\beta = 1\) for real entries and \(\beta = 2\) for complex entries (Tracy and Widom, 1994, 1996; Johnstone, 2001). The joint eigenvalue density is the Laguerre ensemble density, which is exact at finite \(N\), so a test built on the eigenvalues can be calibrated in finite samples and not only asymptotically.

A cointegrating relation generates a separated eigenvalue, which is the spiked covariance model. For a population covariance \(I + \ell\, v v^\top\) and aspect ratio \(\gamma\), the Baik-Ben Arous-Peche transition is sharp at

\[ \ell_c = \sqrt{\gamma} : \]

for \(\ell < \ell_c\) the spiked eigenvalue does not separate from the bulk and the sample eigenvector is asymptotically orthogonal to \(v\), whereas for \(\ell > \ell_c\) the outlier sits at \((1+\ell)(1+\gamma/\ell)\) and the squared overlap of the sample and population eigenvectors converges to \((1 - \gamma/\ell^2)/(1 + \gamma/\ell)\) (Baik, Ben Arous and Peche, 2005; Benaych-Georges and Nadakuditi, 2011). The threshold is the detection floor of the rank problem: below it no consistent eigenvector is available, however large \(T\) is.

The sample covariance is ill-conditioned at these aspect ratios, so it is replaced by a shrinkage estimator before any eigen-decomposition. The optimal linear shrinkage toward a scaled identity,

\[ \hat\Sigma = (1 - \delta) S + \delta\, \frac{\operatorname{tr}(S)}{N} I, \]

has a closed-form, Marchenko-Pastur-consistent intensity \(\delta\) (Ledoit and Wolf, 2004).

A structural caveat conditions the transfer of all of the above to cointegration. The spectral results assume independent or weakly mixing entries with a growing sample. The matrices of the rank problem are built from the levels \(p_{t-1}\), which are \(I(1)\) and therefore non-mixing: they are functionals of matrix Brownian motion, not a Wishart of independent rows. Replacing the Wishart thresholds by Marchenko-Pastur or Baik-Ben Arous-Peche thresholds in the rank problem is therefore not justified by the results as stated. What the correct threshold is in that case is open and is stated as such in the plan.

Optimal stopping of a mean-reverting spread

Once \(\beta\) is fixed, trading a detected relation is an optimal stopping problem for the spread. Centre the residual, \(x = z - \mu\), so that with a constant cost \(c\) per round trip the payoff of an entry at level \(x\) and exit at the mean is the deviation measured net of cost. The state follows

\[ d x_t = -\theta x_t\, dt + \sigma\, dW_t, \]

with infinitesimal generator \(\mathcal{L} = \tfrac{1}{2} \sigma^2 \partial_{xx} - \theta x \partial_x\). On a continuation region the value function satisfies the homogeneous equation \(\mathcal{L} u = 0\), which integrates directly. Writing \(\mathcal{L} u = 0\) as \(u'' / u' = 2\theta x / \sigma^2\) and integrating once gives \(u'(x) = c_1 \exp(\theta x^2 / \sigma^2)\), hence

\[ u(x) = c_0 + c_1 \int_0^{x} \exp\!\left( \frac{\theta s^2}{\sigma^2} \right) ds . \]

The integral is the Dawson function, equivalently expressible through the parabolic cylinder function or the confluent hypergeometric \(M(1/2, 3/2, \theta x^2/\sigma^2)\). The optimal thresholds therefore appear as the unknowns of a system of value-matching and smooth-pasting equations with right-hand sides in this family, and the thresholds are characterized by transcendental equations rather than an elementary closed form. The full transition density of the free process is the exact Gaussian kernel

\[ p(\Delta, x, y) = \frac{1}{\sqrt{2\pi v(\Delta)}} \exp\!\left( -\frac{(y - m(\Delta, x))^2}{2 v(\Delta)} \right), \qquad m = x e^{-\theta \Delta}, \quad v = \frac{\sigma^2}{2\theta} \left( 1 - e^{-2\theta \Delta} \right), \]

which is the same discrete law used in the previous section. With killing at two barriers the transition density has a spectral representation in the Ornstein-Uhlenbeck eigenfunctions, which are the Hermite functions with eigenvalues \(-n\theta\); this representation is what would make a tick-quantized observation model exactly solvable, and is stated here only as the structure of the known result.

Assumptions and rigor

The results above are classified by the project rigor levels. The model and technical assumptions are stated at the top of the page.

Validation targets

The results on this page define the checks that an implementation must pass, and the reductions that show each check is the right one. No implementation exists yet; these are targets, not passed tests.

• With \(r = 0\), the system is \(I(1)\) with no relation, and the detector must reject at its nominal size. This is the null against which the spectral thresholds, including the Ledoit-Wolf regularized ones, are calibrated.
• With \(N = 2\) and \(r = 1\), the reduced-rank regression has a hand-computable generalized eigenvalue problem, and the estimator can be checked in closed form.
• Setting \(s_{j,t} \equiv 1\) reduces the switching model to the contiguous model, and the filter likelihood reduces to the Gaussian likelihood; the two must agree.
• As \(\theta_j \Delta \to 0\) the exact AR(1) law converges to a random walk, which is the resolvability failure and must appear as a loss of power.
• As \(\sigma_j \to 0\) the residual follows the deterministic path \(\mu_j + (z_0 - \mu_j)e^{-\theta_j t}\), which fixes the half-life reduction.
• The stopping thresholds reduce to the known mean-reversion band as \(c \to 0\) and reproduce the exact Gaussian transition kernel above; the two-barrier spectral density must match the free kernel when the barriers are removed.