---
title: "Workbook 5 — Learning Covariances from One Path"
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.4
    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
set.seed(542)
```

Workbook 4 explored a Bayesian approach to learning about a process law by
placing a prior over several candidate process laws. We now return to the
setting where we consider a collection of fixed process laws corresponding
to different parameter values. In particular, we will focus on weakly
stationary laws where the mean $\mu$ and autocovariance function $\gamma$ are
fixed features of that law, even when we do not know their values.

Workbook 2 showed when the sample mean can converge to $\mu$. Consistency says
that the error eventually disappears. It does not tell us how large the error
is at the observed sample size $T$, or whether one path can reveal the
covariances needed to calculate that uncertainty.

## A Mathematical Preliminary: Consistency

Let $\theta\in\R$ be a fixed feature of the process law, and let
$\widehat\theta_T$ be a sequence of random variables, typically a statistic
computed from $X_1,\ldots,X_T$. The sequence is *consistent* for $\theta$ if
it converges in probability to $\theta$: for every $\varepsilon>0$,
$$
\P\bigl(|\widehat\theta_T-\theta|>\varepsilon\bigr)\longrightarrow0.
$$
Equivalently, $\widehat\theta_T\xrightarrow{P}\theta$.

Workbook 2 used a stronger mode. We say that $\widehat\theta_T$ *converges in
mean square* to $\theta$ if
$$
\E[(\widehat\theta_T-\theta)^2]\longrightarrow0.
$$
Chebyshev's inequality turns a vanishing second moment into the probability
bound above, so mean-square convergence implies consistency. The converse is
false: a consistent sequence need not have vanishing mean squared error.

## Consistency Leaves the Error Uncalibrated

Consider a weakly stationary process $(X_t)$ with mean $\mu$ and ACVF $\gamma$ under the law $\P \equiv \P_{\mu, \gamma}$, and write
$$
\bar X_T:=\frac1T\sum_{t=1}^T X_t.
$$
Mean-square consistency is the statement
$$
\E[(\bar X_T-\mu)^2]\longrightarrow0.
$$

Equivalently, $\Var(\bar X_T) \to 0$. If we wish make uncertainty statements about how close $\bar X_T$ is to $\mu$, we need to analyze the variance of $\bar X_T$ more carefully.

**Proposition 5.1 (Exact variance of the sample mean).** If $(X_t)$ is weakly
stationary with ACVF $\gamma$, then
$$
\begin{aligned}
v_T(\gamma)
:=\Var(\bar X_T)
&=\frac1{T^2}\sum_{s=1}^T\sum_{t=1}^T\gamma(t-s)\\
&=\frac1T\left\{
\gamma(0)+2\sum_{h=1}^{T-1}
\left(1-\frac hT\right)\gamma(h)
\right\}.
\end{aligned}
$$

*Proof.* This is the variance expansion already used in the proof of
Proposition 2.18; we record the grouping. The variance of a sum is the
double sum of its covariances, and weak stationarity writes each of them
as $\gamma(t-s)$:
$$
\Var(\bar X_T)
=\frac1{T^2}\sum_{s=1}^T\sum_{t=1}^T\gamma(t-s).
$$
For each lag $h\in\{-(T-1),\ldots,T-1\}$ there are $T-|h|$ ordered pairs
$(s,t)$ with $t-s=h$. Since $\gamma(-h)=\gamma(h)$,
$$
\sum_{s=1}^T\sum_{t=1}^T\gamma(t-s)
=T\gamma(0)+2\sum_{h=1}^{T-1}(T-h)\gamma(h).
$$
Dividing by $T^2$ produces the weights $1-h/T$. $\square$

Weak stationarity turns $T^2$ pairwise covariances into one weighted sum over
lags. Proposition 5.1 is an exact finite-sample identity. It requires neither
Gaussianity nor any law of large numbers.

## If the ACVF Were Known

If $\gamma$ were a known function, Proposition 5.1 would make $v_T(\gamma)$ known.
Chebyshev's inequality would then provide a confidence interval without any
distributional assumption.

**Proposition 5.2 (A known-ACVF Chebyshev interval).** Suppose the conditions of
Proposition 5.1 hold and $\gamma$ is known. For $0<\alpha<1$, define
$$
C_T^{\mathrm{known}}
:=\left[
\bar X_T-\sqrt{\frac{v_T(\gamma)}{\alpha}},
\ \bar X_T+\sqrt{\frac{v_T(\gamma)}{\alpha}}
\right].
$$
Then
$$
\P_{\mu,\gamma}\{\mu\in C_T^{\mathrm{known}}\} \equiv \P\{\mu\in C_T^{\mathrm{known}}\}\ge1-\alpha.
$$

*Proof.* If $v_T(\gamma)>0$, Chebyshev's inequality gives
$$
\P\left\{
|\bar X_T-\mu|>
\sqrt{\frac{v_T(\gamma)}{\alpha}}
\right\}
\le\alpha.
$$
The complement is the displayed coverage event. If $v_T(\gamma)=0$, then
$\bar X_T=\mu$ with probability one. $\square$

The process mean $\mu$ is fixed and the interval varies across possible paths: The probability in Proposition 5.2 is coverage for $\mu$ under one fixed
process law with $\mu$ as true mean. 

### Forecasting with an Estimated Drift

The error for the mean estimate can be used to improve the coverage of a forecast. The estimation error for the mean then
contributes to the error of predicting a future observation.

**Example 5.3 (A conservative prediction interval).** 
We consider a random walk with a fixed but unknown drift $\delta$:
$$
X_t=X_{t-1}+\delta+Z_t,
\qquad t\ge1,
\qquad
(Z_t)\sim\mathrm{WN}(0,\sigma^2).
$$
The initial value $X_0=x_0$ is fixed, the drift $\delta\in\R$ is fixed but
unknown, and $0<\sigma^2<\infty$ is known. Write $\P_\delta$ for any fixed
process law satisfying these conditions. The white-noise assumption specifies
zero means, common variance, and zero covariances at distinct times;
Gaussianity and independence are not assumed.

Observe $X_0,\ldots,X_T$, where $T\ge1$, and define the increments and their
sample mean by
$$
D_t:=X_t-X_{t-1},
\qquad
\widehat\delta_T:=\frac1T\sum_{t=1}^T D_t
=\delta+\frac1T\sum_{t=1}^T Z_t.
$$
The levels are nonstationary, but the increments are weakly stationary with
mean $\delta$ and ACVF zero away from lag zero. Proposition 5.1 applied to
$(D_t)$ gives
$\Var_\delta(\widehat\delta_T)=\sigma^2/T$.

For an integer horizon $h\ge1$, use the fitted forecast
$$
\widehat X_{T+h\mid T}:=X_T+h\widehat\delta_T.
$$
We evaluate this rule under the stated moment assumptions; those assumptions
alone do not identify it as a conditional-mean forecast.

Its forecast error is
$$
\begin{aligned}
e_{T,h}
&:=X_{T+h}-\widehat X_{T+h\mid T}\\
&=\underbrace{\sum_{j=1}^h Z_{T+j}}_{\text{future noise}}
-\underbrace{\frac hT\sum_{t=1}^T Z_t}_{\text{drift estimation error}}.
\end{aligned}
$$
The two sums involve disjoint times, so their covariance is zero. Hence
$$
\E_\delta[e_{T,h}]=0,
\qquad
\Var_\delta(e_{T,h})
=h\sigma^2+\frac{h^2\sigma^2}{T}.
$$

For $0<\alpha<1$, define the prediction interval
$$
\begin{aligned}
c_{T,h}&:=\sqrt{\frac{h\sigma^2+h^2\sigma^2/T}{\alpha}},\\
I_{T,h}&:=[\widehat X_{T+h\mid T}-c_{T,h},
\ \widehat X_{T+h\mid T}+c_{T,h}].
\end{aligned}
$$
Chebyshev's inequality applied to $e_{T,h}$ gives
$$
\P_\delta\{X_{T+h}\notin I_{T,h}\}
=\P_\delta\{|e_{T,h}|>c_{T,h}\}
\le\alpha,
$$
so $\P_\delta\{X_{T+h}\in I_{T,h}\}\ge1-\alpha$ for every $T,h\ge1$
and every fixed $\delta$ under the stated assumptions.

The coverage probability repeats both the observed history and its future
while holding $\delta$ fixed. Drift uncertainty enters through the variation
of $\widehat\delta_T$ across those histories. The bound is marginal at one
chosen horizon; it does not give conditional coverage for every observed
history or simultaneous coverage for a future path.

Treating the fitted drift as known would omit $h^2\sigma^2/T$ from the
forecast-error variance. The ratio of this estimation contribution to the
future-noise contribution is $h/T$, so the two contributions are equal when
$h=T$. The corresponding Chebyshev half-width is larger by a factor
$\sqrt{1+h/T}$ than the half-width computed from future noise alone.

This interval is computable because $\sigma^2$ is known. A known upper bound
on $\sigma^2$ can replace it and preserve the coverage guarantee. Replacing
it by a sample variance requires a separate coverage argument. $\Diamond$

### Dependence Changes Mean Uncertainty

We now have the tools to describe how dependence between observations can change the uncertainty about the mean of a weakly stationary process. The following example shows how dependence can increase the uncertainty about the mean.

**Example 5.4 (Equal marginals, different uncertainty).** Compare two laws
having the same fixed mean $\mu$. Under the first law,
$$
X_t\overset{\text{iid}}{\sim}N(\mu,1).
$$
Under the second, the selected stationary causal solution satisfies
$$
X_t-\mu=0.9(X_{t-1}-\mu)+Z_t,
\qquad
Z_t\overset{\text{iid}}{\sim}N(0,1-0.9^2).
$$
Every $X_t$ has the $N(\mu,1)$ distribution under either law. The first ACVF
vanishes away from zero; the second is $\gamma(h)=0.9^{|h|}$.

At $T=200$, Proposition 5.1 gives
$$
200v_{200}=1
\quad\text{under independence},
\qquad
200v_{200}
=1+2\sum_{h=1}^{199}\left(1-\frac h{200}\right)0.9^h
\approx18.1
$$
under the AR(1) law. The 95% Chebyshev half-widths are $0.316$ and $1.345$.
Equal marginal distributions do not give equal information about their common
mean. $\Diamond$

```{r}
#| label: fig-wb5-mean-distributions
#| code-fold: true
#| code-summary: "R code"
#| fig-cap: "Distributions of the scaled sample mean under i.i.d. N(0,1) data and a stationary Gaussian AR(1) with the same marginal law and phi = 0.9. Thin curves are Monte Carlo densities from 5,000 paths of length 200; thick curves are the exact Gaussian densities. Simulated data."
# Simulate one mean-zero AR(1) path of length n and return its sample mean.
simulate_ar_mean <- function(n, phi) {
  x <- numeric(n)
  # Start in the stationary N(0,1) distribution.
  x[1] <- rnorm(1)
  # This innovation variance keeps the marginal variance equal to one.
  for (t in 2:n) {
    x[t] <- phi * x[t - 1] + rnorm(1, sd = sqrt(1 - phi^2))
  }
  mean(x)
}

