Workbook 9 — Assessing a Stationary Approximation

STA 542 · Introduction to Time Series Analysis

\[\def\E{\mathbb{E}} \def\P{\mathbb{P}} \def\R{\mathbb{R}} \def\Z{\mathbb{Z}} \def\Var{\operatorname{Var}} \def\Cov{\operatorname{Cov}}\]

Stationary Fluctuations or Accumulated Shocks?

Workbook 8 considered how to transform a series so that a stationary model might be appropriate. We now examine statistical evidence for that choice.

We begin with observations \(X_1,\ldots,X_T\) from a process whose mean is known to be zero. Even in this setting, sustained movements can be difficult to interpret: both a stationary process and a random walk can produce them.

We compare a stationary Gaussian AR(1) with coefficient \(0.98\) and a random walk. Both use innovation variance one and start at the same value, drawn from the AR(1)’s stationary distribution independently of both sequences of subsequent innovations.

set.seed(542)
T_path <- 200
ar_path <- numeric(T_path + 1)
ar_path[1] <- rnorm(1, sd = 1 / sqrt(1 - 0.98^2))
for (j in 2:length(ar_path)) ar_path[j] <- 0.98 * ar_path[j - 1] + rnorm(1)
rw_path <- ar_path[1] + c(0, cumsum(rnorm(T_path)))
path_ylim <- range(ar_path, rw_path)
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), bg = "gray92")
plot(0:T_path, ar_path, type = "l", col = "#00539B", ylim = path_ylim,
     xlab = "Time index", ylab = "Level", main = "Stationary AR(1)")
plot(0:T_path, rw_path, type = "l", col = "#AD5B00", ylim = path_ylim,
     xlab = "Time index", ylab = "Level", main = "Random walk")

Figure 1: A stationary AR(1) with coefficient 0.98 and a zero-drift random walk, each with 200 transitions. Both start at the same value, drawn from the AR(1)’s stationary distribution. The panels share a vertical scale. Simulated data, seed 542.

Earlier, we have established that a stationary AR(1) with coefficient \(0.98\) is weakly stationary. In particular, its iid Gaussian innovations here make it strictly stationary. The random walk is not stationary. However, merely looking at their paths does not make it clear which one is stationary, versus which one is not.

lag_max <- 40
acvf_ar <- as.numeric(acf(ar_path, type = "covariance", lag.max = lag_max,
                          plot = FALSE)$acf)
acvf_rw <- as.numeric(acf(rw_path, type = "covariance", lag.max = lag_max,
                          plot = FALSE)$acf)
acvf_ylim <- range(acvf_ar, acvf_rw)
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), bg = "gray92")
acf(ar_path, type = "covariance", lag.max = lag_max, ci = 0, ylim = acvf_ylim,
    col = "#00539B", main = "Stationary AR(1)",
    xlab = "Lag", ylab = "Sample ACVF")
acf(rw_path, type = "covariance", lag.max = lag_max, ci = 0, ylim = acvf_ylim,
    col = "#AD5B00", main = "Random walk",
    xlab = "Lag", ylab = "Sample ACVF")

Figure 2: Sample ACVF of the two records above, using divisor T. The panels share a vertical scale. Simulated data, seed 542.

Their sample ACVFs also decay slowly. A single lag plot does not settle the distinction either.

How Large Should Partial Sums Be Under Stationarity?

To investigate these movements, consider the running sums of the observations: \[ S_t:=\sum_{j=1}^{t}X_j. \] A stretch of positive observations makes this sum rise; a stretch of negative observations makes it fall. How large should these sums be under a stationary model?

For this calculation, we use a new pair of 200-observation records: a stationary Gaussian AR(1) with coefficient \(0.5\) and a random walk, both with innovation variance one. The weaker AR dependence makes the contrast in their partial sums easier to see.

set.seed(54296)
T_kpss <- 200
kpss_records <- list("AR(1), coefficient 0.5" = as.numeric(arima.sim(
  model = list(ar = 0.5), n = T_kpss)),
  "Random walk" = cumsum(rnorm(T_kpss)))
kpss_sums <- lapply(kpss_records, cumsum)
par(mfrow = c(2, 2), mar = c(3.5, 4, 3, 1), mgp = c(2.2, 0.7, 0), bg = "gray92")
for (j in seq_along(kpss_records)) {
  plot(kpss_records[[j]], type = "l", col = "#00539B",
       xlab = "Observation", ylab = "Level", main = names(kpss_records)[j])
}
for (j in seq_along(kpss_records)) {
  plot(kpss_sums[[j]], type = "l", ylim = range(unlist(kpss_sums)),
       col = "#AD5B00", xlab = "Observation",
       ylab = "Partial sum", main = names(kpss_records)[j])
  abline(h = 0, col = "gray60")
}

Figure 3: A stationary AR(1) with coefficient 0.5 and a random walk, with their partial sums. Both generating processes have mean zero; no sample mean is subtracted. The lower panels share a vertical scale. Simulated data, seed 54296.

The random walk has much larger partial sums in these two records. We measure their total squared size by \(\sum_{t=1}^{T}S_t^2\). To judge whether this quantity is unusually large, we first calculate its expected size under a stationary model.

Proposition 9.1 (The scale of centered partial sums). Suppose \((X_t)\) is mean-zero and weakly stationary, with \(\sum_{h\in\mathbb Z}|\gamma(h)|<\infty\) and positive long-run variance \(v:=\sum_{h\in\mathbb Z}\gamma(h)>0\). For \(S_t=\sum_{j=1}^{t}X_j\), \[ \E[S_t^2] =t\gamma(0)+2\sum_{h=1}^{t-1}(t-h)\gamma(h), \qquad \frac{\E[S_t^2]}{t}\longrightarrow v. \] Moreover, \[ \E\left[\frac{1}{T^2v}\sum_{t=1}^{T}S_t^2\right] \longrightarrow\frac12. \]

Proof. Expanding the variance of the sum gives the covariance identity. Since \(S_t=t\bar X_t\), Workbook 6, Proposition 6.2, gives \(\E[S_t^2]/t=t\Var(\bar X_t)\to v\). Summing the covariance identity over \(t\) gives \[ \frac{1}{T^2}\sum_{t=1}^{T}\E[S_t^2] =\frac{T+1}{2T}\gamma(0) +\sum_{h=1}^{T-1} \frac{(T-h)(T-h+1)}{T^2}\gamma(h). \] For each fixed \(h\), the coefficient in the last sum tends to one and is bounded by one. Absolute summability therefore gives the limit \(\gamma(0)/2+\sum_{h\ge1}\gamma(h)=v/2\). Dividing by \(v\) proves the final assertion. \(\square\)

