---
title: "Workbook 9 — Assessing a Stationary Approximation"
subtitle: "STA 542 · Introduction to Time Series Analysis"
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 2
    number-sections: false
    code-copy: true
    fig-width: 9
    fig-height: 3.8
    fig-dpi: 150
execute:
  warning: false
  message: false
editor:
  markdown:
    wrap: 72
---

::: hidden
$$\def\E{\mathbb{E}}
\def\P{\mathbb{P}}
\def\R{\mathbb{R}}
\def\Z{\mathbb{Z}}
\def\Var{\operatorname{Var}}
\def\Cov{\operatorname{Cov}}$$
:::

```{=html}
<style>
.exercise { border-left: 3px solid #00539B; padding: 0.7em 1.1em; background: #f6f8fa; margin: 1.4em 0; border-radius: 0 4px 4px 0; }
.exercise p:first-child { margin-top: 0; }
figure figcaption { color: #555; font-size: 0.9em; }
</style>
```

```{r}
#| label: setup
#| include: false
knitr::opts_chunk$set(comment = "#>")
```

## 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.

```{r}
#| label: fig-wb9-persistent-paths
#| code-fold: true
#| code-summary: "Simulation code"
#| fig-cap: "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."
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")
```

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.

```{r}
#| label: fig-wb9-persistent-acvf
#| code-fold: true
#| code-summary: "Sample ACVF of the two records"
#| fig-cap: "Sample ACVF of the two records above, using divisor T. The panels share a vertical scale. Simulated data, seed 542."
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")
```

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.

```{r}
#| label: fig-wb9-kpss-paths
#| code-fold: true
#| code-summary: "Simulate centered series and their partial sums"
#| fig-height: 5.4
#| fig-cap: "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."
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")
}
```

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
\emph{sample} 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 \emph{population} 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.[^kpss-zero]

::: {.callout-note collapse="true" #kpss-reference}
## More 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.
:::

[^kpss-zero]: 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](https://doi.org/10.1177/1536867X221106371).

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.

```{r}
#| label: wb9-test-functions
#| code-fold: true
#| code-summary: "R functions for the tests"
# 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)
}
```

```{r}
#| label: wb9-centered-kpss
#| code-fold: true
#| code-summary: "Compute KPSS with the mean known to be zero"
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.")
```

For the random walk, the squared sums divided by $T^2$ equal
`r sprintf("%.3f", sum(kpss_centered[[2]]$partial_sums^2) / T_kpss^2)`,
but the Bartlett estimate is also large, at
`r sprintf("%.3f", kpss_centered[[2]]$lrv)`.
Their ratio is only `r sprintf("%.3f", kpss_centered[[2]]$statistic)`,
below 1.66. The statistic is larger than the stationary record's
`r sprintf("%.3f", kpss_centered[[1]]$statistic)`, 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:

```{r}
#| label: wb9-kpss-critical-values
#| echo: false
knitr::kable(data.frame(
  `Mean specification` = c("Known zero", "Unknown constant", "Unknown linear trend"),
  `10%` = c(1.20, 0.347, 0.119), `5%` = c(1.66, 0.463, 0.146),
  `1%` = c(2.79, 0.739, 0.216), check.names = FALSE
), caption = "Approximate asymptotic KPSS critical values. Known zero: Kagalwala (2022), rounded. Fitted constant and trend: Kwiatkowski et al. (1992).")
```

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
\[
\widehat U_t=Y_t-\widehat a-\widehat bX_t.
\]Second, examine the persistence of those residuals. Fit
\[
\nabla\widehat U_t
=c\,\widehat U_{t-1}
+\sum_{j=1}^{k}d_j\nabla\widehat U_{t-j}
+e_t.
\]This is the same regression form used by ADF:
- The lagged level \(\widehat 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 \(\widehat c/\operatorname{se}(\widehat c)\) provides evidence against no cointegration, under the test’s assumptions.
For intuition, if a spread followed \(U_t=\phi U_{t-1}+Z_t\), its changes would satisfy
\[
\nabla U_t=(\phi-1)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.

```{r}
#| label: wb9-airline-setup
#| code-fold: true
#| code-summary: "Annual log growth and its fitted line"
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)
```

```{r}
#| label: fig-wb9-airline-question
#| code-fold: true
#| code-summary: "Plot annual log growth with its fitted line"
#| fig-cap: "Annual log growth in passenger totals and an affine trend fitted by OLS. Real data: datasets::AirPassengers."
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")
```

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.

```{r}
#| label: fig-wb9-airline-residuals
#| code-fold: true
#| code-summary: "Plot the fitted residuals and sample ACVF"
#| fig-cap: "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."
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")
```

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.

```{r}
#| label: wb9-airline-kpss
#| code-fold: true
#| code-summary: "Compute trend KPSS with bandwidth eight"
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.")
```

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.
$$

```{r}
#| label: fig-wb9-airline-kpss-path
#| code-fold: true
#| code-summary: "Plot the scaled cumulative residuals"
#| fig-cap: "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."
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)
```

With bandwidth eight, the statistic is
`r sprintf("%.3f", air_kpss$statistic)`, 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.[^adf-calibration]
These dependence assumptions are maintained under both hypotheses;
ADF does not test whether the innovations are independent.

[^adf-calibration]: 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](https://www.ruf.rice.edu/~econ/papers/2001papers/02Chang.pdf).

### 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)](https://www.ruf.rice.edu/~econ/papers/2001papers/02Chang.pdf)
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](https://doi.org/10.1177/1536867X221106371)
also gives the known-zero-mean reference used here. The Stata `dfuller`
manual, page 3, and [Zivot's unit-root notes](https://faculty.washington.edu/ezivot/econ584/notes/unitrootlecture.pdf)
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`](https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/AirPassengers.html)
identifies the monthly international passenger totals, their units, and
the Box--Jenkins source.

## Exercises

::: {.exercise}
**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).

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

(b) 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?

(c) 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?

(d) 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?

(e) 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}
**Exercise 9.2 (Was differencing necessary?).** Let $(U_t)_{t\in\Z}$
be mean-zero and $I(0)$, with ACVF $\gamma_U$.

(a) 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.

(b) 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.

(c) 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}
**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$.

(a) 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.

(b) 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.

(c) 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}
**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.

(a) 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?

(b) 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?

(c) 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?

(d) 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.
:::