---
title: "Workbook 3 — Forecasting When the Process Law Is Known"
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\N{\mathbb{N}}
\def\cH{\mathcal{H}}
\def\cP{\mathcal{P}}
\def\Var{\operatorname{Var}}
\def\Cov{\operatorname{Cov}}
\def\Corr{\operatorname{Corr}}$$
:::

```{=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 2 supplied process laws and generating equations. It hinted at conditions under which
patterns recur over time, in which case there is some hope of estimation. The second thing for which recurrence can be useful is forecasting. We postpone the question of fitting a model for now, and assume (for now) that
we know the law that governs the process exactly. The
problem is to use an observed window to describe unobserved (future or past) values. 

## A Mathematical Preliminary: Chebyshev's Inequality

Before constructing forecasts, we recall a simple way to turn two moments into
a probability statement. Let $U$ be a real-valued random variable with mean
$\mu$. *Markov's inequality* states that if $V\ge0$ and $c>0$,
then
$$
\P(V\ge c)\le\frac{\E[V]}{c}.
$$
If in addition $U$ has finite variance $\sigma^2$, we obtain *Chebyshev's inequality*, which states that
states that
$$
\P(|U-\mu|\ge c)\le\frac{\sigma^2}{c^2}.
$$
(Apply the bound to
$V=(U-\mu)^2$ with threshold $c^2$.)

Taking $c=\sigma/\sqrt{\alpha}$ for $0<\alpha<1$,
$$
\P\!\left(
  \mu-\frac{\sigma}{\sqrt{\alpha}}
  \leq U\leq
  \mu+\frac{\sigma}{\sqrt{\alpha}}
\right)
\geq1-\alpha.
$$
For 95% coverage, $1/\sqrt{0.05}=\sqrt{20}\approx4.47$. Thus a mean and
variance alone guarantee a 95% interval extending 4.47 standard deviations on
either side of the mean. The interval is often conservative because the two
moments say nothing about the shape of the distribution: If
$U\sim N(\mu,\sigma^2)$, the same coverage uses $1.96\sigma$ in place of
$4.47\sigma$.

The same argument given $W$ produces a conditional bound. If
$\E[U\mid W]=\mu(W)$ and $\Var(U\mid W)=\sigma^2(W)$, then for every $c>0$,
$$
\P\bigl(|U-\mu(W)|\ge c\mid W\bigr)\le\frac{\sigma^2(W)}{c^2}.
$$


## The Target and the Information

What is a good predictor? We want a function of the observed window to lie
'close' to the unobserved value of the target. 

What counts as close depends on the task. A *loss function* $L(a,y)$ assigns
a numerical cost to predicting $a$ when the target equals $y$. 
The default in this course is *squared-error loss*,
$$
L(a,y)=(a-y)^2,
$$
the squared Euclidean distance between prediction and target. It penalizes
large misses more than small ones and treats overprediction and
underprediction symmetrically.

Other losses match other costs. Absolute loss $|a-y|$ downweights isolated
large misses relative to squared error. Sometimes it makes have asymetric losses: Undershooting a capital reserve
typically costs more than overshooting it, so a trading-desk would be wise to pick a loss function that is
asymmetric in the sign of the error.





**Definition 3.1 (Population prediction problem).** Consider an $\R^d$-valued time series process $(X_t)$. Let
$$
W:=(X_{t_1},\ldots,X_{t_m})^\top
$$
be the observed window, and let $Y:=X_s$ be a target not in
that window. A *predictor* is a function $g(W)$ taking values in $\R^d$, and a *loss function* is a function $L : \R^d \times \R^d \to \R$ that measures the loss of a prediction through its *risk*:
$$
\E L(Y, g(W)).
$$
The risk is the expected loss of the prediction. Why do we consider the expected loss? The target has not been
realized/observed when the predictor is formed, so closeness is not a distance
between two already-recorded numbers. Another benefit is that bounding expected loss gives us a probability statement about the prediction for free (recall Markov's inequality above!)

When $s>t_m$, we call $g(W)$ a forecast. When $s<t_1$, it is often called a
*backcast*; a missing target between observed times gives an interpolation or
smoothing problem. 

The word *known* includes the noise distribution and its parameters. Replacing
them by estimates produces a different procedure whose uncertainty includes
the randomness of the fit.

## The Conditional Mean


In this course, the default forecasting risk is *risk under squared-error loss*, which is
$$
\E[(Y-g(W))^2],
$$-
where the expectation is computed under the process law. It turns out that the conditional mean is the optimal predictor under squared-error loss.

**Proposition 3.2 (Conditional mean under squared-error loss).** Suppose
$\E[Y^2]<\infty$. Among all predictors $g(W)$ with finite second moment,
$$
m(W):=\E[Y\mid W]
$$
minimizes $\E[(Y-g(W))^2]$.

*Proof.* Conditional on $W$,
$$
\E[(Y-g(W))^2\mid W]
=\Var(Y\mid W)+(m(W)-g(W))^2.
$$
The first term does not depend on $g$. The second is minimized by
$g(W)=m(W)$. Taking expectations gives the claim. $\square$

Note that to compute the conditional mean, we need to know the process law.

The proposition compares point predictors under squared loss. Different losses select different optimal predictors. Absolute loss
instead selects the conditional median, and other losses can select other
features of the conditional distribution.

We consider two simple examples to illustrate how given a known process law, we can compute the conditional mean.

**Example 3.3 (Same marginals, different forecasts).** In Example 1.2, we
saw two Gaussian processes with the same one-dimensional
marginal $N(0,1)$, yet opposite dynamics: Process A is i.i.d.,
$X_t \overset{\text{iid}}{\sim} N(0,1)$. Process B draws a single
$X_0\sim N(0,1)$ once and repeats it: $X_t:=X_0$ for every $t$. Under A,
$\Cov(X_s,X_t)=0$ for $s\neq t$; under B, $\Cov(X_s,X_t)=1$ for every
pair.

Take the window $W=(X_1,\ldots,X_T)^\top$ and the target $Y=X_{T+1}$.
Under Process A, $X_{T+1}$ is independent of $W$, so
$$
\E[X_{T+1}\mid X_1,\ldots,X_T]=0,
\qquad
\Var(X_{T+1}\mid X_1,\ldots,X_T)=1.
$$
Under Process B, $X_{T+1}=X_T$ almost surely, so
$$
\E[X_{T+1}\mid X_1,\ldots,X_T]=X_T,
\qquad
\Var(X_{T+1}\mid X_1,\ldots,X_T)=0.
$$
The one-dimensional marginals agree, while the optimal forecasts and
their errors differ. ◊


**Example 3.4 (A precipitation indicator).** Let $X_t=1$ if a rain gauge
records precipitation on day $t$, and $X_t=0$ otherwise. Suppose $(X_t)$ is
a time-homogeneous Markov chain on $\{0,1\}$ with
$$
\P(X_{t+1}=1\mid X_t=1)=0.7,
\qquad
\P(X_{t+1}=1\mid X_t=0)=0.2.
$$
The one-step conditional mean is the conditional probability of
precipitation:
$$
\E[X_{T+1}\mid X_T]
=0.7X_T+0.2(1-X_T)
=0.2+0.5X_T.
$$
A wet day is followed by the forecast $0.7$; a dry day, by $0.2$. The
conditional forecast-error variance is $p(1-p)$ at that same $p$, hence
$0.21$ after a wet day and $0.16$ after a dry day.

Because the chain is Markov of order one,
$$
\E[X_{T+1}\mid X_1,\ldots,X_T]=\E[X_{T+1}\mid X_T].
$$
The earlier observations in the window do not change the squared-error
forecast. ◊

A useful result for Gaussian processes is the following Gaussian conditioning formula.

**Proposition 3.5 (Gaussian conditioning).** Suppose the target vector $Y$ and
observed vector $W$ have the known joint Gaussian law
$$
\begin{pmatrix}Y\\W\end{pmatrix}
\sim N\!\left[
\begin{pmatrix}\mu_Y\\\mu_W\end{pmatrix},
\begin{pmatrix}
\Sigma_{YY}&\Sigma_{YW}\\
\Sigma_{WY}&\Sigma_{WW}
\end{pmatrix}
\right],
$$
where $\Sigma_{WW}$ is nonsingular. Then
$$
Y\mid W=w
\sim N\!\left(
\mu_Y+\Sigma_{YW}\Sigma_{WW}^{-1}(w-\mu_W),
\ \Sigma_{YY}-\Sigma_{YW}\Sigma_{WW}^{-1}\Sigma_{WY}
\right).
$$

The proof is straightforward but is best followed along with pen and paper (therefor hidden in a foldable box); go through it upon second reading -- it builds character.

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

Write $A:=\Sigma_{YW}\Sigma_{WW}^{-1}$, which is well defined because
$\Sigma_{WW}$ is nonsingular, and subtract from the target the part of it
carried by the window:
$$
V:=(Y-\mu_Y)-A(W-\mu_W).
$$
The pair $(V,W)$ is an affine function of $(Y,W)$, so it is again jointly
Gaussian, and $\E[V]=0$. Its covariance with the window is
$$
\Cov(V,W)
=\Sigma_{YW}-A\Sigma_{WW}
=\Sigma_{YW}-\Sigma_{YW}\Sigma_{WW}^{-1}\Sigma_{WW}
=0.
$$
For jointly Gaussian vectors, zero covariance gives independence, so $V$ is
independent of $W$. Joint Gaussianity does real work here: uncorrelatedness
alone would not supply it. Using $\Sigma_{WY}=\Sigma_{YW}^\top$ and the
symmetry of $\Sigma_{WW}$,
\begin{align*}
\Var(V)
&=\Sigma_{YY}-A\Sigma_{WY}-\Sigma_{YW}A^\top+A\Sigma_{WW}A^\top\\
&=\Sigma_{YY}-\Sigma_{YW}\Sigma_{WW}^{-1}\Sigma_{WY},
\end{align*}
since the last three terms each equal
$\Sigma_{YW}\Sigma_{WW}^{-1}\Sigma_{WY}$ up to sign.

Now $Y=\mu_Y+A(W-\mu_W)+V$. Conditioning on $W=w$ fixes the second term at
the constant $A(w-\mu_W)$ and leaves the law of $V$ unchanged, so $Y\mid W=w$
is Gaussian with mean $\mu_Y+A(w-\mu_W)$ and covariance $\Var(V)$. $\square$
:::


Example 3.6 is the case $Y=X_s$, with $\Sigma_{WW}=\Sigma_W$ assembled from
the covariance kernel.


**Example 3.6 (A Gaussian process forecast).** Let
$$
(X_t) \sim\operatorname{GP}(\mu,\Sigma)
$$
be a Gaussian process as in Example 2.2, and take the window $W$ and
target $Y=X_s$ from Definition 3.1. Write
$$
\mu_W:=(\mu(t_1),\ldots,\mu(t_m))^\top,
\qquad
\Sigma_W:=\bigl(\Sigma(t_i,t_j)\bigr)_{i,j=1}^m,
\qquad
\sigma_{sW}:=\bigl(\Sigma(s,t_1),\ldots,\Sigma(s,t_m)\bigr)^\top.
$$
If $\Sigma_W$ is nonsingular, Gaussian conditioning gives
$$
X_s\mid W
\sim N\!\left(
\mu(s)+\sigma_{sW}^\top\Sigma_W^{-1}(W-\mu_W),
\ \Sigma(s,s)-\sigma_{sW}^\top\Sigma_W^{-1}\sigma_{sW}
\right).
$$
The conditional mean is an affine linear function of $W$:
$$
\E[X_s\mid W]= \mu(s)+\sigma_{sW}^\top\Sigma_W^{-1}(W-\mu_W).
$$
 The formula uses only
the values of $\mu$ and $\Sigma$ at the times $\{s,t_1,\ldots,t_m\}$;
those times need not be ordered, and $(\mu,\Sigma)$ need not be stationary. ◊

Note that in the above example, $X_s$ could lie in the future, the past, or even somewhere in the middle of the observed window. If it lies in the past, we speak of a *backcast*, in the middle, we speak of an *interpolation*, and the word *forecast* typically refers to the case where $X_s$ lies in the future; although literature often uses it for the `past and middle' cases too.