Under stationarity, the variance of \(S_t\) is approximately \(tv\) for large \(t\). The long-run variance therefore sets the scale against which the partial sums are assessed.

What Changes When the Observations Accumulate?

Consider a random walk starting at zero, \[ X_0=0,\qquad X_t=\sum_{j=1}^{t}Z_j, \qquad (Z_j)\sim\mathrm{WN}(0,\sigma^2),\quad \sigma^2>0. \] Here each observation \(X_t\) is already a sum of shocks. Adding the observations to form \(S_t\) counts each early shock repeatedly: \[ S_t=\sum_{j=1}^{t}X_j =\sum_{j=1}^{t}(t-j+1)Z_j. \] The white-noise covariances give \[ \Var(S_t)=\sigma^2\sum_{j=1}^{t}j^2 =\frac{\sigma^2t(t+1)(2t+1)}{6} \sim\frac{\sigma^2t^3}{3}. \] Thus the variance of these sums grows as \(t^3\), compared with \(tv\) under the stationary model. Summing over the record gives \[ \E\left[\sum_{t=1}^{T}S_t^2\right] =\frac{\sigma^2T(T+1)^2(T+2)}{12} \sim\frac{\sigma^2T^4}{12}. \] This suggests comparing the squared partial sums with their stationary scale \(T^2v\).

Constructing the KPSS Statistic

The Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test uses this comparison. Its null is the centered stationary model in Proposition 9.1; the conditions needed to calibrate the test are stated below. We estimate \(v\) with the Bartlett estimator from Workbook 6.

Definition 9.2 (KPSS statistic with known mean zero). Given \(X_1,\ldots,X_T\) from a process with known mean zero, put \[ S_t:=\sum_{j=1}^{t}X_j, \qquad \hat\gamma_0(h):=\frac1T\sum_{t=1}^{T-h}X_tX_{t+h}. \] For a specified nonnegative integer bandwidth \(H<T\), estimate the long-run variance by \[ \widehat v_H:=\hat\gamma_0(0)+ 2\sum_{h=1}^{H}\left(1-\frac{h}{H+1}\right)\hat\gamma_0(h). \] When \(\widehat v_H>0\), the KPSS statistic is the random variable \[ K_T:=\frac{\sum_{t=1}^{T}S_t^2}{T^2\widehat v_H}. \] Large values of \(K_T\) count against this stationary model. To decide how large is unusual, we need its null distribution.

Remark (A sample ACVF does not require a population ACVF). If \((X_t)\) is not weakly stationary, \(\Cov(X_t,X_{t+h})\) may depend on \(t\) as well as on \(h\). There is then no population ACVF \(\gamma(h)\). The products \(\hat\gamma_0(h)\) remain well defined: they are functions of the observed record. The whole process of computing \(K_T\), does not require that a ACVF exists.

Because the mean is known to be zero, \(\hat\gamma_0(h)\) averages the products \(X_tX_{t+h}\) without subtracting a sample mean. The Bartlett weights are the same as in the workbook on long-run variance. They account for the contribution of serial dependence to the size of the partial sums.

Why Does Accumulation Make the Ratio Large?

For the stationary model, the expected numerator is asymptotic to \(T^2v/2\), and a consistent Bartlett estimate makes the denominator approximately \(T^2v\). For the random walk above, the expected numerator instead grows as \(\sigma^2T^4/12\).

The Bartlett estimate can also become large: it is being computed from nonstationary levels, which have no stationary long-run variance. Could this growth cancel the growth in the numerator? For every record, Cauchy–Schwarz gives \(|\hat\gamma_0(h)|\le\hat\gamma_0(0)\), so \[ 0\le\widehat v_H\le(H+1)\hat\gamma_0(0). \] For this random walk, \(\E[\hat\gamma_0(0)]=\sigma^2(T+1)/2\). Hence \[ \E[T^2\widehat v_H] \le\frac{\sigma^2}{2}(H+1)T^2(T+1). \] The numerator’s expectation grows as \(T^4\), whereas the denominator’s grows at most as \((H+1)T^3\). These rates suggest growth of the statistic on the scale \(T/(H+1)\) when the bandwidth grows more slowly than the sample size.

Comparing expectations alone does not prove that a random ratio grows. If the white-noise innovations are also independent and identically distributed, a functional CLT supplies the stronger conclusion: for this random walk, if \(H_T/T\to0\), then for every fixed \(c>0\), \[ \P(K_T>c)\longrightarrow1. \] We take this result without proof. It also holds for \(X_t=\sum_{j=1}^{t}V_j\) when the dependent increments \((V_t)\) satisfy the linear-process conditions in the box below. Thus accumulation can make KPSS reject with probability tending to one even though its estimated denominator grows.

Choosing a Critical Value

An expectation does not determine a test’s critical value. We need the null distribution of \(K_T\). Under the conditions in the box below, the known-zero-mean KPSS test rejects above the approximate asymptotic 5% critical value 1.66.¹

NoteMore Details ♠: The KPSS Null Distribution

Suppose \((X_t)\) satisfies the stationary assumptions of Proposition 9.1. Put \(S_0=0\) and write \(\widehat v_T:=\widehat v_{H_T}\) for the Bartlett estimate at the chosen bandwidth. The usual asymptotic calibration follows from two additional conditions:

1. A functional CLT for the partial sums: \[ \frac{S_{\lfloor Tr\rfloor}}{\sqrt{Tv}} \Rightarrow W(r),\qquad 0\le r\le1, \] as a process, where \(W\) is standard Brownian motion: a mean-zero Gaussian process with continuous paths and \(\Cov(W(r),W(s))=\min\{r,s\}\). The convergence concerns the whole partial-sum path.

2. A consistent long-run variance estimate: \[ \widehat v_T\xrightarrow{P}v. \]

Rewrite the statistic as \[ K_T =\frac{v}{\widehat v_T}\frac1T \sum_{t=1}^{T}\left(\frac{S_t}{\sqrt{Tv}}\right)^2 \xrightarrow{d}\int_0^1W(r)^2\,dr. \] The functional CLT gives the limit of the Riemann sum of squared path values. Consistency makes the factor \(v/\widehat v_T\) converge to one. The limit distribution therefore does not depend on the unknown \(v\).

