Unit root tests in Python: ADF, KPSS, Phillips–Perron and DF-GLS

Before you model a time series you need to know what kind of series it is. A stationary series returns to its mean, or to its trend, after a shock. A series with a unit root does not: every shock shifts its path permanently, and its variance grows without bound. The difference decides whether you should difference the data or detrend it, which standard errors are valid, and whether a regression of one trending series on another means anything at all. Unit root tests turn this question into a formal test.

What unit root tests test

The augmented Dickey–Fuller (ADF) test runs the regression

Δy_t = a + d·t + b·y_{t−1} + c₁Δy_{t−1} + … + c_kΔy_{t−k} + e_t

and tests H0: b = 0 (a unit root) against b < 0 (stationarity). Under the null the t ratio of b does not follow a t distribution; Dickey and Fuller (1979) derived its distribution, whose critical values lie well below the usual −1.96. The lagged differences absorb serial correlation in the errors. Three choices shape the test:

• Deterministic terms. Include a constant for a series that fluctuates around a level and a constant and trend for a trending series. The alternative should be a plausible description of the data: a trending series tested without a trend will rarely reject.
• Lag length. Too few lags distort the test's size; too many cost power. Information criteria or a sequential t rule are common choices.
• Power. Unit-root tests have little power against stationary series with roots close to one. Failing to reject is weak evidence for a unit root.

The other tests vary this recipe:

• Phillips–Perron runs the Dickey–Fuller regression without lagged differences and corrects the statistics for serial correlation with a Newey–West estimate of the long-run variance.
• DF-GLS (Elliott, Rothenberg and Stock 1996) first removes the constant and trend by generalized least squares and then runs the Dickey–Fuller regression on the detrended series. It has better power than the ADF test when the series is close to a unit root.
• KPSS (Kwiatkowski, Phillips, Schmidt and Shin 1992) reverses the hypotheses: H0 is stationarity around a level or a trend. Using it with an ADF test helps. If the ADF test rejects and KPSS does not, the series looks stationary; if KPSS rejects and the ADF test does not, it looks integrated; if both or neither reject, the data cannot tell.
• Zivot–Andrews (Zivot and Andrews 1992) allows one break in the intercept, the trend or both at a date chosen from the data. A one-time shift in a stationary series can make it look like a unit-root process to tests that ignore the break.

In OpenEconometrics

The five tests are functions that return result tables, not model fits. Each table prints and exports with .to_latex(), and the scalar results are in .attrs: statistic, p_value, critical_values ("1%", "5%", "10%"), lags, trend, nobs and notes on the sources used. The call shapes, with the Stata commands they correspond to:

kpss and zandrews are user-written Stata commands by Baum. Note the different default of oe.kpss: trend=False tests level stationarity, whereas the Stata command includes a trend unless notrend is given.