## Best Linear Prediction

Recall from Definition 2.14 that a process $(X_t)$ is *weakly stationary*
when $\E[X_t]=\mu$ for a single $\mu$ and
$\Cov(X_t,X_{t+h})=\gamma(h)$ for a function $\gamma$ of the lag $h$
alone. The function $\gamma$ is the *autocovariance function* (ACVF).

A full process law determines the conditional mean. If we are unwilling speculate on a full process law, but still want a reasonable forecast because we assume the process is weakly stationary, a principled way is to determine the best predictor within the class of affine functions
of our finite window of observed data.

**Definition 3.7 (Linear predictor).** Let
$$
W:=(X_{t_1},\ldots,X_{t_m})^\top
$$
be an observed window and let $Y = X_{t_0}$ be a square-integrable target. A *linear
predictor* of $Y$ from $W$ is an (affine) linear function of the observed window:
$$
g(W)=a_0+a^\top W,
$$
for some $a_0\in\R$ and $a\in\R^m$.

A best linear predictor exists for weak stationary processes with invertible covariance matrices.

**Proposition 3.8 (Best linear predictor).** Let $(X_t)$ be weakly stationary
with mean $\mu$ and ACVF $\gamma$. Fix a horizon $h\ge1$ and use the last $p$
observations as the window. Define
$$
W_{T,p}:=
\begin{pmatrix}
X_T-\mu\\
X_{T-1}-\mu\\
\vdots\\
X_{T-p+1}-\mu
\end{pmatrix},
\qquad
\Gamma_p:=\bigl(\gamma(i-j)\bigr)_{i,j=1}^p,
\qquad
\gamma_{p,h}:=
\begin{pmatrix}
\gamma(h)\\
\gamma(h+1)\\
\vdots\\
\gamma(h+p-1)
\end{pmatrix}.
$$
If $\Gamma_p$ is nonsingular, then among all linear predictors of $X_{T+h}$
from $(X_T,\ldots,X_{T-p+1})$, the predictor with the smallest mean squared
error is
$$
\widehat X^{\mathrm{lin}}_{T+h\mid T}
=\mu+a_h^\top W_{T,p}
$$
with
$$
a_h=\Gamma_p^{-1}\gamma_{p,h}.
$$
Its mean squared forecast error is
$$
\gamma(0)-\gamma_{p,h}^\top\Gamma_p^{-1}\gamma_{p,h}.
$$