Let \(q_{1-\alpha}\) be the \((1-\alpha)\) quantile of \(\int_0^1W(r)^2\,dr\). Then \[ \P_{H_0}(K_T>q_{1-\alpha})\longrightarrow\alpha. \] For \(\alpha=0.05\), this quantile is approximately 1.66.

One sufficient setting for both conditions is \[ X_t=\sum_{j=0}^{\infty}\psi_jZ_{t-j}, \qquad \sum_{j\ge0}(1+j)|\psi_j|<\infty, \qquad \sum_{j\ge0}\psi_j\ne0, \] where \((Z_t)\) is i.i.d. with mean zero, positive variance, and finite fourth moment, and \(H_T\to\infty\) with \(H_T=o(\sqrt T)\). Here \(v=\Var(Z_0)(\sum_{j\ge0}\psi_j)^2>0\). We take the functional CLT and Bartlett consistency for this class without proof.

For the two simulated records above, using Bartlett bandwidth \(H=8\) gives the following results. The table shows \(T^{-2}\sum_{t=1}^T S_t^2\) and \(\widehat v_H\) separately before taking their ratio.

Run the following code once to define the test functions used below.

# Fixed-specification test helpers for STA 542 Workbooks 9 and 10.
# Base R only. Inputs must be finite observations at equally spaced times;
# missing observations are rejected, never silently removed.
# These are approximate tests under their stated model assumptions, not
# universal tests of weak stationarity. Choose deterministic terms, ADF lags,
# and KPSS bandwidth before inspecting the resulting rejection decisions.
# ADF assumes the chosen lag augmentation adequately describes short-run
# dynamics under a unit-root null. EG additionally presumes two I(1) series
# and a specified deterministic part in the cointegrating regression.
# KPSS uses its conventional short-memory stationary null, positive long-run
# variance, and large-sample calibration. A fixed classroom bandwidth is a
# sensitivity choice, not a proof that a consistent LRV estimate was obtained.
# No automatic lag selection and no interpolated or clipped p-values.
#
# MacKinnon (2010), "Critical Values for Cointegration Tests", QED WP 1227:
# https://www.econ.queensu.ca/research/working-papers/1227
# Coefficient arrays and response-surface convention below are transcribed
# from statsmodels 0.14.5, statsmodels/tsa/adfvalues.py (N = 1 and N = 2 only):
# https://github.com/statsmodels/statsmodels/blob/v0.14.5/statsmodels/tsa/adfvalues.py
# Exception: statsmodels attributes the no-constant ADF array to MacKinnon
# (1996); that case was not updated in 2010. All constant/trend arrays used
# here are the 2010 response surfaces. Their values are approximations.
# Its BSD 3-clause notices apply to the transcribed arrays:
# Copyright (C) 2006, Jonathan E. Taylor. All rights reserved.
# Copyright (c) 2006-2008 Scipy Developers. All rights reserved.
# Copyright (c) 2009-2018 statsmodels Developers. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# a. Redistributions of source code must retain the above copyright notice,
#    this list of conditions and the following disclaimer.
# b. Redistributions in binary form must reproduce the above copyright
#    notice, this list of conditions and the following disclaimer in the
#    documentation and/or other materials provided with the distribution.
# c. Neither the name of statsmodels nor the names of its contributors
#    may be used to endorse or promote products derived from this software
#    without specific prior written permission.
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL STATSMODELS OR CONTRIBUTORS BE LIABLE FOR
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
# DAMAGE.

.wb_series <- function(x, name = "x") {
  if (!is.numeric(x) || is.complex(x) || !is.null(dim(x)) ||
      length(x) < 4L || any(!is.finite(x))) {
    stop(name, " must be a finite numeric vector with at least four observations.")
  }
  as.numeric(x)
}

.wb_integer <- function(x, name, upper) {
  if (length(x) != 1L || !is.numeric(x) || !is.finite(x) ||
      x < 0 || x != floor(x) || x > upper) {
    stop(name, " must be a nonnegative integer no larger than ", upper, ".")
  }
  as.integer(x)
}

.wb_ols <- function(y, design) {
  fit <- lm.fit(design, y)
  k <- ncol(design)
  df <- length(y) - k
  if (fit$rank < k || df <= 0L) stop("Regression is rank deficient or has no residual degrees of freedom.")
  sigma2 <- sum(fit$residuals^2) / df
  if (!is.finite(sigma2) || sigma2 <= 0) stop("Regression has zero or invalid residual variance.")
  se_pivot <- sqrt(diag(chol2inv(qr.R(fit$qr))) * sigma2)
  se <- numeric(k)
  se[fit$qr$pivot] <- se_pivot
  names(se) <- colnames(design)
  list(coefficients = fit$coefficients, se = se,
       residuals = fit$residuals, fitted = fit$fitted.values,
       df = df, sigma2 = sigma2)
}

.wb_mackinnon <- function(nobs, deterministic, N = 1L) {
  # Rows: 1%, 5%, 10%. Columns: a0, a1, a2, a3.
  # Critical value = a0 + a1/nobs + a2/nobs^2 + a3/nobs^3.
  key <- paste(N, deterministic, sep = "_")
  coefficients <- switch(key,
    "1_none" = c(-2.56574, -2.2358, -3.627, 0,
                 -1.94100, -0.2686, -3.365, 31.223,
                 -1.61682, 0.2656, -2.714, 25.364),
    "1_constant" = c(-3.43035, -6.5393, -16.786, -79.433,
                     -2.86154, -2.8903, -4.234, -40.040,
                     -2.56677, -1.5384, -2.809, 0),
    "1_trend" = c(-3.95877, -9.0531, -28.428, -134.155,
                  -3.41049, -4.3904, -9.036, -45.374,
                  -3.12705, -2.5856, -3.925, -22.380),
    "2_constant" = c(-3.89644, -10.9519, -33.527, 0,
                     -3.33613, -6.1101, -6.823, 0,
                     -3.04445, -4.2412, -2.720, 0),
    "2_trend" = c(-4.32762, -15.4387, -35.679, 0,
                  -3.78057, -9.5106, -12.074, 0,
                  -3.49631, -7.0815, -7.538, 21.892),
    stop("Unsupported MacKinnon specification."))
  a <- matrix(coefficients, nrow = 3, byrow = TRUE)
  critical <- as.numeric(a %*% (1 / nobs)^(0:3))
  setNames(critical, c("1%", "5%", "10%"))
}