# In R, n denotes the sample size T; the seed makes the simulation reproducible.
set.seed(542)
n_mean <- 200
phi_mean <- 0.9
# Each replication generates a new path and records one sample mean.
iid_means <- replicate(5000, mean(rnorm(n_mean)))
ar_means <- replicate(5000, simulate_ar_mean(n_mean, phi_mean))
# Put both mean errors on the same sqrt(T) scale; their population means are zero.
scaled_iid_means <- sqrt(n_mean) * iid_means
scaled_ar_means <- sqrt(n_mean) * ar_means
# Compute the exact variance of sqrt(T) times the AR(1) sample mean.
ar_variance_factor <- 1 + 2 * sum(
  (1 - (1:(n_mean - 1)) / n_mean) * phi_mean^(1:(n_mean - 1))
)

# Thick curves show exact Gaussian densities; thin curves estimate them from the paths.
density_grid <- seq(-14, 14, length.out = 500)
par(bg = "gray92")
plot(density_grid, dnorm(density_grid), type = "l", col = 4, lwd = 3,
     xlab = expression(sqrt(T) * (bar(X)[T] - mu)), ylab = "Density",
     main = "Dependence changes the error scale", ylim = c(0, 0.45))
lines(density(scaled_iid_means), col = 4, lwd = 1)
lines(density_grid, dnorm(density_grid, sd = sqrt(ar_variance_factor)),
      col = 2, lwd = 3, lty = 2)
