*! version 2.0.0 03may2026

Title

pocoint -- Phillips-Ouliaris cointegration test

Syntax

pocoint varlist [if] [in] [, options]

varlist contains 2 to 13 numeric time-series variables. The first
variable is the dependent variable in the cointegrating regression; the
remaining variables are the regressors. The data must be tsset as
a pure time series (no panel).

test(string) choose the test statistic:
Za, Zt, Pu, or Pz; default is Zt

defaults across implementations: pocoint = Zt;
Python arch = Zt; R urca::ca.po = Pu; R
tseries::po.test reports Za only; EViews reports both
Za and Zt side-by-side

trend(string) deterministic terms in the cointegrating
regression: n (none), c (constant, default),
ct (constant + trend), or ctt (constant + trend + trend^2)

kernel(string) long-run-variance kernel:
bartlett (default), parzen, or quadratic-spectral
(synonym qs)

lags(#) user-specified kernel bandwidth (truncation lag); overrides auto, lshort, and llong

auto data-driven bandwidth via the Andrews/Newey-West rule
(matches Python arch's default)

lshort use the short bandwidth trunc((n-1)/100) (default)

llong use the long bandwidth trunc((n-1)/30)

resid(name) save the cointegrating-regression residuals
in a new variable name

tseries force backwards-compatible mode that
reproduces R's tseries::po.test exactly: only Za with
trend(c) or trend(n); Bartlett kernel only; p-value via
R's built-in tabulated grid

urca force backwards-compatible mode that reproduces R's
urca::ca.po exactly: only Pu and Pz, Bartlett kernel,
trend(n), trend(c), or trend(ct); up to 6 variables;
fixed urca critical-value table; no p-value

eviews force backwards-compatible mode that reproduces
the EViews default Phillips-Ouliaris output: both Z(alpha) and
Z(t) are computed; Bartlett kernel; Newey-West fixed bandwidth
trunc(4*(T/100)^(2/9)); demeaned cointegrating regression; no
d.f. adjustment; no other option may be specified

Description

pocoint computes the four Phillips-Ouliaris (1990) residual-based
tests of the null hypothesis of no cointegration:

- Za: Phillips Z(alpha) coefficient test

- Zt: Phillips Z(t) t-ratio test (default)

- Pu: variance ratio test introduced in the paper

- Pz: multivariate trace test (invariant to the
normalization of the cointegrating regression)

The first stage runs a static OLS cointegrating regression of the first
variable on the remaining variables (and the chosen deterministic
terms). The residual-based statistic is then constructed from a
Bartlett-kernel long-run-variance estimator with the truncation lag
controlled by lags(), lshort, or llong.

Critical values at the 1%, 5%, and 10% levels and asymptotic p-values
are obtained from a response-surface fit (Hansen-Haug type), simulated
by Kevin Sheppard for the Python arch package, with up to 25
million Monte Carlo replications per cell. The same response surface
allows up to 13 variables in the cointegrating regression.

The tseries option reproduces the older R routine
tseries::po.test (A. Trapletti) exactly to within machine
precision. In that mode only Za with trend(c) or
trend(n) is computed. The p-value is obtained by linear
interpolation on R's built-in critical-value table (returning one of
0.01, 0.025, 0.05, 0.075, 0.10, 0.125, 0.15 for statistics outside the
table boundaries) with R's "smaller/greater than printed p-value"
warning.

The urca option reproduces the R routine urca::ca.po
(B. Pfaff) exactly. Only Pu and Pz are available with
trend(n), trend(c), or trend(ct), the bandwidth is
fixed to urca's lshort/llong formulas, and critical values
come from urca's internal table; no p-value is reported. Note that
urca's Pz uses a different moment matrix from the
arch/Hansen-Haug formulation, so the two implementations
intentionally produce different statistics.

The eviews option reproduces the EViews default
Phillips-Ouliaris output exactly. Both Za and Zt are
reported in a single table, with Bartlett kernel, Newey-West fixed
bandwidth trunc(4*(T/100)^(2/9)), demeaned cointegrating
regression, and no degrees-of-freedom adjustment. No other option
may be specified in this mode. Test statistics agree with EViews to
machine precision; p-values use the Hansen-Haug response surface
(EViews itself reports MacKinnon (1996) p-values, which differ
slightly).

Options

test(string) selects which Phillips-Ouliaris statistic to
compute: Zt (default), Za, Pu, or Pz.
The Z-tests are lower-tail (reject if statistic is below the
critical value). The P-tests are upper-tail (reject if statistic
is above the critical value).

trend(string) controls the deterministic terms included in the
cointegrating regression and in the residual-based step. The four
choices and their arch equivalents are: n/none,
c/constant, ct, and ctt. Different trend
choices have different asymptotic distributions, and pocoint
selects the matching response-surface coefficients automatically.

kernel(string) selects the kernel used to estimate the long-run
variance of the residual-based AR(1) innovations. All
Phillips-Ouliaris statistics depend on a HAC long-run-variance
estimate; different kernels give different finite-sample size and
power but converge to the same asymptotic distribution. Choices
are bartlett (default; the original Newey-West kernel), parzen
(smoother weights, lower bias), and quadratic-spectral (or its
synonym qs; optimal in a mean-squared-error sense for many
data-generating processes). Python arch, R tseries, and R
urca all use Bartlett by default.

lags(#) fixes the kernel truncation lag (the bandwidth) at a
user-supplied value. The bandwidth controls how many sample
autocovariances enter the long-run-variance estimate: lags(0)
gives the White covariance with no autocovariance correction (the
test statistic then reduces to a simple Phillips Z(alpha)/Z(t)),
while larger values capture more serial correlation but increase
finite-sample variance. lags() overrides auto,
lshort, and llong.

auto requests an automatic data-driven bandwidth via the
Andrews (1991) / Newey-West (1994) plug-in rule. This is the
default behaviour in Python arch, and using auto together
with kernel(bartlett) reproduces arch's default to machine
precision. The chosen bandwidth depends on the kernel and on a
preliminary AR(1) fit to the residual-based innovations; it is
displayed in the output along with the rule used.

lshort (default) sets the truncation lag to the simple
short-bandwidth rule trunc((n-1)/100) (Schwert (1989), used by
R tseries::po.test). The lag grows slowly with the sample size
and is appropriate for residuals with mild serial correlation.

llong sets the truncation lag to the long-bandwidth rule
trunc((n-1)/30) (also from Schwert (1989) and used by R
tseries::po.test when lshort=FALSE). Use this when the
data show strong persistence or large negative moving-average
components, since a longer window captures more of the spectral mass
near the origin.

resid(name) saves the residuals from the cointegrating
regression as a new variable name. These are the OLS residuals
u_t from the first-stage long-run regression of the dependent
variable on the regressors and deterministic terms. Useful for
post-estimation work such as plotting (e.g. tsline name),
running an additional unit-root test on the residuals
(e.g. dfuller name), or supplying the error-correction term to a
subsequent VECM step.

tseries forces the backwards-compatible mode that mirrors
tseries::po.test from R. Only Za with trend(c) is
computed, with the older Newey-West correction (slightly different
from the arch/Hansen-Haug formulation). In this mode no asymptotic
p-value is reported; only the test statistic and the original
Phillips-Ouliaris (1990) Tables Ia/Ib critical values would apply.

urca forces the backwards-compatible mode that mirrors R's
urca::ca.po. Supports Pu and Pz with
trend(n), trend(c), or trend(ct) only, up to 6
variables. The bandwidth defaults to trunc(4*((n-1)/100)^0.25)
(lshort) or trunc(12*((n-1)/100)^0.25) (llong).
Critical values come from the fixed table embedded in urca. The
Pz formulation in urca uses M_zz computed from raw lagged
levels (without prior trend demeaning) and therefore differs from the
arch/Hansen-Haug formulation that follows the article's remark
(f) on p. 173. No asymptotic p-value is reported.

eviews forces the backwards-compatible mode that reproduces the
EViews default Phillips-Ouliaris cointegration test. Both
Z(alpha) and Z(t) are computed and reported in the same
table. All settings are fixed to the EViews defaults and no other
option (test(), trend(), kernel(), lags(),
auto, lshort, llong, resid()) may be
specified: Bartlett kernel, Newey-West fixed bandwidth
trunc(4*(T/100)^(2/9)), demeaned cointegrating regression
(trend(c)), no degrees-of-freedom adjustment for the long-run
variance. The algorithm uses the bias-corrected autoregressive
coefficient rho* - 1 = (rho - 1) - T_eff*lambda/sum(u^2_lag)
together with the EViews convention T_eff = T - 1 (the number of
residuals after the AR(1) regression on the cointegrating
residuals). This differs from pocoint's default formula, which
follows Phillips and Ouliaris (1990) and uses the original sample
size T as the normalization factor (see Examples below for the
algorithmic difference). Critical values and p-values come from the
Hansen-Haug response surface; EViews itself reports MacKinnon (1996)
p-values. Test statistics agree with EViews to machine precision;
p-values may differ slightly because of the different response
surface.

Stored results

pocoint stores different results in r() depending on which
mode is used. All four modes share r(statistic), r(n),
r(lag), and r(method).

Default mode (Hansen-Haug response surface):

Scalars

r(statistic) test statistic

r(p) asymptotic p-value (Hansen-Haug response surface)

r(cv_1pct) 1% critical value

r(cv_5pct) 5% critical value

r(cv_10pct) 10% critical value

r(minN) smallest N used to fit the response surface

r(lag) kernel truncation lag (bandwidth) used

r(n) effective sample size

r(nvars) number of variables in the cointegrating regression

Macros

r(method) "Phillips-Ouliaris (1990) cointegration test"

r(test) test type used: Za / Zt / Pu / Pz

r(trend) deterministic specification: n / c / ct / ctt

r(kernel) kernel used: bartlett / parzen / quadratic-spectral

r(bw_method) bandwidth selection method (e.g. "lshort", "auto: Andrews/Newey-West")

r(depvar) dependent variable in the cointegrating regression

r(indeps) regressors in the cointegrating regression

tseries mode (tseries option):

Scalars

r(statistic) Z(alpha) statistic

r(p) p-value from R's tabulated grid (one of 0.01, 0.025, 0.05, 0.075, 0.10, 0.125, 0.15)

r(alpha) OLS estimate of the AR(1) coefficient on the residuals

r(lag) Bartlett bandwidth used

r(n) effective sample size

r(tseries) 1 (mode flag)

Macros

r(method) "Phillips-Ouliaris (tseries::po.test compatible)"

r(trend) deterministic specification: n or c

urca mode (urca option):

Scalars

r(statistic) test statistic (Pu or Pz)

r(cv_1pct) 1% critical value (urca's internal table)

r(cv_5pct) 5% critical value

r(cv_10pct) 10% critical value

r(lag) Bartlett bandwidth used

r(n) effective sample size

r(nvars) number of variables in the cointegrating regression

r(urca) 1 (mode flag)

Macros

r(method) "Phillips-Ouliaris (urca::ca.po compatible)"

r(test) test type used: Pu or Pz

r(trend) deterministic specification: n / c / ct

eviews mode (eviews option):

Scalars

r(z) Z(alpha) test statistic

r(tau) Z(t) test statistic

r(p_z) p-value for Z(alpha) (Hansen-Haug response surface)

r(p_tau) p-value for Z(t) (Hansen-Haug response surface)

r(rho) OLS estimate of the AR(1) coefficient

r(rho_star) bias-corrected AR(1) coefficient

r(rho_se) standard error of the bias-corrected coefficient

r(sigma2) residual variance (no d.f. adjustment)

r(omega2) long-run residual variance

r(lambda) one-sided long-run residual autocovariance

r(lag) Newey-West fixed bandwidth

r(n) effective sample size

r(nvars) number of variables in the cointegrating regression

r(eviews) 1 (mode flag)

Macros

r(method) "Phillips-Ouliaris (EViews compatible)"

r(kernel) "bartlett"

r(trend) "c"

Examples

The standard Z(t) demeaned test (the default, matching Python arch):

. webuse lutkepohl2

. pocoint consump inc

All four Phillips-Ouliaris tests with a constant + linear trend in the
cointegrating regression:

. pocoint consump inc, test(Za) trend(ct)

. pocoint consump inc, test(Zt) trend(ct)

. pocoint consump inc, test(Pu) trend(ct)

. pocoint consump inc, test(Pz) trend(ct)

A test with a manual bandwidth of 6 and the multivariate-trace
statistic Pz:

. pocoint consump inc, test(Pz) trend(c) lags(6)

Use a Parzen kernel with the Andrews/Newey-West data-driven bandwidth
(matches Python arch's default behaviour):

. pocoint consump inc, kernel(parzen) auto

Quadratic-spectral kernel with auto bandwidth:

. pocoint consump inc, kernel(qs) auto

Save the cointegrating-regression residuals for further inspection
(e.g. for use as the error-correction term in a VECM step):

. pocoint consump inc, test(Za) trend(c) resid(uhat)

. tsline uhat

. dfuller uhat, noconstant

Backwards-compatible mode reproducing R tseries::po.test:

. pocoint consump inc, tseries

Backwards-compatible mode reproducing R urca::ca.po:

. pocoint consump inc, urca test(Pu) trend(c)

. pocoint consump inc, urca test(Pz) trend(ct) llong

Backwards-compatible mode reproducing the EViews default
Phillips-Ouliaris output (both Z(alpha) and Z(t) reported in one
table; no other option is allowed):

. pocoint consump inc, eviews

Examples adapted from other softwares

Each source package (R tseries, R urca, Python arch) ships
its own example data set. We have re-saved each one as a Stata-ready
file and a CSV. The Stata files are loaded directly from the web with
use (no manual download required):

tseries_example_no_cointegration.dta 1001-obs bivariate random walk

tseries_example_cointegrated.dta 1001-obs cointegrated pair y = 2 - 3*x + N(0,5^2)

urca_ecb_example.dta quarterly Euro-area data: real M3, real GDP, long bond rate

arch_crude_example.dta 393 monthly Brent and WTI log crude-oil prices, 1987m5-2020m1

Equivalent CSVs are available at the same URL by replacing
.dta with .csv (load with import delimited). The
ECB example is also provided as an EViews workfile at
https://www.eruygurakademi.com/datasets/urca_ecb_example.wf1.

R: tseries::po.test (A. Trapletti)

Bivariate random-walk example (no cointegration), 1001 observations.
R and Stata both read the same data from the web (R as CSV, Stata as
DTA), so the test statistics agree to machine precision:

R

install.packages("tseries") # only once

library(tseries)

d <- read.csv("https://www.eruygurakademi.com/datasets/tseries_example_no_cointegration.csv")

po.test(ts(cbind(d$x1, d$x2)))

Stata equivalent:

. use https://www.eruygurakademi.com/datasets/tseries_example_no_cointegration.dta, clear

. pocoint x1 x2, tseries

Cointegrated case (1001 observations, y = 2 - 3*x + N(0,5^2)):

R

library(tseries)

d <- read.csv("https://www.eruygurakademi.com/datasets/tseries_example_cointegrated.csv")

po.test(ts(cbind(d$x, d$y)))

Stata equivalent:

. use https://www.eruygurakademi.com/datasets/tseries_example_cointegrated.dta, clear

. pocoint x y, tseries

Same data with the modern arch-style test (continuous p-value via the
Hansen-Haug response surface):

. pocoint x y, test(Za) trend(c)

. pocoint x y, test(Pz) trend(c) auto

R: urca::ca.po (B. Pfaff)

Quarterly Euro-area data (1997Q3 to 2003Q4, 26 observations): real M3,
real GDP, and the long bond rate. R reads the CSV; Stata reads the
DTA (already tsset on quarterly variable quarter):

R

install.packages("urca") # only once

library(urca)

d <- read.csv("https://www.eruygurakademi.com/datasets/urca_ecb_example.csv")

summary(ca.po(cbind(d$m3_real, d$gdp_real, d$rl), type="Pz"))

Stata equivalent:

. use https://www.eruygurakademi.com/datasets/urca_ecb_example.dta, clear

. pocoint m3_real gdp_real rl, urca test(Pz) trend(n)

The same multivariate-trace test using the article's own Hansen-Haug
formulation (which differs from urca's; see the urca option
above):

. pocoint m3_real gdp_real rl, test(Pz) trend(n)

. pocoint m3_real gdp_real rl, test(Pz) trend(n) auto

Python: arch.unitroot.cointegration.phillips_ouliaris (K. Sheppard)

Monthly WTI and Brent log crude-oil prices, 1987m5 to 2020m1, 393
observations. Python reads the CSV; Stata reads the DTA (already
tsset on monthly variable month):

Python

# pip install arch pandas # only once, in your shell

import numpy as np, pandas as pd

from arch.unitroot.cointegration import phillips_ouliaris

d = pd.read_csv("https://www.eruygurakademi.com/datasets/arch_crude_example.csv")

phillips_ouliaris(d["WTI"], d[["Brent"]], trend="c", test_type="Zt")

phillips_ouliaris(d["WTI"], d[["Brent"]], trend="c", test_type="Pz")

Stata equivalent. With the auto option pocoint
reproduces arch's default behaviour to machine precision (data-driven
bandwidth via the Andrews/Newey-West rule):

. use https://www.eruygurakademi.com/datasets/arch_crude_example.dta, clear

. pocoint WTI Brent, test(Zt) trend(c) auto

. pocoint WTI Brent, test(Za) trend(c) auto

. pocoint WTI Brent, test(Pu) trend(c) auto

. pocoint WTI Brent, test(Pz) trend(c) auto

arch defaults to the Bartlett kernel. Try the other two:

. pocoint WTI Brent, test(Zt) trend(c) kernel(parzen) auto

. pocoint WTI Brent, test(Zt) trend(c) kernel(qs) auto

With pocoint's own fixed-bandwidth default (lshort) the
statistics are close to arch's defaults but not identical, because
arch picks an automatic data-driven bandwidth:

. pocoint WTI Brent, test(Zt) trend(c)

. pocoint WTI Brent, test(Pz) trend(c) llong

EViews: Phillips-Ouliaris cointegration test

The eviews option reproduces the EViews default
Phillips-Ouliaris output to machine precision. Using the same Euro-area
ECB data set as the urca example (1997Q3-2003Q4, 26 observations),
the EViews dialog (Group -> View -> Cointegration Test, Test method =
Phillips-Ouliaris, Trend specification = Constant (Level), HAC Options
= Bartlett kernel + Newey-West Fixed bandwidth, d.f. adjustment off)
gives:

EViews output (excerpt)

Dependent tau-statistic z-statistic

GDP_REAL -2.380265 -6.749947

M3_REAL -1.927034 -6.393164

RL -2.930983 -11.513503

The matching EViews workfile is available at the same URL with the
.wf1 extension:

https://www.eruygurakademi.com/datasets/urca_ecb_example.wf1

Stata equivalent:

. use https://www.eruygurakademi.com/datasets/urca_ecb_example.dta, clear

. pocoint gdp_real m3_real rl, eviews

. pocoint m3_real gdp_real rl, eviews

. pocoint rl m3_real gdp_real, eviews

The eviews mode also prints the intermediate quantities shown in
the EViews "Intermediate Results" panel (Rho - 1, bias-corrected
Rho* - 1, Rho* S.E., residual variance, long-run residual variance,
long-run residual autocovariance), all of which agree with EViews to
machine precision. Note that these test statistics differ from the
Za and Zt produced by pocoint's default mode (which
follows the original Phillips-Ouliaris (1990) formula and uses the
sample size T as the normalization factor): EViews uses
T_eff = T - 1, the number of AR(1) residuals. Use eviews
when reproducing EViews output and the default mode (with auto)
when reproducing Python arch.

References

Phillips, P. C. B., and S. Ouliaris. 1990. Asymptotic Properties of
Residual Based Tests for Cointegration. Econometrica 58: 165-193.

Hansen, B. E. 1992. Efficient estimation and testing of cointegrating
vectors in the presence of deterministic trends. Journal of
Econometrics 53: 87-121.

Sheppard, K. arch: Autoregressive Conditional Heteroskedasticity,
unit-root and cointegration tests for Python.

Trapletti, A. tseries: Time Series Analysis and Computational
Finance. R package, function po.test.

Author

H. Ozan Eruygur

AHBV University, Ankara, Turkiye.

Department of Economics

https://www.ozaneruygur.com

eruygur@gmail.com

Eruygur Academy and Consulting (Eruygur Akademi ve Danismanlik), Ankara, Turkiye.

https://www.eruygurakademi.com

eruygurakademi@gmail.com

pocoint v3.2.0 - May 2026

Please cite as:

Eruygur, H. O. 2026. pocoint: Phillips-Ouliaris cointegration test.
Stata package version 3.2.0. Available from: https://www.eruygurakademi.com

Critical-value and p-value tables embedded in this package are
sourced from K. Sheppard's arch Python package and reflect very
large Monte Carlo simulations of the Hansen-Haug response surface.