wb_adf <- function(x, lags = 0, deterministic = c("none", "constant", "trend")) {
  deterministic <- match.arg(deterministic)
  x <- .wb_series(x)
  n <- length(x)
  ndet <- switch(deterministic, none = 0L, constant = 1L, trend = 2L)
  lags <- .wb_integer(lags, "lags", floor(n / 2) - ndet - 1L)
  if (max(x) == min(x)) stop("ADF is undefined for a constant series.")
  dx <- diff(x)
  t <- seq.int(lags + 2L, n)
  design <- matrix(x[t - 1L], ncol = 1L, dimnames = list(NULL, "level"))
  if (lags > 0L) {
    lagged <- vapply(seq_len(lags), function(j) dx[t - j - 1L], numeric(length(t)))
    colnames(lagged) <- paste0("difference_lag", seq_len(lags))
    design <- cbind(design, lagged)
  }
  if (ndet >= 1L) design <- cbind(design, constant = 1)
  if (ndet == 2L) design <- cbind(design, trend = seq_along(t))
  fit <- .wb_ols(dx[t - 1L], design)
  if (sum(fit$residuals^2) <= 100 * .Machine$double.eps^2 * sum(dx[t - 1L]^2)) {
    stop("The ADF regression is essentially exact; its statistic is unreliable.")
  }
  statistic <- unname(fit$coefficients["level"] / fit$se["level"])
  critical <- .wb_mackinnon(length(t), deterministic, N = 1L)
  list(statistic = statistic, stat = statistic, nobs = length(t), n = n,
       critical_nobs = length(t), lags = lags, deterministic = deterministic,
       critical = critical, reject = statistic < critical,
       null = "Unit root", alternative = "Stationary relative to the specified deterministic terms",
       coefficients = fit$coefficients, se = fit$se,
       residuals = fit$residuals, fitted = fit$fitted, df = fit$df)
}

wb_eg <- function(y, x, lags = 0, deterministic = c("constant", "trend")) {
  deterministic <- match.arg(deterministic)
  y <- .wb_series(y, "y")
  x <- .wb_series(x, "x")
  if (length(y) != length(x)) stop("x and y must have the same length and aligned dates.")
  n <- length(y)
  design <- cbind(x = x, constant = 1)
  if (deterministic == "trend") design <- cbind(design, trend = seq_len(n))
  levels <- .wb_ols(y, design)
  if (sum(levels$residuals^2) < 100 * .Machine$double.eps * sum((y - mean(y))^2)) {
    stop("The levels relation is nearly exact; the residual test is numerically unreliable.")
  }
  adf <- wb_adf(levels$residuals, lags = lags, deterministic = "none")
  # EG uses N = 2 and the levels regression's deterministic specification.
  # Match statsmodels.coint: critical-value nobs is n - 1 regardless of lags;
  # adf$nobs = n - 1 - lags records the actual second-stage regression size.
  critical <- .wb_mackinnon(n - 1L, deterministic, N = 2L)
  list(statistic = adf$statistic, stat = adf$statistic,
       nobs = adf$nobs, n = n, critical_nobs = n - 1L,
       lags = adf$lags, deterministic = deterministic,
       critical = critical, reject = adf$statistic < critical,
       null = "No cointegration between the two I(1) series",
       alternative = "Cointegration with the specified deterministic terms",
       coefficients = levels$coefficients, residuals = levels$residuals,
       fitted = levels$fitted, adf = adf)
}

wb_kpss <- function(x, bandwidth, deterministic = c("constant", "trend", "none")) {
  deterministic <- match.arg(deterministic)
  x <- .wb_series(x)
  n <- length(x)
  if (deterministic != "none" && max(x) == min(x)) stop("KPSS is undefined for a constant series.")
  bandwidth <- .wb_integer(bandwidth, "bandwidth", n - 1L)
  if (deterministic == "none") {
    # Known population mean zero: do not estimate or subtract a sample mean.
    e <- x
    coefficients <- numeric(0)
  } else {
    design <- matrix(1, nrow = n, ncol = 1L, dimnames = list(NULL, "constant"))
    if (deterministic == "trend") design <- cbind(design, trend = seq_len(n))
    fit <- .wb_ols(x, design)
    e <- fit$residuals
    coefficients <- fit$coefficients
    if (sum(e^2) <= 100 * .Machine$double.eps^2 * sum(x^2)) {
      stop("The fitted deterministic part is essentially exact; KPSS is undefined.")
    }
  }
  partial_sums <- cumsum(e)
  gamma <- numeric(bandwidth + 1L)
  gamma[1L] <- sum(e^2) / n
  if (bandwidth > 0L) {
    for (h in seq_len(bandwidth)) gamma[h + 1L] <- sum(e[(h + 1L):n] * e[1L:(n - h)]) / n
  }
  # bandwidth is the maximum included lag L, with weight 1 - h/(L+1).
  lrv <- gamma[1L]
  if (bandwidth > 0L) lrv <- lrv + 2 * sum((1 - seq_len(bandwidth) / (bandwidth + 1)) * gamma[-1L])
  if (!is.finite(lrv) || lrv <= 0) stop("KPSS requires a positive estimated long-run variance.")
  statistic <- sum(partial_sums^2) / (n^2 * lrv)
  # Fitted cases: KPSS (1992) Table 1, also used by tseries::kpss.test.
  # Known mean zero: Kagalwala (2022), Table 2, Bartlett KPSS, rounded.
  # https://doi.org/10.1177/1536867X221106371
  # Its limit is integral_0^1 W(r)^2 dr, not the demeaned Brownian bridge law.
  critical <- switch(deterministic, none = c(2.79, 1.66, 1.20),
                     constant = c(0.739, 0.463, 0.347), trend = c(0.216, 0.146, 0.119))
  names(critical) <- c("1%", "5%", "10%")
  list(statistic = statistic, stat = statistic, nobs = n, n = n,
       bandwidth = bandwidth, deterministic = deterministic,
       critical = critical, reject = statistic > critical,
       null = switch(deterministic, none = "Stationarity with known mean zero",
                     constant = "Level stationarity", trend = "Trend stationarity"),
       alternative = "A stochastic trend under the KPSS model",
       lrv = lrv, residuals = e, partial_sums = partial_sums,
       autocovariances = gamma, coefficients = coefficients)
}

