Cointegration

In-depth Articles


Cointegration

The variables of a model are said to be cointegrated if there exists a STATIONARY linear combination (remembering that the sum of stochastic processes is a stochastic process) such that the series, taken individually, are non-stationary. If this relationship exists, then it is said that there is a long-term relationship between the variables of the model. (In fact, if I consider a stationary process and an integrated one (non-stationary), the linear combination of these can result stationary if they are integrated of the same order. This can happen if the considered phenomena share the set of factors responsible for the non-stationarity.)

A multivariate time series Xt = [X1t, . , Xkt] is said to be cointegrated if each individual Xit is an I(1) time series and there exists at least one weight vector a' = [a1, . , ak] such that the time series obtained as a linear combination of a'Xt is stationary.
a = is called the cointegrating vector
The cointegrating vector is not unique, in fact every multiple of the vector a (b*a with b being any scalar) is also a cointegrating vector; in particular for k=2 it is possible to obtain a single cointegrating vector while for k>2 it is possible to obtain h cointegrating vectors (such that h>k) usually collected in a matrix A called the Cointegration matrix.
Theorem: if the matrix A' is made up of h linearly independent row vectors, A'Xt is a stationary process and c' is a vector linearly independent from the rows of A', then c' is not stationary.

To explain the cointegrated system, the metaphor of the two drunks is often used:
Two drunks leave the bar and we decide to observe the path they take to reach their homes. Let "t=1,2,." indicate a discrete measure of the time that passes from the moment the two drunks leave the bar, and let zt and xt be the distances of the two from the bar. Each step taken at time t is described by the quantities ut and vt. The first drunk has more difficulty moving around because he has one leg shorter than the other and therefore asks the other to be able to walk arm in arm.
Given the circumstances, it can be assumed that the first drunk walks following a path described by a random walk with drift: zt =z(t-1) + Φ + ut and that he pulls his friend towards him before each step, so that the starting point of the second drunk at time t is proportional to that of the first drunk, i.e., a process written as xt = γz(t-1) + vt with γ being a real constant.
The path followed by the two can be described as a VAR(1):
We have thus shown that Xt and Zt share the same non-stationary component gt (RWD) that is in I(1), known as the common trend of the system. Therefore, the system [Zt Xt]T is said to be cointegrated. ``` n <- 100 eps <- rnorm(n);eps2 <- rnorm(n) z0<- x0 <- rep(0, 100) for (i in seq.int(2, 100)){ # Drunken person with a limp z0[i] <- 1*z0[i-1] + 0.1 + eps[i] # Drunken person helping the one with a limp x0[i] <- 0.6*z0[i-1] + eps2[i] } plot(ts(z0), col=3) lines(ts(x0),lty=2, col=2) legend("topleft",c("Drunk with a limp","Drunk helping the one with a limp"), lty=c(1,2), col = c(3,2)) ```
plot of chunk unnamed-chunk-9

The goal of the Error Correction Representation is to decompose the cointegrated process into a component that describes long-term behavior (non-stationary) and one that describes short-term behavior (stationary).
Let us consider as an example a bivariate VAR(1) which derives from two RWs without drift:
X1t = γ X2t-1 + u1t + γ u2t
X2t = X2t-1 + u2t

Subtracting X1t-1 and X2t-1 respectively from these models, we obtain:

ΔX1t = -X1t-1 + γ X2t-1 + u1t + γ u2t
ΔX2t = u2t

Π = the matrix of the coefficients of the model constructed with first differences.
Performing the calculations using matrix forms, we obtain:

[ΔX1t ΔX2t]T = [-1 0]T (X1t-1-γ X2t-1) + [u1t + γ u2t u2t]T <- ERROR CORRECTION REPRESENTATION

where:
In this case, the rank of matrix Π is easy to calculate and turns out to be 1, but what does this mean?
In our specific case, the series is therefore cointegrated.
The matrix Π can also be thought of as the product of the equilibrium correction matrix Β with the Cointegration matrix AT (amplitude of the deviation from the cointegration relationship at time t-1).

The Johansen technique is used to estimate the matrix Π through the estimation of a Cointegration matrix AT and the matrix of equilibrium correction weights. This method can be used in the case where the Cointegration rank is known. This method is divided into two steps: in the first step, the equilibrium correction matrix Β is estimated using the least squares method; in the second, the Cointegration matrix AT is estimated using maximum likelihood estimation.
Specifically, we can distinguish three cases:

CONCENTRATED REGRESSION - both Cointegration rank and Cointegration matrix are known

Assuming h and A are known, first verify using the ADF test that the k variables of Xt are all I(1); then, if they are, calculate A'Xt and apply the same tests to verify the stationarity of the h univariate processes that compose it; if the h tests lead to the rejection of the null hypothesis of the existence of a unit root, then the hypothesis that A is the Cointegration matrix can be accepted.
Starting the estimation, we consider the 2 auxiliary regression models to estimate the link between ΔXt and Xt-1 net of other variables using the Frisch-Waugh theorem:

ΔXt = Π0 + Π1 ΔXt-1 + ... + Πp-1 ΔXt-p+1 + ut

Xt-1 = θ0 + θ1 ΔXt-1 + ... + θp-1 ΔXt-p+1 + vt

ut contains all the information on ΔXt net of the short-term component, while vt contains everything that Xt-1 not explained by the short-term component.
By performing a regression between the residuals, we find the estimate we are looking for; this model is known as the concentrated regression model:
ut = Π vt + εt =
For convenience, we consider:
(ut)T = (AT vt)T ΒT + (εt)T

then we can estimate ΒT using the least squares method, resulting in:

estimated ΒT = [(AT vt)(AT vt)T]-1 (AT vt) ut

which can also be rewritten as:
= (AT vt vtT A)-1 AT vt ut =
where:
Σvv = vt vtT
Σvu = vt ut


Therefore:
estimated Β = (AT Σvv A)-1 A Σvu

ESTIMATION OF A - the Cointegration rank is known, but the Cointegration matrix is not

In the case where matrix A is unknown, we must distinguish between two scenarios: the first occurs when the maximum Cointegration rank k=2 (i.e., a bivariate process), in which case the Engle-Granger technique can be used; whereas if the maximum rank is greater than 2, we are forced to use the Johansen technique.


Maximizing this function, we find the estimated matrix A; now that matrix A is known, we obtain the estimate of Β from the least squares formula derived in step 1, and thus the estimate of Π0=Β AT is equal to:

Π0 = Σvu' A AT



COINTEGRATION RANK VERIFICATION - neither the Cointegration rank nor the Cointegration matrix is known

The estimation reasoning to follow in this situation is the same as in the previous point, except that we must initially estimate the Cointegration rank using tests based on the likelihood ratio, which can be:

5 tabulated cases for Cointegration rank verification

In the tables created for the two tests used for verifying the Cointegration rank, 5 cases are distinguished based on the presence of a deterministic component. Let us hypothesize a generic model:

ΔXt = Β AT Xt + μ0 + μ1 t + εt
s.t. μ0 = a0Β + γ0 ; μ1 = a1 Β + γ1


Then the 5 cases that we can distinguish are:
CASE 1:
CASE 2:
CASE 3:
CASE 4:
CASE 5:

### VECM model and Cointegration
    # Let's revisit the metaphor of the two drunks
    n <- 100
    eps <- rnorm(n);eps2 <- rnorm(n)
    z0<- x0 <- rep(0, 100)
    for (i in seq.int(2, 100)){
      # Drunken person with a limp
      z0[i] <- 1*z0[i-1] + 0.1 + eps[i]
      # Drunken person helping the one with a limp
      x0[i] <- 0.6*z0[i-1] + eps2[i]
    }
    dati<-as.data.frame(cbind(z0,x0))

    # Johansen Test (trace)
    vecm <- ca.jo(dati, type = "trace", ecdet = "none", K =2)
    vec2var(vecm,r=1)

##
## Coefficient matrix of lagged endogenous variables:
##
## A1:
##        z0.l1        x0.l1
## z0 0.9799309 -0.093289080
## x0 0.5410979  0.003056476
##
##
## A2:
##        z0.l2       x0.l2
## z0 0.1916015 -0.19661290
## x0 0.0766579 -0.04710763
##
##
## Coefficient matrix of deterministic regressor(s).
##
##     constant
## z0 0.1953953
## x0 0.0612081

    summary(vecm)

