Cointegration in Python: ARDL bounds test, Engle–Granger and Johansen

Two series can each wander like a random walk and still be tied together, as consumption and income are. If some linear combination of them is stationary, they are cointegrated. A regression in levels then estimates a genuine long-run relation, and deviations from that relation help predict how the series move next. Without cointegration the same regression is spurious: it can show a large R-squared and significant t statistics for series that are unrelated. Cointegration tests tell the two cases apart, and autoregressive distributed lag (ARDL) models estimate the long-run relation and the short-run dynamics in one equation.

What cointegration tests and ARDL models estimate

Suppose y_t and x_t are integrated of order one, I(1), so their first differences are stationary (see unit-root tests). They are cointegrated if y_t − θx_t is stationary for some θ. Three approaches are common.

• Engle–Granger. Engle and Granger (1987) proposed two steps: regress y on x by OLS, then test the residuals for a unit root with a Dickey–Fuller regression. Because the residuals are estimated, the usual Dickey–Fuller critical values do not apply; the test has its own tables. It handles one cointegrating relation and depends on which variable is put on the left.
• Johansen. Johansen (1988) works with the whole system as a vector autoregression and tests the number r of independent cointegrating relations with trace and maximum-eigenvalue statistics. Their distributions are non-standard; Osterwald-Lenum (1992) tabulates the critical values. Johansen's method can find more than one relation and leads directly to a vector error-correction model.
• ARDL and the bounds test. An ARDL(p, q) model regresses y on p of its own lags and on x with lags 0 to q. Rewritten in error-correction form, it separates the long run from the short run:

Δy_t = c − a (y_{t−1} − θ'x_t) + Σ ψ_i Δy_{t−i} + Σ ω_l Δx_{t−l} + e_t

a is the speed of adjustment, the share of last period's deviation from the long-run relation that is corrected this period, and θ holds the long-run coefficients. Pesaran, Shin and Smith (2001) test for a level relationship with an F test that the level terms are zero, and a t test of a. The critical values come as bounds: the lower bound applies when all regressors are I(0), the upper bound when all are I(1). A statistic beyond the upper bound rejects "no level relationship" whatever the order of integration of the regressors; one inside the lower bound does not reject; in between, the test is inconclusive. The approach is a single-equation, conditional model: Pesaran, Shin and Smith assume that the level of y does not feed back into x.

As a rule of thumb, Engle–Granger suits a pair of I(1) series, Johansen a system with possibly several relations, and ARDL a single equation whose regressors may mix I(0) and I(1) series.

In OpenEconometrics

All three functions take a time= column (optional; without it the row order is the time order). Rows are sorted by it, integer periods must be consecutive and datetime values are taken as consecutive periods. Lags need consecutive periods, so interior gaps are refused.

• oe.egranger(df, "y", ["x"], time="t", lags=1) corresponds to egranger y x, lags(1), Schaffer's user-written Stata command. trend= is "none", "constant" (default), "trend" or "quadratic", and lags= is an integer or "aic", "bic" or "t". The result holds the test, with MacKinnon's approximate p-value and critical values, and the cointegrating regression (attrs["cointegrating_vector"]); ecm=True adds an error-correction regression.
• oe.vecrank(data=df, y=["y", "x"], time="t", lags=2) corresponds to vecrank y x, lags(2). lags counts lags of the VAR in levels and trend= takes the five cases "none", "rconstant", "constant" (default), "rtrend" and "trend". Each row of the table is a rank, with the trace and maximum-eigenvalue statistics and their 5% and 1% critical values. attrs["selected_rank"] is the first rank whose trace statistic does not exceed its 5% critical value.
• oe.ardl(data=df, y="y", x=["x"], time="t", maxlags=4, ec=True) corresponds to ardl y x, maxlags(4) aic ec followed by estat ectest, using Kripfganz and Schneider's user-written ardl command. Fixed orders go in lags=[p, q_1, ...]; otherwise every combination up to maxlags is fitted on a common sample and the one with the lowest ic="aic" (the default; the Stata command defaults to BIC) or "bic" is reported. ec=True reports the error-correction form, with the speed of adjustment as ADJ:L.y, the long-run coefficients as LR: terms with delta-method standard errors, and the short-run terms as SR:. trend= and restricted= select the Pesaran–Shin–Smith case I to V, and exog= adds regressors that enter only contemporaneously and stay out of the long run. fit.tests["bounds_f"] and fit.tests["bounds_t"] carry the I(0) and I(1) bounds at 10%, 5% and 1% and a decision at each level.

The ARDL guide, the Engle–Granger section of the unit-root guide and the Johansen section of the VAR guide give the formulas and tables.

Supported options

Options of oe.ardl

oe.egranger

Stata command: egranger.

oe.vecrank

Stata command: vecrank.

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 200 quarters. The regressor x is a random walk with drift 0.1 and unit shocks. The outcome follows an ARDL(1,1), y_t = 0.3 + 0.7y_{t−1} + 0.4x_t − 0.16x_{t−1} + e_t with e_t ~ N(0, 0.5²). Its long-run relation is y = 1 + 0.8x, deviations shrink by a = 0.3 per quarter, and in error-correction form the coefficient on Δx_t is 0.16 and the constant 0.3. The script runs the Engle–Granger and Johansen tests and fits an ARDL model with automatic lag selection.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n=200, seed=51):
    """x is a random walk with drift; y follows an ARDL(1,1) whose long-run relation is y = 1 + 0.8 x."""
    rng = np.random.default_rng(seed)
    x = np.cumsum(0.1 + rng.normal(size=n))
    y = np.zeros(n)
    y[0] = 1.0 + 0.8 * x[0]
    for t in range(1, n):
        y[t] = 0.3 + 0.7 * y[t - 1] + 0.4 * x[t] - 0.16 * x[t - 1] + rng.normal(scale=0.5)
    return pd.DataFrame({"quarter": np.arange(1, n + 1), "y": y, "x": x})