kpss_centered <- lapply(kpss_records, wb_kpss, bandwidth = 8,
                        deterministic = "none")
knitr::kable(data.frame(
  Series = names(kpss_records),
  Squared_sums_over_T2 = vapply(kpss_centered, function(z)
    sum(z$partial_sums^2) / T_kpss^2, numeric(1)),
  Bartlett_estimate = vapply(kpss_centered, function(z) z$lrv, numeric(1)),
  Statistic = vapply(kpss_centered, function(z) z$statistic, numeric(1)),
  Critical_5pct = 1.66
), digits = 3, row.names = FALSE,
col.names = c("Series", "Squared sums / T²", "Bartlett estimate", "KPSS", "5% cutoff"),
caption = "Known mean zero; Bartlett bandwidth eight. KPSS = (squared sums / T²) / Bartlett estimate.")

Known mean zero; Bartlett bandwidth eight. KPSS = (squared sums / T²) / Bartlett estimate.

Series	Squared sums / T²	Bartlett estimate	KPSS	5% cutoff
AR(1), coefficient 0.5	0.619	2.660	0.233	1.66
Random walk	107.658	82.764	1.301	1.66

For the random walk, the squared sums divided by \(T^2\) equal 107.658, but the Bartlett estimate is also large, at 82.764. Their ratio is only 1.301, below 1.66. The statistic is larger than the stationary record’s 0.233, but neither record leads to rejection. Growth with sample size does not guarantee rejection for every path of length 200.

Which Processes Does KPSS Cover?

The stationary model used above has absolutely summable autocovariances and positive long-run variance. We now give this class a name and define its accumulated counterpart.

Definition 9.3 (Integration of order zero and one). In this workbook, a process is integrated of order zero, written \(I(0)\), if it is weakly stationary, its ACVF \(\gamma\) satisfies \(\sum_{h\in\mathbb Z}|\gamma(h)|<\infty\), and its long-run variance is \[ v:=\sum_{h\in\mathbb Z}\gamma(h)>0. \] In the mean-zero setting, a process is integrated of order one, written \(I(1)\), if \[ X_t=X_0+\sum_{j=1}^{t}V_j,\qquad t\geq1, \] where \((V_t)\) is mean-zero and \(I(0)\), \(\E[X_0]=0\), and \(\E[X_0^2]<\infty\).

Remark (Higher integration orders). For an integer \(d\ge 1\), we call \((X_t)\) \(I(d)\) if \((\nabla X_t)_{t\ge 1}\) is \(I(d-1)\), taking \(I(0)\) from Definition 9.3. In particular, \(\nabla X_t=V_t\) is \(I(0)\) when \((X_t)\) is \(I(1)\), and two differences of an \(I(2)\) process are \(I(0)\). The order records how many differences reach this stationary class.

Proposition 9.4 (Variance of an accumulated level). Under the \(I(1)\) representation in Definition 9.3, with increment long-run variance \(v\), \[ \frac{\Var(X_t)}{t}\longrightarrow v. \]

Proof. Let \(A_t=\sum_{j=1}^{t}V_j\). The partial-sum variance result from the workbook on long-run variance gives \(\Var(A_t)/t\to v\). The initial value contributes \(\Var(X_0)/t\to0\), and Cauchy–Schwarz gives \[ \frac{|\Cov(X_0,A_t)|}{t} \leq \sqrt{\Var(X_0)} \sqrt{\frac{\Var(A_t)}{t}}\,t^{-1/2} \longrightarrow0. \] Expanding \(\Var(X_0+A_t)\) proves the result. \(\square\)

What Changes When the Mean Is Unknown?

Suppose instead that we observe \(Y_t=\mu+X_t\), where \(\mu\) is unknown and \((X_t)\) is mean-zero and \(I(0)\) under the null. The usual constant-mean KPSS test estimates \(\mu\) by \(\bar Y_T\) and applies the same statistic to the residuals \(e_t=Y_t-\bar Y_T\).

Write \(S_t=\sum_{j=1}^{t}X_j\) for the partial sums of the unobserved centered process. The observable residual sums are \[ R_t:=\sum_{j=1}^{t}e_j=S_t-\frac{t}{T}S_T. \] Unlike \(S_T\), the final residual sum \(R_T\) is always zero. Estimating the mean changes the path used in the numerator and its null reference distribution. We need a different critical value. The Bartlett denominator must also be computed from the residuals \(e_t\).

For an unknown linear mean, consider \(Y_t=a+bt+X_t\). If \((X_t)\) is mean-zero and \(I(0)\), we call \((Y_t)\) trend-stationary. Subtracting the population mean \(a+bt\) then recovers stationary fluctuations. If instead \((X_t)\) is \(I(1)\), subtracting the line leaves a process whose variance grows as in Proposition 9.4.

Estimate \(a,b\) by OLS and use \(e_t=Y_t-\hat a-\hat bt\) in the numerator and Bartlett denominator. The fitted residuals satisfy \[ \sum_{t=1}^T e_t=0,\qquad \sum_{t=1}^T t e_t=0. \] Fitting the slope imposes a second constraint, so the fitted-constant reference no longer applies. Under the sufficient linear-process and bandwidth conditions in the KPSS reference box, the critical values are:

Approximate asymptotic KPSS critical values. Known zero: Kagalwala (2022), rounded. Fitted constant and trend: Kwiatkowski et al. (1992).

Mean specification	10%	5%	1%
Known zero	1.200	1.660	2.790
Unknown constant	0.347	0.463	0.739
Unknown linear trend	0.119	0.146	0.216

In wb_kpss, one can use deterministic = "none" when the population mean is known' (suspected) to be zero. Versus“constant”or“trend”` when estimating a mean or line.

Does Annual Passenger Growth Need Another Difference?

