Stationarity Pitfalls in Climate Data#

The unit-root tests in Testing for stationarity and cointegration were developed for macroeconomic series: a few hundred quarterly observations, no seasonality left after adjustment, a plausible constant or linear trend. Climate series break most of those assumptions at once. This page collects the failure modes and what to do about each.

Why climate series stress the tests#

Four properties of climate data push the Augmented Dickey-Fuller (ADF) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests into their weak regimes:

Strong deterministic trends. Forced warming produces a series with an unmistakable upward slope. ADF with only a constant cannot represent “trending but mean reverting”, so it reports a unit root regardless of whether the trend is stochastic or deterministic.

Low power against persistent alternatives. Ocean heat content, soil moisture, and multi-decadal oscillation indices are highly persistent but not necessarily integrated. ADF fails to reject a unit root for an AR(1) with a coefficient of 0.95 far more often than not, at any sample size a climate record realistically offers.

Structural breaks. Volcanic eruptions, regime shifts in observing systems, and changes in instrumentation all inject level or variance shifts. A level shift makes a stationary series look integrated to KPSS while ADF may still reject — the textbook route to a conflicting verdict.

Long memory. Several climate series are better described as fractionally integrated than as cleanly I(0) or I(1). The tests have no way to say so; they answer the binary question they were asked.

The practical consequence: conflicting and inconclusive entries in the joint_status column are the normal outcome for climate data, not a sign that something went wrong. Read the table, do not read only the verdict.

Trend stationarity versus difference stationarity#

A series with an upward slope can be trend stationary (deterministic trend plus stationary noise) or difference stationary (a stochastic trend, with or without drift). The tests distinguish these poorly, and the choice changes what a VAR fitted on the data actually means.

Test both ways before deciding:

from impulso import adf_test, kpss_test

adf_test(data, regression="ct")     # allow a deterministic trend
kpss_test(data, regression="ct")    # test stationarity around that trend

Three routes are available once you have decided, and they are not equivalent:

  1. Fit in levels with regression="ct" reasoning behind it. Appropriate if you believe the trend is deterministic. The VAR sees the trending series; impulse responses are responses of the level, and they need not die out.

  2. Difference the series. Appropriate if you believe the trend is stochastic. Impulse responses are now responses of the growth rate, and they do die out. This is a different question, answered about a different object — not a preprocessing detail.

  3. Keep levels and pass the trend as an exogenous regressor. VARData accepts exog, so a deterministic time trend (or a forcing series) can enter the model directly, leaving the endogenous block to carry the dynamics.

There is a good argument that anthropogenically forced warming is closer to trend stationary than to difference stationary — the trend has a physical driver rather than being an accumulation of shocks. That argument, not the p-value, is what should decide the specification. State which route you took and why.

Anomalies do not make a series stationary#

Converting to anomalies (subtracting a climatological baseline, usually per calendar month) removes the seasonal cycle. It is easy to assume this makes the series stationary. It does not.

Anomalies remove the seasonal unit roots. They leave the zero-frequency stochastic trend untouched, because subtracting a fixed monthly baseline is a level shift per month, not a difference. Deseasonalised Mauna Loa carbon dioxide is the clean illustration: the sawtooth is gone and the series is still visibly I(1).

Test the anomalies, not the raw series and not your intuition:

from impulso import integration_order

integration_order(anomalies, max_order=2).summary()

Spurious regression#

Two independent series that both trend upward will look strongly related. Regress one on the other and you get a large coefficient, a high goodness-of-fit, and residuals that are themselves integrated. Nothing about the relationship is real.

This is the reason the pretests matter for a VAR and not just for a single regression. A VAR fitted in levels on integrated series inherits the same problem: the estimated dynamics can reflect shared trending rather than any transmission mechanism, and standard inference on those coefficients is not valid. Establishing the integration order first, and the cointegration rank second, is what separates a real long-run relationship from a shared drift.

Cointegration and what to do about it#

If two climate series share a stochastic trend — cumulative emissions and global mean surface temperature is the obvious pair, given the near-linear relationship between cumulative carbon dioxide and warming — they are cointegrated, and the cointegrating relationship is the physically meaningful part.

from impulso import select_lag_order, johansen_test

p = select_lag_order(data, max_lags=8).bic
result = johansen_test(data, det_order=0, k_ar_diff=p - 1)
result.rank

A rank of 1 or more is a warning about differencing. Differencing every series removes the levels information, and the long-run relationship lives entirely in the levels. You would be throwing away the part of the data that answers the question.

Impulso does not implement a vector error-correction model (VECM). That is the textbook answer to a cointegrated system, and it is out of scope. Two routes remain:

Fit the VAR in levels. This is the standard Bayesian position (Sims, Stock and Watson): a VAR in levels is consistent whether or not the series are cointegrated, and it neither imposes nor discards the cointegrating relationships. The classical objection is about the validity of hypothesis tests on the coefficients, which does not bite for posterior inference in the same way. The Minnesota prior already shrinks each equation toward a random walk, which is the right prior mean for integrated data. This is the recommended default.

Difference and accept the cost. Legitimate if the long-run relationship is not what you are asking about — if the question is genuinely about short-run dynamics between growth rates. Say explicitly that the long-run relationship was discarded, so that nobody reads a levels interpretation into the results.

What to record#

Whatever you decide, the pretests are part of the specification, so record them the way you would record the lag order:

  • The full tables, not the verdicts. adf.summary(), kpss.summary(), integration_order(...).summary(), johansen_test(...).summary().

  • Every place the two tests conflicted, and what you did about it. The inconclusive list on IntegrationOrderResult is the short version.

  • The significance level and the deterministic terms used. A verdict without its alpha and its regression is not reproducible.

  • d_max, the highest integration order in the system. It is the augmentation term a Toda-Yamamoto style procedure needs, and it is far easier to carry forward now than to reconstruct later. See Granger causality and Toda-Yamamoto for the procedure that consumes it — including why it refuses to run when inconclusive is non-empty.

  • The cointegration rank and the lag order it was conditioned on. The rank is not invariant to k_ar_diff.