##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.4194067418 0.0006979154
##
## Values of teststatistic and critical values of test:
##
##           test 10pct  5pct  1pct
## r <= 1 |  0.07  6.50  8.18 11.65
## r = 0  | 53.35 15.66 17.95 23.52
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
##           z0.l2     x0.l2
## z0.l2  1.000000 1.0000000
## x0.l2 -1.690071 0.2712602
##
## Weights W:
## (This is the loading matrix)
##
##          z0.l2        x0.l2
## z0.d 0.1715324 -0.003432119
## x0.d 0.6177558  0.001062155

    # Johansen Test (max eigenvalue)
    vecm2 <- ca.jo(dati, type = "eigen", ecdet = "none", K =2)
    summary(vecm2)

##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.4194067418 0.0006979154
##
## Values of teststatistic and critical values of test:
##
##           test 10pct  5pct  1pct
## r <= 1 |  0.07  6.50  8.18 11.65
## r = 0  | 53.28 12.91 14.90 19.19
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
##           z0.l2     x0.l2
## z0.l2  1.000000 1.0000000
## x0.l2 -1.690071 0.2712602
##
## Weights W:
## (This is the loading matrix)
##
##          z0.l2        x0.l2
## z0.d 0.1715324 -0.003432119
## x0.d 0.6177558  0.001062155

    s = 1.000*z0 - 1.663959*x0 #Cointegration link
    plot(s, type="l")
plot of chunk unnamed-chunk-10
    adfTest(s,type="nc",lags=1)

##
## Title:
##  Augmented Dickey-Fuller Test
##
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -6.0945
##   P VALUE:
##     0.01
##
## Description:
##  Fri Sep 07 15:44:38 2018 by user: gieck

    vecm.r1 <- cajorls(vecm, r = 1)
    alpha <- coef(vecm.r1$rlm)[1,]
    beta <- vecm.r1$beta


    # 5 Cointegration cases
    data(Canada)
    # Case 1: No intercept/trend (not directly supported in R's ca.jo default)
    # Case 2: No intercept in the cointegrating equation
    vecm2 <- ca.jo(Canada, type = "trace", ecdet = "none", K =3)
    # Case 3: Constant restricted to the cointegrating equation
    vecm3 <- ca.jo(Canada, type = "trace", ecdet = "const", K =3)
    # Case 4: Trend restricted to the cointegrating equation
    vecm4 <- ca.jo(Canada, type = "trace", ecdet = "trend", K =3)
    # Case 5: Unrestricted constant and trend (not directly supported in R)
    summary(vecm2)

##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.4178095906 0.1828654502 0.1239725052 0.0007523871
##
## Values of teststatistic and critical values of test:
##
##           test 10pct  5pct  1pct
## r <= 3 |  0.06  6.50  8.18 11.65
## r <= 2 | 10.78 15.66 17.95 23.52
## r <= 1 | 27.14 28.71 31.52 37.22
## r = 0  | 70.96 45.23 48.28 55.43
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
##               e.l3    prod.l3      rw.l3       U.l3
## e.l3     1.0000000  1.0000000  1.0000000  1.0000000
## prod.l3  0.6409112  0.2993239 -0.5765703 -0.1674879
## rw.l3   -0.6996053 -0.3775646 -0.3098725 -0.3210560
## U.l3     0.5506937  2.7029914  1.7800952  1.5046076
##
## Weights W:
## (This is the loading matrix)
##
##               e.l3    prod.l3       rw.l3         U.l3
## e.d    0.025577740  0.1542901 -0.01923271 -0.004492451
## prod.d 0.001322905  0.1470512  0.08043698  0.015693650
## rw.d   0.095275595 -0.1653593  0.13744785 -0.001148642
## U.d    0.007709430 -0.1220844 -0.02034877  0.005043952

    summary(vecm3)

##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1]  5.619084e-01  2.067893e-01  1.241824e-01  5.502064e-02 -8.115665e-15
##
## Values of teststatistic and critical values of test:
##
##            test 10pct  5pct  1pct
## r <= 3 |   4.58  7.52  9.24 12.97
## r <= 2 |  15.32 17.85 19.96 24.60
## r <= 1 |  34.09 32.00 34.91 41.07
## r = 0  | 100.94 49.65 53.12 60.16
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
##                  e.l3      prod.l3        rw.l3         U.l3     constant
## e.l3        1.0000000    1.0000000    1.0000000    1.0000000    1.0000000
## prod.l3     0.5175539    0.2784319   -0.5024247    0.5680653   -0.1357442
## rw.l3      -0.5809564   -0.3514234   -0.3185323   -0.4769518   -0.3283532
## U.l3        1.3287436    2.8794520    1.8360495    2.5158096    1.5437606
## constant -900.1428894 -930.1029338 -616.0750504 -989.6170094 -758.5196674
##
## Weights W:
## (This is the loading matrix)
##
##              e.l3     prod.l3       rw.l3         U.l3      constant
## e.d    0.04212272  0.10050067 -0.01583762  0.029356910  4.788044e-11
## prod.d 0.04600305  0.17137682  0.08507953 -0.057954672  7.044160e-11
## rw.d   0.11985585 -0.20631119  0.14834704  0.004323796  6.291878e-13
## U.d    0.01106869 -0.08943698 -0.02578978 -0.025521676 -4.083406e-11

    summary(vecm4)