The monthly international airline passenger totals in datasets::AirPassengers cover 1949–1960, in thousands of passengers. Write \(P_t\) for the It is a two-stage procedure: estimate a relationship between the series, then test whether its deviations still contain a unit root. Suppose (X_t) and (Y_t) are individually (I(1)). Cointegration means that, for some constants (a,b), [ U_t=Y_t-a-bX_t ]is (I(0)). The individual series can wander while their difference from this relationship remains stationary. First, estimate the relationship. Regress (Y_t) on an intercept and (X_t), giving fitted residuals [ U_t=Y_t-a-bX_t. ]Second, examine the persistence of those residuals. Fit [ U_t =c,U_{t-1} +{j=1}^{k}d_jU{t-j} +e_t. ]This is the same regression form used by ADF: - The lagged level (U_{t-1}) lets us assess whether deviations predict movement back toward zero. - The lagged changes accommodate short-run dependence. - A sufficiently negative value of (c/(c)) provides evidence against no cointegration, under the test’s assumptions. For intuition, if a spread followed (U_t=U_{t-1}+Z_t), its changes would satisfy [ U_t=()U_{t-1}+Z_t. ]A unit root gives coefficient zero. A stationary AR(1) gives a negative coefficient. The augmented regression extends this idea beyond that simple dependence model. The crucial difference is the critical value. We estimated the spread using the same observations. Even unrelated random walks can be fitted together surprisingly well, so the null calibration must account for that fitting step. Engle–Granger therefore uses cointegration critical values; simply running an ordinary ADF test on the fitted residuals and accepting its usual p-value gives the wrong calibration. Engle–Granger implementation documentationpassenger total in month \(t\). Annual log growth is \[ W_t=(1-B^{12})\log P_t=\log\frac{P_t}{P_{t-12}}. \] This gives 132 monthly observations, from January 1950 through December 1960. Workbook 8 considered taking another ordinary difference; the plots left its necessity unresolved.

Can a linear trend with stationary fluctuations describe annual growth? The proposed model is \[ W_t=a+bt+U_t, \] with \((U_t)\) mean-zero and \(I(0)\) under the trend-KPSS null. We estimate \(a,b\) by OLS and plot the fitted line alongside annual growth.

air_growth <- diff(log(datasets::AirPassengers), lag = 12)
air_t <- seq_along(air_growth)
air_line <- lm(as.numeric(air_growth) ~ air_t)
air_residual <- ts(residuals(air_line), start = start(air_growth),
                   frequency = 12)

par(mar = c(4, 4.5, 2, 1), bg = "white")
plot(air_growth, col = "#00539B", xlab = "Year",
     ylab = "Annual log growth", main = "Can a trend describe the changing mean?")
lines(ts(fitted(air_line), start = start(air_growth), frequency = 12),
      col = "#AD5B00", lwd = 2)
legend("topright", c("Annual log growth", "Fitted affine trend"),
       col = c("#00539B", "#AD5B00"), lty = 1, lwd = c(1, 2), bty = "n")

Figure 4: Annual log growth in passenger totals and an affine trend fitted by OLS. Real data: datasets::AirPassengers.

Subtracting the fitted line gives residuals \(\widehat U_t=W_t-\widehat a-\widehat bt\). These measure deviations from the fitted trend in annual growth.

par(mfrow = c(1, 2), mar = c(4, 4.5, 3, 1), bg = "white")
plot(air_residual, col = "#00539B", xlab = "Year",
     ylab = "Residual annual log growth", main = "After subtracting the fitted line")
abline(h = 0, col = "gray60", lty = 2)
acf(as.numeric(air_residual), type = "covariance", lag.max = 24,
    col = "#00539B", xlab = "Lag (months)", ylab = "Sample autocovariance",
    main = "Dependence remains")

Figure 5: Residuals after fitting an affine trend to annual log growth, with their sample ACVF at monthly lags. Correlation remains after the mean adjustment. Real data: datasets::AirPassengers.

The residuals remain above or below zero for stretches of months. Their sample ACVF is positive at short lags and negative around the annual lag. A stationary process can have these features. The Bartlett estimate must account for dependence when setting the scale of cumulative residuals.

For an initial calculation, use Bartlett bandwidth \(H=8\). This includes short-lag covariances but omits the negative covariances around twelve months. Apply trend KPSS to \(W_t\); it fits the line and computes both the partial sums and the Bartlett estimate from the residuals.

air_kpss <- wb_kpss(air_growth, bandwidth = 8, deterministic = "trend")
air_T <- length(air_growth)
air_scaled_sums <- air_kpss$partial_sums / sqrt(air_T * air_kpss$lrv)
knitr::kable(data.frame(
  `Squared sums / T²` = sum(air_kpss$partial_sums^2) / air_T^2,
  `Bartlett estimate` = air_kpss$lrv,
  KPSS = air_kpss$statistic,
  `5% cutoff` = unname(air_kpss$critical["5%"]), check.names = FALSE),
  digits = c(6, 6, 3, 3), format.args = list(scientific = FALSE),
  caption = "An affine trend is fitted; Bartlett bandwidth H = 8.")

An affine trend is fitted; Bartlett bandwidth H = 8.

Squared sums / T²	Bartlett estimate	KPSS	5% cutoff
0.000892	0.016525	0.054	0.146

The following plot shows \(R_t/\sqrt{T\widehat v_8}\), where \(R_t=\sum_{j=1}^t\widehat U_j\) and \(T=132\). The KPSS statistic is its average squared height: \[ K_T=\frac1T\sum_{t=1}^T \left(\frac{R_t}{\sqrt{T\widehat v_8}}\right)^2. \]

par(mar = c(4, 4.5, 2, 1), bg = "white")
plot(ts(air_scaled_sums, start = start(air_growth), frequency = 12),
     col = "#00539B", xlab = "Year", ylab = "Scaled cumulative residual",
     main = "How large are the cumulative residuals after scaling?")
abline(h = 0, col = "gray60", lty = 2)

Figure 6: Cumulative residuals divided by the estimated partial-sum scale. Their average squared height is the KPSS statistic. The endpoint is zero because OLS fitted an intercept. Real data: datasets::AirPassengers.

With bandwidth eight, the statistic is 0.054, below the fitted-trend 5% cutoff \(0.146\). KPSS therefore does not reject stationary fluctuations around a linear trend in annual growth. This calculation supplies no evidence that an additional difference is needed.

The fitted residuals describe departures from the trend in annual growth. An additional difference, \((1-B)W_t\), would instead describe month-to-month changes in annual growth.

Does the Series Have a Unit Root?