*Proof.* Any linear predictor of $X_{T+h}$ from $(X_T,\ldots,X_{T-p+1})$ can
be written as $\mu+a^\top W_{T,p}$. Expanding the squared error gives
$$
\E[(X_{T+h}-\mu-a^\top W_{T,p})^2]
=\gamma(0)-2a^\top\gamma_{p,h}+a^\top\Gamma_pa.
$$
Completing the square, or differentiating with respect to $a$, gives the
stated coefficient and minimum. $\square$

Weak stationarity makes $\Gamma_p$ and $\gamma_{p,h}$ functions of lags. The
same coefficient vector therefore applies at every forecast origin $T$.
Without stationarity, the corresponding covariance matrices can depend on
$T$. 

**Remark 3.9 (The Gaussian case, and what it adds).** When the target and
window are jointly Gaussian, the conditional mean is affine, so the best
linear predictor of Proposition 3.8 coincides with the conditional mean of
Proposition 3.5, and its mean squared error equals the conditional variance.



**Example 3.10 (One AR(1), two noise specifications).** Consider the stationary
causal solution of the (shifted) AR(1) 
$$
X_t-\mu=\phi(X_{t-1}-\mu)+Z_t,
\qquad
|\phi|<1.
$$
Iterating forward (check) gives
$$
X_{T+h}-\mu
=\phi^h(X_T-\mu)
+\sum_{j=0}^{h-1}\phi^jZ_{T+h-j}.
$$