lines(density(scaled_ar_means), col = 2, lwd = 1, lty = 2)
legend("topright", c("i.i.d.", "AR(1), phi = 0.9"),
       col = c(4, 2), lty = c(1, 2), lwd = 2, bty = "n")
```

Gaussianity makes the two exact density curves in the figure available. It
was not used in Proposition 5.1 or Proposition 5.2.

## The ACVF Is Usually Unknown

The interval in Proposition 5.2 is computable only because the entire ACVF was
supplied. In data analysis, $\gamma$ is usually another unknown feature of the
process law. Hence, we have to estimate it from the data. A natural estimator is the sample ACVF.

**Definition 5.5 (Sample ACVF).** For $0\le h<T$, define
$$
\hat\gamma_T(h)
:=\frac1T\sum_{t=1}^{T-h}
(X_t-\bar X_T)(X_{t+h}-\bar X_T),
\qquad
\hat\gamma_T(-h):=\hat\gamma_T(h).
$$
The divisor is $T$, rather than $T-h$. This convention makes the complete
sample ACVF nonnegative definite.

The sample ACVF is itself a collection of time averages. Learning the mean
does not automatically make these new averages consistent.

**Example 5.6 (A persistent random scale).** Let $\Lambda$ be independent of an
i.i.d. sequence $(R_t)_{t\in\Z}$, with
$$
\P(\Lambda=1)=\P(\Lambda=2)=\frac12,
\qquad
\P(R_t=-1)=\P(R_t=1)=\frac12,
$$
and set $X_t:=\Lambda R_t$. The random scale $\Lambda$ is part of this one
process law; it does not label different candidate laws.

A common time shift changes only the indices of the i.i.d. sequence, so
$(X_t)$ is strictly stationary. Its mean is zero and
$$
\gamma(0)=\E[\Lambda^2]=\frac52,
\qquad
\gamma(h)=\E[\Lambda^2R_tR_{t+h}]=0
\quad(h\ne0).
$$
The sample mean therefore converges to zero in mean square.

Since $X_t^2=\Lambda^2$ for every $t$,
$$
\hat\gamma_T(0)
=\frac1T\sum_{t=1}^T(X_t-\bar X_T)^2
=\Lambda^2-\bar X_T^2
\xrightarrow{P}\Lambda^2.
$$
The population value is $5/2$, while the random limit is either $1$ or $4$.
One path learns its realized scale rather than the variance obtained by
averaging over both scales. $\Diamond$

```{r}
#| label: fig-wb5-random-scale
#| fig-height: 4.2
#| fig-cap: "Running means and second moments for two paths from each scale in Example 5.6. Every mean approaches zero; each second moment remains at its realized scale squared rather than the population value 5/2 (dashed). Simulated data."
set.seed(542)
n_scale <- 800
scale_values <- c(1, 1, 2, 2)
scale_paths <- sapply(scale_values, function(scale) {
  scale * sample(c(-1, 1), n_scale, replace = TRUE)
})
running_means <- apply(scale_paths, 2, cumsum) / seq_len(n_scale)
running_seconds <- apply(scale_paths^2, 2, cumsum) / seq_len(n_scale)
scale_colors <- c("#0072B2", "#56B4E9", "#D55E00", "#E69F00")

