Workbook 5 — Learning Covariances from One Path

STA 542 · Introduction to Time Series Analysis

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 code 
# 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")
Figure 1: 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.

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\)

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)
Figure 2: 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.

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

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

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

  3. 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?

  4. Which conclusions in parts (a)–(c) require Gaussianity? What additional conclusion becomes available if the innovations are Gaussian?

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.

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

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

  3. 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)\).

  4. Carry out the requested simulation.

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

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