If we suppose only that
$(Z_t)\sim\mathrm{WN}(0,\sigma^2)$, the causal representation is
$$
X_t-\mu=\sum_{k=0}^{\infty}\phi^kZ_{t-k}.
$$
By Proposition 3.8, the best linear predictor of $X_{T+h}$ from $(X_1,\ldots,X_T)$ is (check!)
$$
\widehat X^{\mathrm{lin}}_{T+h\mid T}
=\mu+\phi^h(X_T-\mu)
$$
with mean squared error (check this too!)
$$
s_h^2:=\sigma^2\sum_{j=0}^{h-1}\phi^{2j}
=\sigma^2\frac{1-\phi^{2h}}{1-\phi^2}.
$$
The white noise specification does not identify the conditional mean. Indeed, by linearity of (conditional) expectation,
$$
\begin{aligned}
\E[X_{T+h}\mid X_1,\ldots,X_T]
&=\mu+\phi^h(X_T-\mu)\\
&\quad+\sum_{j=0}^{h-1}\phi^j
\E[Z_{T+h-j}\mid X_1,\ldots,X_T].
\end{aligned}
$$
The white-noise conditions make the future innovations uncorrelated with the
observed window, but they do not force the conditional expectations in this
display to be zero ($(Z_t)$ being iid would suffice in this case -- check!). Thus the conditional mean is not determined by this
second-order specification, and can differ from the best linear predictor.

Now strengthen the assumption to
$$
Z_t\overset{\mathrm{iid}}{\sim}N(0,\sigma^2).
$$
The future innovations are then independent of $(X_1,\ldots,X_T)$, so every
conditional expectation in the preceding display is zero. Their weighted sum
is also Gaussian. Consequently,
$$
X_{T+h}\mid X_1,\ldots,X_T
\sim N\!\left(
\mu+\phi^h(X_T-\mu),
\ s_h^2
\right).
$$

Under white noise, $\mu+\phi^h(X_T-\mu)$ and $s_h^2$ describe a linear
projection and its marginal mean squared error. Under i.i.d. Gaussian noise,
they are the conditional mean and conditional variance. In either case, the linear forecast
approaches $\mu$ and $s_h^2$ approaches $\sigma^2/(1-\phi^2)$ as our forecast window grows (as $h\to\infty$) due to stationarity. ◊


## From Predictive Moments to Predictive Intervals

**Definition 3.11 (Prediction set and coverage).** Let $Y=X_{T+h}\in\R^d$
be the future quantity of interest and let $W$ denote the information used to
predict it. A *prediction set* is a set $I_h(W)\subseteq\R^d$ computed from
$W$. When $Y$ is scalar, $I_h(W)$ is often an interval.

The set has *conditional coverage* at least $1-\alpha$ given $W$ if
$$
\P\{Y\in I_h(W)\mid W\}\geq 1-\alpha
$$
for every relevant observed window. It has *marginal coverage* at least
$1-\alpha$ if
$$
\P\{Y\in I_h(W)\}\geq 1-\alpha.
$$
Marginal coverage does not mean that the prediction set is nonrandom:
$I_h(W)$ may depend on $W$. It means that coverage is averaged over the
possible observed windows and their possible futures. Conditional coverage
for every window implies marginal coverage by the law of total probability,
$$
\P\{Y\in I_h(W)\}
=\E\!\left[\P\{Y\in I_h(W)\mid W\}\right]
\geq 1-\alpha,
$$
but marginal coverage alone need not imply conditional coverage. ◊

Write
$$
m_h(W):=\E[X_{T+h}\mid W],
\qquad
v_h(W):=\Var(X_{T+h}\mid W).
$$
Suppose these two conditional moments are available, while the shape of the
conditional distribution is left unspecified.

The predictive mean describes the center of the forecast, while the predictive
variance describes its scale. By the conditional form of Chebyshev's
inequality,
$$
m_h(W)\ \pm\ \sqrt{\frac{v_h(W)}{\alpha}}
$$
has conditional coverage at least $1-\alpha$ for each fixed horizon $h$. At
95%, this moment-only interval extends approximately
$4.47\sqrt{v_h(W)}$ on either side of the predictive mean. It is valid under
any known law with finite conditional variance, but it is usually wide.

There is a related guarantee when only the process mean and ACVF are known.
Let $s_h^2$ be the mean squared error of the best linear predictor from
Proposition 3.8. Ordinary Chebyshev's inequality gives
$$
\P\!\left(
  |X_{T+h}-\widehat X^{\mathrm{lin}}_{T+h\mid T}|
  \leq\frac{s_h}{\sqrt{\alpha}}
\right)
\geq1-\alpha.
$$
This gives only marginal coverage. The ACVF alone does not generally determine
the conditional variance at the particular window we observed, so it does not
provide the preceding conditional-coverage statement. Example 3.10 shows the
distinction: its white-noise specification supports this marginal bound, while
its i.i.d. Gaussian specification determines a conditional law. Conditional coverage with the Chebyshev bound requires us to specify the conditional mean and variance.

Can we do better if we know the shape of the conditional distribution? Let
$$
G_{h,W}(x):=\P(X_{T+h}\leq x\mid W)
$$
denote its conditional CDF. When $G_{h,W}$ is continuous, the equal-tail
interval
$$
\left[
G_{h,W}^{-1}(\alpha/2),
G_{h,W}^{-1}(1-\alpha/2)
\right]
$$
has conditional coverage exactly $1-\alpha$. In the special case
$$
X_{T+h}\mid W\sim N(m_h(W),v_h(W)),
$$
the equal-tail interval becomes
$$
m_h(W)\pm z_{1-\alpha/2}\sqrt{v_h(W)}
$$
because $z_q$ is the $q$th quantile of a standard Gaussian distribution. At
95%, $z_{0.975}\approx1.96$, much smaller than the moment-only multiplier
$4.47$. The sharper interval comes from knowing the Gaussian shape; it does not
follow from the mean and variance alone.