par(mfrow = c(1, 2), mar = c(4, 4, 2.5, 1), bg = "gray92")
matplot(seq_len(n_scale), running_means, type = "l", lty = 1,
        col = scale_colors, xlab = "T", ylab = "Running mean",
        main = "The mean")
abline(h = 0, lty = 2)
matplot(seq_len(n_scale), running_seconds, type = "l", lty = 1,
        col = scale_colors, xlab = "T", ylab = "Running second moment",
        main = "The second moment")
abline(h = 2.5, lty = 2)
```

Just like with the mean, we need further assumptions beyond weak stationarity to ensure that the sample ACVF converges to the true ACVF. 

## Beyond mean ergodicity

In Example 5.6, averaging $X_t$ recovers its population mean, while averaging
$X_t^2$ does not recover its population second moment. To estimate a lagged
covariance, we also need to average products $X_tX_{t+h}$ at a fixed lag $h$.
Each calculation uses a different function of the observations, so each
requires a convergence statement for that function's average.

Let $g$ specify what we compute from a window of consecutive observations
before averaging across time. A window may contain one observation, as when
we square it, or several observations, as when we multiply its first and last
entries. Its length stays fixed while the observed record grows. The following
definition asks whether the average of these computed values converges in
probability to their common expectation.

**Definition 5.7 ($g$-ergodicity).** Let $(X_t)$ be strictly stationary. Fix a
window length $r+1$, where $r$ is a nonnegative integer, and a function
$g:\R^{r+1}\to\R$ for which
$$
\E\!\left[|g(X_0,X_1,\ldots,X_r)|\right]<\infty.
$$
The process is *$g$-ergodic* if
$$
\frac1{T-r}\sum_{t=1}^{T-r}g(X_t,X_{t+1},\ldots,X_{t+r})
\xrightarrow{P}
\E[g(X_0,X_1,\ldots,X_r)].
$$

Strict stationarity makes the target unambiguous: every shifted window has the
same distribution and hence the same expectation. The definition asks only
for the one time average determined by $g$. It makes no claim about every
possible feature of the path.

For covariance estimation, the identity
$$
\gamma(h)=\E[X_0X_h]-\mu^2
$$
identifies the required averages: the sample mean must recover $\mu$, and
each raw product average must recover $\E[X_0X_h]$. Weak stationarity makes
these particular population moments independent of the time origin, so we
can state the required convergence without assuming strict stationarity.

**Definition 5.8 (ACVF-learnable through lag $m$).** A weakly stationary
process is *ACVF-learnable through lag $m$* if
$$
\bar X_T\xrightarrow{P}\mu
$$
and, for every $h=0,1,\ldots,m$,
$$
\frac1T\sum_{t=1}^{T-h}X_tX_{t+h}
\xrightarrow{P}\E[X_0X_h].
$$

Under strict stationarity, these are equivalent to the $g$-ergodicity
requirements for $g(x)=x$ and the finitely many functions
$g_h(x_0,\ldots,x_h)=x_0x_h$, since $(T-h)/T\to1$ for every fixed $h$.
Definition 5.8 itself needs only the weakly stationary population moments
displayed in it.

**Proposition 5.9 (Fixed-lag sample-ACVF consistency).** If $(X_t)$ is
ACVF-learnable through lag $m$, then, for each fixed
$h=0,1,\ldots,m$,
$$
\hat\gamma_T(h)\xrightarrow{P}\gamma(h).
$$

*Proof.* Expanding Definition 5.5 gives
$$
\begin{aligned}
\hat\gamma_T(h)
={}&\frac1T\sum_{t=1}^{T-h}X_tX_{t+h}
-\bar X_T\frac1T\sum_{t=1}^{T-h}X_t\\
&-\bar X_T\frac1T\sum_{t=1}^{T-h}X_{t+h}
+\left(1-\frac hT\right)\bar X_T^2.
\end{aligned}
$$

Each partial average differs from $\bar X_T$ by $h$ omitted observations
divided by $T$. For the omitted terms at the right endpoint, weak stationarity
and Cauchy--Schwarz give
$$
\E\!\left[\left|\frac1T\sum_{t=T-h+1}^{T}X_t\right|\right]
\le\frac hT\sqrt{\mu^2+\gamma(0)}
\longrightarrow0.
$$
The same bound applies to the $h$ omitted terms at the left endpoint. Markov's
inequality therefore makes both contributions converge to zero in probability.
Both partial averages converge to $\mu$ because $h$ is fixed.

The four terms in the expansion consequently converge to
$\E[X_0X_h]-\mu^2-\mu^2+\mu^2=\gamma(h)$. $\square$

We will use one standard fact without proof: a well-defined causal function of
an i.i.d. sequence is $g$-ergodic for every integrable fixed-window function
used here. This covers finite moving averages and the selected causal ARMA and
GARCH solutions encountered in the course. We use the fact only to obtain the
required averages.

## Fixed Lags Leave a Growing-Sum Problem

Proposition 5.9 controls any lag fixed before $T$ grows. Proposition 5.1 uses
lags $0,1,\ldots,T-1$. Pointwise consistency at each fixed lag does not control
a sum whose number of estimated terms grows with the sample size.

We can learn any fixed collection of autocovariances, but the variance of the
sample mean involves a collection that grows with $T$. Fixed-lag consistency
supplies neither a common rate nor control of that growing sum.
Workbook 6, on estimating uncertainty, identifies a stable population scale
and develops covariance estimators that yield asymptotically conservative
intervals. Workbook 7, on central limit theorems, asks what additional
conditions justify a Gaussian approximation and narrower intervals.

## Exercises

::: {.exercise}
**Exercise 5.1 (Known dependence and uncertainty about a mean).** Let
$$
X_t-\mu=\phi(X_{t-1}-\mu)+Z_t,
\qquad
\E[Z_t]=0,
\qquad
\Var(Z_t)=1-\phi^2,
\qquad |\phi|<1,
$$
where $(Z_t)$ is i.i.d. Use the selected stationary causal solution. Treat
$\phi$ as known and $\mu$ as fixed but unknown.

*Method and report.* Derive all population quantities by hand. Use R or a
calculator for the numerical sums. Report the three 95% Chebyshev half-widths
to three decimal places and interpret their ordering.

(a) Show that $\gamma(h)=\phi^{|h|}$.

(b) Derive the exact finite-sample variance of $\bar X_T$.

(c) At $T=50$, calculate the 95% Chebyshev half-width for each
$\phi\in\{-0.6,0,0.6\}$. All three laws have the same marginal variance.
Explain the ordering of the half-widths using the signs of the lagged
covariances. Does dependence always make estimation of a mean less precise?

(d) Which conclusions in parts (a)--(c) require Gaussianity? What additional
conclusion becomes available if the innovations are Gaussian?
:::

::: {.exercise}
**Exercise 5.2 (What one MA(1) path can learn).** Let
$$
X_t=Z_t+\theta Z_{t-1},
$$
where $(Z_t)_{t\in\Z}$ are i.i.d. with mean zero and variance
$\sigma^2<\infty$. Treat $\theta$ and $\sigma^2$ as fixed features of the law.

*Method and report.* Derive the ACVF and identify each function being averaged.
For the computation, use independent $N(0,1)$ innovations, $\theta=0.6$, and
one path of length 2,000. Plot running estimates of $\gamma(0)$ and
$\gamma(1)$, and report their final values.

(a) Derive $\gamma(0)$, $\gamma(1)$, and $\gamma(h)$ for $|h|>1$.

(b) Identify the fixed-window functions $g$ needed to learn the mean and the
raw products through lag two.

(c) Use the causal-i.i.d. fact and Proposition 5.9 to justify consistency of
$\hat\gamma_T(0)$, $\hat\gamma_T(1)$, and $\hat\gamma_T(2)$.

(d) Carry out the requested simulation.

(e) Explain why the argument does not justify summing all sample
autocovariances through lag $T-1$.

```{r}
#| eval: false
# Student starter code
set.seed(542)
n <- 2000
theta <- 0.6
z <- rnorm(n + 1)
x <- z[2:(n + 1)] + theta * z[1:n]

# Complete the functions and plot the running estimates.
sample_acvf <- function(x, h) {
  # ...
}
sizes <- seq(50, n, by = 25)
```
:::

## Sources

- van der Vaart, A. W. (2010). *Time Series*. Theorem 8.10 supplies the
  causal-function route to fixed-window ergodic averages.
- Brockwell, P. J., and Davis, R. A. (2016). *Introduction to Time Series and
  Forecasting*, 3rd ed. Sections 2.3 and 7.3 discuss sample autocovariances and
  large-sample inference for a process mean.