Stationarity and Unit Roots

"They keep telling me to return to the mean. I have looked everywhere and I cannot find one. Every shock I absorb, I keep forever, and I am quietly proud of my long memory."

A Unit Root That Refuses to Revert to the Mean

The previous section, Section 3.1, treated a time series as a single realization of a stochastic process and built the vocabulary of expectation, autocovariance, and the ergodic hope that one long sample can stand in for the ensemble. That hope has a price, and stationarity is the price. We can only estimate a mean from one trajectory if the mean is the same at every point in time; we can only estimate an autocovariance from one trajectory if the autocovariance depends on the gap between two times and not on where in the series we happen to look. This section makes those requirements exact, shows what goes wrong when they fail, and gives the standard recipe for restoring them. The thread running through it is a warning that Chapter 2 first raised at evaluation time and that returns here at modeling time: non-stationarity is the quiet way that a temporal analysis fools itself, and the cure is to detect it deliberately rather than to discover it in production.

1. Strict Versus Weak (Covariance) Stationarity Intermediate

A stochastic process $\{X_t\}$ is strictly stationary when its entire probabilistic law is invariant under a shift of the time origin. Formally, for every finite set of indices $t_1 < t_2 < \dots < t_k$ and every lag $\tau$, the joint distribution satisfies

Nothing about the joint behavior of any group of observations depends on when they were observed, only on their spacing. This is a beautiful and almost useless condition, because verifying it requires knowing every finite-dimensional distribution, and estimating those from one finite sample is hopeless. Strict stationarity is the ideal we rarely test and the assumption we sometimes prove holds for a model we have constructed.

The workhorse weakens this to its first two moments. A process is weakly stationary (equivalently covariance stationary or second-order stationary) when three conditions hold for all $t$ and all lags $h$:

The mean does not drift, the variance neither explodes nor vanishes, and the autocovariance between two observations depends only on the lag $h$ that separates them, never on the absolute time $t$. That third condition is the one that makes estimation possible: because $\gamma(h)$ is the same wherever we measure it, we can average over the whole series to estimate it, which is exactly the ergodic bargain of Section 3.1. The function $\gamma(h)$ is the autocovariance function, its normalized form $\rho(h) = \gamma(h)/\gamma(0)$ is the autocorrelation function (the subject of Section 3.3), and weak stationarity is precisely the condition under which these two functions are well defined as properties of the process rather than artifacts of a window. The notation $\mu$, $\sigma^2$, $\gamma(h)$ follows the unified table in Appendix A.