• Time. time= is optional; without it the row order is the time order. With it, rows are sorted, integer periods must be consecutive and distinct, and datetimes are taken as consecutive in sorted order. Missing values are refused, because dropping an interior observation would change every lag.
• oe.dfuller. trend= is "none", "constant" (default), "trend" or "drift", and lags= (default 0) sets the lagged differences. The p-value is MacKinnon's approximate asymptotic p-value; critical_values come from MacKinnon's finite-sample response surfaces and critical_values_fuller from Fuller's table. regress=True also returns the test regression.
• oe.pperron. Reports Z(rho) and Z(t); the default truncation is int(4(n/100)^(2/9)) Newey–West lags. Only Z(t) has a p-value.
• oe.dfgls. Reports one row per lag up to maxlag (default Schwert's rule) and selects the lag by the modified AIC (MAIC); the sequential-t and SC choices are in attrs. With a trend the critical values are those of Elliott, Rothenberg and Stock and there is no p-value.
• oe.kpss. Reports every truncation from 0 to lags (default Schwert's int(12(T/100)^(1/4))); auto=True uses an automatic Newey–West bandwidth.
• oe.zandrews. break_= is "intercept", "trend" or "both", trim=0.15 excludes the ends of the sample, and lags= is an integer or "aic", "bic" or "t". attrs["break_period"] is the last observation before the break.

The unit-root family guide gives the formulas and the sources of every table, and also covers the panel unit-root and structural-break tests.

Supported options

oe.dfuller

Stata command: dfuller.

oe.pperron

Stata command: pperron.

oe.kpss

Stata command: kpss.

oe.dfgls

Stata command: dfgls.

oe.zandrews

Stata command: zandrews.

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

Example

The script simulates 300 observations of three series. walk is a random walk with drift 0.1 and unit shocks, so it has a unit root. trend is a linear trend with slope 0.1 plus AR(1) noise with coefficient 0.7, so it is trend stationary. shift is AR(1) noise with coefficient 0.8 around a level that jumps by 6 after period 150, so it is stationary with a break. Every test includes a constant and a trend.

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


def ar1(rng, n, phi):
    """Stationary AR(1) noise with unit innovation variance."""
    u = np.zeros(n)
    for t in range(1, n):
        u[t] = phi * u[t - 1] + rng.normal()
    return u


def simulate(n=300, seed=41):
    rng = np.random.default_rng(seed)
    t = np.arange(1, n + 1)
    return pd.DataFrame({
        "t": t,
        "walk": np.cumsum(0.1 + rng.normal(size=n)),          # unit root with drift: I(1)
        "trend": 0.1 * t + ar1(rng, n, 0.7),                  # stationary around a linear trend
        "shift": 6.0 * (t > 150) + ar1(rng, n, 0.8),          # stationary around a level that jumps at t = 150
    })


if __name__ == "__main__":
    df = simulate()
    rows = []
    for series in ["walk", "trend", "shift"]:
        tests = {
            "ADF": oe.dfuller(df, series, time="t", lags=2, trend="trend"),
            "Phillips-Perron": oe.pperron(df, series, time="t", trend="trend"),
            "DF-GLS": oe.dfgls(df, series, time="t"),
            "Zivot-Andrews": oe.zandrews(df, series, time="t", break_="intercept", lags=2),
            "KPSS": oe.kpss(df, series, time="t", trend=True),
        }
        for name, res in tests.items():
            a = res.attrs
            rows.append({"series": series, "test": name, "H0": "stationary" if name == "KPSS" else "unit root",
                         "lags": a["lags"], "statistic": a["statistic"], "5% cv": a["critical_values"]["5%"],
                         "p_value": a["p_value"]})
            if name == "Zivot-Andrews":
                rows[-1]["test"] += f" (break at t={a['break_period']})"
    print(pd.DataFrame(rows).round(3).to_string(index=False))

With openecon 0.3.20a1 from PyPI the script prints:

series                           test         H0  lags  statistic  5% cv  p_value
  walk                            ADF  unit root     2     -2.013 -3.425    0.594
  walk                Phillips-Perron  unit root     5     -2.288 -3.425    0.441
  walk                         DF-GLS  unit root     2     -1.991 -2.890      NaN
  walk  Zivot-Andrews (break at t=95)  unit root     2     -3.854 -4.800      NaN
  walk                           KPSS stationary    15      0.226  0.146    0.010
 trend                            ADF  unit root     2     -7.609 -3.425    0.000
 trend                Phillips-Perron  unit root     5     -7.830 -3.425    0.000
 trend                         DF-GLS  unit root     9     -3.728 -2.890      NaN
 trend Zivot-Andrews (break at t=225)  unit root     2     -8.097 -4.800      NaN
 trend                           KPSS stationary    15      0.070  0.146    0.100
 shift                            ADF  unit root     2     -4.803 -3.425    0.000
 shift                Phillips-Perron  unit root     5     -5.111 -3.425    0.000
 shift                         DF-GLS  unit root     3     -3.921 -2.890      NaN
 shift Zivot-Andrews (break at t=150)  unit root     2     -6.967 -4.800      NaN
 shift                           KPSS stationary    15      0.152  0.146    0.045

Reading the output

• The random walk keeps its unit root. None of the four unit-root tests rejects: the ADF statistic is −2.013 against a 5% critical value of −3.425 (p = 0.594), Phillips–Perron gives p = 0.441, DF-GLS −1.991 against −2.890 and Zivot–Andrews −3.854 against −4.800. KPSS rejects stationarity: 0.226 exceeds its 5% critical value of 0.146. The printed p-value of 0.010 is the lower bound of the KPSS table, so the true p-value is below 0.01.
• The trend-stationary series is recognized. ADF (−7.609) and Phillips–Perron (−7.830) reject the unit root with p-values below 0.001, and DF-GLS (−3.728) is beyond its critical value. KPSS (0.070) does not reject; its 0.100 is the upper bound of the table. This series should be detrended, not differenced.
• A break misleads KPSS here. For shift the unit-root tests reject even without modelling the break (ADF −4.803, p = 0.0005). KPSS, however, rejects stationarity at the 5% level (0.152 against 0.146, p = 0.045), because it reads the jump in the level as non-stationarity. The ADF and KPSS verdicts conflict, and the break is the reason.
• Zivot–Andrews finds the break. For shift it dates the break at t = 150, the last observation before the jump, and its statistic of −6.967 is far below the 5% critical value of −4.800. For the other two series the reported break dates have no meaning; the test always reports the date of its minimum statistic.
• The lag columns differ by design. ADF and Zivot–Andrews use the two lags requested, Phillips–Perron its default of 5 Newey–West lags and KPSS Schwert's default of 15 for 300 observations. DF-GLS chooses 2, 9 and 3 lags by MAIC. DF-GLS and Zivot–Andrews have no tabulated p-value, so compare their statistics with the critical values.

Limitations

• There is no missing="drop": series must be complete and gap free.
• oe.dfgls uses the Elliott–Rothenberg–Stock critical values only (Stata's default Cheung–Lai response surface is not reproduced) and reports no p-value with a trend.
• oe.kpss uses the Bartlett kernel only, and its p-value is interpolated in the published table and bounded to [0.01, 0.10].
• oe.zandrews allows one break, uses one lag order for all candidate break dates and reports no p-value.
• oe.pperron reports no p-value for Z(rho), only critical values.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Dickey, D. A. and Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74(366a), 427–431. doi:10.1080/01621459.1979.10482531
• Elliott, G., Rothenberg, T. J. and Stock, J. H. (1996). Efficient tests for an autoregressive unit root. Econometrica, 64(4), 813–836. doi:10.2307/2171846
• Kwiatkowski, D., Phillips, P. C. B., Schmidt, P. and Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root. Journal of Econometrics, 54(1–3), 159–178. doi:10.1016/0304-4076(92)90104-Y
• Zivot, E. and Andrews, D. W. K. (1992). Further evidence on the great crash, the oil-price shock, and the unit-root hypothesis. Journal of Business & Economic Statistics, 10(3), 251–270. doi:10.1080/07350015.1992.10509904