Stationarity: ADF, KPSS, differencing and fractional differencing.

What stationary means

A series is stationary (in the weak sense) if its mean, its variance and the covariance between two observations $k$ periods apart do not depend on the date. Put another way: a stretch from 2005 and one from 2023 are two samples from the same distribution. That is what makes it meaningful to estimate a mean, a volatility or a regression with past data and use it in the future.

Prices are not stationary. An index that traded at 1,000 in 2003 and at 5,000 in 2024 does not have “a mean” to return to; its variance grows with the horizon. The simplest model that captures this is the random walk: today's price is yesterday's plus a shock, $p_t = p_{t-1} + \varepsilon_t$. Its difference, the return, is stationary. In the jargon, the price is integrated of order 1, $I(1)$, and the return, $I(0)$.

flowchart LR
  P[Price I1] --> L[Logarithm]
  L --> D1[Difference d=1]
  L --> DF[Fractional diff.]
  D1 --> R[Return I0]
  DF --> M[I0 with memory]

Why it matters in a trading system:

• Spurious regressions. Two independent integrated series give huge correlations and t's by chance (lesson 0.7.2). A regression of one price on another only makes sense if they are cointegrated (lesson 2.1.5).
• Mean reversion. A rule that buys when the price is “below its mean” assumes there is a mean to return to. In an $I(1)$ price there is none; in a stationary spread, there is (lesson 2.3.3).
• Model features. A model trained on the level of a price learns the relationship of a period that will not repeat (level 3).

The data in this lesson are synthetic (level 2 README): an index future that is a random walk with drift, commodity futures that return very slowly to an anchor (half-lives of years) and pairs whose true spread is stationary or not, depending on the pair. The answer is known in advance; that makes it possible to see where the tests get it right and where they do not.

Dickey-Fuller: testing the unit root

The idea of the Dickey-Fuller test is to estimate how much the series “comes back”. In the regression

a negative $\gamma$ means that, when $y$ is high, it tends to go down: reversion. If $\gamma = 0$, the series is a random walk (it has a unit root). The ADF test (augmented, because of the lags $\Delta y_{t-i}$, which absorb short-term autocorrelation) uses the null $H_0: \gamma = 0$ against $\gamma < 0$.

The catch is that, under the null, $y_{t-1}$ is not stationary and the t-statistic of $\gamma$ does not follow a Student's t. It follows the Dickey-Fuller distribution, shifted further to the left. You can see it by simulating seeded random walks and computing the statistic on each one:

import numpy as np
import pandas as pd
from scipy import stats

rng = np.random.default_rng(212)
sims, n = 4000, 1000
y = rng.standard_normal((sims, n)).cumsum(axis=1)  # walks
dy, ylag = np.diff(y, axis=1), y[:, :-1]

def t_gamma(dy, ylag):
    # OLS of dy = a + g·ylag, row by row
    xm = ylag - ylag.mean(1, keepdims=True)
    g = (xm * dy).sum(1) / (xm**2).sum(1)
    a = dy.mean(1) - g * ylag.mean(1)
    e = dy - a[:, None] - g[:, None] * ylag
    s2 = (e**2).sum(1) / (dy.shape[1] - 2)
    return g / np.sqrt(s2 / (xm**2).sum(1))

t = t_gamma(dy, ylag)
print(f"5% percentile of the Dickey-Fuller t: "
      f"{np.percentile(t, 5):.2f}  (normal: "
      f"{stats.norm.ppf(0.05):.2f})")
print(f"“stationary” walks with the normal threshold: "
      f"{np.mean(t < stats.norm.ppf(0.05)):.1%}")

5% percentile of the Dickey-Fuller t: -2.88  (normal: -1.64)
“stationary” walks with the normal threshold: 45.3%