Strict stationarity implies weak stationarity whenever the second moment is finite, but the converse fails: a process can have a constant mean, constant variance, and lag-only autocovariance while its higher moments drift, so it is weakly but not strictly stationary. The two notions coincide for one important family, the Gaussian process, because a multivariate normal law is determined entirely by its mean vector and covariance matrix, so fixing the first two moments fixes everything. Most classical time-series theory lives in the weakly stationary world and leans on Gaussian or near-Gaussian behavior to make that world the strictly stationary one. The models in the chapters ahead each assume weak stationarity of something: an AR($p$) process assumes the series itself is stationary (which constrains its coefficients), a moving-average MA($q$) process is stationary by construction, an ARIMA model assumes the series becomes stationary after differencing (this section's subject), and a GARCH model in Chapter 5 assumes the conditional variance is stationary even when the level is not.

2. Why Stationarity Matters Beginner

Stationarity is what makes the past a guide to the future. A forecast is an extrapolation of an estimated law beyond the data, and the extrapolation is only valid if the law that held during estimation still holds during the forecast horizon. When a series is stationary, the mean we measured is the mean that will persist, the variance bounds the spread we will see, and the autocovariance structure we fitted is the dependence the future will obey. When a series is non-stationary, every one of those estimates is attached to a moving target, and the model confidently predicts a world that has already moved on.

The failure is not merely a loss of accuracy; it corrupts the inference machinery itself. Ordinary least squares, the workhorse behind regression-based forecasting, derives its standard errors, t-statistics, and p-values under assumptions that include a stable second-moment structure. Apply that machinery to a non-stationary regressor and the standard errors are computed against a variance that grows with the sample. The t-statistics inflate, the p-values shrink far too small, and the model declares significant relationships that are pure accidents of trending. Subsection four makes this concrete with two independent random walks that a regression insists are strongly related. This is the same family of self-deception that Section 2.6 warned about for evaluation: there, a leaky split let the future grade the past; here, a non-stationary series lets a spurious trend masquerade as a signal. Both are ways a temporal analysis reports a number it cannot deliver, and both are defended by the same discipline of checking the assumption before trusting the result. The connection runs deeper, because a backtest over a non-stationary series can look honest fold by fold and still mislead, since the regime each fold learns is not the regime the next fold faces.

3. The Random Walk and Unit Roots Intermediate

The canonical non-stationary process is the random walk,

where each step adds a fresh zero-mean shock $\varepsilon_t$ to the previous value. It is the simplest possible model of a quantity with no anchor: a price, an accumulated position, a drifting sensor bias. As Figure 3.2.1 contrasts, such a walk has no spring pulling it home, unlike a stationary series tethered to its mean. Unrolling the recursion from a fixed start $X_0$ reveals what the walk really is,

a running sum of every shock the process has ever received. No shock is ever discounted or forgotten; each one is added to the level and stays there forever. This is the defining feature of a unit-root process: shocks are permanent. The variance follows immediately, since the shocks are independent,

which grows without bound as $t$ increases. A process whose variance depends on $t$ cannot be weakly stationary, and the violation is structural, not a matter of degree: the spread of the random walk widens like $\sqrt{t}$ with no equilibrium it tends toward.

The general lens for this is the lag operator $L$, defined by $L X_t = X_{t-1}$. The random walk written with $L$ is $(1 - L)X_t = \varepsilon_t$, and more generally a first-order autoregression is

The behavior of the process is governed entirely by the root of the characteristic polynomial $1 - \phi z = 0$, which sits at $z = 1/\phi$. The stationarity condition is that this root lies strictly outside the unit circle, $|z| > 1$, equivalently $|\phi| < 1$. When $|\phi| < 1$ the process is a stationary AR(1): it has a fixed mean of zero, a finite variance $\sigma^2/(1-\phi^2)$, and a memory that decays geometrically, so a shock today fades as $\phi^h$ after $h$ steps and the series reverts to its mean. When $\phi = 1$ the root sits exactly on the unit circle at $z = 1$, the polynomial $1 - \phi z$ has a unit root, and the process is the random walk, with permanent shocks and no mean reversion. The knife edge between a forgetful, mean-reverting world and a perfect-memory, wandering one is the single value $\phi = 1$.

Code 3.2.1 simulates the three regimes and confirms numerically that the unit-root variance grows linearly with $t$ while the stationary variances settle to their theoretical ceilings. We average over many independent paths to estimate the variance at each time, which is the ensemble view that subsection one said we usually cannot afford but can always simulate.

import numpy as np

def ar1_paths(phi, T, n_paths, sigma=1.0, seed=0):
    """Simulate n_paths of an AR(1): X_t = phi * X_{t-1} + eps_t, starting at 0."""
    rng = np.random.default_rng(seed)
    eps = rng.normal(0.0, sigma, size=(n_paths, T))   # one shock per step per path
    X = np.zeros((n_paths, T))
    for t in range(1, T):
        X[:, t] = phi * X[:, t - 1] + eps[:, t]        # the AR(1) recursion
    return X

T, N = 200, 20000
for phi in (0.5, 0.9, 1.0):
    X = ar1_paths(phi, T, N)
    var_T = X[:, -1].var()                             # empirical Var(X_T) across paths
    if phi < 1:
        theory = 1.0 / (1.0 - phi**2)                  # stationary AR(1) variance ceiling
        print(f"phi={phi}: Var(X_T)~{var_T:6.2f}  theory={theory:6.2f}  (stationary)")
    else:
        print(f"phi={phi}: Var(X_T)~{var_T:6.2f}  theory={T:6.2f}  (= T*sigma^2, grows)")

phi=0.5: Var(X_T)~  1.34  theory=  1.33  (stationary)
phi=0.9: Var(X_T)~  5.27  theory=  5.26  (stationary)
phi=1.0: Var(X_T)~199.4   theory=200.00  (= T*sigma^2, grows)

4. Testing for Unit Roots and Stationarity Advanced

Eyeballing a plot is not enough, because a persistent stationary AR(1) with $\phi = 0.97$ and a true random walk look almost identical over a few hundred points. Two formal tests dominate practice, and they are most powerful when read together because their null hypotheses are opposite. The Augmented Dickey-Fuller (ADF) test takes the null hypothesis that a unit root is present (the series is non-stationary). It estimates the regression

where the lagged differences ($p$ of them) are the "augmentation" that soaks up short-run autocorrelation. The coefficient of interest is $\gamma$: a unit root corresponds to $\gamma = 0$ (rewriting $X_t = \phi X_{t-1} + \dots$ gives $\gamma = \phi - 1$, so $\phi = 1 \Leftrightarrow \gamma = 0$). The test statistic is the t-ratio on $\hat\gamma$, but it does not follow a normal distribution under the null; it follows the non-standard Dickey-Fuller distribution, so the critical values are more negative than the usual ones. A large negative ADF statistic (below the critical value, small p-value) rejects the unit root, evidence the series is stationary.

The Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test reverses the burden of proof: its null hypothesis is that the series is stationary (around a level or a deterministic trend), and the alternative is a unit root. It decomposes the series into a deterministic trend, a random walk, and a stationary error, and tests whether the random-walk component has zero variance. A large KPSS statistic (small p-value) rejects stationarity. Because ADF and KPSS frame the question from opposite sides, crossing them gives a four-way diagnosis, summarized in Table 3.2.1, and the agreement of two tests with opposite nulls is far more convincing than either alone.

The most vivid reason to take all of this seriously is spurious regression. Regress one independent random walk on another, two series with no causal or statistical link whatsoever, and ordinary least squares will routinely report a large coefficient, a high $R^2$, and a t-statistic so large that the relationship looks overwhelmingly significant. The trends wander, sometimes in loose visual sympathy, and OLS mistakes shared wandering for shared meaning. Code 3.2.2 demonstrates this from scratch over many trials and shows how often the naive regression is fooled.

import numpy as np

def spurious_trial(T, rng):
    """Regress one random walk on an INDEPENDENT random walk; return |t-stat| and R^2."""
    x = np.cumsum(rng.normal(size=T))          # random walk 1
    y = np.cumsum(rng.normal(size=T))          # random walk 2, fully independent of x
    X = np.column_stack([np.ones(T), x])       # design matrix [intercept, x]
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    resid = y - X @ beta
    s2 = resid @ resid / (T - 2)               # residual variance
    se = np.sqrt(s2 * np.linalg.inv(X.T @ X)[1, 1])   # std error of the slope
    t_stat = beta[1] / se
    ss_tot = ((y - y.mean())**2).sum()
    r2 = 1.0 - (resid @ resid) / ss_tot
    return abs(t_stat), r2

rng = np.random.default_rng(1)
trials = 2000
ts, r2s = zip(*(spurious_trial(200, rng) for _ in range(trials)))
ts, r2s = np.array(ts), np.array(r2s)

# Under a TRUE null with stationary data, |t| > 1.96 should happen ~5% of the time.
print(f"fraction with |t| > 1.96 : {(ts > 1.96).mean():.1%}   (should be ~5%)")
print(f"median R^2 of nonsense  : {np.median(r2s):.3f}")

fraction with |t| > 1.96 : 74.8%   (should be ~5%)
median R^2 of nonsense  : 0.235

5. Achieving Stationarity: Differencing and the Order of Integration Intermediate

If a unit root is the disease, differencing is the cure. The first difference operator is $\Delta X_t = X_t - X_{t-1} = (1-L)X_t$. Applying it to the random walk gives $\Delta X_t = \varepsilon_t$, white noise, which is stationary. A series that becomes stationary after differencing $d$ times, but not after $d-1$, is called integrated of order $d$, written $X_t \sim I(d)$. A stationary series is $I(0)$; a random walk is $I(1)$; a series whose changes themselves wander (an accelerating trend) may be $I(2)$ and need differencing twice. The integration order $d$ is exactly the $d$ in ARIMA($p,d,q$): it counts how many differences turn the series into something the stationary AR and MA machinery of Chapter 5 can model.

Two refinements matter in practice. Seasonal differencing removes a periodic unit root with the operator $\Delta_m X_t = X_t - X_{t-m}$, subtracting the value one full season ago (for monthly data with yearly seasonality, $m = 12$), and is the right tool when the non-stationarity repeats on a fixed calendar cycle rather than wandering freely. Variance-stabilizing transforms address a different failure: when the spread of a series grows with its level (multiplicative seasonality, exponential-looking growth), no amount of differencing fixes the changing variance, but a logarithm or a Box-Cox transform does, turning multiplicative structure into additive structure that differencing can then handle. The standard recipe is therefore transform first to stabilize variance, then difference to stabilize the mean.

Code 3.2.3 implements first differencing and the integration-order search from scratch, then re-tests with an ADF call from statsmodels as the library half of the pair. The from-scratch routine repeatedly differences until a hand-rolled ADF-style check stops rejecting stationarity, and the library version collapses the whole detect-and-difference cycle into a couple of calls.

import numpy as np
from statsmodels.tsa.stattools import adfuller

# ---- from scratch: difference until stationary, counting the integration order ----
def difference(series):
    """First difference: returns X_t - X_{t-1}, one element shorter."""
    return series[1:] - series[:-1]

def integration_order(series, max_d=3, alpha=0.05):
    """Find the smallest d such that the d-times-differenced series tests stationary."""
    s = np.asarray(series, dtype=float)
    for d in range(max_d + 1):
        pvalue = adfuller(s, autolag="AIC")[1]   # ADF p-value on the current series
        if pvalue < alpha:                       # reject unit root -> stationary
            return d, pvalue
        s = difference(s)                        # not stationary: difference once more
    return max_d, pvalue

rng = np.random.default_rng(7)
rw = np.cumsum(rng.normal(size=400))             # an I(1) random walk

d, p = integration_order(rw)
print(f"detected integration order d = {d}  (ADF p after differencing = {p:.4f})")

# sanity: the random walk itself is non-stationary, its first difference is not
print(f"ADF p-value, levels        : {adfuller(rw)[1]:.4f}  (>0.05 -> unit root)")
print(f"ADF p-value, 1st difference: {adfuller(difference(rw))[1]:.4f}  (<0.05 -> stationary)")

detected integration order d = 1  (ADF p after differencing = 0.0000)
ADF p-value, levels        : 0.6312  (>0.05 -> unit root)
ADF p-value, 1st difference: 0.0000  (<0.05 -> stationary)