For intervals $I_1(W),\ldots,I_H(W)$, pointwise conditional coverage means
$$
\P\{X_{T+h}\in I_h(W)\mid W\}\geq1-\alpha
$$
for each fixed $h$. Simultaneous coverage instead concerns
$$
\P\!\left\{
X_{T+1}\in I_1(W),\ldots,X_{T+H}\in I_H(W)
\mathrel{\Big|} W
\right\}\geq1-\alpha.
$$
The first collection of statements does not imply the second. Thus, intervals
plotted over several horizons form a pointwise prediction fan unless stated
otherwise.



**Example 3.12 (What Gaussianity buys for a random walk).** Starting from a
fixed $X_0$, suppose
$$
X_t=X_{t-1}+\delta+Z_t,
\qquad
(Z_t)\sim\mathrm{WN}(0,\sigma^2).
$$
At horizon $h$,
$$
X_{T+h}
=X_T+h\delta+\sum_{j=1}^h Z_{T+j}.
$$
White noise fixes the mean and variance of the future sum:
$$
\E\!\left[\sum_{j=1}^h Z_{T+j}\right]=0,
\qquad
\Var\!\left(\sum_{j=1}^h Z_{T+j}\right)=h\sigma^2.
$$
It also makes the future sum uncorrelated with the observed past. Therefore
$X_T+h\delta$ is the best linear forecast from the observed window, with mean
squared error $h\sigma^2$. Chebyshev's inequality gives the pointwise marginal
guarantee
$$
\P\!\left\{
  \left|X_{T+h}-(X_T+h\delta)\right|
  \leq\sigma\sqrt{\frac{h}{\alpha}}
\right\}
\geq1-\alpha.
$$

Now strengthen the noise specification to
$$
Z_t\overset{\mathrm{iid}}{\sim}N(0,\sigma^2).
$$
The future innovations are then independent of the observed window and their
sum is Gaussian, so
$$
X_{T+h}\mid X_0,\ldots,X_T
\sim N(X_T+h\delta,h\sigma^2).
$$
The exact pointwise conditional interval is
$$
X_T+h\delta
\ \pm\
z_{1-\alpha/2}\sigma\sqrt h.
$$
At 95%, the Chebyshev and Gaussian half-widths are respectively
$$
4.47\sigma\sqrt h
\qquad\text{and}\qquad
1.96\sigma\sqrt h.
$$
Both fans widen at rate $\sqrt h$, but the moment-only fan is about
$4.47/1.96\approx2.28$ times as wide at every horizon.

```{r}
#| label: fig-random-walk-fans
#| fig-height: 4.2
#| code-fold: true
#| code-summary: "R code"
#| fig-cap: "Two pointwise forecast fans for the same simulated random-walk history, forecast center, and vertical scale. Left: the exact 95% conditional fan under i.i.d. Gaussian noise. Right: the conservative Chebyshev band with at least 95% marginal coverage under every white-noise law. Simulated data."
set.seed(542)
rw_T <- 60
rw_H <- 30
rw_delta <- 0.12
rw_sigma <- 1
rw_alpha <- 0.05

rw_time_observed <- 0:rw_T
rw_x <- c(0, cumsum(rw_delta + rnorm(rw_T, sd = rw_sigma)))
rw_h <- 0:rw_H
rw_time_future <- rw_T + rw_h
rw_center <- tail(rw_x, 1) + rw_delta * rw_h

rw_gaussian_half <- qnorm(1 - rw_alpha / 2) * rw_sigma * sqrt(rw_h)
rw_chebyshev_half <- rw_sigma * sqrt(rw_h / rw_alpha)
rw_xlim <- c(rw_T - 30, rw_T + rw_H)
rw_recent <- rw_time_observed >= rw_xlim[1]
rw_ylim <- range(
  rw_x[rw_recent],
  rw_center - rw_chebyshev_half,
  rw_center + rw_chebyshev_half
)

draw_rw_fan <- function(half_width, fill, title, ylab) {
  plot(
    NA,
    xlim = rw_xlim,
    ylim = rw_ylim,
    xlab = "time",
    ylab = ylab,
    main = title
  )
  polygon(
    c(rw_time_future, rev(rw_time_future)),
    c(rw_center - half_width, rev(rw_center + half_width)),
    border = NA,
    col = fill
  )
  lines(rw_time_observed, rw_x, col = "#303030", lwd = 1.2)
  lines(rw_time_future, rw_center, col = "#00539B", lwd = 2)
  lines(rw_time_future, rw_center - half_width,
        col = "#555555", lty = 2)
  lines(rw_time_future, rw_center + half_width,
        col = "#555555", lty = 2)
  abline(v = rw_T, col = "#777777", lty = 3)
}

old_par <- par(no.readonly = TRUE)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
draw_rw_fan(
  rw_gaussian_half,
  adjustcolor("#00539B", alpha.f = 0.22),
  "Gaussian: 1.96 SD",
  "level"
)
draw_rw_fan(
  rw_chebyshev_half,
  adjustcolor("#E5A823", alpha.f = 0.28),
  "White noise: 4.47 SD",
  ""
)
par(old_par)
```

The common vertical scale makes the width difference visible. At each fixed
horizon, the Gaussian interval has 95% coverage conditional on the observed
history. White noise alone does not determine the conditional distribution
given that history, so the Chebyshev guarantee is marginal. If the future
innovations are additionally independent of the observed past, the Chebyshev
guarantee is conditional as well. Neither collection of intervals is a
simultaneous region for the whole future path. The level process is
nonstationary in either case; stationarity is not needed for this forecast
calculation. ◊




## Backcasting a Balloon from Radar Readings