The lower 5% of the simulated statistic is at −2.88 (the tabulated critical value for the regression with a constant and this length is −2.86; the difference is the noise of 4,000 simulations), and not at the normal's −1.64. Anyone who uses the Student's t threshold will declare 45.3% of random walks stationary: almost half. statsmodels.tsa.stattools.adfuller already uses the correct critical values (MacKinnon's) and chooses the number of lags $p$ by information criterion (autolag="AIC").

Two opposite nulls: ADF and KPSS

The ADF has a design problem: the null is the unit root, so “not rejecting” can mean that the series is a random walk or that the sample is too short to see it. The KPSS test (Kwiatkowski, Phillips, Schmidt and Shin, 1992) swaps the roles: its null is that the series is stationary around a constant, and it rejects it when the cumulative sum of the deviations grows too much. Using both together gives more information than either one separately.

import warnings
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tools.sm_exceptions import InterpolationWarning

# kpss warns when the statistic falls outside its table of
# p-values (it clips them to 0.01 or 0.10): the warning is expected
warnings.filterwarnings("ignore", category=InterpolationWarning)

def contrastes(x):
    x = pd.Series(x).dropna()
    a = adfuller(x, regression="c", autolag="AIC",
                 result_object=True)
    k = kpss(x, regression="c", nlags="auto", result_object=True)
    return {"n": len(x), "adf": a.statistic, "p_adf": a.pvalue,
            "kpss": k.statistic, "p_kpss": k.pvalue}

cont = pd.read_parquet("futuros2_continuos.parquet")
ida = np.log(cont[cont["mercado"] == "IDA"]
             .set_index("fecha")["cierre_ajustado"])
pd.DataFrame({"IDA log price": contrastes(ida),
              "IDA return": contrastes(ida.diff())}).T.round(3)

                    n     adf  p_adf    kpss  p_kpss
IDA log price  6539.0   0.176  0.971  10.817    0.01
IDA return     6538.0 -26.296  0.000   0.264    0.10

With result_object=True, both functions return an object with names (statistic, pvalue, critical_values) instead of a tuple whose length changes with the arguments; statsmodels warns with a FutureWarning that this will be the default behavior, and asking for it now silences it. kpss returns p-values interpolated in a table that runs from 0.01 to 0.10; outside that range it clips them and emits an InterpolationWarning. The warning does not signal an error, but that the real p-value is below 0.01 or above 0.10; only that category is silenced so that the output is clean. The KPSS critical values with a constant are 0.347 (10%), 0.463 (5%) and 0.739 (1%).

The log price of IDA, a fictional equity index built as a random walk with drift, gives an ADF of 0.18 (p = 0.97): the unit root is not rejected. The KPSS is 10.82, far above 0.463: stationarity is rejected. The two say the same thing in opposite languages. The daily return gives the opposite: ADF of −26.3 (p ≈ 0) and KPSS of 0.26, below the threshold. It is the textbook case.

Spreads with a known answer

The pairs in pares2.parquet make it possible to go further. All their relationships are in logarithms: in pairs A, B and C, the spread $u = \ln y - \alpha - \beta \ln x$ is a stationary Ornstein-Uhlenbeck process with a half-life of 5, 15 and 40 sessions (a percentage, with a standard deviation of 2% to 4.5%); in E, the same spread is a random walk, even though the two assets have a return correlation of 0.86. Here the true $\alpha$ and $\beta$ from pares2_verdad.csv are used: estimating them changes the critical values (Engle-Granger, lesson 2.1.5).

pares = pd.read_parquet("pares2.parquet").pivot(
    index="fecha", columns="activo", values="precio")
verdad = pd.read_csv("pares2_verdad.csv").set_index("par")

def spread(par):
    # all relationships are in logarithms (formulacion = "log")
    v = verdad.loc[par]
    y, x = pares[v["activo_y"]], pares[v["activo_x"]]
    return np.log(y) - v["alfa"] - v["beta"] * np.log(x)

tabla = {}
for par in ["A", "B", "C", "E"]:
    vm = verdad.loc[par, "vida_media_sesiones"]
    s = spread(par)
    tabla[f"{par}, half-life {vm:g}"] = contrastes(s)
    tabla[f"{par}, last year"] = contrastes(s.iloc[-252:])
pd.DataFrame(tabla).T.round(3)

                       n     adf  p_adf   kpss  p_kpss
A, half-life 5    3018.0 -14.929  0.000  0.108   0.100
A, last year       252.0  -4.251  0.001  0.287   0.100
B, half-life 15   3018.0  -8.975  0.000  0.385   0.084
B, last year       252.0  -2.825  0.055  0.133   0.100
C, half-life 40   3018.0  -5.562  0.000  0.826   0.010
C, last year       252.0  -1.643  0.461  0.333   0.100
E, half-life nan  3018.0  -2.012  0.281  4.775   0.010
E, last year       252.0  -2.371  0.150  0.586   0.024

With the full twelve years (3,018 sessions), the ADF gets all four right: it rejects the unit root in A, B and C (statistics of −14.9, −9.0 and −5.6) and does not reject it in E (p = 0.28). The KPSS gets A and E right, but rejects the stationarity of C (0.826 > 0.739, p ≤ 0.01), which is stationary. With a single year of data, the picture changes: the ADF no longer rejects either in B (p = 0.055) or in C (p = 0.46). With a half-life of 40 sessions, a year contains about six reversion cycles, and that is not enough to tell them apart from a walk.

The four combinations

The top two rows are the clear answer. The bottom two say that the test cannot answer with these data, and they are more frequent than it seems with short or persistent series, which are exactly the ones of interest for mean reversion.

flowchart TB
  A{ADF rejects?} -- yes --> K1{KPSS rejects?}
  A -- no --> K2{KPSS rejects?}
  K1 -- no --> E[Stationary]
  K1 -- yes --> P[Very persistent]
  K2 -- yes --> U[Unit root]
  K2 -- no --> M[Not enough sample]

h = spread("H")
ruptura = verdad.loc["H", "fecha_ruptura"]
pd.DataFrame({"full": contrastes(h),
              "before the break": contrastes(h[:ruptura].iloc[:-1]),
              }).T.round(3)

                       n    adf  p_adf   kpss  p_kpss
full              3018.0 -1.448  0.559  6.643   0.010
before the break  1875.0 -5.702  0.000  0.407   0.074

Differencing: what you gain and what you lose

Differencing once ($d = 1$) turns the price into a return and solves stationarity. The cost is memory: today's return says nothing about the level the price is at. It is measured with the correlation between the transformed series and the original:

r = ida.diff().dropna()
print(f"correlation of the return with the log price: "
      f"{np.corrcoef(r, ida.loc[r.index])[0, 1]:.3f}")

# over-differencing a spread that was already stationary
ds = spread("A").diff().dropna()
print(f"lag-1 autocorrelation of spread A: "
      f"{spread('A').autocorr(1):.3f}; of its difference: "
      f"{ds.autocorr(1):.3f}")

correlation of the return with the log price: 0.025
lag-1 autocorrelation of spread A: 0.865; of its difference: -0.055

The return has a correlation of 0.025 with the price level: all the information about where the market is has disappeared. And differencing what was already stationary damages the series: spread A has an autocorrelation of 0.865 at lag 1 (the theoretical Ornstein-Uhlenbeck value is $0.5^{1/5} = 0.871$); its difference has −0.055, close to the value of an over-differenced AR(1) ($-(1-\phi)/2 = -0.065$): the positive persistence that revealed the reversion has disappeared and what is left is a series with almost no memory. The rule: difference just enough, not by reflex.

Fractional differencing

Is there anything between $d = 0$ (the price, with all its memory and non-stationary) and $d = 1$ (the return, stationary and without memory)? Yes: the difference operator $(1-B)^d$, with $B$ the lag operator ($B\,y_t = y_{t-1}$), admits any real $d$. Its binomial expansion gives an infinite weighted sum of the present and the past,

With $d = 1$ the weights are $1, -1, 0, 0, \ldots$: the usual difference. With $0 < d < 1$ all the past weights are negative and decay slowly, like a power: the series remembers the level for a long time.

def pesos(d, umbral=1e-4, kmax=100_000):
    w = [1.0]
    for k in range(1, kmax):
        wk = -w[-1] * (d - k + 1) / k
        if abs(wk) < umbral:
            break
        w.append(wk)
    return np.array(w)

for d in [0.1, 0.3, 0.5, 1.0]:
    w = pesos(d)
    print(f"d={d:.1f}: {len(w):4d} weights; first "
          f"{np.round(w[:5], 3).tolist()}")

d=0.1:  503 weights; first [1.0, -0.1, -0.045, -0.028, -0.021]
d=0.3:  388 weights; first [1.0, -0.3, -0.105, -0.06, -0.04]
d=0.5:  200 weights; first [1.0, -0.5, -0.125, -0.062, -0.039]
d=1.0:    2 weights; first [1.0, -1.0]

The fixed-window version (fixed-width window fractional differencing, López de Prado, 2018, chapter 5) cuts the sum where the weights fall below a threshold (here $10^{-4}$) and applies the same window to every date. That way all values carry the same amount of history and none uses future data: the value on date $t$ only combines $y_t, y_{t-1}, \ldots$

def ffd(x, d, umbral=1e-4):
    w = pesos(d, umbral)
    # “valid” convolution: the value at t uses y_t … y_(t-L+1)
    v = np.convolve(x.to_numpy(), w, mode="valid")
    return pd.Series(v, index=x.index[len(w) - 1:])

# check: d = 1 is the log return
print(np.allclose(ffd(ida, 1.0), ida.diff().dropna()))
x3 = ffd(ida, 0.3)
print(f"d = 0.3: starts on {x3.index[0]:%Y-%m-%d}, "
      f"{len(ida) - len(x3)} sessions after the price")

True
d = 0.3: starts on 2001-07-17, 387 sessions after the price

Now $d$ is swept and, for each value, the ADF of the transformed series and its correlation with the log price are computed. Following López de Prado, the ADF uses a single lag (maxlag=1, with no automatic selection), so that all the $d$ values are compared with the same regression:

filas = []
for d in np.round(np.arange(0, 1.01, 0.05), 2):
    x = ffd(ida, d)
    a = adfuller(x, maxlag=1, regression="c", autolag=None,
                 result_object=True)
    crit = a.critical_values
    filas.append({"d": d, "adf": a.statistic, "p": a.pvalue,
                  "corr": np.corrcoef(x, ida.loc[x.index])[0, 1]})
barrido = pd.DataFrame(filas).set_index("d")
d_min = barrido.index[barrido["adf"] < crit["5%"]].min()
print(f"5% critical value: {crit['5%']:.2f}; minimum d: {d_min}")
barrido.loc[[0.0, 0.1, 0.2, 0.25, 0.3, 0.5, 1.0]].round(3)

5% critical value: -2.86; minimum d: 0.25
         adf      p   corr
d
0.00   0.393  0.981  1.000
0.10  -0.877  0.795  0.998
0.20  -2.166  0.219  0.987
0.25  -3.152  0.023  0.977
0.30  -4.173  0.001  0.963
0.50 -12.092  0.000  0.845
1.00 -57.513  0.000  0.025

The figure at the top draws the full sweep. The ADF falls continuously as $d$ rises and crosses the critical value (−2.86) at $d = 0.25$ (ADF of −3.15; with $d = 0.20$ it is still −2.17); with that $d$, the series is stationary according to the ADF and keeps a correlation of 0.977 with the log price. The return ($d = 1$) is far more stationary than necessary (ADF of −57.5) and has lost all its memory (0.025). Between the two there is a range of series that preserve information on the level without the problems of a non-stationary regressor.

What this changes in a system

Three warnings. The first, that “not rejecting” is not “accepting”: an ADF that does not reject with one year of data does not prove that the spread is a random walk. The second, that the tests in this lesson measure stationarity of the past: pair H was stationary until it stopped being so. The third, that sweeping $d$, windows or pairs until one passes the ADF is a multiple search like any other (lesson 2.1.8).

Exercises

1. Repeat the sweep of $d$ with the gold (MEA) and light crude (ENA) futures. What minimum $d$ comes out for each one and what correlation does it keep?

def d_minima(serie, crit=-2.86):
    for d in np.round(np.arange(0, 1.01, 0.05), 2):
        x = ffd(serie, d)
        adf = adfuller(x, maxlag=1, regression="c", autolag=None,
                       result_object=True).statistic
        if adf < crit:
            return d, adf, np.corrcoef(x, serie.loc[x.index])[0, 1]

for m in ["IDA", "MEA", "ENA"]:
    s = np.log(cont[cont["mercado"] == m]
               .set_index("fecha")["cierre_ajustado"])
    d, adf, c = d_minima(s)
    print(f"{m}: minimum d {d:.2f}, ADF {adf:.2f}, correlation {c:.3f}")

IDA: minimum d 0.25, ADF -3.15, correlation 0.977
MEA: minimum d 0.15, ADF -3.05, correlation 0.981
ENA: minimum d 0.10, ADF -3.40, correlation 0.988

2. The ADF has more power the longer the sample. For spread B (half-life of 15 sessions), compute the ADF p-value with the last 126, 252, 504 and 1,008 days. From what length does it reject the unit root at 5%?

sb = spread("B")
for n in [126, 252, 504, 1008]:
    p = adfuller(sb.iloc[-n:], regression="c", autolag="AIC",
                 result_object=True).pvalue
    print(f"{n:5d} sessions: p = {p:.4f}")

  126 sessions: p = 0.2449
  252 sessions: p = 0.0548
  504 sessions: p = 0.0090
 1008 sessions: p = 0.0000

3. A colleague proposes a reversion rule: buy IDA when its price is 10% below its 200-session mean. What do the tests in this lesson say about that rule, and what would you test instead?