---
title: "Workbook 10 — Multivariate Stationarity and Cointegration"
subtitle: "STA 542 · Introduction to Time Series Analysis"
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 2
    toc-location: right
    number-sections: false
    code-copy: true
    fig-width: 8
    fig-height: 3.2
    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>
:root { --wb-ink: #173753; --wb-accent: #a35f32; }
#quarto-document-content { max-width: 820px; }
body { color: #25313b; line-height: 1.62; }
h1, h2, h3 { color: var(--wb-ink); }
h1.title { font-weight: 650; letter-spacing: -.025em; }
h2 { margin-top: 2.15rem; padding-bottom: .4rem; border-bottom: 1px solid #dce4ea; }
.subtitle { color: #607080; }
a { color: #145784; }
.exercise { border-left: 3px solid var(--wb-ink); padding: .9em 1.2em; background: #f3f6f8; margin: 1.5em 0; border-radius: 0 5px 5px 0; }
.exercise p:first-child { margin-top: 0; }
figure { margin: 1.25rem 0 1.6rem; }
figure figcaption { color: #586572; font-size: .88em; line-height: 1.45; }
details.code-fold > summary { color: #586572; font-size: .88em; }
.callout-tip { --bs-callout-color: #173753; }
table { font-size: .94em; }
@media print {
  #TOC, .quarto-title-meta, .code-copy-button, details.code-fold { display: none !important; }
  body { font-size: 10pt; }
  h2 { break-after: avoid; }
  figure, table, .exercise { break-inside: avoid; }
}
</style>
```

```{r}
#| label: setup
#| include: false
knitr::opts_chunk$set(comment = "#>")
wb_blue <- "#1C5276"
wb_copper <- "#AA5F32"
wb_gray <- "#657582"
wb_par <- function(rows = 1, columns = 1, simulated = FALSE) {
  par(mfrow = c(rows, columns), mar = c(3.6, 4.0, 2.0, 0.8),
      mgp = c(2.45, 0.7, 0), las = 1, tcl = -0.25,
      col.axis = "#43515D", col.lab = "#25313B", col.main = "#173753",
      cex.main = 1.0, font.main = 2,
      bg = if (simulated) "gray94" else "white")
}
```

## Vector Means and Weak Stationarity

Several quantities may be recorded at each time: prices of related goods,
temperatures at nearby stations, or measurements from several sensors.
We model such observations by vectors $X_t\in\R^d$. Their mean is a vector;
their covariance records both variation within coordinates and variation
between them. For random vectors $U,V$ with finite second moments, define
$$
\Cov(U,V):=\E[(U-\E[U])(V-\E[V])^\top].
$$
Its $(i,j)$ entry is $\Cov(U_i,V_j)$.

The definition of weak stationarity for vectors is directly analogous to that for scalars.

**Definition 10.1 (Vector weak stationarity).** A process $(X_t)$ in
$\R^d$ is *weakly stationary* if it has finite second moments, a constant
mean vector $\mu$, and
$$
\Cov(X_{t+h},X_t)=\Gamma(h)
$$
for every admissible $t,h$, where the matrix $\Gamma(h)$ depends only on
the lag. The function $\Gamma$ is its *matrix autocovariance function*.

The entry $\Gamma_{ij}(h)$ compares coordinate $i$ at time $t+h$ with
coordinate $j$ at time $t$. An important difference between the scalar and vector definitions:Transposing the covariance gives
$\Gamma(-h)=\Gamma(h)^\top$. Thus $\Gamma(0)$ is symmetric and positive
semidefinite; $\Gamma(h)$ need not be symmetric at a nonzero lag (note that this means that it generalizes the scalar definition regardless).
Stationarity of each coordinate checks only the diagonal entries. Joint
stationarity also requires the cross-covariances to depend only on lag.

## Filters Across Time and Coordinates

The backshift operator $B$ moves an observation one step backward in time for vectors all the same.
For a vector series, it shifts every coordinate together:
$$
BX_t:=X_{t-1}=(X_{1,t-1},\ldots,X_{d,t-1})^\top,
\qquad B^jX_t=X_{t-j}.
$$
A scalar filter combines observations at different times. With vectors,
matrix coefficients can also combine coordinates.

**Definition 10.2 (Finite matrix filter).** For fixed real $r\times d$
matrices $C_0,\ldots,C_m$, the filtered process is
$$
W_t:=\sum_{j=0}^{m}C_jB^jX_t
=\sum_{j=0}^{m}C_jX_{t-j}\in\R^r.
$$
The output dimension $r$ can differ from the input dimension $d$.

The difference filter uses $C_0=I_d$, $C_1=-I_d$ and produces
$\nabla X_t=(1-B)X_t=X_t-X_{t-1}$, where $I_d$ is the identity matrix.
It differences each coordinate separately.

**Proposition 10.3 (Weak stationarity under finite matrix filtering).**
Let $(X_t)$ be a weakly stationary process in $\R^d$ with mean $\mu$
and matrix ACVF $\Gamma$. For any fixed real $r\times d$ matrices
$C_0,\ldots,C_m$, the process $(W_t)$ in Definition 10.2 is weakly
stationary, with mean
$$
\mu_W=\sum_{j=0}^{m}C_j\mu
$$
and matrix ACVF
$$
\Gamma_W(h)=\sum_{j=0}^{m}\sum_{k=0}^{m}
C_j\Gamma(h-j+k)C_k^\top.
$$
These formulas hold wherever all required observations are defined.

::: {.callout-note collapse="true"}
## Proof

Each coordinate of $W_t$ is a finite linear combination of random
variables with finite second moments, so it also has a finite second
moment. By linearity,
$$
\E[W_t]=\sum_{j=0}^{m}C_j\E[X_{t-j}]
=\sum_{j=0}^{m}C_j\mu=\mu_W,
$$
which does not depend on $t$.

For the covariance, the lag between $t+h-j$ and $t-k$ is $h-j+k$.
Thus
$$
\begin{aligned}
\Cov(W_{t+h},W_t)
&=\sum_{j=0}^{m}\sum_{k=0}^{m}
C_j\Cov(X_{t+h-j},X_{t-k})C_k^\top\\
&=\sum_{j=0}^{m}\sum_{k=0}^{m}
C_j\Gamma(h-j+k)C_k^\top.
\end{aligned}
$$
This depends only on $h$, so $(W_t)$ is weakly stationary. $\square$
:::

In a similar way as for scalars, we can define $I(1)$ and $I(0)$ for vectors.

**Definition 10.4 (Vector integration orders).** A process $(X_t)$ in
$\R^d$ is *$I(0)$* if it is jointly weakly stationary, its matrix ACVF
satisfies $\sum_{h\in\Z}|\Gamma_{ij}(h)|<\infty$ for every $i,j$, and
its long-run covariance matrix
$$
\Omega:=\sum_{h\in\Z}\Gamma(h)
$$
has positive diagonal entries. We allow $\Omega$ to be singular.
It is *$I(1)$* if
$$
X_t=X_0+\sum_{j=1}^{t}V_j,\qquad t\ge1,
$$
where $(V_t)$ is a mean-zero $I(0)$ vector process and $X_0$ has finite
second moments.

Just like with scalars, first differencing makes $I(1)$ vector processes into $I(0)$ vector processes. However, in the multivariate setting, there is a further possibility.
Even when every coordinate of $X_t$ is $I(1)$, a fixed linear combination
$\beta^\top X_t$ can be $I(0)$. This filter combines coordinates at the
same time: $m=0$ and $C_0=\beta^\top$. Before we proceed with differencing and developing the theory of vector AR and MA models, we will explore this additional possibility to obtain stationary processes.

**Definition 10.5 (Cointegration).** Suppose every coordinate of
$X_t\in\R^d$ is $I(1)$. The coordinates are *cointegrated* if there is a
fixed nonzero vector $\beta\in\R^d$ such that $\beta^\top X_t$ is $I(0)$.
We call $\beta$ a *cointegrating vector*.

The scalar combination is called a *spread*. Multiplying $\beta$ by a
nonzero constant describes the same relation. Subtracting a constant from
the spread only changes its mean.

**Example 10.6 (Cancelling a stochastic trend).** Let $(Z_t)$, $(U_t)$,
and $(V_t)$ be mutually independent $I(0)$ processes, indexed by
$t\in\Z$, with $\E[Z_t]=0$. Write their positive long-run variances as
$\tau_Z^2$, $\tau_U^2$, and $\tau_V^2$, respectively. Define
$$
S_0:=0,\qquad S_t:=S_{t-1}+Z_t,\quad t\geq1,
$$
and consider
$$
X_t:=S_t+U_t,\qquad Y_t:=2S_t+V_t.
$$
Their differences are
$$
\nabla X_t=Z_t+\nabla U_t,\qquad
\nabla Y_t=2Z_t+\nabla V_t.
$$
These are mean-zero, weakly stationary processes with absolutely summable
autocovariances. Their long-run variances are $\tau_Z^2$ and
$4\tau_Z^2$, respectively, because $\nabla U_t$ and $\nabla V_t$ have
zero long-run variance. Both differences are therefore $I(0)$, so both
levels are $I(1)$. Yet
$$
Y_t-2X_t=V_t-2U_t
$$
is also $I(0)$, with long-run variance $\tau_V^2+4\tau_U^2>0$. The common
accumulated component cancels, so the coordinates are cointegrated with
cointegrating vector $(-2,1)^\top$.

The gap can be serially correlated, so its history may help predict it...
 $\diamond$

The requirement for $(Y_t,X_t)$ to be $I(1)$ in each coordinate is crucial for cointegration.
There are many processes that posses a common trend that are not cointegrated, and their behavior
is fundamentally different.

**Example 10.7 (Cancelling a deterministic trend).** Now replace $S_t$
by $t$:
$$
X_t:=t+U_t,\qquad Y_t:=2t+V_t.
$$
Neither process is stationary, because its mean changes with time. Their
combination is again
$$
Y_t-2X_t=V_t-2U_t,
$$
which is $I(0)$. Nevertheless, the coordinates are not cointegrated.
Neither is $I(1)$: their variances remain constant, whereas an $I(1)$
process has variance growing asymptotically linearly. A stationary
combination of nonstationary processes is therefore not sufficient for
cointegration. $\diamond$

The mathematical difference between the two examples is subtle: both processes result in a stationary.
However, in the first case, we can draw a specific conclusion about a *statistical relationship* between the two processes.
In Example 10.6, the relation between the levels reflects a *shared
stochastic component that is (potentially) unobserved*: both processes contain $S_t$. They are $I(1)$,
and subtracting $2X_t$ from $Y_t$ cancels that component, leaving an
$I(0)$ gap.

In Example 10.7, a regression can fit the levels closely because both
means increase with time, even though the processes are independent.
Interpreting this fit as evidence of some kind of dependence (a statistical relationship) would be wrong!
More generally, these so called *spurious regressions* refers to misleading evidence
from regression when series happen to have a common trend.
It is not defined by whether the gap is stationary: here the gap is
$I(0)$, but the coordinates are not $I(1)$ and hence are not cointegrated: there was no *common stochastic component* to cancel,
the processes just happened to travel in "at the same rate" over time.

Given observations, how might we find a candidate cointegrating
coefficient? A *levels regression* fits the undifferenced observations
$y_t$ on an intercept and $x_t$:
$$
y_t\approx a+cx_t.
$$
Least squares chooses the coefficients to make the observed squared gaps
small. This criterion establishes neither that the original processes
are $I(1)$ nor that $Y_t-a-cX_t$ is $I(0)$. Trends and persistence can produce high $R^2$ values and apparently
significant regression coefficients even for independent processes. This misleading regression evidence is called spurious regression. 

## Spurious Regression

Australian air passengers and rice production in Guinea make an unusual
regression pair. Hyndman and Athanasopoulos (2021, Section 7.3) use these
42 annual World Bank observations, 1970–2011, to illustrate the effect of
trends on a *levels regression*, which uses the undifferenced values of
the series. Here we regress each year's passenger numbers on an intercept
and that year's rice production. Regressing their year-to-year changes on
each other would instead be a regression in differences. Passengers are
measured in millions and rice production in million metric tons.

```{r}
#| label: air-rice-data
#| code-fold: true
#| code-summary: "Load the 42 observations and fit the line"
air_rice <- read.csv(text = '"Year","Passengers","Production"
1970,7.3187,0.311405
1971,7.3266,0.32518
1972,7.7956,0.339565
1973,9.3846,0.35458
1974,10.6647,0.370265
1975,11.0551,0.38664
1976,10.8643,0.40374
1977,11.3065,0.4216
1978,12.1223,0.44025
1979,13.0225,0.45972
1980,13.6488,0.480055
1981,13.2195,0.50129
1982,13.1879,0.52346
1983,12.6015,0.54661
1984,13.2368,0.57079
1985,14.4121,0.596035
1986,15.4973,0.622395
1987,16.8802,0.649925
1988,18.8163,0.67867
1989,15.1143,0.70869
1990,17.5534,0.740035
1991,21.8601,0.772765
1992,23.8866,0.806945
1993,26.9293,0.842635
1994,26.8885,0.879905
1995,28.8314,0.918825
1996,30.0751,0.95946
1997,30.9535,1.0019
1998,30.1857,1.04622
1999,31.5797,1.09249
2000,32.577569,1.14081
2001,33.477398,1.03352
2002,39.021581,1.08867
2003,41.386432,1.14676
2004,41.596552,1.20796
2005,44.657324,1.27242
2006,46.951775,1.34031
2007,48.728837,1.40159
2008,51.488427,1.53409
2009,50.026967,1.45593
2010,60.640913,1.49896
2011,63.3603103378,1.5')
fit_air_rice <- lm(Passengers ~ Production, data = air_rice)
```

```{r}
#| label: fig-air-rice
#| fig-height: 2.9
#| code-fold: true
#| code-summary: "R code for the two views of the same observations"
#| fig-cap: "Real data: Australian air passengers and Guinea rice production, 1970–2011, World Bank via Hyndman–Athanasopoulos. Both series rise over time, producing a high-fit levels regression between the two measurements."
wb_par(columns = 2)
indexed <- 100 * sweep(as.matrix(air_rice[, c("Passengers", "Production")]),
                      2, as.numeric(air_rice[1, c("Passengers", "Production")]), "/")
matplot(air_rice$Year, indexed, type = "l", lty = 1, lwd = 1.7,
        col = c(wb_blue, wb_copper), xlab = "Year", ylab = "Index (1970 = 100)",
        main = "Growth over time")
legend("topleft", c("Air passengers", "Rice production"),
       col = c(wb_blue, wb_copper), lty = 1, lwd = 1.7, bty = "n", cex = .8)
plot(air_rice$Production, air_rice$Passengers, pch = 19, cex = .65,
     col = wb_blue, xlab = "Rice production (million metric tons)",
     ylab = "Air passengers (millions)", main = "The levels regression")
abline(fit_air_rice, col = wb_copper, lwd = 1.8)
par(mfrow = c(1, 1))
```

The fitted line has $R^2=$ `r round(summary(fit_air_rice)$r.squared, 3)`:
it accounts for about `r round(100 * summary(fit_air_rice)$r.squared, 1)`%
of the variation in passenger numbers in this sample. Both series rise
over time, so later years tend to have both more rice production and more
air passengers. A line can therefore fit the observations closely even
when the gap around that line is nonstationary. A high $R^2$ also does not
establish that rice production causes air travel.

*Spurious regression* occurs when a regression gives misleading evidence
of a relationship because it ignores trends or persistence. Even
independent persistent series can yield a high $R^2$ and apparently
significant coefficients: conventional regression tests can reject far
more often than their stated significance levels suggest. The following
example uses independent random walks, for which we can check directly
whether any fixed nonzero linear combination is stationary.

**Example 10.8 (Independent trends).** Consider two independent random
walks starting at $A_0=B_0=0$:
$$
A_t=A_{t-1}+\varepsilon_t,\qquad
B_t=B_{t-1}+\eta_t,\qquad t\ge1,
$$
where $(\varepsilon_t)$ and $(\eta_t)$ are independent sequences of
independent $N(0,1)$ variables. First differencing recovers the shocks:
$\nabla A_t=\varepsilon_t$ and $\nabla B_t=\eta_t$, so both differences
are stationary white noise.

Each level is the sum of $t$ independent shocks of variance one, so
$\Var(A_t)=\Var(B_t)=t$. Independence also gives $\Cov(A_t,B_t)=0$.
Thus, for any fixed coefficients $(c_1,c_2)\ne(0,0)$,
$$
\Var(c_1A_t+c_2B_t)
=c_1^2\Var(A_t)+c_2^2\Var(B_t)
=t(c_1^2+c_2^2).
$$
This variance grows with $t$, whereas weak stationarity requires a
constant variance. No nonzero fixed linear combination of the levels is
therefore stationary. In particular, there are no fixed constants $a,b$
for which the gap $B_t-a-bA_t$ is stationary; subtracting an intercept
does not change the variance.

Over a finite record, the independent paths can nevertheless happen to
move together, giving a high $R^2$ when $b_t$ is regressed on an intercept
and $a_t$ (a levels regression). The usual regression $t$-test for a zero
slope can then reject
far more often than its stated significance level suggests: the
nonstationary regression error violates the test's assumptions.
Exercise 10.2 investigates how often independent random walks produce a
high $R^2$ over repeated records. $\diamond$

## Cointegration and Common Trends


Recall that an $I(0)$ process is weakly stationary, has absolutely
summable autocovariances, and has positive long-run variance. An $I(1)$
process accumulates mean-zero $I(0)$ increments:
$$
X_t=X_0+\sum_{j=1}^{t}V_j,\qquad t\ge1,
$$
where $(V_t)$ is mean-zero and $I(0)$, and $X_0$ has finite variance.
First differencing recovers these increments: $\nabla X_t=V_t$.

Here we allow the initial value $X_0$ to have a nonzero mean. Adding a
constant $c$ to $X_0$ shifts every $X_t$ by the same amount; it changes
neither the differences nor the variances.

Let $v_V:=\sum_{h\in\Z}\gamma_V(h)>0$ be the long-run variance of the
increments, where $\gamma_V$ is their ACVF. The partial-sum variance
result gives
$$
\frac{\Var(X_t)}{t}\longrightarrow v_V.
$$
Thus $\Var(X_t)$ is approximately $t v_V$ for large $t$. Its growth rules
out weak stationarity, which requires a constant variance. Under this
convention, an $I(1)$ process is therefore not itself $I(0)$.

::: {.exercise}
**Exercise 10.3 (Can stationary series be cointegrated?).** Let $(U_t)$
and $(V_t)$ be independent sequences of independent $N(0,1)$ variables,
indexed by $t\in\mathbb Z$. Define
$$
X_t=U_t,\qquad Y_t=2U_t+V_t.
$$
A student claims that $X_t$ and $Y_t$ are cointegrated because
$Y_t-2X_t=V_t$ is stationary.

(a) Determine whether $(X_t)$, $(Y_t)$, and $(Y_t-2X_t)$ are $I(0)$.
    Justify each answer using the definition.

(b) Assess the student's claim under Definition 10.5.

(c) Show that $(Y_t-bX_t)$ is $I(0)$ for every fixed $b\in\mathbb R$.
    Can stationarity identify the coefficient $2$? What distinguishes
    $2$ from the other coefficients?
:::

**Example 10.9 (Modeling the cointegrating spread).** We specialize the
construction of Example 10.6. Let $(\eta_t)$, $(\varepsilon_t)$, and
$(\xi_t)$ be mutually independent sequences, each i.i.d. with mean zero
and positive finite variance $\sigma_\eta^2$, $\sigma_\varepsilon^2$, and
$\sigma_\xi^2$, respectively. Set $S_0=0$ and
$S_t=S_{t-1}+\eta_t$ for $t\ge1$. Independently of $S$, let $U_t$ and $V_t$ be the
stationary causal solutions of
$$
U_t=\phi_1U_{t-1}+\varepsilon_t,\qquad
V_t=\phi_2V_{t-1}+\xi_t,\qquad |\phi_1|,|\phi_2|<1.
$$
For constants $a$ and $b\ne0$, define
$$
X_t=S_t+U_t,\qquad Y_t=a+bS_t+V_t.
$$
We observe $X_t$ and $Y_t$; the common trend and the stationary
components are not separately recorded. Neither coordinate equals $S_t$.

Write
$$
q_U:=\Var(U_t)=\frac{\sigma_\varepsilon^2}{1-\phi_1^2},\qquad
q_V:=\Var(V_t)=\frac{\sigma_\xi^2}{1-\phi_2^2}.
$$
The level variances are $\Var(X_t)=t\sigma_\eta^2+q_U$ and
$\Var(Y_t)=b^2t\sigma_\eta^2+q_V$. Their differences are $I(0)$ with
long-run variances $\sigma_\eta^2$ and $b^2\sigma_\eta^2$, respectively:
the differenced stationary terms $\nabla U_t$ and $\nabla V_t$ contribute
zero long-run variance. Both levels are therefore $I(1)$.

For known $a,b$, the observations give the mean-zero spread
$$
E_t:=Y_t-a-bX_t=V_t-bU_t.
$$
Independence of $U$ and $V$ gives
$$
\gamma_E(h)=q_V\phi_2^{|h|}+b^2q_U\phi_1^{|h|},\qquad
\tau_E^2=\frac{\sigma_\xi^2}{(1-\phi_2)^2}
  +\frac{b^2\sigma_\varepsilon^2}{(1-\phi_1)^2}>0.
$$
Thus $E_t$ is $I(0)$, and $(-b,1)^\top$ is a cointegrating vector for
$(X_t,Y_t)^\top$.

Compare a fixed regression line in this model with one between the
independent random walks of Example 10.8. For a trial slope $c$, the
residual here is $Y_t-a-cX_t=(b-c)S_t+V_t-cU_t$. The residual variances are
$$
\begin{aligned}
\Var(B_t-a-cA_t)&=t(1+c^2),\\
\Var(Y_t-a-cX_t)&=t(b-c)^2\sigma_\eta^2+q_V+c^2q_U.
\end{aligned}
$$
For the independent walks, every fixed slope leaves accumulating error.
In the shared-trend model, choosing $c=b$ cancels that accumulation and
leaves the stationary error $E_t$. A high sample $R^2$ can occur in either
example; cointegration requires a fixed population relation with an
$I(0)$ error.

The simulation below uses normal innovations and independent stationary
initial values $U_0\sim N(0,q_U)$ and $V_0\sim N(0,q_V)$, independent of
future innovations.

```{r}
#| label: shared-trend-data
#| code-fold: true
#| code-summary: "R code to simulate the common trend and stationary components"
set.seed(5421005)
n_shared <- 200L # n_shared is the sample size T in the text.
b_shared <- 1.3
a_shared <- 2
phi1_shared <- .35
phi2_shared <- .8
sd_eta <- .7
sd_epsilon <- .3
sd_xi <- .3
q_u_shared <- sd_epsilon^2 / (1 - phi1_shared^2)
q_v_shared <- sd_xi^2 / (1 - phi2_shared^2)
s_shared <- cumsum(rnorm(n_shared, sd = sd_eta))
u_shared <- v_shared <- numeric(n_shared)
u_previous <- rnorm(1, sd = sqrt(q_u_shared))
v_previous <- rnorm(1, sd = sqrt(q_v_shared))
for (j in seq_len(n_shared)) {
  u_shared[j] <- phi1_shared * u_previous + rnorm(1, sd = sd_epsilon)
  v_shared[j] <- phi2_shared * v_previous + rnorm(1, sd = sd_xi)
  u_previous <- u_shared[j]
  v_previous <- v_shared[j]
}
x_shared <- s_shared + u_shared
y_shared <- a_shared + b_shared * s_shared + v_shared
e_shared <- y_shared - a_shared - b_shared * x_shared
```

```{r}
#| label: fig-common-trend
#| fig-height: 2.9
#| code-fold: true
#| code-summary: "R code for the observed levels and their cancelling spread"
#| fig-cap: "Simulated data: 200 observations, a = 2, b = 1.3, phi_1 = 0.35, phi_2 = 0.8, normal innovation SDs 0.7, 0.3, and 0.3, seed 5421005. Both levels contain a random walk and a stationary component. Their known spread is V minus bU."
wb_par(columns = 2, simulated = TRUE)
matplot(cbind(y_shared, a_shared + b_shared * x_shared), type = "l",
        col = c(wb_blue, wb_copper), lty = 1, lwd = c(1.1, 1.4),
        xlab = "Time index", ylab = "Level (simulated units)",
        main = "The common trend in the levels")
legend("topleft", c("Y", "a + bX"), col = c(wb_blue, wb_copper),
       lty = 1, lwd = 1.5, bty = "n", cex = .85)
plot(e_shared, type = "l", col = wb_blue, lwd = 1.3,
     xlab = "Time index", ylab = "Spread (simulated units)",
     main = "The known spread: Y − a − bX")
abline(h = 0, col = "gray60", lty = 3)
par(mfrow = c(1, 1))
```

This model has several filters leading to stationary outputs. For example, both
$$
\begin{pmatrix}\nabla X_t\\ \nabla Y_t\end{pmatrix}
=\begin{pmatrix}\eta_t+\nabla U_t\\b\eta_t+\nabla V_t\end{pmatrix}
\quad\text{and}\quad
\begin{pmatrix}\nabla X_t\\Y_t-a-bX_t\end{pmatrix}
=\begin{pmatrix}\eta_t+\nabla U_t\\E_t\end{pmatrix}
$$
are jointly weakly stationary. The second records the current gap as one
coordinate. Further differencing that gap gives $\nabla E_t$, which is
stationary but has zero long-run variance and hence is outside our $I(0)$
class. The choice of a stationarity-producing filter affects which
quantities are directly available for modeling. Exercise 10.1 gives a
three-coordinate example with two independent stationary relations.
$\diamond$

## Prediction from the Stationary Gap

Retain the known model and parameters of Example 10.9. We can compute
$E_t=Y_t-a-bX_t$ from the observed levels. To forecast $E_{t+1}$ using
the current gap alone, consider linear predictors $rE_t$. Since the
spread has mean zero, an intercept cannot reduce the mean squared error.

Expanding the squared error gives
$$
\E[(E_{t+1}-rE_t)^2]
=\gamma_E(0)-2r\gamma_E(1)+r^2\gamma_E(0).
$$
The minimizing coefficient and forecast are
$$
\kappa:=\frac{\gamma_E(1)}{\gamma_E(0)}
=\frac{q_V\phi_2+b^2q_U\phi_1}{q_V+b^2q_U},\qquad
\widehat E^{\mathrm{lin}}_{t+1\mid t}=\kappa E_t,
$$
with mean squared error
$$
\E[(E_{t+1}-\widehat E^{\mathrm{lin}}_{t+1\mid t})^2]
=\gamma_E(0)-\frac{\gamma_E(1)^2}{\gamma_E(0)}.
$$
This is the best linear forecast based on $E_t$ alone. It is the
one-observation case of Proposition 3.8 in the workbook on prediction.

The coefficient $\kappa$ is a weighted average of $\phi_1$ and $\phi_2$,
so $|\kappa|<1$. The predicted change is $(\kappa-1)E_t$: a positive gap
has a negative predicted change, and a negative gap has a positive one.
Compared with repeating the current gap, the mean squared error is
smaller by $(1-\kappa)^2\gamma_E(0)$.

```{r}
#| label: spread-forecast-risk
#| include: false
gamma_e_shared <- function(h) {
  q_v_shared * phi2_shared^abs(h) +
    b_shared^2 * q_u_shared * phi1_shared^abs(h)
}
kappa_shared <- gamma_e_shared(1) / gamma_e_shared(0)
mse_gap_shared <- gamma_e_shared(0) - gamma_e_shared(1)^2 / gamma_e_shared(0)
mse_repeat_shared <- 2 * (gamma_e_shared(0) - gamma_e_shared(1))
```

For the parameters used in the simulation, $\kappa=$
`r round(kappa_shared, 3)`. The model mean squared error is
`r round(mse_gap_shared, 3)`, compared with `r round(mse_repeat_shared, 3)`
for repeating the current gap.

When $\phi_1\ne\phi_2$, older gaps improve the linear forecast. Indeed,
$$
\begin{aligned}
\Cov(E_{t+1}-\kappa E_t,E_{t-1})
&=\gamma_E(2)-\frac{\gamma_E(1)^2}{\gamma_E(0)}\\
&=\frac{b^2q_Uq_V(\phi_1-\phi_2)^2}{q_V+b^2q_U}>0.
\end{aligned}
$$
The forecast error is still correlated with an earlier observation.
In particular, the spread is not an AR(1) process: its ACVF does not
satisfy $\gamma_E(2)\gamma_E(0)=\gamma_E(1)^2$. If
$\phi_1=\phi_2=\phi$, however,
$E_t=\phi E_{t-1}+\xi_t-b\varepsilon_t$, and the current gap suffices.

For a forecast using several past gaps, apply Proposition 3.8 with the
ACVF in Example 10.9. Joint models for $(\nabla X_t,E_t)^\top$ can also
use changes in $X_t$; the workbooks on vector autoregression and matrix
filters develop such models.

To forecast $X_{t+1}$ and $Y_{t+1}$ separately, we also need a forecast
of one coordinate's change. We predict $\nabla X_{t+1}$ linearly from
the same current gap $E_t$. The model gives
$$
\nabla X_{t+1}=(\phi_1-1)U_t+\eta_{t+1}+\varepsilon_{t+1}.
$$
The new innovations are independent of $E_t$, and
$\Cov(U_t,E_t)=-bq_U$. Hence the best linear forecast of this change
from $E_t$ is
$$
\widehat{\nabla X}^{\mathrm{lin}}_{t+1\mid t}
=\frac{\Cov(\nabla X_{t+1},E_t)}{\Var(E_t)}E_t
=\frac{b(1-\phi_1)q_U}{q_V+b^2q_U}E_t.
$$

Adding the predicted change to $X_t$ gives a forecast of $X_{t+1}$.
We then use $Y_{t+1}=a+bX_{t+1}+E_{t+1}$ to obtain
$$
\begin{aligned}
\widehat X^{\mathrm{lin}}_{t+1\mid t}
&=X_t+\frac{b(1-\phi_1)q_U}{q_V+b^2q_U}E_t,\\
\widehat Y^{\mathrm{lin}}_{t+1\mid t}
&=a+b\widehat X^{\mathrm{lin}}_{t+1\mid t}+\kappa E_t\\
&=Y_t-\frac{(1-\phi_2)q_V}{q_V+b^2q_U}E_t.
\end{aligned}
$$
These minimize mean squared error among forecasts of the forms
$X_t+rE_t$ and $Y_t+sE_t$, respectively: each uses only the current
gap to predict the next change.

For $b>0$, a positive gap predicts an increase in $X$ and a decrease
in $Y$; the signs reverse for a negative gap. The forecasted levels
satisfy $\widehat Y^{\mathrm{lin}}_{t+1\mid t}-a-b\widehat X^{\mathrm{lin}}_{t+1\mid t}
=\kappa E_t$, agreeing with the gap forecast above.

## Estimating and Assessing a Relation

The population spread uses unknown coefficients in an application. Given
observations $(x_t,y_t)$, choose a response and estimate its relation by
ordinary least squares:
$$
(\widehat a,\widehat b)
:=\operatorname*{arg\,min}_{a,b}\sum_{t=1}^T(y_t-a-bx_t)^2,
\qquad \widehat u_t:=y_t-\widehat a-\widehat b x_t.
$$
This fixes the coefficient of $Y$ at one; rescaling all coefficients
would describe the same population direction. In R, the calculation is
`fit <- lm(y ~ x)`, followed by `resid(fit)`.

```{r}
#| label: shared-trend-fit
#| code-fold: true
#| code-summary: "Fit the relation in the simulated record"
fit_shared <- lm(y_shared ~ x_shared)
```

In the model of Example 10.9, OLS consistently estimates $a$ and $b$.
A fixed incorrect slope leaves a multiple of $S_t$ in the residual; the
correct slope cancels it. The estimated slope in the simulated record is
`r round(unname(coef(fit_shared)[2]), 3)`, compared with the true value
`r b_shared`. Consistency does not give conventional regression standard
errors their usual interpretation.

```{r}
#| label: fig-estimated-spread
#| fig-height: 2.3
#| code-fold: true
#| code-summary: "R code to compare the fitted spread with the known spread"
#| fig-cap: "The same simulated record as the common-trend figure. The thick line is the known stationary spread; the thin line uses the estimated intercept and slope. Fitting the relation changes the spread."
wb_par(simulated = TRUE)
plot(e_shared, type = "l", lwd = 2.5, col = "gray55",
     xlab = "Time index", ylab = "Spread (simulated units)",
     ylim = range(e_shared, resid(fit_shared)))
lines(resid(fit_shared), lwd = .9, col = wb_blue)
abline(h = 0, col = "gray70", lty = 3)
legend("top", c("Known spread", "Fitted spread"),
       col = c("gray55", wb_blue), lwd = c(2.5, .9), bty = "n", cex = .85,
       horiz = TRUE, inset = c(0, -.25), xpd = NA)
```

For a coefficient $b_0$ specified before looking at the observations, the
spread $Y_t-b_0X_t$ is an observed scalar series. The level-KPSS procedure
tests whether it is $I(0)$ around an unknown constant mean, subject to the
same dependence and bandwidth conditions as before. Rejection rules
against this particular coefficient choice; it does not rule out every
possible cointegrating relation.

The OLS spread has instead been chosen from the observations. Even under
Example 10.9, its fitted residuals satisfy
$$
\widehat u_t
=e_t-(\widehat a-a)-(\widehat b-b)x_t,\qquad e_t:=y_t-a-bx_t.
$$
Its sample mean is zero because an intercept was fitted. Neither that
constraint nor a small plotted spread establishes a stationary population
relation. Ordinary KPSS critical values do not account for estimating the
slope from the same record. We use OLS to estimate a candidate relation
and examine its later behavior; formal testing of the estimated relation
requires a separate calibration.

## Oil and Gasoline

The `astsa::oil` series records weekly West Texas Intermediate crude-oil
spot prices in dollars per barrel; `astsa::gas` records New York Harbor
regular gasoline spot prices in cents per gallon. These 545 matched EIA
observations run from 2000 to mid-2010. Crude oil is an input to gasoline
production, which motivates examining their relative prices. Write $O_t$
and $G_t$ for the recorded oil and gasoline prices. We specify
the unit-slope log relation first:
$$
R_t:=\log G_t-\log O_t.
$$
Changing the recorded units adds a constant to $R_t$, leaving its
level-stationarity property unchanged.

Use weeks 1–350 for assessment and fitting, leaving weeks 351–545 for a
later check. On the first period, apply level KPSS to $R_t$, using the
Bartlett estimate with bandwidth $H=8$ and the supplied 5% cutoff $0.463$.
The null is that this fixed log ratio is $I(0)$ around a constant mean.
If both log prices are $I(1)$, this null specifies a particular cointegrating
relation. 

```{r}
#| label: oil-gas-calculations
#| code-fold: true
#| code-summary: "Load the prices, test the fixed ratio, and fit the general relation"
# Level KPSS: estimate a constant mean, use divisor n at every lag, and
# include lags 1,...,H with Bartlett weights 1-h/(H+1). Here n corresponds
# to the mathematical sample size T.
kpss_level <- function(x, H = 8L) {
  stopifnot(is.numeric(x), is.null(dim(x)), length(x) >= 4L,
            all(is.finite(x)), length(H) == 1L, is.finite(H),
            H >= 0L, H == floor(H), H < length(x))
  n <- length(x)
  e <- as.numeric(x) - mean(x)
  gamma <- vapply(0:H, function(h) {
    sum(e[seq_len(n - h)] * e[h + seq_len(n - h)]) / n
  }, numeric(1))
  weights <- 1 - (0:H) / (H + 1)
  lrv <- 2 * sum(weights * gamma) - gamma[1L]
  if (!is.finite(lrv) || lrv <= 0) {
    stop("KPSS requires a positive estimated long-run variance.")
  }
  partial_sums <- cumsum(e)
  statistic <- sum(partial_sums^2) / (n^2 * lrv)
  list(statistic = statistic, lrv = lrv, bandwidth = H, n = n,
       critical_5pct = 0.463, reject_5pct = statistic > 0.463,
       residuals = e, partial_sums = partial_sums,
       autocovariances = gamma)
}
data("oil", package = "astsa")
data("gas", package = "astsa")
stopifnot(length(oil) == 545L, length(gas) == length(oil),
          identical(tsp(oil), tsp(gas)), all(oil > 0), all(gas > 0))
log_oil <- log(as.numeric(oil))
log_gas <- log(as.numeric(gas))
fuel_time <- as.numeric(time(oil))
formation_end <- 350L
formation <- seq_len(formation_end)
later <- (formation_end + 1L):length(log_oil)

# The coefficient one is prespecified; KPSS is applied only to this ratio.
fixed_spread <- log_gas - log_oil
fixed_kpss <- kpss_level(fixed_spread[formation], H = 8L)

# Fit on weeks 1-350 and retain the same coefficients for every later week.
fuel_fit <- lm(log_gas ~ log_oil, subset = formation)
fuel_a <- unname(coef(fuel_fit)[1])
fuel_b <- unname(coef(fuel_fit)[2])
fuel_spread <- log_gas - fuel_a - fuel_b * log_oil
spread_summary <- function(z) {
  c(Mean = mean(z), SD = sd(z),
    `Lag-one ACF` = as.numeric(acf(z, lag.max = 1, plot = FALSE)$acf[2]))
}
fuel_windows <- rbind(
  `Weeks 1-350: fitting` = spread_summary(fuel_spread[formation]),
  `Weeks 351-545: later` = spread_summary(fuel_spread[later])
)

```

The KPSS statistic is `r round(fixed_kpss$statistic, 3)`, below $0.463$.
We do not reject stationarity of the prespecified log ratio at 5%.
This is evidence compatible with the proposed relation over these weeks,
subject to the KPSS assumptions; nonrejection does ofcourse not establish the null.

Now fit $\log G_t=a+b\log O_t+U_t$ by OLS on the same 350 weeks. The
estimates are $\widehat a=$ `r round(fuel_a, 3)` and
$\widehat b=$ `r round(fuel_b, 3)`. The estimated slope is close to the
prespecified value one. For every later observation, keep these estimates
fixed and compute $\log G_t-\widehat a-\widehat b\log O_t$.


```{r}
#| label: fig-oil-gas
#| fig-height: 4.4
#| code-fold: true
#| code-summary: "R code for the frozen relation and its later spread"
#| fig-cap: "Real data: weekly WTI crude oil and New York Harbor gasoline spot prices, EIA via astsa. The vertical line ends the 350-week fitting period. The lower panel compares the frozen OLS spread with the prespecified log ratio centered using its first-period mean; the two curves are close because the fitted slope is close to one."
wb_par(rows = 2)
matplot(fuel_time, cbind(log_gas, fuel_a + fuel_b * log_oil),
        type = "l", col = c(wb_blue, wb_copper), lty = 1, lwd = 1.3,
        xlab = "", ylab = "Log price (cents/gallon)",
        main = "Gasoline and the relation fitted to oil")
abline(v = fuel_time[formation_end], lty = 3, col = wb_gray)
legend("topleft", c("Observed log gasoline", "Frozen oil-based relation"),
       col = c(wb_blue, wb_copper), lty = 1, bty = "n", cex = .8)
fixed_centered <- fixed_spread - mean(fixed_spread[formation])
matplot(fuel_time, cbind(fuel_spread, fixed_centered), type = "l",
        col = c(wb_blue, wb_copper), lty = c(1, 2), lwd = c(1.5, 1.0),
        xlab = "Year (dataset index)", ylab = "Log spread",
        main = "The same coefficients in the later period")
abline(v = fuel_time[formation_end], lty = 3, col = wb_gray)
abline(h = 0, lty = 3, col = "gray65")
legend("topleft", c("Frozen OLS spread", "Fixed ratio, centered"),
       col = c(wb_blue, wb_copper), lty = c(1, 2), bty = "n", cex = .8,
       horiz = TRUE)
par(mfrow = c(1, 1))
```

```{r}
#| label: oil-gas-window-table
#| echo: false
knitr::kable(fuel_windows, digits = 3, row.names = TRUE,
             col.names = c("Mean", "SD", "Lag-one ACF"))
```

The first-period fitted spread has mean zero by construction. Its later
mean is `r round(fuel_windows[2, "Mean"], 3)`, and its later SD is
`r round(fuel_windows[2, "SD"] / fuel_windows[1, "SD"], 2)` times its
fitting-period SD. These changes can reflect sampling variation, coefficient
error, or a changing relation. The comparison describes how the frozen
relation behaves beyond its fitting period; it does not separate those
explanations.

The fixed-ratio test, the fitted slope, and the later spread answer
different questions. The first checks one proposed stationary relation;
the second estimates a more general relation; the third examines its
behavior on later observations. Prediction additionally requires a model
for how the gap evolves. A trading strategy would also require positions,
timing, and transaction costs.

## Exercises

::: {.exercise}
**Exercise 10.1 (Independent stationary relations).** Let $Q_0=0$ and
$Q_t=Q_{t-1}+\eta_t$ for $t\ge1$. Suppose $(\eta_t)$, $(U_t)$, and
$(V_t)$ are mutually independent sequences of independent $N(0,1)$
variables. Define
$$
X_t=(Q_t+U_t,\ Q_t+V_t,\ Q_t)^\top.
$$

(a) Find every nonzero $c\in\mathbb R^3$ for which $c^\top X_t$ is
    weakly stationary. Give its variance and ACVF.

(b) Find two linearly independent coefficient vectors with this property
    and express every stationary combination using their two spreads.
    Does $X_{1t}-X_{2t}$ supply a third independent stationary relation?
:::


::: {.exercise}
**Exercise 10.2 (What does a high $R^2$ establish?).** Compare two models
for a pair $(X_t,Y_t)$. In both, $X_0=0$ and
$X_t=X_{t-1}+\eta_t$. In Model A, $Y_0=0$ and
$Y_t=Y_{t-1}+\zeta_t$. In Model B,
$Y_t=1+1.5X_t+U_t$. The sequences $(\eta_t)$, $(\zeta_t)$, and
$(U_t)$ are mutually independent, with independent $N(0,1)$ variables.

(a) Using seed `542102`, simulate 1,000 independent records of length
    $T=200$. Within each record, use the same $X_t$ for both models
    and fit each $y_t$ on an intercept and $x_t$. For the first record,
    report both fitted slopes and $R^2$ values, and plot both residual
    series; overlay $U_t$ on the Model B plot.

(b) For each model, report the median fitted slope, median $R^2$, and
    fraction $\hat p$ of records with $R^2>1/2$. Give the Monte Carlo
    standard error $\sqrt{\hat p(1-\hat p)/1000}$ of each fraction.
    Explain why a high $R^2$ does not establish cointegration. In Model B,
    distinguish the known stationary spread from the fitted residuals.
:::


::: {.exercise}
**Exercise 10.4 (Can a small extra trend destroy cointegration?).** Let
$S_0=R_0=0$ and
$$
S_t=S_{t-1}+\eta_t,\qquad R_t=R_{t-1}+\xi_t,\qquad t\ge1.
$$
Suppose $(\eta_t)$, $(\xi_t)$, and $(U_t)$ are mutually independent
sequences of independent $N(0,1)$ variables, with $U_t$ defined also at
$t=0$. For a fixed $\delta\ge0$, define
$$
X_t=S_t,\qquad Y_t=2S_t+U_t+\delta R_t.
$$

(a) Set $\delta=0$. Find a cointegrating vector for $(X_t,Y_t)^\top$.
    Explain which accumulating component cancels in its spread.

(b) Now suppose $\delta>0$. For an arbitrary fixed slope $b$, calculate
    $\Var(Y_t-bX_t)$. Can any slope make this gap stationary? Are the
    coordinates cointegrated?

(c) Using seed `542104`, simulate one record of length $T=5000$ for
    $\delta=0$ and $\delta=0.02$, using the same realizations of $S_t$,
    $R_t$, and $U_t$ in both cases. Plot the gaps $Y_t-2X_t$ over the
    first 200 observations and over the full record, using the same
    vertical scale. Describe what changes when the record is extended.
    Does a stationary-looking gap establish cointegration? Use (b) to
    explain why a small nonzero $\delta$ can be difficult to detect in
    a short record.
:::

## Sources

[Hyndman and Athanasopoulos (2021, Section 7.3)](https://otexts.com/fpp3/regression-evaluation.html#spurious-regression)
provide the air-passenger/rice example and attribute the observations to
the World Bank. [Granger and Newbold (1974)](https://doi.org/10.1016/0304-4076(74)90034-7)
study spurious time-series regressions.
[Engle and Granger (1987)](https://doi.org/10.2307/1913236) develop
cointegration and error correction. The OLS consistency statement here is
restricted to the explicitly specified common-trend model.

[Kwiatkowski, Phillips, Schmidt, and Shin (1992)](https://doi.org/10.1016/0304-4076(92)90104-Y)
give the level-KPSS reference used for the fixed spread. The folded R code
implements the course's divisor-$T$ sample covariances and Bartlett weights
$1-h/(H+1)$. The oil and gasoline observations are distributed with
[astsa](https://cran.r-project.org/package=astsa), with original spot-price
data from the [U.S. Energy Information Administration](https://www.eia.gov/dnav/pet/pet_pri_spt_s1_w.htm).
Shumway and Stoffer use the same pair in their time-series regression examples.