library(forecast); library(ggplot2)
set.seed(123)
y_stat <- arima.sim(n = 1000, list(ar = 0.5))
autoplot(y_stat) + theme_bw(base_size = 14) +
labs(title = "Stationary AR(1), phi = 0.5", y = "")Module 2: Testing for Stationarity
Econ 6376 — Applied Time Series Econometrics
How to Use These Notes
This chapter has one main reading path and three optional layers.
- Core material is the main text. It contains the concepts, notation, interpretations, and applied skills expected of everyone.
- Deeper Dive sections explain why a result works or develop it more fully. They are useful, but they can be skipped on a first reading.
- Technical Note sections state qualifications or formal details that matter for precise reasoning.
- Looking Ahead sections introduce an idea that will be taught formally in a later module.
If you missed lecture, read the main text, run the Core code, and complete the Core Practice problems. Then return to the optional sections that address your questions or interests.
Core Learning Objectives
By the end of this module, you should be able to:
- Apply the difference operator \(\Delta\) and compute first and second differences both algebraically and in R.
- Conduct visual inspection of a time series and its ACF to form a preliminary judgment about stationarity.
- Build the Dickey-Fuller test regression from the AR(1) equation and explain why it takes that specific form.
- Distinguish the three ADF specifications (no constant, drift only, drift + trend) and choose the appropriate one for a given series.
- Explain why the DF \(t\)-statistic does not follow a \(t\)-distribution under the null and demonstrate this via Monte Carlo.
- Use MacKinnon critical values to make a test decision and recognize the role of sample size.
- Recognize KPSS as a reverse-null test and use it as a complement to ADF (confirmatory analysis).
- Apply the rule-of-thumb ratio \(\sigma_{\Delta y} / \sigma_y\) and understand its underlying logic.
- Define the order of integration and recognize the over-differencing problem when a trend-stationary series is differenced.
- Synthesize multiple pieces of evidence into a stationarity judgment, erring toward differencing when in doubt.
2.1 Setup: Why We Are Doing This
Last module we ended on spurious regression — two completely independent random walks producing a beautiful \(R^2\) and a fake significant coefficient. Granger and Newbold’s (1974) result is the entire reason this module exists. If you ignore non-stationarity, your regression output will lie to you with confidence. The first defense against that is detecting non-stationarity before you run any regression that depends on the time series properties of your data.
Recall the boundary from Module 1:
- \(|\phi| < 1\): stationary, finite variance, well-behaved inference.
- \(\phi = 1\): random walk, time-dependent and unbounded variance, spurious regression with anything else.
- \(\phi = 0.95\) vs \(\phi = 1.0\): nearly indistinguishable in the small samples we typically work with.
Where we are on the mixing board. Module 2 turns on no new sliders in the master equation:
\[y_t = \alpha + \delta t + \sum_{j=1}^{p} \phi_j y_{t-j} + \sum_{l=1}^{q} \theta_l \epsilon_{t-l} + \epsilon_t.\]
The DGP is still Module 1’s AR(1) — \(\alpha\) and \(\phi_1\) are the only active channels. What changes is the question we ask of it: is \(\phi_1 = 1\)? Along the way, the deterministic-trend slider \(\delta t\) makes its first cameo, because deciding whether a trending series is trend-stationary or difference-stationary turns out to be part of the same testing problem.
Today’s question is concrete: given a time series, which side of the boundary are we on? We will develop several different sources of evidence — visual inspection, the autocorrelation function, formal statistical tests, and a rule of thumb — and combine them into a workflow that produces a defensible verdict. None of these tools is conclusive on its own, especially in the borderline cases. The discipline of stationarity analysis is using them together and being honest about what you don’t know.
2.2 The Difference Operator
Before we test for stationarity, we need the tool that will appear repeatedly in this module and throughout the course: the difference operator.
First Difference
The first difference of \(y_t\) is the change from one period to the next: \[\Delta y_t = y_t - y_{t-1}\]
In the lag operator notation from Module 1, \(\Delta = (1 - L)\). The difference operator is the discrete-time analog of a derivative — it measures the instantaneous rate of change of the series. When we say a series is “growing at 2% per year,” we are making a statement about its first difference (or, more precisely, its log difference).
Second Difference
The second difference is the difference of the first difference: \[\Delta^2 y_t = \Delta(\Delta y_t)\]
Apply the inner \(\Delta\) first: \[\Delta^2 y_t = \Delta(y_t - y_{t-1}) = (y_t - y_{t-1}) - (y_{t-1} - y_{t-2})\]
Combining terms: \[\Delta^2 y_t = y_t - 2y_{t-1} + y_{t-2}\]
The lag operator gives the same result more cleanly: \[\Delta^2 y_t = (1 - L)^2 y_t = (1 - 2L + L^2) y_t = y_t - 2y_{t-1} + y_{t-2}\]
The analogy to calculus is exact: the first difference is like a first derivative (rate of change), and the second difference is like a second derivative (acceleration — change in the rate of change). For most economic time series you encounter, you almost never need beyond a second difference. We will return to this in Section 2.10 when we define the order of integration.
Why This Matters Right Now
The difference operator does double duty in this module:
- It is the dependent variable in the Dickey-Fuller test regression we will build in Section 2.5.
- It is the tool we use to make a non-stationary series stationary, once we have detected non-stationarity (Section 2.10).
The algebra needs to be comfortable before we go further.
2.3 Building Evidence I: Visual Inspection
The first source of evidence is the cheapest: look at the data. Visual inspection cannot give you a definitive answer about stationarity, but it can rule out the easy cases and flag the hard ones.
What Stationarity Looks Like
What should a stationary series look like? Recall the three conditions from Module 1:
- Constant mean across time.
- Constant variance across time.
- Autocovariance that depends only on the lag, not on when you measure it.
Visually, this means: if you take any two non-overlapping windows of the series, the means should be similar, and the spreads (standard deviations) should be similar. The series will have its ups and downs — a stationary series can be highly volatile — but the envelope of its movement should be roughly the same throughout the sample.
Compare this to a random walk:
set.seed(123)
y_rw <- cumsum(rnorm(1000))
autoplot(ts(y_rw)) + theme_bw(base_size = 14) +
labs(title = "Random Walk", y = "")The stationary series oscillates around a fixed mean (zero) with a roughly consistent envelope. The random walk wanders — there is no level the series returns to, and the spread of values you see in the second half of the sample looks different from the first half. This is the visual fingerprint of a unit root: the series goes places and stays there.
The Trap: Trend-Stationary Series
Visual inspection alone is not enough. Consider a series with a deterministic trend on top of stationary AR(1) deviations: \[y_t = 0.05\, t + u_t, \qquad u_t = 0.5\, u_{t-1} + \epsilon_t\]
set.seed(456)
n <- 300
# Build the DGP literally: stationary AR(1) deviations u_t ...
u <- numeric(n)
u[1] <- 0
for (t in 2:n) {
u[t] <- 0.5 * u[t-1] + rnorm(1) # u_t = 0.5 u_{t-1} + epsilon_t
}
# ... riding on a deterministic trend: y_t = 0.05 t + u_t
y_trend <- 0.05 * (1:n) + u
autoplot(ts(y_trend)) + theme_bw(base_size = 14)This series clearly trends upward. A naive viewer would say “non-stationary” and reach for the difference operator. But this series is trend-stationary (the concept introduced in Module 1, Section 1.4): the deviations from the trend are stationary. The right treatment is to detrend (subtract \(\hat{\alpha} + \hat{\delta} t\)), not to difference. Visual inspection cannot distinguish trend-stationary from difference-stationary — both will look “trended” in a level plot.
This is why we need more tools.
Three Things to Look For
When you plot a series for the first time, ask yourself:
- Does the series wander without returning to a level? This suggests possible non-stationarity (a unit root or a structural break).
- Is the variance changing over time? Periods of high volatility followed by calm periods suggest variance non-stationarity. For now we will set this aside; we revisit it in Module 14 when we discuss ARCH/GARCH.
- Is there an obvious trend? If yes, the question becomes “trend-stationary or difference-stationary?” — and visual inspection alone cannot answer it.
2.4 Building Evidence II: The Autocorrelation Function
The second source of evidence is the ACF, which we introduced in Module 1. Recall that for a stationary AR(1), \(\rho_k = \phi^k\) — the autocorrelation decays geometrically at rate \(\phi\). This decay pattern is what we look for in the correlogram.
Patterns to Recognize
- Fast geometric decay: the ACF drops quickly toward zero. The bars are inside the confidence bands within a handful of lags. Suggests stationarity with short-to-moderate memory.
- Slow but visible decay: the bars decline gradually but are still inside the bands by 15-20 lags. Suggests stationarity with high persistence (large \(\phi\), perhaps 0.85-0.95).
- Persistent / no decay: the bars stay near 1.0 for many lags, often dozens. The decay looks linear rather than exponential, or there is barely any decay at all. Strong evidence of a unit root.
The visual difference between “very slow decay” and “no decay” is exactly the \(\phi = 0.95\) vs \(\phi = 1.0\) problem from Module 1. In small samples, the ACF often cannot tell them apart.
(A note on packages: patchwork sits outside the core course package stack. It does no modeling — it only arranges several ggplots into one figure, here and in the Module 2 application below — so we use it as a stated display-only convenience.)
library(forecast); library(patchwork)
set.seed(123)
phis <- c(0.3, 0.7, 0.95, 1.0)
plots <- lapply(phis, function(p) {
if (p < 1) {
y <- arima.sim(n = 500, list(ar = p))
} else {
y <- cumsum(rnorm(500))
}
ggAcf(y, lag.max = 30) + theme_bw(base_size = 12) +
ggtitle(paste("phi =", p))
})
wrap_plots(plots, ncol = 2)For \(\phi = 0.3\) the ACF dies out by lag 5. For \(\phi = 0.7\) it is gone by lag 12 or so. For \(\phi = 0.95\) it stretches over 25+ lags but is gradually declining. For \(\phi = 1\) (random walk) it barely moves.
“Build a Case” Framing
We now have two sources of evidence — visual inspection of the series and visual inspection of the ACF. Neither is conclusive. Both can be misled by trend-stationary series, near-unit-root processes, and small samples.
The framing for the rest of the module is this: stationarity is a verdict, not a measurement. No single test gives you the answer. We collect multiple lines of evidence and weigh them together. When the evidence is consistent, we have a clear verdict. When it conflicts, we have to make a judgment call, and we will lean on the cost asymmetry (Section 2.10): spurious regression is more dangerous than over-differencing.
2.5 The Dickey-Fuller and Augmented Dickey-Fuller Tests
The third source of evidence is the formal statistical test. The most widely used test for a unit root is the Dickey-Fuller (DF) test and its augmented version, the ADF.
What We Are Testing
Return to the AR(1) data generating process: \[y_t = \alpha + \phi y_{t-1} + \epsilon_t\]
We want to test:
- \(H_0: \phi = 1\) (unit root, non-stationary)
- \(H_1: \phi < 1\) (stationary)
This is a one-sided test. We do not consider \(\phi > 1\) because explosive economic series essentially never occur (a country whose output literally doubled every year would not be a useful object of study).
Why Standard OLS Inference Fails
The naive approach is: estimate the AR(1) by OLS, get \(\hat{\phi}\), compute its standard error, form a \(t\)-statistic for \(H_0: \phi = 1\), and compare to standard \(t\)-tables. This fails. Under \(H_0\), the regressor \(y_{t-1}\) is itself a random walk — non-stationary. Standard OLS asymptotics require the regressors to be well-behaved (loosely: they have finite second moments and the regressor moment matrix converges to a positive-definite limit). A random walk’s variance grows linearly with \(t\), so this assumption is violated. As a result:
- \(\hat{\phi}\) is still consistent for \(\phi\), but the limit distribution is no longer normal.
- The standard error formula does not give a valid Wald-type statistic.
- The \(t\)-statistic exists, but the standard critical values do not apply.
This is the same conceptual problem as spurious regression in Module 1. Whenever non-stationary regressors enter an OLS regression, standard inference breaks down. The DF test is essentially a fix for this in the unit-root testing context, and the MacKinnon critical values we will introduce below are the empirical quantiles of the resulting non-standard distribution.
Building the DF Regression
The trick is to rearrange the AR(1) so that the dependent variable is stationary under the null. Subtract \(y_{t-1}\) from both sides: \[y_t - y_{t-1} = \alpha + \phi y_{t-1} - y_{t-1} + \epsilon_t\] \[\Delta y_t = \alpha + (\phi - 1) y_{t-1} + \epsilon_t\]
Define \(\gamma = \phi - 1\). The DF regression is: \[\Delta y_t = \alpha + \gamma y_{t-1} + \epsilon_t\]
The hypothesis test becomes:
- \(H_0: \gamma = 0\), equivalently \(\phi = 1\) (unit root)
- \(H_1: \gamma < 0\), equivalently \(\phi < 1\) (stationary)
Technical Note — Two different \(\gamma\)’s
The \(\gamma\) in the DF regression is a regression coefficient, \(\gamma = \phi - 1\). It is unrelated to the autocovariances \(\gamma_k\) from Module 1 (which always carry a lag subscript). The unit-root literature’s use of \(\gamma\) here is standard, so we keep it — just read the symbol in context.
Why this rearrangement? Two reasons:
The dependent variable is stationary under \(H_0\). Under the null, \(\Delta y_t = \epsilon_t\) — pure white noise. The left-hand side of the regression is stationary regardless of whether the null or alternative holds, which is at least somewhat better-behaved than having a non-stationary dependent variable. The regression is still subject to the non-standard asymptotics caused by the non-stationary regressor \(y_{t-1}\) on the right, but the left side is not the source of the problem.
The hypothesis maps directly to the parameter. The test statistic of interest is the \(t\)-stat on \(\hat{\gamma}\). If \(\hat{\gamma}\) is significantly less than zero, we reject the unit root. The interpretation is direct.
This \(t\)-statistic has its own conventional symbol: \(\tau\) (tau), the Dickey-Fuller test statistic. It is computed like an ordinary \(t\)-ratio but compared against DF critical values, not \(t\)-tables (the next two subsections show why). The name survives in software: urca::ur.df() labels its critical values tau2 for the drift specification and tau3 for the trend specification. Throughout this module, “the ADF statistic” and \(\tau\) are the same object.
The Dickey-Fuller test was designed in 1979 and the choice of regression form was deliberate — it gave Dickey and Fuller a tractable way to derive the limiting distribution of the test statistic, even though that distribution turned out to be non-standard.
The Augmented Version
The DF test as stated above is the simple version, and it has a problem: it assumes the errors \(\epsilon_t\) are white noise. If the true DGP is an AR(\(p\)) with \(p > 1\), the errors of the simple DF regression will have serial correlation, biasing the test. The fix is to include lagged differences as additional regressors: \[\Delta y_t = \alpha + \gamma y_{t-1} + \sum_{i=1}^{p} \beta_i \Delta y_{t-i} + \epsilon_t\]
These lagged differences soak up any short-run dynamics so that what is left is white noise. The test of interest is still \(H_0: \gamma = 0\). The augmented version is what people actually use in practice; the “ADF test” is the default. We will discuss how to choose \(p\) later in this section.
The Three ADF Specifications
This is where most applied mistakes happen. There are three versions of the ADF regression depending on what deterministic terms you include. Choosing the wrong specification can flip your conclusion.
Specification 1: No constant, no trend (“none”) \[\Delta y_t = \gamma y_{t-1} + \sum_{i=1}^{p} \beta_i \Delta y_{t-i} + \epsilon_t\] Use when the series has zero mean and no trend. This is rare in practice — almost all economic series have a non-zero mean. Even if they do have a non-zero mean the inclusion of a mean term is a generalization of this (though subject to estimation error).
Specification 2: Constant, no trend (“drift”) \[\Delta y_t = \alpha + \gamma y_{t-1} + \sum_{i=1}^{p} \beta_i \Delta y_{t-i} + \epsilon_t\] Use when the series has a non-zero mean but no obvious trend. This is the default for most stationary economic series — interest rates, unemployment rates, inflation rates, debt-to-GDP ratios.
Specification 3: Constant and trend (“trend”) \[\Delta y_t = \alpha + \delta t + \gamma y_{t-1} + \sum_{i=1}^{p} \beta_i \Delta y_{t-i} + \epsilon_t\] Use when the series has a clear deterministic trend. This specification tests whether the series is stationary around the trend (trend-stationary) or has a unit root with drift.
The critical values differ across specifications because including more deterministic terms changes the asymptotic distribution of the test statistic. Using critical values from the wrong specification is a real and common error.
The general guidance:
- Look at a plot of the series. Does it have a non-zero mean? A trend?
- Match the deterministic terms in the test to what you see.
- If unsure, err toward including more deterministic terms — but be aware that this reduces test power.
If you want a formal, defensible procedure for making this choice — and an explanation of the extra phi1/phi2/phi3 statistics that urca::ur.df() prints in its output — the two Deeper Dive sections below develop both. On a first reading, you can move straight to the Monte Carlo demonstration.
Deeper Dive — A Systematic Approach: The Pantula Principle
The “match the specification to the plot” guidance above is informal and practical. For students who want a more systematic procedure, Enders (Ch. 4, pp. 215-218) describes the Pantula principle, which is essentially a general-to-specific algorithm for choosing the right ADF specification.
The motivation: including unnecessary deterministic terms reduces test power, while omitting necessary ones can bias the test toward non-rejection. The Pantula principle resolves this trade-off by starting with the most general specification, testing for the unit root, and then testing whether the deterministic terms can be dropped — if they can, you re-test with a simpler (and more powerful) specification.
The procedure:
Start with the trend specification: \[\Delta y_t = \alpha + \delta t + \gamma y_{t-1} + \sum_{i=1}^{p} \beta_i \Delta y_{t-i} + \epsilon_t\] Run the ADF test for \(H_0: \gamma = 0\).
If you reject the unit root at this stage, the procedure can stop — the series is stationary (possibly around a trend), and you have your answer. But you should also check whether the trend term is needed: if \(\hat{\delta}\) is statistically insignificant (using a standard \(t\)-test, since under \(H_1\) standard inference is back in force), the deterministic trend may be redundant. Re-fit without the trend if you want a tidier model, but the unit root verdict already stands.
If you fail to reject the unit root at the trend specification, you cannot yet conclude there is a unit root. Low power may be the issue — and the trend specification is the lowest-power version. You need to test whether you should drop the trend term to gain power. Use the joint \(F\)-test \(\Phi_3\) (defined in the next Deeper Dive) for the null \(H_0: \gamma = 0\) and \(\delta = 0\) jointly.
If \(\Phi_3\) rejects but the marginal \(t\)-test on \(\gamma\) did not, the rejection is being driven by the trend, not the unit root. Conclude: trend-stationary.
If \(\Phi_3\) does not reject, drop the trend and re-run with the drift specification: \[\Delta y_t = \alpha + \gamma y_{t-1} + \sum_{i=1}^{p} \beta_i \Delta y_{t-i} + \epsilon_t\] This has more power. Test for the unit root again.
If the drift specification rejects, conclude: stationary around a non-zero mean. Done.
If the drift specification fails to reject, use the joint \(F\)-test \(\Phi_1\) for \(H_0: \gamma = 0\) and \(\alpha = 0\) jointly. If \(\Phi_1\) also fails to reject, drop the constant and run the no-constant specification. (In practice this last step rarely matters for economic data, since most series have non-zero means.)
The Pantula principle is a defensible way to choose specifications without ad hoc judgment calls. It is more robust than eyeballing the plot, especially for students who are not yet comfortable trusting their visual intuition. For most applied work, however, looking at the plot and matching the specification gives the same answer faster.
Reference: Enders, Ch. 4, pp. 215-218; Pantula, S.G. (1989), “Testing for Unit Roots in Time Series Data,” Econometric Theory, 5, 256-271.
Deeper Dive — The Joint F-Tests: \(\Phi_1\), \(\Phi_2\), \(\Phi_3\)
When you run urca::ur.df(...) and look at the summary output, you will see joint test statistics in addition to the main DF \(t\)-statistic on z.lag.1. With type = "drift", you get phi1; with type = "trend", you get phi2 and phi3. These are joint \(F\)-tests from Dickey & Fuller’s original (1981) paper that test combined hypotheses about the unit root and the deterministic terms together.
The three statistics are defined relative to the trend specification: \[\Delta y_t = \alpha + \delta t + \gamma y_{t-1} + \sum_{i=1}^{p} \beta_i \Delta y_{t-i} + \epsilon_t\]
| Statistic | Joint null hypothesis | Interpretation |
|---|---|---|
| \(\Phi_1\) | \(\alpha = 0, \gamma = 0\) (in drift specification) | Random walk with no drift vs. stationary AR(1) with intercept |
| \(\Phi_2\) | \(\alpha = 0, \delta = 0, \gamma = 0\) (in trend spec) | Random walk with no drift vs. trend-stationary process |
| \(\Phi_3\) | \(\delta = 0, \gamma = 0\) (in trend specification) | Random walk with drift vs. trend-stationary process |
These are computed as standard \(F\)-statistics from the regression sum of squares, but — like the marginal DF \(t\)-statistic — they do not follow standard \(F\)-distributions under their nulls. The non-stationarity of \(y_{t-1}\) under the unit root null means the joint distribution of \((\hat{\alpha}, \hat{\delta}, \hat{\gamma})\) is non-standard. Dickey and Fuller (1981) tabulated the critical values; modern software (including urca) reports both the test statistic and the critical values automatically.
How they fit into the Pantula principle. The joint \(F\)-tests are how you formally check whether the deterministic terms can be dropped in step 3 of the procedure above. If \(\Phi_3\) rejects but the marginal DF \(t\)-statistic does not, the rejection is being driven by the trend term — implying trend stationarity rather than a unit root. If \(\Phi_3\) does not reject, you can drop the trend and gain power by moving to the simpler specification.
Practical advice. For most applied work, you do not need to compute \(\Phi_1\), \(\Phi_2\), \(\Phi_3\) by hand or even interpret them in detail. Just be aware that when you run summary(ur.df(y, type = "trend")), the additional test statistics in the output are these joint tests. The critical values reported alongside them tell you whether to reject the joint null. If you are following the Pantula principle, \(\Phi_3\) is the one you care about most.
References: Dickey, D.A. and Fuller, W.A. (1981), “Likelihood Ratio Statistics for Autoregressive Time Series with a Unit Root,” Econometrica, 49, 1057-1072; Enders, Ch. 4, p. 213.
The Non-Standard Distribution Under H₀: A Monte Carlo Demonstration
We have asserted that the DF \(t\)-statistic does not follow a \(t\)-distribution under the null. Rather than derive this (the derivation involves functional Brownian motion arguments — see Phillips, 1987, or Hamilton Ch. 17 if you are curious), let us show it with a Monte Carlo simulation. This is this module’s central assumption-break: we take the machinery a standard \(t\)-test trusts and watch it fail under the unit-root null.
The strategy: generate many independent random walks (so \(H_0\) is true by construction), run the DF regression on each, and collect the resulting \(t\)-statistics. The empirical distribution of these statistics is the Dickey-Fuller distribution.
set.seed(42)
n_sims <- 5000
T <- 200
df_stats <- replicate(n_sims, {
# Generate a random walk under H_0: phi = 1
y <- cumsum(rnorm(T))
dy <- diff(y)
y_lag <- y[-T]
# Run the DF regression (drift / intercept specification)
reg <- lm(dy ~ y_lag)
# Extract the t-statistic on y_lag
summary(reg)$coefficients["y_lag", "t value"]
})
# Summary statistics
mean(df_stats) # Should be around -1.5, not 0[1] -1.546135
sd(df_stats) # Should be around 0.86, not 1[1] 0.8525777
quantile(df_stats, 0.05) # Should be around -2.86, not -1.65 5%
-2.892739
Now plot the resulting distribution alongside the standard normal density:
library(ggplot2)
ggplot(data.frame(stat = df_stats), aes(x = stat)) +
geom_histogram(aes(y = after_stat(density)), bins = 60,
fill = "steelblue", alpha = 0.6) +
stat_function(fun = dnorm, color = "red", linewidth = 1.2) +
geom_vline(xintercept = -1.645, linetype = "dashed",
color = "red", linewidth = 0.8) +
geom_vline(xintercept = quantile(df_stats, 0.05),
linetype = "dashed", color = "steelblue", linewidth = 0.8) +
theme_bw(base_size = 14) +
labs(title = "DF distribution (blue) vs standard normal (red)",
subtitle = paste0("5% critical values: DF = ",
round(quantile(df_stats, 0.05), 2),
" | Normal = -1.65"),
x = "t-statistic under H_0: phi = 1",
y = "Density")What you see: The DF distribution is shifted to the left — its mean is about \(-1.5\), not zero. It is not symmetric. The 5% left-tail critical value is about \(-2.86\) for a sample of size \(T = 200\) with the drift specification, not the \(-1.65\) you would use for a standard normal one-sided test.
The implication. If you used the wrong critical value — say, \(-1.65\) from a standard normal — you would reject the unit root null roughly 46% of the time when it is actually true, instead of the nominal 5%. Your test would be massively oversized. The MacKinnon critical values we present below are not arbitrary numbers; they are the empirical quantiles of exactly this distribution, computed once carefully so we do not have to redo the simulation every time.
This is the same problem as spurious regression: standard OLS inference breaks down when regressors are non-stationary. The DF test does not solve that problem; it works around it by using the correct critical values for the resulting non-standard distribution.
MacKinnon Critical Values
MacKinnon (1996) computed the critical values of the DF distribution for many sample sizes and specifications, and provided response surface formulas so you do not need lookup tables. For the drift specification (Specification 2), the critical values are approximately:
\[ \begin{aligned} \text{CV}_{0.01} &= -3.43035 - \frac{6.5393}{T} - \frac{16.786}{T^2} - \frac{79.433}{T^3} \\ \text{CV}_{0.05} &= -2.86154 - \frac{2.8903}{T} - \frac{4.234}{T^2} - \frac{40.04}{T^3} \\ \text{CV}_{0.10} &= -2.56677 - \frac{1.5384}{T} - \frac{2.809}{T^2} \end{aligned} \]
(The subscript is the significance level; we write \(\text{CV}\) rather than \(\alpha\) to avoid a collision with the intercept \(\alpha\) in the test regression.) Notice that the critical values depend on \(T\). For small samples they are more extreme; as \(T \to \infty\) they converge to their asymptotic values. The corrections shrink rapidly as \(T\) grows — for \(T > 500\) they are negligible.
Reference: MacKinnon, J.G. (1996), “Numerical Distribution Functions for Unit Root and Cointegration Tests,” Journal of Applied Econometrics, 11, 601-618. Enders (Ch. 4, p. 215) also tabulates these values.
A simple R function:
mackinnon_cv <- function(T, level = 0.05, spec = "drift") {
# Returns the MacKinnon critical value for the ADF test
# spec = "drift" gives the constant-only specification
if (spec != "drift") {
stop("Only drift specification implemented in this teaching example.
Use urca::ur.df() for other specifications.")
}
if (level == 0.01) {
return(-3.43035 - 6.5393/T - 16.786/T^2 - 79.433/T^3)
}
if (level == 0.05) {
return(-2.86154 - 2.8903/T - 4.234/T^2 - 40.04/T^3)
}
if (level == 0.10) {
return(-2.56677 - 1.5384/T - 2.809/T^2)
}
stop("Level must be 0.01, 0.05, or 0.10")
}
# Sanity check against the Monte Carlo result from above
mackinnon_cv(T = 200) # Returns ~-2.876[1] -2.876102
quantile(df_stats, 0.05) # The empirical 5% quantile from the MC 5%
-2.892739
The two values should be very close — MacKinnon’s formula is a smoothed and refined version of exactly the kind of Monte Carlo we just ran, with vastly more replications and across many sample sizes.
For the trend specification (Spec 3), the critical values are more negative still (around \(-3.41\) at the 5% level for large \(T\)). For the no-constant specification (Spec 1), they are less negative (around \(-1.95\) at the 5% level). The full set of formulas is in MacKinnon’s paper. In practice, the urca::ur.df() function returns the critical values automatically for whichever specification you chose.
(This function also lives in helpers/testing.R as the canonical course copy, so later modules can source() it instead of re-defining it.)
Lag Selection in the ADF
How many lagged differences (\(p\)) should you include in the ADF regression?
- Too few: residuals will still have serial correlation, biasing the test.
- Too many: you lose degrees of freedom and reduce power.
Common approaches:
Information criteria: minimize AIC or BIC over a range of candidate \(p\) values. BIC tends to choose smaller models and is the more conservative choice. In
urca::ur.df(), setselectlags = "BIC"or"AIC". We will talk more about information criteria such as AIC and BIC later.Sequential testing (top-down): start with a maximum lag \(p_{\max}\) and drop the highest insignificant lag. Schwert (1989) suggested: \[p_{\max} = \left\lfloor 12 \cdot \left(\frac{T}{100}\right)^{1/4} \right\rfloor\] For \(T = 100\), this gives \(p_{\max} = 12\). For \(T = 200\), \(p_{\max} \approx 14\). Then test the lagged-difference coefficients sequentially and drop the top one if insignificant, refitting until the highest remaining lag is significant.
Default: For most applied work, BIC selection over a range like \(0\) to \(\lfloor T^{1/3} \rfloor\) is fine. We will use BIC as the default in the problem sets.
Building the ADF in R: From Scratch and From the Package
Live coding the test “from scratch” demystifies what urca::ur.df() is doing under the hood. Recall the embed() function from Module 1.
library(urca); library(forecast)
set.seed(8675309)
y <- arima.sim(n = 200, list(ar = 0.7)) # Truly stationary
T <- length(y)
p_lags <- 4
# Construct the regression matrix manually
dy <- diff(y)
n <- length(dy)
X <- embed(dy, p_lags + 1) # [dy_t, dy_{t-1}, ..., dy_{t-p}]
dy_t <- X[, 1]
dy_lags <- X[, -1, drop = FALSE]
y_lag <- y[(p_lags + 1):(T - 1)] # y_{t-1} aligned with the truncated dy
# Run the regression
adf_reg <- lm(dy_t ~ y_lag + dy_lags)
summary(adf_reg)
Call:
lm(formula = dy_t ~ y_lag + dy_lags)
Residuals:
Min 1Q Median 3Q Max
-2.76518 -0.58199 0.06084 0.58090 2.93186
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.002613 0.069252 -0.038 0.970
y_lag -0.324343 0.065915 -4.921 1.87e-06 ***
dy_lags1 0.066426 0.080404 0.826 0.410
dy_lags2 0.010452 0.078874 0.133 0.895
dy_lags3 0.086795 0.074705 1.162 0.247
dy_lags4 0.080070 0.072493 1.105 0.271
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9666 on 189 degrees of freedom
Multiple R-squared: 0.1501, Adjusted R-squared: 0.1276
F-statistic: 6.677 on 5 and 189 DF, p-value: 9.432e-06
# The test statistic is the t-stat on y_lag
adf_stat <- summary(adf_reg)$coefficients["y_lag", "t value"]
critical <- mackinnon_cv(T = T, level = 0.05)
cat("ADF statistic:", round(adf_stat, 3), "\n")ADF statistic: -4.921
cat("5% critical value:", round(critical, 3), "\n")5% critical value: -2.876
cat("Reject H_0:", adf_stat < critical, "\n")Reject H_0: TRUE
For this stationary AR(1) with \(T = 200\), you get an ADF statistic around \(-5\) (the exact value depends on the simulated draw), well below the critical value of \(-2.88\). We comfortably reject the unit root null in favor of stationarity.
The packaged version:
adf_test <- ur.df(y, type = "drift", lags = 4)
summary(adf_test)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression drift
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-2.76518 -0.58199 0.06084 0.58090 2.93186
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.002613 0.069252 -0.038 0.970
z.lag.1 -0.324343 0.065915 -4.921 1.87e-06 ***
z.diff.lag1 0.066426 0.080404 0.826 0.410
z.diff.lag2 0.010452 0.078874 0.133 0.895
z.diff.lag3 0.086795 0.074705 1.162 0.247
z.diff.lag4 0.080070 0.072493 1.105 0.271
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9666 on 189 degrees of freedom
Multiple R-squared: 0.1501, Adjusted R-squared: 0.1276
F-statistic: 6.677 on 5 and 189 DF, p-value: 9.432e-06
Value of test-statistic is: -4.9206 12.1109
Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
The output will show the same regression, the same test statistic (modulo small numerical differences in how the package handles intercepts), and the critical values for 1%, 5%, and 10% levels for the drift specification at the appropriate sample size.
For a real workflow, you would use the packaged version with selectlags = "BIC":
adf_test <- ur.df(y, type = "drift", lags = 12, selectlags = "BIC")
summary(adf_test)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression drift
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-2.57956 -0.64418 0.03028 0.67596 2.88927
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.00902 0.07172 -0.126 0.900
z.lag.1 -0.29040 0.05569 -5.215 4.93e-07 ***
z.diff.lag 0.02940 0.07417 0.396 0.692
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9802 on 184 degrees of freedom
Multiple R-squared: 0.1411, Adjusted R-squared: 0.1317
F-statistic: 15.11 on 2 and 184 DF, p-value: 8.396e-07
Value of test-statistic is: -5.2145 13.5965
Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
This searches for the best lag length up to \(p = 12\) using BIC, runs the test, and reports the results.
2.6 KPSS: The Reverse-Null Test
The Dickey-Fuller test has a particular structural feature: its null hypothesis is a unit root. In statistical testing, we can only reject the null, not accept it. So if ADF fails to reject the unit root, the correct interpretation is “we lack evidence against the unit root,” not “the series has a unit root.” Combined with the low power problem (Section 2.8), this creates an asymmetry — the test is biased toward concluding non-stationarity when the data are uninformative.
The KPSS test (Kwiatkowski, Phillips, Schmidt, and Shin, 1992) is the conceptual mirror: its null hypothesis is stationarity, and the alternative is a unit root.
| Test | \(H_0\) | \(H_1\) |
|---|---|---|
| ADF | unit root | stationary |
| KPSS | stationary | unit root |
Why Use Both
Because both tests have low power against close alternatives, neither one alone is fully informative. Used together they enable confirmatory analysis:
| ADF says | KPSS says | Conclusion |
|---|---|---|
| Reject unit root | Fail to reject stationarity | Strong evidence of stationarity |
| Fail to reject unit root | Reject stationarity | Strong evidence of a unit root |
| Reject unit root | Reject stationarity | Contradiction — possible structural break or misspecification |
| Fail to reject unit root | Fail to reject stationarity | Inconclusive — borderline case, defer to other evidence |
The diagonal cases (when both tests agree) give you high confidence. The off-diagonal cases tell you something is unusual and you need to think harder. The “inconclusive” case (both fail to reject) is itself useful information: it tells you honestly that the data are too noisy or the sample is too small to discriminate between stationarity and a unit root. In that case, you fall back on the cost asymmetry from Section 2.10.
Using KPSS: What Matters
You can use KPSS correctly knowing just three things:
- The null is stationarity, opposite of ADF.
- A large test statistic is evidence against stationarity (the opposite of how ADF works — this is a common point of confusion).
- Critical values come from the KPSS asymptotic distribution, tabulated in their paper and built into R packages.
The mechanics behind the statistic are in the Deeper Dive below; they are not essential for using the test.
Deeper Dive — How KPSS Works
The KPSS test decomposes the series into a deterministic component, a random walk component, and a stationary error: \[y_t = \xi_t + r_t + u_t, \quad r_t = r_{t-1} + v_t, \quad v_t \sim N(0, \sigma_v^2)\]
Under \(H_0\) (stationarity), \(\sigma_v^2 = 0\) — the random walk component does not exist. The KPSS statistic measures how much the partial sums of the residuals from regressing \(y_t\) on the deterministic component “wander.” If the residuals are stationary, the partial sums stay bounded; if they wander, the test rejects.
KPSS in R
library(urca)
# Level stationarity (no deterministic trend)
kpss_level <- ur.kpss(y, type = "mu")
summary(kpss_level)
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 4 lags.
Value of test-statistic is: 0.08
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739
# Trend stationarity (allow for a deterministic trend)
kpss_trend <- ur.kpss(y, type = "tau")
summary(kpss_trend)
#######################
# KPSS Unit Root Test #
#######################
Test is of type: tau with 4 lags.
Value of test-statistic is: 0.0795
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.119 0.146 0.176 0.216
The output shows the test statistic and the critical values at 10%, 5%, 2.5%, and 1% levels. If the statistic exceeds the critical value, you reject stationarity. Note again: this is the opposite direction of comparison from ADF, where you reject the null when the statistic is below the critical value. Get this confused at your peril.
Choosing Between mu and tau
- Use
type = "mu"when testing whether the series is stationary around a constant level (no trend in the deterministic component). This pairs naturally with the drift specification of ADF. - Use
type = "tau"when testing whether the series is stationary around a deterministic trend. This pairs naturally with the trend specification of ADF.
For the FRED unemployment series in Section 2.12, we will use type = "mu" because there is no obvious deterministic trend in unemployment.
2.7 Other Tests: Phillips-Perron and ERS
There is a small zoo of unit root tests. Two more deserve a brief mention — you need to recognize the names, not master the mechanics.
Phillips-Perron (PP): Like ADF, but uses a non-parametric correction for serial correlation in the errors instead of including lagged differences. Same null, similar critical values, generally similar conclusions. Available as urca::ur.pp(). Sometimes more powerful than ADF when the error structure is too complex for a small number of lags to capture, but in most cases the two give similar answers.
Elliott-Rothenberg-Stock (ERS / DF-GLS): A more powerful variant of the ADF test that uses generalized least squares (GLS) detrending instead of OLS detrending. ERS has higher power than ADF, particularly when the series has a non-zero mean or trend. Available as urca::ur.ers(). If you only want to remember one alternative to ADF as a “more powerful version,” this is the one.
There are other tests in the literature — Ng-Perron, Bierens, Im-Pesaran-Shin (for panels) — and active research continues. The point is not that you need to know them all. The point is that these are different tools for the same job, none of them solves the fundamental low-power problem, and they can be used as complements to ADF and KPSS when you want additional evidence.
2.8 Power, the Borderline Case, and Why We Use Multiple Tests
The single most important practical fact about unit root tests is that they have low power against highly persistent stationary alternatives.
To see what this means concretely, consider a stationary AR(1) with \(\phi = 0.95\) — clearly stationary, but very close to the boundary. How often does the ADF test correctly reject the unit root null?
set.seed(2024)
n_sims <- 1000
# Power at phi = 0.95 for various sample sizes
power_results <- sapply(c(50, 100, 200, 500), function(T) {
rejections <- replicate(n_sims, {
y <- arima.sim(n = T, list(ar = 0.95))
test <- ur.df(y, type = "drift", lags = 4)
test_stat <- test@teststat[1]
crit <- mackinnon_cv(T = T, level = 0.05)
test_stat < crit # TRUE = reject H_0 = correct decision
})
mean(rejections)
})
names(power_results) <- c("T=50", "T=100", "T=200", "T=500")
power_results T=50 T=100 T=200 T=500
0.059 0.085 0.267 0.922
Typical results (from the run above; your exact numbers will wobble with the seed):
- \(T = 50\): power ≈ 0.06 (barely above the nominal 5% level)
- \(T = 100\): power ≈ 0.09
- \(T = 200\): power ≈ 0.27
- \(T = 500\): power ≈ 0.92
Read these numbers carefully. With 100 observations and \(\phi = 0.95\), the ADF test correctly rejects the unit root less than 10% of the time. The truly stationary series is misclassified as a unit root roughly 90% of the time. This is not a defect of the ADF test specifically; KPSS, PP, ERS all suffer from analogous power problems near the boundary.
This is why we cannot rely on a single test. The borderline cases are exactly the cases where stationarity testing is most useful and most fragile. The cure is the “build a case” framework: collect multiple lines of evidence and make a defensible judgment. It is also why the cost asymmetry in Section 2.10 matters: when the evidence is genuinely uninformative, we fall back on the cost of being wrong in each direction.
2.9 The Rule of Thumb
The rule of thumb is the fourth and last source of evidence — fast, informal, and surprisingly informative.
The Rule
Taught to me by Jeff Mills:
- Compute \(\sigma_y\), the standard deviation of the series in levels.
- Compute \(\sigma_{\Delta y}\), the standard deviation of the first differences.
- Compute the ratio: \[R = \frac{\sigma_{\Delta y}}{\sigma_y}\]
- If \(R < 0.5\): differencing reduces the variance sharply. The series is probably non-stationary. Difference it.
- If \(R \geq 0.5\): differencing does not reduce variance much. The series is probably already stationary. Keep the levels.
Why It Works
The rule captures a simple intuition. For a stationary AR(1) with \(|\phi| < 1\), the variance of the levels is \(\sigma_\epsilon^2 / (1 - \phi^2)\), while the variance of the first differences is \(\text{Var}(\Delta y_t) = 2(\gamma_0 - \gamma_1) = 2\gamma_0(1 - \phi) = 2\sigma_\epsilon^2 / (1 + \phi)\). Take the ratio: \[\frac{\text{Var}(\Delta y_t)}{\text{Var}(y_t)} = \frac{2\sigma_\epsilon^2 / (1 + \phi)}{\sigma_\epsilon^2 / (1 - \phi^2)} = 2(1 - \phi)\]
Take square roots to get the ratio of standard deviations and evaluate at a few values:
- \(\phi = 0\) (white noise): \(R = \sqrt{2} \approx 1.41\)
- \(\phi = 0.5\): \(R = \sqrt{1.0} = 1.00\)
- \(\phi = 0.7\): \(R = \sqrt{0.6} \approx 0.77\)
- \(\phi = 0.9\): \(R = \sqrt{0.2} \approx 0.45\)
- \(\phi = 0.95\): \(R = \sqrt{0.1} \approx 0.32\)
- \(\phi = 1.0\) (random walk): \(R \to 0\)
Two observations:
- The cutoff \(R < 0.5\) corresponds roughly to \(\phi > 0.875\). Series with very high persistence get flagged for differencing.
- The rule cannot tell the difference between \(\phi = 0.95\) and \(\phi = 1.0\) — both produce small \(R\). This is the same low-power problem as ADF, viewed from a different angle.
When to Use It
The rule of thumb is:
- Fast: one line of R code.
- Easy to interpret: no critical values, no specifications, no choices.
- Not a formal test: you should not put it in a paper as your main evidence. Use it as a sanity check, especially when ADF and KPSS disagree.
ratio <- sd(diff(y)) / sd(y)
ratio[1] 0.7599642
Below 0.5? Difference. Above 0.5? Keep the levels. Treat the verdict as one piece of evidence among several.
2.10 Differencing, Order of Integration, and the Over-Differencing Trap
We have spent most of this module detecting non-stationarity. Now we address the resolution.
Differencing as the Resolution
If a series is non-stationary because of a unit root, the most common remedy is differencing. Recall that for a random walk \(y_t = y_{t-1} + \epsilon_t\): \[\Delta y_t = y_t - y_{t-1} = \epsilon_t\]
The first difference of a random walk is white noise — manifestly stationary. More generally, if a series has exactly one unit root, taking the first difference produces a stationary series. If it has two unit roots (rare), you need to difference twice. And so on.
Order of Integration
The order of integration of a series is the number of times you need to difference it to make it stationary. Notation: \(y_t \sim I(d)\) if the series requires \(d\) differences to become stationary.
- \(I(0)\): Already stationary. No differencing needed.
- \(I(1)\): Stationary after one difference. The most common case for non-stationary economic series — GDP, prices, asset prices, exchange rates in levels are all typically \(I(1)\).
- \(I(2)\): Stationary after two differences. Rare but occasionally encountered. The aggregate price level can be \(I(2)\) when both the level and the inflation rate are non-stationary — that is, when inflation itself has a unit root and needs to be differenced.
- \(I(d)\) for \(d > 2\): Almost never seen in practice. If your testing suggests you need \(d > 2\), something is probably wrong with your specification — perhaps you have a structural break being mistaken for high-order integration.
A useful analogy: how many licks does it take to get to the center of a Tootsie Pop? The order of integration is how many differences it takes to reach stationarity. We almost never need more than two licks.
The Practical Workflow
- Test the levels with ADF (and KPSS if you want confirmatory analysis).
- If you fail to reject the unit root, take a first difference and test again.
- In most applied economics work, we stop there unless we have reason to think the series could be \(I(2)\).
- If \(I(2)\) is plausible, do not just keep differencing upward. Start from the highest order you think is possible and test downward, which is the logic behind Dickey-Pantula style procedures.
# Simple helper for the common I(0)/I(1) case
find_order_i01 <- function(y, level = 0.05) {
cv_col <- if (level <= 0.01) {
1
} else if (level <= 0.05) {
2
} else {
3
}
test_level <- ur.df(y, type = "drift", lags = 4)
if (test_level@teststat[1] < test_level@cval[1, cv_col]) {
return(0)
}
test_diff <- ur.df(diff(y), type = "drift", lags = 4)
if (test_diff@teststat[1] < test_diff@cval[1, cv_col]) {
return(1)
}
return(NA) # Inconclusive; if I(2) is plausible, use downward Dickey-Pantula-style testing
}For the kind of series we usually see in applied economics, that simple two-step \(I(0)\)/\(I(1)\) check is often enough. If you think an \(I(2)\) process is on the table, use the downward-testing logic instead. (This helper also lives in helpers/testing.R, alongside the more general find_order() that tests up to a chosen maximum order.)
The Over-Differencing Trap
What happens if you difference a series that was actually trend-stationary, not difference-stationary? Suppose the true DGP is: \[y_t = \alpha + \delta t + u_t, \quad u_t \sim N(0, \sigma^2) \text{ i.i.d.}\]
This is a deterministic trend with white noise around it — the simplest possible trend-stationary process. The correct treatment is to detrend by regressing \(y_t\) on a constant and \(t\), then working with the residuals (which are stationary by construction).
Suppose instead we difference. Compute \(\Delta y_t\): \[\Delta y_t = y_t - y_{t-1}\] \[\Delta y_t = (\alpha + \delta t + u_t) - (\alpha + \delta (t-1) + u_{t-1})\] \[\Delta y_t = \delta + (u_t - u_{t-1})\]
The constant \(\alpha\) cancels (good — that’s not the problem). The linear trend \(\delta t - \delta(t-1) = \delta\) collapses to a constant (also good — the differenced series has the constant drift \(\delta\)). But the noise component is now \(u_t - u_{t-1}\), and this is not white noise.
To see what it is, write it in the master-equation form: \[\Delta y_t = \delta + u_t - u_{t-1}\]
If we relabel \(u_t\) as the new innovation \(\epsilon_t\), this is exactly an MA(1) process: \[\Delta y_t = \delta + \epsilon_t + \theta \epsilon_{t-1}, \quad \theta = -1\]
In lag operator form, the MA polynomial is: \[\Theta(L) = 1 + \theta L = 1 - L\]
The root of \(\Theta(L) = 0\) is \(L = 1\) — exactly on the unit circle. The MA polynomial has a unit root.
Why this is bad. A moving average process is invertible if all roots of its MA polynomial lie outside the unit circle. Invertibility is what allows us to recover the innovations from past observations of the series — it is essential for forecasting and for maximum likelihood estimation. When the MA polynomial has a unit root (it is non-invertible):
- The innovations cannot be recovered uniquely from the data.
- Maximum likelihood estimation will struggle. The likelihood surface is flat or ill-conditioned near the boundary, and estimation routines often pin \(\hat{\theta}\) exactly at \(-1\) (or report convergence failures).
- The model is not identified — different parameter combinations give the same likelihood.
- Forecasts behave poorly because the implicit forecasting weights do not converge.
In short: differencing a trend-stationary series replaces a deterministic trend with a non-invertible MA(1). You have not fixed a stationarity problem; you have substituted one pathology for another. The differenced series is technically stationary (its mean and variance are constant), but the model is malformed.
For now, the lesson is: trend-stationary and difference-stationary series require different treatments. Use the trend-augmented ADF specification to distinguish them when you suspect a deterministic trend might be present.
Looking Ahead — Invertibility and diagnostics
You have just met the vocabulary “unit root in the MA polynomial” ahead of schedule. Module 3 develops invertibility formally as the MA counterpart of AR stationarity — and demonstrates experimentally what a nearly non-invertible MA(1) does to a sample ACF and to arima() estimates. Module 5 (diagnostics) then returns to the practical symptom: a fitted ARMA model whose MA coefficient sits suspiciously close to \(-1\) is often telling you the series was over-differenced.
Difference vs. Detrend: Choosing Correctly
The correct treatment depends on what the series actually is:
- Difference-stationary (unit root, possibly with drift): difference. Removes the stochastic trend without introducing artifacts.
- Trend-stationary (deterministic trend with stationary deviations): detrend. Removes the deterministic trend without introducing a unit root in the MA polynomial.
How do you tell which? Use the trend-augmented ADF specification (Specification 3 from Section 2.5). It explicitly tests for a unit root in the presence of a deterministic trend:
- If you reject the unit root with the trend specification: the series is likely trend-stationary. Detrend.
- If you fail to reject the unit root even with the trend included: the series probably has a true unit root. Difference.
This is one of the few cases where the choice of ADF specification really determines your modeling strategy. Get it right.
The Cost Asymmetry, Revisited
We told you earlier that when in doubt, err toward differencing. That is still good advice — but now you know there is a real cost to over-differencing, too. The right framing is:
- Spurious regression (the cost of failing to difference a unit-root series) gives you confidently wrong answers with significant t-stats and high \(R^2\). You publish results that do not replicate.
- Over-differencing (the cost of differencing a trend-stationary series) introduces a non-invertible MA(1). Your coefficient estimates are not systematically biased, but the moving-average structure you just created in the errors means the usual standard errors can be substantially wrong until you account for it. Less catastrophic than spurious regression — but not harmless.
The asymmetry remains: spurious regression is the bigger inference risk. But you are now aware of the over-differencing tax, and you can use the trend-augmented ADF to make the choice deliberately rather than defaulting to differencing in all cases.
2.11 Synthesis: Building a Case
We now have four sources of evidence:
- Visual inspection of the series: Does it wander? Is the variance changing? Is there a trend?
- Visual inspection of the ACF: Does it decay quickly (stationary) or slowly/not at all (unit root)?
- Statistical tests: ADF (and PP, ERS) test $H_0 = $ unit root. KPSS tests $H_0 = $ stationary. Use them together when possible.
- Rule of thumb: \(\sigma_{\Delta y} / \sigma_y\). Quick informal check.
Combine them. When all four point the same way, the answer is clear. When they conflict, you have to make a judgment call.
The Decision Framework
- Plot the series. What is the gross structure? Trend? Drift? Volatility changes?
- Plot the ACF. Fast decay or slow decay?
- Choose the appropriate ADF specification based on what you saw in step 1. Run it.
- Run KPSS with the matching specification (
muif you ran ADF with drift,tauif you ran ADF with trend). - Compute the rule-of-thumb ratio as a sanity check.
- Synthesize. Do the four sources agree? If yes, you have your answer. If no, fall back on the cost asymmetry: difference unless you have strong reason not to.
- Remember the over-differencing trap. If you suspect a deterministic trend, detrend rather than difference. The trend-augmented ADF is your tool for distinguishing.
Why Multiple Pieces of Evidence?
The recurring theme: every test we have, every visual diagnostic, every rule of thumb has limited power against the cases that matter most — the borderline ones where \(\phi\) is close to 1 but not equal. No single tool can give you certainty in those cases. The discipline is to use multiple tools and accept that sometimes the honest answer is “I don’t know, but I’m going to difference anyway because the cost asymmetry favors that choice.”
This is not a weakness of time series analysis. It is the truth: extracting causal claims from passive observation of one realization of a stochastic process is hard, and pretending otherwise is how you get spurious regressions.
Practicing the Workflow
The decision framework above only becomes useful once it comes naturally — and that takes repetition on series where you don’t already know the answer. Two ways to get that practice before the problem set:
- Generate your own mystery series. Use the simulators you built in Module 1: draw \(\phi\) somewhere in \(\{0.5, 0.8, 0.9, 0.95, 1.0\}\) without looking, simulate, and run the full workflow — plot, ACF, the right ADF specification, matching KPSS, rule of thumb, verdict. Then check against the \(\phi\) you drew. The borderline draws are the ones that teach you the most.
- Work the Core Practice problems at the end of this module. They are built around exactly this workflow, and the power study in particular previews the kind of reasoning the problem sets reward.
2.12 Core Application: US Unemployment Rate (continued from Module 1)
This application completes the Module 2 sequence: Concept \(\rightarrow\) Math \(\rightarrow\) Simulate \(\rightarrow\) Real Data.
In Module 1 we pulled the US civilian unemployment rate from FRED and plotted it. We noted that the ACF decays slowly and that the series has high persistence. Now let us run the full stationarity workflow.
To pull the series live, store your FRED key in the environment (never hardcode it in course files) and run:
library(fredr); library(forecast); library(urca); library(patchwork)
fredr_set_key(Sys.getenv("FRED_KEY"))
# Pull the full series (1948-present)
unrate <- fredr(series_id = "UNRATE",
observation_start = as.Date("1948-01-01"))
y <- ts(unrate$value, start = c(1948, 1), frequency = 12)
T <- length(y)So that these notes render identically for everyone — with or without an API key or an internet connection — the code below loads the same series from the course’s cached copy in data/UNRATE.csv. The workflow from here on is identical either way.
library(forecast); library(urca); library(patchwork)
unrate <- read.csv("../data/UNRATE.csv") # cached FRED pull; see data/README.md
# The cache has one missing month (2025-10, a federal data-release gap).
# Keep the complete run from 1948-01 through the last month before the gap
# so the monthly ts index stays aligned.
first_na <- which(is.na(unrate$UNRATE))
if (length(first_na) > 0) unrate <- unrate[seq_len(min(first_na) - 1), ]
y <- ts(unrate$UNRATE, start = c(1948, 1), frequency = 12)
T <- length(y)
T[1] 933
Step 1: Visual Inspection
p1 <- autoplot(y) + theme_bw(base_size = 12) +
labs(title = "US Unemployment Rate, 1948-present", y = "Percent")
p2 <- ggAcf(y, lag.max = 60) + theme_bw(base_size = 12) +
labs(title = "ACF of unemployment rate")
p1 / p2What you see:
- Levels: The series oscillates roughly between 3% and 11% (with a spike near 15% at the 2020 pandemic shock), with no obvious deterministic trend. There are clear cyclical movements driven by recessions. No obvious variance changes (perhaps slightly more volatile in the post-1970 period, but not dramatically).
- ACF: Slow decay. The autocorrelation at lag 12 is still around 0.7. At lag 36 it is still positive and likely significant. This is the visual fingerprint of high persistence — possibly a unit root, or possibly \(\phi\) very close to 1.
Verdict from visual inspection: ambiguous. Looks persistent. Could go either way.
Step 2: ADF
The series has a non-zero mean but no obvious trend, so we use the drift specification.
adf_result <- ur.df(y, type = "drift", lags = 12, selectlags = "BIC")
summary(adf_result)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression drift
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-1.8826 -0.1323 -0.0195 0.1015 10.3134
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.179413 0.047927 3.743 0.000193 ***
z.lag.1 -0.031484 0.008067 -3.903 0.000102 ***
z.diff.lag 0.050759 0.032968 1.540 0.123995
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4137 on 917 degrees of freedom
Multiple R-squared: 0.01755, Adjusted R-squared: 0.01541
F-statistic: 8.191 on 2 and 917 DF, p-value: 0.0002979
Value of test-statistic is: -3.9027 7.6157
Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.43 -2.86 -2.57
phi1 6.43 4.59 3.78
mackinnon_cv(T = T, level = 0.05)[1] -2.864643
The statistic printed above is around \(-3.9\), below the 5% critical value of approximately \(-2.86\): taken alone, the ADF test rejects the unit root for this sample. Hold that rejection loosely. It is nowhere near the comfortable \(-5\) we got for a cleanly stationary simulated AR(1); it is sensitive to the sample window (the enormous, fast-reverting 2020 pandemic spike does real work pulling the test toward rejection); and — as the next step shows — the other formal test is about to disagree with it. One test statistic is one piece of evidence, not a verdict.
Step 3: KPSS
kpss_result <- ur.kpss(y, type = "mu")
summary(kpss_result)
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 6 lags.
Value of test-statistic is: 0.9703
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739
The KPSS statistic for the unemployment rate is typically large enough to reject stationarity at standard levels.
Step 4: Rule of Thumb
sd(diff(y)) / sd(y)[1] 0.2431225
The ratio is about \(0.24\) in this sample (the 2020 spike inflates the variance of the differences; pre-2020 samples give values closer to \(0.13\)) — either way, well below the 0.5 threshold. The rule of thumb says: difference.
Step 5: Synthesis
Lining up the evidence:
- Visual / ACF: Suggests high persistence. Inconclusive between \(\phi = 0.95\) and \(\phi = 1.0\).
- ADF: Rejects the unit root — but the rejection is sample-sensitive and leans on the extraordinary 2020 episode.
- KPSS: Rejects stationarity.
- Rule of thumb: Well below 0.5 — difference.
Notice where this lands us in the confirmatory table from Section 2.6: ADF and KPSS both reject — the contradiction row, which flags a possible structural break or misspecification. That is a plausible reading here: the 2020 pandemic spike is exactly the kind of extraordinary one-off event that unit-root and stationarity tests were not built for. So the formal evidence genuinely conflicts, the visual evidence says “very persistent,” and the rule of thumb says “difference.” When the evidence conflicts, we fall back on the cost asymmetry. The verdict for our purposes: treat the unemployment rate as \(I(1)\). We will work with \(\Delta y_t\) when we get to modeling.
Deeper Dive — A Note on the Macro Debate
This is real data being honest about its borderline nature. There is genuine and unresolved disagreement among macroeconomists about whether the unemployment rate is stationary. The arguments:
Stationary camp: The unemployment rate is bounded — it cannot go below 0% or much above 25% without the economy ceasing to function. A bounded series cannot have a true unit root in the strict sense, because a random walk has unbounded variance. So the unemployment rate must, at long horizons, be stationary. The persistence we see in 75 years of data is just \(\phi\) very close to 1, not \(\phi = 1\) exactly.
Unit root camp: The persistence is so strong that for any practical sample size, the unemployment rate is operationally indistinguishable from a unit root process. Whether the “true” \(\phi\) is 0.998 or 1.000 makes no difference for inference at horizons we care about. For modeling purposes, treat it as \(I(1)\).
Both arguments have merit. The data themselves cannot resolve the question — even with \(T \approx 930\) observations, the test power is not enough to distinguish \(\phi = 0.998\) from \(\phi = 1.000\). This is exactly the borderline case the module has been preparing you for. The cost asymmetry resolves the dilemma: when in doubt, treat as \(I(1)\), because spurious regression is more costly than over-differencing.
This sets up Problem Set 1, where you will run the same workflow on a different FRED series and write up your verdict and reasoning.
Common Pitfalls and Misconceptions
“ADF failing to reject means the series has a unit root.” No. Failing to reject means the data are insufficient to rule out the unit root. This could be because the series truly has a unit root, or because the test has low power against a near-unit-root stationary alternative. Don’t confuse “I cannot reject” with “I have proven.”
“The DF \(t\)-statistic follows a \(t\)-distribution.” No. It follows the Dickey-Fuller distribution, which is not symmetric and has a heavier left tail. Use MacKinnon critical values, not standard \(t\)-tables. The reason is the same as for spurious regression: standard OLS asymptotics break down with non-stationary regressors.
“Always use the trend specification to be safe.” No. Including unnecessary deterministic terms reduces test power. Match the specification to what you see in the plot. Use the trend specification only when there is an obvious trend.
“ADF and KPSS test the same thing.” No. They have opposite nulls. Use them as complements: confirmatory evidence when they agree, honest “don’t know” when both fail to reject.
“Over-differencing is harmless.” No. Differencing a trend-stationary series introduces a non-invertible MA(1) with \(\theta = -1\). The differenced series is technically stationary but the implied model is malformed.
“The rule of thumb is just folklore.” It is a rule of thumb, but it has a real basis in the variance ratio of stationary AR(1) processes. It is not a substitute for formal testing, but it captures the same intuition more quickly.
“More tests are always better.” Up to a point. Running ADF, KPSS, PP, and ERS on the same series gives you four pieces of evidence, but if they all use the same underlying data and similar testing principles, they are not fully independent. The best complement to ADF is KPSS (because the nulls are reversed); adding PP or ERS gives you incremental information at most.
Connection to Enders
- The difference operator and difference equations: Enders Ch. 1, pp. 1-46
- Stochastic difference equations: Enders Ch. 2, pp. 47-52
- The Dickey-Fuller test, all three specifications: Enders Ch. 4, pp. 206-215
- The Pantula principle and choosing among specifications: Enders Ch. 4, pp. 215-218
- Phillips-Perron and other tests: Enders Ch. 4, pp. 222-227
- KPSS and the reverse-null approach: Enders Ch. 4, pp. 228-230
- Order of integration and trends: Enders Ch. 4, pp. 181-189
- Trend-stationary vs difference-stationary: Enders Ch. 4, pp. 189-200
A convention warning when you cross-reference. Enders writes the AR(1) as \(y_t = a_0 + a_1 y_{t-1} + \epsilon_t\) (his \(a_1\) is our \(\phi\)), and he states stationarity conditions in terms of characteristic roots lying inside the unit circle. This course states the equivalent condition as roots of the lag polynomial \(\Phi(L) = 0\) lying outside the unit circle. Both are correct — they describe reciprocal roots — but keep the framing straight when reading Ch. 4. Hamilton (Ch. 17) is the rigorous reference for the asymptotics behind the DF distribution; Hyndman & Athanasopoulos discuss unit-root testing from a forecasting-workflow angle (and use the KPSS test as their default, the reverse of our ADF-first habit).
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, 427-431. (In
literature/Time Series/1979_dickey.pdf) - Dickey, D.A. and Fuller, W.A. (1981), “Likelihood Ratio Statistics for Autoregressive Time Series with a Unit Root,” Econometrica, 49, 1057-1072. (Joint \(F\)-tests \(\Phi_1\), \(\Phi_2\), \(\Phi_3\).)
- Pantula, S.G. (1989), “Testing for Unit Roots in Time Series Data,” Econometric Theory, 5, 256-271.
- MacKinnon, J.G. (1996), “Numerical Distribution Functions for Unit Root and Cointegration Tests,” Journal of Applied Econometrics, 11, 601-618.
- 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, 159-178.
- Phillips, P.C.B. (1987), “Time Series Regression with a Unit Root,” Econometrica, 55, 277-301.
- Schwert, G.W. (1989), “Tests for Unit Roots: A Monte Carlo Investigation,” Journal of Business and Economic Statistics, 7, 147-159.
- Elliott, G., Rothenberg, T.J., and Stock, J.H. (1996), “Efficient Tests for an Autoregressive Unit Root,” Econometrica, 64, 813-836. (In
literature/Time Series/1996_elliott.pdf) - Granger, C.W.J. and Newbold, P. (1974), “Spurious Regressions in Econometrics,” Journal of Econometrics, 2, 111-120. (In
literature/Time Series/1974_granger.pdf)
Practice Problems
Core Practice
Variance ratio computation. Compute the theoretical ratio \(\sigma_{\Delta y} / \sigma_y\) analytically for an AR(1) with \(\phi = 0.3, 0.7, 0.95\). Then verify by simulating long series (\(T = 5000\)) and computing the empirical ratio. Do they match?
Power of the ADF test. Conduct a Monte Carlo study to compute the power of the ADF test against a stationary AR(1) with \(\phi = 0.95\), for sample sizes \(T \in \{50, 100, 200, 500, 1000\}\). Plot power as a function of \(T\). At what sample size does the test achieve at least 50% power? At least 80%? Repeat the exercise for \(\phi = 0.90\) and \(\phi = 0.80\) and overlay the curves on a single plot. (This style of power study reappears on the problem sets — it is worth doing carefully now.)
The over-differencing trap in practice. Generate a trend-stationary series with \(T = 300\): \[y_t = 5 + 0.05\, t + u_t, \qquad u_t = 0.5\, u_{t-1} + \epsilon_t\] where \(u_t\) is the stationary AR(1) deviation from the trend (the same construction as Section 2.3) and \(\epsilon_t\) is white noise. Then:
- Detrend correctly by regressing \(y_t\) on a constant and \(t\). Plot the residuals and compute their ACF.
- Difference \(y_t\) to get \(\Delta y_t\). Plot it and compute its ACF.
- Fit an MA(1) to the differenced series using
arima(). What is the estimated \(\hat{\theta}\)? Is it close to \(-1\)?
- Fit an MA(1) to the differenced series using
- Discuss what each treatment tells you.
ADF specification matters. Generate two series, both with \(T = 200\):
- Series A (trend-stationary):
0.05 * (1:200) + arima.sim(n = 200, list(ar = 0.5)) - Series B (random walk with drift):
cumsum(rnorm(200) + 0.05)
Plot both. They will look similar. Now run the ADF test on each, with both the drift specification and the trend specification. What do you conclude in each case? Which specification gives you the right answer for which series?
- Series A (trend-stationary):
Real data workflow. Pull a different FRED series (suggestions: real GDP
GDPC1, industrial productionINDPRO, the federal funds rateFEDFUNDS, the CPICPIAUCSL). Run the full stationarity workflow: visual inspection, ACF, ADF (with the appropriate specification), KPSS, and the rule of thumb. Synthesize the evidence and write up a verdict in 2-3 paragraphs. Include the cost asymmetry in your reasoning.
Key Takeaways
The difference operator \(\Delta = (1 - L)\) is the discrete-time analog of a derivative. It is both the dependent variable in the DF test and the tool we use to make non-stationary series stationary.
Visual inspection of the series and the ACF are the cheapest evidence and are sufficient for the easy cases. They cannot distinguish trend-stationary from difference-stationary, and they fail in small samples.
The Dickey-Fuller test rearranges the AR(1) into a regression of \(\Delta y_t\) on \(y_{t-1}\). The test of \(H_0: \phi = 1\) becomes a test of \(H_0: \gamma = 0\) where \(\gamma = \phi - 1\). The \(t\)-statistic does not follow a standard \(t\)-distribution under the null because non-stationary regressors break OLS asymptotics — you must use MacKinnon critical values.
There are three ADF specifications (none, drift, trend). Choose based on what you see in the plot. Critical values differ across specifications.
KPSS is the reverse-null test: \(H_0\) = stationary, \(H_1\) = unit root. Use ADF and KPSS together for confirmatory analysis.
All unit root tests have low power against highly persistent stationary alternatives. With \(\phi = 0.95\) and \(T = 100\), ADF rejects less than 10% of the time. This is why we use multiple sources of evidence.
The rule of thumb \(\sigma_{\Delta y} / \sigma_y < 0.5\) flags series with persistence above roughly \(\phi > 0.875\) as candidates for differencing.
Order of integration \(I(d)\): the number of differences needed to achieve stationarity. Most economic series are \(I(0)\) or \(I(1)\); \(I(2)\) is rare; higher orders are essentially never seen.
The over-differencing trap: differencing a trend-stationary series introduces a non-invertible MA(1) with \(\theta = -1\). Use the trend-augmented ADF to distinguish trend-stationary from difference-stationary; detrend the former and difference the latter.
The cost asymmetry: when in doubt, difference. Spurious regression is a more dangerous inference error than over-differencing. But be aware of the over-differencing tax — it shows up in your ARMA model fits later.
Looking Ahead — Module 3
Two loose ends from this module become Module 3’s opening material:
- The over-differencing trap named a boundary. We showed that differencing a trend-stationary series manufactures an MA(1) with a unit root in \(\Theta(L)\) — a violation of invertibility. Module 3 defines invertibility properly as the MA counterpart of AR stationarity and probes that boundary experimentally with a nearly non-invertible MA(1).
- The power problem needs a complementary tool. Formal tests cannot reliably separate \(\phi = 0.95\) from \(\phi = 1.0\) in realistic samples. Module 3 adds the correlogram fingerprint — reading ACF/PACF pairs to identify AR(\(p\)) and MA(\(q\)) structure — which works alongside testing rather than replacing it. The mixing board gets its first new sliders: higher-order AR memory (\(\phi_1, \ldots, \phi_p\)) and memory of past shocks (\(\theta_1, \ldots, \theta_q\)).