if __name__ == "__main__":
    df = simulate()

    eg = oe.egranger(df, "y", ["x"], time="quarter", lags=1)
    print(f"Engle-Granger tau = {eg.attrs['statistic']:.3f}, 5% cv = {eg.attrs['critical_values']['5%']:.3f}, "
          f"p = {eg.attrs['p_value']:.2g}")
    print("Cointegrating regression:", {k: round(v, 3) for k, v in eg.attrs["cointegrating_vector"].items()})

    johansen = oe.vecrank(data=df, y=["y", "x"], time="quarter", lags=2)
    print(johansen[["rank", "eigenvalue", "trace", "trace_cv5", "max", "max_cv5"]].round(3).to_string(index=False))
    print("Johansen rank estimate:", johansen.attrs["selected_rank"])

    fit = oe.ardl(data=df, y="y", x=["x"], time="quarter", maxlags=4, ec=True)
    print(fit.summary())
    for name in ["bounds_f", "bounds_t"]:
        b = fit.tests[name]
        print(f"{name}: 5% bounds I(0) = {b['critical_values']['5%']['I0']}, I(1) = {b['critical_values']['5%']['I1']}"
              f" -> {b['decision']['5%']}")

With openecon 0.3.20a1 from PyPI the script prints:

Engle-Granger tau = -5.427, 5% cv = -3.367, p = 2.5e-05
Cointegrating regression: {'x': 0.736, 'Intercept': 1.473}
 rank  eigenvalue  trace  trace_cv5    max  max_cv5
    0         NaN 40.671      15.41 37.681    14.07
    1       0.173  2.990       3.76  2.990     3.76
    2       0.015    NaN        NaN    NaN      NaN
Johansen rank estimate: 1
ARDL(2,1) regression, error-correction form — y
Observations: 196  |  Covariance: nonrobust  |  Confidence: 95%

Term            Estimate  Std. error         t     P>|stat|   CI lower   CI upper
------------  ----------  ----------  --------  -----------  ---------  ---------
[ADJ]
ADJ:L.y        -0.281633    0.042865  -6.57024  4.64929e-10  -0.366183  -0.197084
[LR]
LR:x            0.746403    0.023851   31.2944  4.11203e-77   0.699357   0.793448
[SR]
SR:LD.y       -0.0877921   0.0504625  -1.73975    0.0835139  -0.187327  0.0117433
SR:D.x          0.227508   0.0457286   4.97519  1.44883e-06    0.13731   0.317706
SR:Intercept    0.412574   0.0898711   4.59073  7.99398e-06   0.235307   0.589841

Excluded observations: 4