##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 4.505013e-01 1.962777e-01 1.676668e-01 4.647108e-02 1.896912e-16
##
## Values of teststatistic and critical values of test:
##
##           test 10pct  5pct  1pct
## r <= 3 |  3.85 10.49 12.25 16.26
## r <= 2 | 18.72 22.76 25.32 30.45
## r <= 1 | 36.42 39.06 42.44 48.45
## r = 0  | 84.92 59.14 62.99 70.05
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
##                e.l3    prod.l3      rw.l3        U.l3    trend.l3
## e.l3        1.00000  1.0000000  1.0000000  1.00000000  1.00000000
## prod.l3   -41.92624  0.7723987 -0.3468422  0.23574849 -0.12062153
## rw.l3     -76.94645 -0.2572818 -0.4848433 -0.03286331 -0.29511866
## U.l3     -132.85358  2.6289928  2.5757023  1.61283935  1.51474627
## trend.l3   54.56956 -0.1778682  0.1766543 -0.37545112 -0.03414891
##
## Weights W:
## (This is the loading matrix)
##
##                e.l3     prod.l3       rw.l3        U.l3      trend.l3
## e.d    0.0002028169  0.14777357  0.02351692  0.01691655 -3.920432e-12
## prod.d 0.0001558757 -0.03577745  0.20463284 -0.06258514  2.173585e-11
## rw.d   0.0011022537 -0.18875801  0.04885538  0.05046075  2.731205e-11
## U.d    0.0001125446 -0.07967784 -0.05973735 -0.02781459  1.643948e-12

Python in Practice

Below we demonstrate cointegration concepts with simulated data in Python.

1. Simulating Cointegrated Series

import numpy as np
from statsmodels.tsa.stattools import adfuller, coint

np.random.seed(42)
n = 300

# Common stochastic trend (non-stationary)
trend = np.cumsum(np.random.normal(0, 1, n))

# Two series sharing the same trend (cointegrated)
y1 = trend + np.random.normal(0, 0.5, n)
y2 = 0.7 * trend + np.random.normal(0, 0.5, n)

# Each is I(1) individually
print("ADF on Y1 (p-value):", round(adfuller(y1)[1], 4))
print("ADF on Y2 (p-value):", round(adfuller(y2)[1], 4))

# Output:
# ADF on Y1 (p-value): 0.7832
# ADF on Y2 (p-value): 0.6914
# Both non-stationary (cannot reject unit root)

2. Engle-Granger Cointegration Test

# Test for cointegration using Engle-Granger method
t_stat, p_value, crit_values = coint(y1, y2)
print(f"Engle-Granger test statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Cointegrated: {'Yes' if p_value < 0.05 else 'No'}")

# The cointegrating residual should be stationary
spread = y1 - (1/0.7) * y2
adf_spread = adfuller(spread)
print(f"\nADF on spread (p-value): {adf_spread[1]:.4f}")
print(f"Spread is stationary: {'Yes' if adf_spread[1] < 0.05 else 'No'}")

# Output:
# Engle-Granger test statistic: -14.2176
# p-value: 0.0000
# Cointegrated: Yes
#
# ADF on spread (p-value): 0.0000
# Spread is stationary: Yes

3. Spurious Regression (Non-Cointegrated)

# Two independent random walks
rw1 = np.cumsum(np.random.normal(0, 1, n))
rw2 = np.cumsum(np.random.normal(0, 1, n))

t_stat2, p_value2, _ = coint(rw1, rw2)
print(f"Independent RWs - Engle-Granger p-value: {p_value2:.4f}")
print(f"Cointegrated: {'Yes' if p_value2 < 0.05 else 'No'}")

# Output:
# Independent RWs - Engle-Granger p-value: 0.4521
# Cointegrated: No

Results

>