**Example 3.13 (Where was the balloon before radar contact?).** A radar first
detects a suspected surveillance balloon at time $t=0$. One time step is ten
minutes, and the radar reports a position at $t=0,1,\ldots,8$. Write
$$
R_t=(R_{t,E},R_{t,N})^\top\in\R^2
\quad\text{and}\quad
O_t=(O_{t,E},O_{t,N})^\top\in\R^2
$$
for the balloon's true position and the radar reading, respectively. Both
coordinates are measured in kilometers, east and north of a fixed origin. The
observed data are $O_0,\ldots,O_8$. We want to infer $R_{-6}$, the unobserved
true position one hour before first contact.

The radar does not observe $R_t$ exactly. We use the observation model
$$
O_t=R_t+\varepsilon_t,
\qquad
\varepsilon_t\overset{\mathrm{iid}}{\sim}N_2(0,rI_2),
$$
where $r$ is the measurement-error variance in each coordinate.

To connect positions before and after first contact, we model the true route as
$$
R_t=b+dt+U_t,
\qquad
U_t=\phi U_{t-1}+Z_t,
\qquad
Z_t\overset{\mathrm{iid}}{\sim}N_2(0,qI_2).
$$
The vector $b\in\R^2$ is the mean position at $t=0$, and $d\in\R^2$ is the
mean displacement during one ten-minute step. The stationary AR(1) process
$(U_t)_{t\in\Z}$ describes departures from that mean route, with
$|\phi|<1$ and innovation variance $q$ in each coordinate. The sequences
$(Z_t)$ and $(\varepsilon_t)$ are independent. The quantities
$b,d,\phi,q$, and $r$ are known.

```{r}
#| label: balloon-data
#| include: false
set.seed(542)
balloon_t <- -10:12
phi <- 0.95
q <- 1
r <- 4
b <- c(100, 50)
drift <- c(6, 2)
tau2 <- q / (1 - phi^2)

trend <- function(t) {
  outer(t, drift) + matrix(b, nrow = length(t), ncol = 2, byrow = TRUE)
}

u <- matrix(NA_real_, nrow = length(balloon_t), ncol = 2)
u[1, ] <- rnorm(2, sd = sqrt(tau2))
for (j in 2:length(balloon_t)) {
  u[j, ] <- phi * u[j - 1, ] + rnorm(2, sd = sqrt(q))
}
route <- trend(balloon_t) + u

radar_t <- 0:8
radar_index <- match(radar_t, balloon_t)
radar <- route[radar_index, ] +
  matrix(rnorm(2 * length(radar_t), sd = sqrt(r)), ncol = 2)

C_radar <- tau2 * phi^abs(outer(radar_t, radar_t, "-")) +
  diag(r, length(radar_t))
back_t <- -6:0
K <- sapply(back_t, function(s) tau2 * phi^abs(s - radar_t))
centered_radar <- radar - trend(radar_t)
back_mean <- trend(back_t) + t(K) %*% solve(C_radar, centered_radar)
back_var <- tau2 - colSums(K * solve(C_radar, K))

origin_mean <- back_mean[1, ]
origin_var <- back_var[1]
origin_radius <- sqrt(qchisq(0.95, df = 2) * origin_var)
origin_truth <- route[match(back_t[1], balloon_t), ]
```

For the simulation, we take
$$
\phi=0.95,\quad q=1,\quad r=4,
\quad b=(100,50)^\top,\quad d=(6,2)^\top.
$$
The code generates both the hidden route $(R_t)$ and the noisy readings
$(O_t)$. The backcast uses only the following nine observed readings; the
simulated route is retained so that we can later compare the inference with
the truth.

```{r}
#| label: tbl-radar
#| echo: false
radar_table <- data.frame(
  `Time t` = radar_t,
  `Minutes after first contact` = 10 * radar_t,
  `Observed east (km)` = round(radar[, 1], 1),
  `Observed north (km)` = round(radar[, 2], 1),
  check.names = FALSE
)
knitr::kable(
  radar_table,
  caption = "The nine observed radar readings from first contact onward."
)
```

For observed times $t_1,\ldots,t_m$, let the $i$th row of the $m\times2$
matrix $O$ be $O_{t_i}^\top$, and let the $i$th row of $M$ be the corresponding
mean position $(b+dt_i)^\top$. Write $\tau^2:=q/(1-\phi^2)$ and define
$$
C_{ij}:=\tau^2\phi^{|t_i-t_j|}+r\,1\{i=j\},
\qquad
c_s:=\tau^2
\begin{pmatrix}
\phi^{|s-t_1|}\\
\vdots\\
\phi^{|s-t_m|}
\end{pmatrix}.
$$
The matrix $C$ is the covariance matrix of either observed coordinate across
the radar window, while $c_s$ contains the covariances between the true
coordinate at time $s$ and those readings. Applying Proposition 3.5 separately
to the east and north coordinates gives
$$
\E[R_s\mid O_{t_1},\ldots,O_{t_m}]
=b+ds+c_s^\top C^{-1}(O-M),
$$
and
$$
\Var(R_s\mid O_{t_1},\ldots,O_{t_m})
=\left(\tau^2-c_s^\top C^{-1}c_s\right)I_2,
$$
where row $i$ of $O-M$ is $O_{t_i}-(b+dt_i)$.

At $s=-6$, one hour before first contact, the conditional mean is
$(`r sprintf("%.1f", origin_mean[1])`,
`r sprintf("%.1f", origin_mean[2])`)^\top$ km and the conditional covariance
is `r sprintf("%.2f", origin_var)` $I_2$ km$^2$. The pointwise 95% location
region is a disk of radius `r sprintf("%.1f", origin_radius)` km around that
mean. The simulated true position was
$(`r sprintf("%.1f", origin_truth[1])`,
`r sprintf("%.1f", origin_truth[2])`)^\top$ km.