KPSS starts from a stationary null. The augmented Dickey–Fuller (ADF) test starts from accumulation. First suppose the population mean is known to be zero. We compare \[ \begin{aligned} H_0:&\quad X_t=X_0+\sum_{j=1}^{t}V_j, &&(V_t)\text{ is mean-zero and }I(0),\\ H_1:&\quad (X_t)\text{ is mean-zero and }I(0). \end{aligned} \] Under \(H_0\) we use the centered \(I(1)\) setting of Definition 9.3. The distinction is whether stationary fluctuations accumulate in the observed levels.

Under \(H_0\), \((1-B)X_t=V_t\): one difference recovers an \(I(0)\) process. The difference polynomial \(1-z\) has a root at \(z=1\). In a linear model for the levels, this appears as an uncancelled factor \(1-B\) after any common factors on the two sides of the model equation have been removed. This is the unit root tested by ADF. Under \(H_1\), the levels are already \(I(0)\). Their differences are also weakly stationary, so finding stationary-looking differences alone does not distinguish the hypotheses.

The increments under \(H_0\) may be serially correlated. ADF regresses changes on the preceding level and on earlier changes. The earlier changes approximate this dependence; the coefficient on the preceding level is used to test for a unit root. The number of earlier changes included is chosen for the regression, without specifying an ARMA order for the underlying process.

Definition 9.5 (Augmented Dickey–Fuller statistic). Given \(X_0,\ldots,X_T\) and a nonnegative integer lag count \(k\), fit by least squares \[ \nabla X_t=bX_{t-1}+\sum_{j=1}^{k}c_j\nabla X_{t-j}+E_{t,k}, \qquad t=k+1,\ldots,T. \] Here \(\nabla X_t=X_t-X_{t-1}\). Let \(\widehat b_{T,k}\) be the fitted coefficient on \(X_{t-1}\) and \(\widehat{\operatorname{se}}(\widehat b_{T,k})\) its conventional OLS standard error. When the regression is identified and this standard error is positive, the ADF statistic is the random variable \[ \operatorname{ADF}_{T,k} :=\frac{\widehat b_{T,k}} {\widehat{\operatorname{se}}(\widehat b_{T,k})}. \]

When the increments admit an autoregressive representation, the unit-root null makes the coefficient on the lagged level zero. The ADF regression uses finitely many lagged changes to approximate that representation. A sufficiently negative standardized coefficient gives evidence against \(H_0\) in favor of the stationary alternative.

Calibration requires more than the moment conditions defining \(I(0)\). One sufficient setting is that the stationary process in either hypothesis is a causal, invertible ARMA process of unknown finite order, driven by i.i.d. mean-zero innovations with positive variance and finite fourth moment, with the number of augmentation lags increasing suitably with the sample size.² These dependence assumptions are maintained under both hypotheses; ADF does not test whether the innovations are independent.

Choosing the Regression Terms and Critical Value

As with KPSS, the critical values come from advanced asymptotic theory; we use the supplied values without deriving them. Although the statistic has the form of a regression \(t\) statistic, ordinary regression critical values do not apply. Reject when the ADF statistic is below its supplied lower-tail critical value.

For an unknown constant mean, include an intercept in the ADF regression. For an unknown linear mean, include an intercept and time. Each choice has its own critical values, just as fitting the mean changes the KPSS reference. The conventional specifications are:

Specification	Terms fitted	Unit-root null	Stationary alternative
none	No intercept or trend	Centered accumulation with no drift	Known mean zero
constant	Intercept	Accumulation with no drift	Unknown constant mean
trend	Intercept and time	Accumulation with possible linear drift	Unknown linear mean

Fit these terms within the ADF regression and use its matching reference; subtracting a fitted mean first does not justify the none critical values.

In a finite record, inspect the regression residual ACF and compare a few choices of \(k\). Additional lags can capture more dependence but require estimating more coefficients. Use wb_adf(x, lags = k, deterministic = "none"), replacing "none" by "constant" or "trend" as appropriate. The function returns the statistic, matching critical values, and residuals.

How Should We Read the Two Tests Together?

The tests begin from different null hypotheses:

Test	Null hypothesis within the specified model	Rejection direction
KPSS	An \(I(0)\) remainder around the specified deterministic mean, known or fitted	Statistic above its critical value
ADF	A unit root in the level	Statistic below its critical value

An ADF rejection is evidence against a unit root under the maintained assumptions. A KPSS rejection is evidence against its stationary null. Neither decision directly compares the covariance structure at different times in the record.

Remark 9.6 (A change in dependence). A process can keep the same \(N(0,1)\) marginal distribution at every time while its lag-one covariance changes. It then fails weak stationarity without accumulating a random walk. Exercise 9.3 constructs such a process and compares the test decisions with the sample ACF in each half of the record. It asks whether an ADF rejection together with a KPSS nonrejection can distinguish this process from a stationary one.

For annual passenger growth, trend KPSS at bandwidth eight did not reject the model \(W_t=a+bt+U_t\) with \((U_t)\) mean-zero and \(I(0)\). Exercise 9.4 repeats trend KPSS at other bandwidths, applies ADF to \(W_t\) with a fitted trend, and compares the residuals with \((1-B)W_t\). These comparisons assess whether another ordinary difference is warranted.

The workbook on cointegration asks a different question: can a linear combination of two accumulating series have a stationary remainder?

Sources

Dickey and Fuller (1979) develop the unit-root reference distribution. Chang and Park (2002) give the ADF calibration for autoregressive approximations to invertible linear processes. MacKinnon supplies the numerical approximations used by the helper: the 1996 coefficients for none and the 2010 coefficients for constant and trend. Kwiatkowski, Phillips, Schmidt, and Shin (1992) introduce KPSS and its fitted-mean critical values. Kagalwala (2022), Table 2 also gives the known-zero-mean reference used here. The Stata dfuller manual, page 3, and Zivot’s unit-root notes distinguish the deterministic reference cases. The folded test-function code records its implementation sources and attribution.

Shumway and Stoffer (2025), Time Series Analysis and Its Applications, 5th ed., Section 5.2, discusses integrated models and unit-root tests. The R documentation for AirPassengers identifies the monthly international passenger totals, their units, and the Box–Jenkins source.

Exercises

Exercise 9.1 (Does more data make the test valid?). Use the known-zero-mean KPSS statistic in Definition 9.2. Consider the stationary AR(1) \[ X_t=0.5X_{t-1}+Z_t,\qquad (Z_t)\text{ i.i.d. }N(0,1). \] Work by hand in (a)–(b) and (e), and use R in (c)–(d).