r_squared: 0.518201  |  adjusted_r_squared: 0.508111  |  rmse: 0.513737  |  log_likelihood: -145.035
aic: 300.07  |  bic: 316.46  |  hqic: 306.705  |  durbin_watson: 1.96662
df_model: 4  |  df_resid: 191  |  ssr: 50.4097
F test that all error-correction coefficients except the constant are zero: F(4, 191) = 51.3578, p = 2.611e-29
Pesaran-Shin-Smith bounds F test of no level relationship (Case III: unrestricted intercept, no trend): PSS bounds F(2, 191) = 21.6795
Pesaran-Shin-Smith bounds t test of no level relationship (Case III: unrestricted intercept, no trend): PSS bounds t = -6.57024
Warning: The first 4 observation(s) provide the initial lags.
bounds_f: 5% bounds I(0) = 4.94, I(1) = 5.73 -> reject H0 (no level relationship)
bounds_t: 5% bounds I(0) = -2.86, I(1) = -3.22 -> reject H0 (no level relationship)

Reading the output

• Engle–Granger rejects no cointegration. The residual statistic is −5.427, well below the 5% critical value of −3.367 (p = 2.5e-05). The static first-step regression gives a slope of 0.736 and an intercept of 1.473 against the true 0.8 and 1: it ignores the short-run dynamics, and its estimates are only a rough guide in a sample of this size.
• Johansen finds one relation. For rank 0 the trace statistic is 40.671 against a 5% critical value of 15.41, so "no cointegration" is rejected. For rank 1 it is 2.990 against 3.76, so "at most one relation" is not. The maximum-eigenvalue statistics agree, and the estimated rank is 1.
• AIC picks one lag too many. The selected model is ARDL(2,1) rather than the true ARDL(1,1); the extra short-run term SR:LD.y is −0.088 with p = 0.084. The four excluded observations are those the lag search held back.
• The adjustment speed is recovered. ADJ:L.y is −0.282 (true −0.3): about 28% of a deviation from the long-run relation is corrected each quarter. SR:D.x (0.228) and the constant (0.413) have intervals that contain the true 0.16 and 0.3.
• The long-run coefficient is slightly low. LR:x is 0.746 with a standard error of 0.024, and its interval from 0.699 to 0.793 narrowly misses the true 0.8. One in twenty 95% intervals misses by chance.
• The bounds tests reject. The F statistic of 21.68 is far above the 5% I(1) bound of 5.73, and the t statistic of −6.570 is below the I(1) bound of −3.22. The data reject "no level relationship" whether x is I(0) or I(1).

Limitations

• oe.ardl accepts no weights.
• The common oe.predict and oe.margins adapters do not cover oe.ardl; use its own result tables and family helpers.
• oe.ardl reports only the asymptotic Pesaran–Shin–Smith bounds at 10%, 5% and 1%: there are no small-sample bounds and no p-values for the bounds tests, which always use the classical covariance whatever covariance is.
• The ARDL lag search keeps p ≥ 1 and reports the selected model on the common search sample; refit with lags= to use the observations that a larger maxlags held back. HAC covariance is not offered.
• oe.egranger accepts at most five regressors, and the standard errors of its cointegrating regression are not valid for inference.
• oe.vecrank reports no p-values. Its critical values stop at K − r = 11, and those of trend="none" for K − r ≥ 7 and trend="trend" for K − r ≥ 6 come from a single published source and are flagged in the result.
• All three need complete, gap-free series; a missing value inside the series is refused rather than dropped.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Engle, R. F. and Granger, C. W. J. (1987). Co-integration and error correction: Representation, estimation, and testing. Econometrica, 55(2), 251–276. doi:10.2307/1913236
• Johansen, S. (1988). Statistical analysis of cointegration vectors. Journal of Economic Dynamics and Control, 12(2–3), 231–254. doi:10.1016/0165-1889(88)90041-3
• Osterwald-Lenum, M. (1992). A note with quantiles of the asymptotic distribution of the maximum likelihood cointegration rank test statistics. Oxford Bulletin of Economics and Statistics, 54(3), 461–472. doi:10.1111/j.1468-0084.1992.tb00013.x
• Pesaran, M. H., Shin, Y. and Smith, R. J. (2001). Bounds testing approaches to the analysis of level relationships. Journal of Applied Econometrics, 16(3), 289–326. doi:10.1002/jae.616