```{r}
#| label: fig-balloon-backcast
#| fig-cap: "A simulated balloon route (gray), observed radar readings from first contact onward (blue), and conditional-mean backcasts for the preceding hour (red). Only the blue readings enter the backcast; the gray route is unobserved truth retained for comparison. The circle is a pointwise 95% region for the position at t = −6. Simulated data."
angle <- seq(0, 2 * pi, length.out = 200)
circle <- cbind(
  origin_mean[1] + origin_radius * cos(angle),
  origin_mean[2] + origin_radius * sin(angle)
)
par(bg = "gray92")
display_index <- balloon_t >= min(back_t) & balloon_t <= max(radar_t)
plot(route[display_index, 1], route[display_index, 2],
     type = "l", lty = 2, col = "gray45",
     xlim = range(c(route[display_index, 1], radar[, 1], back_mean[, 1], circle[, 1])),
     ylim = range(c(route[display_index, 2], radar[, 2], back_mean[, 2], circle[, 2])),
     asp = 1,
     xlab = "East coordinate (km)", ylab = "North coordinate (km)",
     main = "Backcasting a route from a radar window")
polygon(circle, border = 2, col = adjustcolor(2, alpha.f = 0.12))
lines(back_mean[, 1], back_mean[, 2], col = 2, lwd = 2)
points(radar[, 1], radar[, 2], pch = 19, col = 4)
points(origin_mean[1], origin_mean[2], pch = 4, col = 2, lwd = 2, cex = 1.4)
legend("topleft",
       legend = c("Unobserved true route", "Observed radar readings", "Backcast mean", "95% region at t = -6"),
       col = c("gray45", 4, 2, 2), lty = c(2, NA, 1, 1),
       pch = c(NA, 19, NA, NA), lwd = c(1, NA, 2, 1), bty = "n")
```

Radar error prevents the first reading $O_0$ from revealing $R_0$ exactly.
Later readings help locate that latent boundary position, and the correlation in
$(U_t)$ carries the information backward. If $r=0$, $R_0$ is observed exactly;
then the later readings add no information about $R_s$ for $s<0$ under this
Gaussian AR(1) law. ◊

## Looking Ahead

Every mean, covariance, coefficient, and conditional distribution in this
workbook was known. Workbook 4 places a prior distribution on an unknown
parameter and carries its posterior uncertainty into the forecast. The
workbooks on one-path learnability, identification, sampling uncertainty, and
transformation then develop the prerequisites for fitting a model from one
observed path.

## Exercises

::: {.exercise}
**Exercise 3.1 (The loss determines the point forecast).** Let $Y$ be a
square-integrable target and $W$ an observed window.

(a) Reproduce the decomposition in Proposition 3.2 for squared loss.

(b) Suppose $Y\mid W=w$ takes the values $0$ and $10$ with probabilities
$0.9$ and $0.1$. Find the conditional mean and a conditional median.

(c) Compare the two candidates under conditional squared loss and conditional
absolute loss.

(d) Explain why the phrase “best forecast” is incomplete until the loss and
the available predictors are specified.
:::


::: {.exercise}
**Exercise 3.2 (Two-lag linear prediction).** Let $(X_t)$ be mean-zero and
weakly stationary, with $\gamma(0)>0$, and use $(X_T,X_{T-1})$ to predict (using a linear predictor)
$X_{T+1}$.

(a) Write $\Gamma_2$ and $\gamma_{2,1}$ from Proposition 3.8.

(b) Solve for the two coefficients when
$\gamma(0)=1$, $\gamma(1)=0.6$, and $\gamma(2)=0.2$.

(c) Compute the mean squared forecast error.

(d) State which step would fail if $\Gamma_2$ were singular.
:::


::: {.exercise}
**Exercise 3.3 (The Gaussian AR(1) in both directions).** Consider the i.i.d.
Gaussian specification in Example 3.10 with $\mu=2$, $\phi=0.8$, and
$\sigma=1$.

(a) Given $X_T=5$, compute the point forecasts and forecast-error variances for
$h=1,5,20$.

(b) Construct the corresponding pointwise 95% prediction intervals.

(c) Given $X_0=5$, compute the conditional mean and variance of $X_{-5}$.

(d) Verify directly from the ACVF that
$X_{-5}-\mu-\phi^5(X_0-\mu)$ is uncorrelated with $X_t$ for every $t\ge0$.
Explain where Gaussianity is used to turn this calculation into the full-window
backcast.
:::


::: {.exercise}
**Exercise 3.4 (One-observation prediction for an MA(1)).** Let
$$
X_t=Z_t+\theta Z_{t-1},
$$
where $(Z_t)$ is white noise with variance $\sigma^2$. Use only $X_T$ to
predict $X_{T+1}$ linearly.

(a) Compute the coefficient in the best linear predictor.

(b) Compute its mean squared forecast error.

(c) Explain why the white-noise assumptions alone do not determine
$\E[X_{T+1}\mid X_T]$ or an exact 95% prediction interval.

(d) State what changes if the process is jointly Gaussian.
:::


::: {.exercise}
**Exercise 3.5 (Where did the balloon come from?).** This exercise reconstructs
the backcast in Example 3.13 and then measures what the later radar readings
contribute. Use R for parts (a)--(d), show the code that produces each answer,
and report distances in kilometers. Answer parts (e)--(f) in writing.

Run the code in Example 3.13 before starting. It creates the following objects:

- `radar_t`, the observed times $0,\ldots,8$;
- `radar`, the $9\times2$ matrix $O$ of observed east and north readings;
- `phi`, `q`, and `r`, the known scalar parameters;
- `b` and `drift`, where `drift` represents $d$;
- `trend(t)`, which returns the mean route $b+dt$ at the supplied times.

The objects `route` and `origin_truth` contain the simulated true route. They
are unavailable in a real backcasting problem and may not be used below. Set
the target time to $s=-6$ and use
$$
\begin{aligned}
m_s&=b+ds+c_s^\top C^{-1}(O-M),\\
s_s^2&=\tau^2-c_s^\top C^{-1}c_s,
\qquad
\tau^2=\frac{q}{1-\phi^2}.
\end{aligned}
$$
The conditional distribution of the target is
$R_s\mid O_0,\ldots,O_8\sim N_2(m_s,s_s^2I_2)$.
In R, use `trend(s)` for $b+ds$.

(a) **In R:** Set `obs_t <- radar_t`, `O <- radar`, `M <- trend(obs_t)`, and
`s <- -6`. Construct $C$ and $c_s$ from the formulas in Example 3.13. Verify
that $C$ is $9\times9$, $c_s$ has length $9$, and $O-M$ is $9\times2$.

(b) **In R:** Use these objects to compute $m_{-6}$ and $s_{-6}^2$. Report the
two coordinates of the conditional mean in kilometers and the conditional
covariance matrix in square kilometers.

(c) **In R:** Compute
$$
\P(\|R_{-6}-m_{-6}\|\leq5\mid O_0,\ldots,O_8).
$$
Report the probability and interpret it as a statement about the hidden
position one hour before contact. *Hint: conditional on the radar readings,
$\|R_{-6}-m_{-6}\|^2/s_{-6}^2$ has a $\chi^2_2$ distribution; use `pchisq()`.*

(d) **In R:** Repeat part (b) using only $O_0$. In this calculation,
$C=\tau^2+r$ and $c_{-6}=\tau^2\phi^6$. Compare the first-reading and
nine-reading conditional variances in a small table. Report the absolute and
percentage reduction due to the eight later readings.

(e) **By hand:** Now set $r=0$. Let $C_0$ be the $9\times9$ covariance matrix
of either coordinate of $(R_0,\ldots,R_8)$ after removing the mean route, and
let $e_1=(1,0,\ldots,0)^\top$. Show that
$c_{-6}=\phi^6C_0e_1$ and hence that
$c_{-6}^\top C_0^{-1}=\phi^6e_1^\top$. Explain what this says about the weights
assigned to readings after $t=0$.

(f) **In words:** In two or three sentences, explain why the disk plotted in
Example 3.13 is a 95% region for $R_{-6}$ but not a 95% region for the
six-position route $(R_{-6},\ldots,R_{-1})$.
:::


::: {.exercise}
**Exercise 3.6 (Forecasting volatility).** Recall that a causal GARCH(1,1)
process satisfies
$$
X_t=\sqrt{V_t}\,Z_t,
\qquad
V_t=\omega+\alpha X_{t-1}^2+\beta V_{t-1},
$$
where $\omega>0$, $\alpha,\beta\geq0$, and the innovations $(Z_t)$ are i.i.d.
with mean zero and variance one. Each innovation is independent of the return
history before it. Suppose $\alpha+\beta<1$, the parameters are known, and the
observed history $\cH_T=(\ldots,X_{T-1},X_T)$ determines the current state
$V_{T+1}$. Parts (a)--(d) require no further assumption on the innovation
distribution.

At horizon $h=1$, $V_{T+1}$ is known from $\cH_T$. At longer horizons,
$V_{T+h}$ is random because it depends on returns that have not yet been
observed. In this exercise, *forecasting volatility* means computing the
conditional variance of the future return,
$$
v_h:=\Var(X_{T+h}\mid\cH_T).
$$
Its square root is the corresponding conditional standard deviation. Part (a)
connects this forecast to the future GARCH state $V_{T+h}$.

(a) Show that, for every $h\geq1$,
$$
\E[X_{T+h}\mid\cH_T]=0
\quad\text{and}\quad
\Var(X_{T+h}\mid\cH_T)=\E[V_{T+h}\mid\cH_T].
$$

(b) Show that
$$
v_1=V_{T+1},
\qquad
v_h=\omega+(\alpha+\beta)v_{h-1},
\quad h\geq2,
$$
and solve the recursion for $v_h$.

(c) Take $\omega=0.1$, $\alpha=0.1$, $\beta=0.8$, and $V_{T+1}=4$.
Compute $v_h$ for $h=1,2,5,20$.

(d) Find the limiting conditional variance and explain what the limit means.

(e) Now suppose additionally that the innovations are standard Gaussian.
Construct the exact one-step 95% prediction interval for $X_{T+1}$.

(f) Explain why replacing $V_{T+1}$ by $v_h$ inside a Gaussian interval need
not give an exact interval for $h>1$.

(g) Interpret the forecasts: which forecast remains constant across observed
histories, and which forecast changes with the observed history?
:::


------------------------------------------------------------------------

*Sources for this workbook: van der Vaart (2010), Time Series, lecture notes,
VU Amsterdam, Ch. 2, for best prediction and linear projection; Brockwell and
Davis (1991), Time Series: Theory and Methods, 2nd ed., Springer, Ch. 5, for
linear prediction; Shumway and Stoffer (2025), Time Series Analysis and Its
Applications, 5th ed., Springer, §3.4, for forecasting and prediction
intervals.*