1. Show that replacing \(X_t\) by \(cX_t\), where \(c\ne0\), leaves \(K_T\) unchanged when the same bandwidth is used.

2. Derive the marginal variance \(\gamma(0)\) and long-run variance \(v\) of this AR(1). If the denominator \(T^2v\) were replaced by \(T^2\gamma(0)\), by what factor would \(\sum_{t=1}^T S_t^2/(T^2v)\) change?

3. Using seed 54903, simulate 5,000 independent records \(X_1,\ldots,X_{1600}\). Start each at \(X_0\sim N(0,\gamma(0))\), independently of its subsequent innovations. For \(T\in\{200,400,800,1600\}\), use the first \(T\) observations of each record to compare the bandwidth rules \[ H_T=8 \qquad\text{and}\qquad H_T=\left\lceil12(T/200)^{1/3}\right\rceil. \] Use wb_kpss() with deterministic = "none". For each \(T\) and rule, report \(H_T\), the fraction \(\widehat p\) of statistics above 1.66, and its Monte Carlo standard error \(\sqrt{\widehat p(1-\widehat p)/5000}\). These rejections are errors. How close are their frequencies to 5%, relative to the Monte Carlo standard errors?

4. Repeat (c) for 5,000 independent random walks \(X_t=X_{t-1}+Z_t\), with \(X_0=0\) and i.i.d. \(N(0,1)\) increments, using seed 54905. These rejection frequencies estimate the test’s power. How does detection change with \(T\), and what is the cost of using the larger bandwidth?

5. For the stationary AR(1), use fixed-lag covariance consistency to find the probability limit of the Bartlett estimate with \(H_T=8\). Compare this limit with \(v\) and explain the consequence for false rejections as \(T\) grows. Which bandwidth rule satisfies the conditions in The KPSS Null Distribution box, and what do those conditions imply about the limiting false-rejection probability? Explain why improving random-walk detection alone does not establish asymptotic validity.

Exercise 9.2 (Was differencing necessary?). Let \((U_t)_{t\in\Z}\) be mean-zero and \(I(0)\), with ACVF \(\gamma_U\).

1. Find the ACVF and long-run variance of \(D_t:=U_t-U_{t-1}\). Simplify \(\sum_{t=1}^{T}D_t\) and show that its variance stays bounded as \(T\) increases.

2. Now let \((V_t)\) be mean-zero and \(I(0)\), with long-run variance \(v_V>0\), and let \(X_t=\sum_{j=1}^{t}V_j\), with \(X_0=0\). Find the long-run variance of \(\nabla X_t\) and the limiting value of \(\Var(X_T)/T\). Both \(D_t\) and \(\nabla X_t\) are weakly stationary. Use their long-run variances to explain why \((X_t)\) is \(I(1)\) but \((U_t)\) is not, under Definition 9.3.

3. Suppose \(W_t=a+bt+U_t\), with \((U_t)\) as above. An analyst applies the constant-mean KPSS test to \(\nabla W_t\). Does the population model for these differences satisfy the positive-long-run-variance condition needed for the usual KPSS reference? Explain why a small KPSS statistic would not settle whether the additional difference was needed.

Exercise 9.3 (Can the tests miss a change in dependence?). We compare a stationary process with one whose lag-one covariance changes while its marginal distributions stay the same. Let \((Z_t)_{t\in\Z}\) be i.i.d. \(N(0,1)\), \(\theta=0.8\), and \(m=400\). The stationary process is \[ C_t=\frac{Z_t+\theta Z_{t-1}}{\sqrt{1+\theta^2}}. \] The changing process has \(X_t=Z_t\) for \(t\le m\) and \(X_t=C_t\) for \(t>m\).

1. Verify that every \(X_t\) and \(C_t\) has the same marginal distribution. Calculate the lag-one covariance strictly within each half of each process.

2. Simulate 1,000 independent pairs of records of length 800 using seed 54904; within each pair, use the same innovations. For the first pair, plot both records and the sample ACF in each half through lag ten. Compare the plots with the population covariances from (a).

For all 1,000 pairs, use constant-only ADF with four lagged differences and constant-only KPSS with bandwidth eight. Report their rejection fractions at the supplied 5% critical values for each process.

3. An analyst calls a record stationary when ADF rejects and KPSS does not. Would this rule distinguish these two processes? Use the covariance calculation, plots, and rejection fractions to explain what the rule misses.

Exercise 9.4 (Does annual growth need another difference?). Use \(W_t=(1-B^{12})\log P_t\) from datasets::AirPassengers, the 132 annual log-growth observations from January 1950 through December 1960. The worked example fitted \(W_t=a+bt+U_t\) and used trend KPSS at bandwidth eight. Assess the case for taking an additional difference.

1. Plot \(\nabla W_t\) and its sample ACVF through lag 24. Compare them with the fitted residuals and sample ACVF in the worked example. What does each series measure, and how has differencing changed the visible dependence?

2. Apply trend KPSS to \(W_t\) with bandwidths \(H\in\{4,8,12,24\}\). Report the Bartlett estimates, test statistics, and 5% decisions. Explain how changing the bandwidth changes the statistic. Why must these calculations use the trend-KPSS reference even though OLS residuals have sample mean zero?

3. Apply ADF to \(W_t\) with an intercept and time trend, using \(k\in\{0,4,8,12\}\) lagged differences. Report the statistics, 5% cutoffs, and decisions. Plot the regression residual ACFs for \(k=0\) and \(k=12\). What do these plots suggest about dependence left in the regressions?

4. Would you retain annual log growth after subtracting a fitted trend, or take an additional difference? Support your choice using the plots and test results, and identify one remaining concern about stationarity. Use Exercise 9.2 to explain why choosing the output with the smaller KPSS statistic would not settle the question.

Footnotes

1. The known-zero-mean critical values are approximately 1.20, 1.66, and 2.79 at levels 10%, 5%, and 1%. These are rounded values from the zero-mean Bartlett KPSS reference in Kagalwala (2022), Table 2.

2. For the known-zero-mean version, one sufficient rule is \(k_T\to\infty\) and \(k_T=o(\sqrt T)\). For each fixed process in the stated class, the null rejection probability then tends to the nominal level when the matching asymptotic critical value is used. See Chang and Park (2002), Assumptions 1–3 and Theorem 3.6.