---
title: "Workbook 15 — Modeling What Is Hidden"
subtitle: "Forecasting with a Hidden State"
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.8
    fig-dpi: 150
execute:
  warning: false
  message: false
editor:
  markdown:
    wrap: 72
---

::: hidden
$$\def\E{\mathbb{E}}
\def\P{\mathbb{P}}
\def\R{\mathbb{R}}
\def\C{\mathbb{C}}
\def\Z{\mathbb{Z}}
\def\N{\mathbb{N}}
\def\cH{\mathcal{H}}
\def\Var{\operatorname{Var}}
\def\Cov{\operatorname{Cov}}
\def\Corr{\operatorname{Corr}}$$
:::

```{=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: .91em; }
@media print {
  #TOC, .quarto-title-meta, .code-copy-button { display: none !important; }
  body { font-size: 10pt; }
  h2 { break-after: avoid; }
  figure, table { break-inside: avoid; }
}
</style>
```

```{r}
#| label: setup
#| include: false
options(digits = 6)
knitr::opts_chunk$set(comment = "#>")
source("R/volatility.R")
source("R/dlm.R")
source("R/hmm.R")
```

In Workbook 2 we asked which processes satisfy a proposed set of equations for a given choice of parameters. In a state-space model the equations describe a hidden state and an observed series together. For the Dow Jones Industrial Average, we observe daily closing prices and calculate returns from them. We do not observe volatility directly, but we can infer it from the size of those returns. We can model the log variance of the daily return as a hidden state. When it is high, large gains and large losses are both more likely. If high volatility persists, large moves tend to occur close together.

![Daily log returns of the Dow Jones Industrial Average, April 2006–April 2016, from Workbook 1, Example 1.11. Large gains and losses cluster in late 2008. The dashed line marks the mean daily log return.](data/djia-daily-log-returns.png){#fig-djia-returns fig-alt="Daily log returns fluctuate around zero. The widest swings in both directions occur close together in late 2008."}

The return history can then help us estimate a parameter shared across days, such as the persistence of volatility, and infer a random quantity specific to one day, such as today's volatility. We use both to forecast the chance of a large move tomorrow.

## Stochastic volatility and the next return {#how-large-might-the-next-return-be}

**Example 15.1 (What can we forecast when the mean stays fixed?).** In Workbook 1, Example 1.13, we drew a log variance $h_t$ for each day and used it to generate a return $X_t$. Now suppose we keep only the returns. The equations that generated them were
$$
\begin{aligned}
h_t &= \mu+\phi(h_{t-1}-\mu)+\eta_t, & \eta_t &\sim N(0,\sigma_\eta^2),\\
X_t &= \delta+e^{h_t/2}Z_t, & Z_t &\sim N(0,1).
\end{aligned}
$$
The parameters are the mean return $\delta$, the long-run mean log variance $\mu$, the persistence coefficient $\phi$, and the state shock variance $\sigma_\eta^2>0$. Under a given candidate law these numbers stay fixed across days and realizations, even when their values are unknown to us. For the simulation we set them to $0.02$, $-1$, $0.95$, and $0.09$, respectively. The shocks in both equations are mutually independent over time and across the two sequences.

Does the first equation give a stationary log variance process? For $|\phi|<1$, the AR(1) construction in Workbook 2, Example 2.8, gives the stationary solution
$$
h_t=\mu+\sum_{j=0}^{\infty}\phi^j\eta_{t-j}.
$$
The series converges in mean square because $\sigma_\eta^2\sum_{j\geq0}\phi^{2j}<\infty$. Separating the first term leaves $\eta_t+\phi(h_{t-1}-\mu)$ for $h_t-\mu$, as the state equation requires. We drew the initial log variance from its stationary law, independently of subsequent state shocks and all return shocks. The return equation then gives $X_t$ from $h_t$ and an independent $Z_t$. Even with all four parameters specified, both $h_t$ and $X_t$ are random. We see only the realized returns $x_t$.

Given $h_t$, the return has standard deviation $e^{h_t/2}$ and mean $\delta=0.02$. Every possible value of tomorrow's log variance gives the same mean, so the forecast mean is $0.02$ after any history. What can the history tell us about the *size* of tomorrow's move? Call a move large if the return falls more than one unit from $0.02$ in either direction. The figure shows the last ten returns available for two forecasts, one made at the end of day 140 and the other at the end of day 155. After which history would we assign a higher probability to a large move the next day?

```{r}
#| label: fig-sv-histories
#| echo: false
#| fig-height: 3.15
#| fig-alt: "Two ten-day windows on the same return scale. All returns in the first window lie within one unit of the mean. Two in the second window fall below that band, including the drop on day 153."
#| fig-cap: "Simulated daily returns (seed 542). The dashed line marks the mean, and the shaded band extends one unit above and below it. Each panel ends on the last day whose return is available for the forecast."
sv_x <- sv_record()
sv_origins <- c(140L, 155L)
sv_history_plot(sv_x, sv_origins)
```

All ten returns in the left panel lie inside the band. In the right panel two fall below it, including the return of $-3.64$ on day 153. A large negative shock $Z_{153}$ could have produced that drop even at moderate volatility. A fall at least this large is more probable at higher volatility, so it gives us evidence that volatility was high on day 153. With $\phi=0.95$ in the simulation, high log variance on day 153 can persist into day 156, when the return has a fresh shock of its own.

We would therefore expect a higher probability of a large move on day 156 than on day 141. To calculate each probability, we use the returns available through the preceding day to predict tomorrow's volatility. For each possible volatility the return equation gives a chance of a large move. We average those chances over the predicted distribution of volatility. ◊

```{r}
#| label: sv-record-code
#| eval: false
#| code-fold: true
#| code-summary: "Reproduce the simulated returns and figure"
# From the workbook folder:
source("R/volatility.R")
# Start the log variance in its stationary Gaussian distribution.
sv_x <- sv_record(seed = 542, n = 400, delta = 0.02, mu = -1,
                  phi = 0.95, sigma_eta = 0.3)
sv_origins <- c(140L, 155L)
sv_history_plot(sv_x, sv_origins)
```

## Markov states and observation laws {#what-must-the-state-remember}

How does the model predict tomorrow's volatility? First suppose we knew today's log variance $h_t=s$. The first equation in Example 15.1 gives tomorrow's log variance the distribution $N(\mu+\phi(s-\mu),\sigma_\eta^2)$. Knowing earlier log variances as well would leave this forecast unchanged, because tomorrow's shock is independent of them. This is the *Markov property*.

The rainfall model in Workbook 3, Example 3.4, had the same property. A rain gauge recorded the realized rain indicator, $1$ on a wet day and $0$ on a dry day. The model supplied the two transition probabilities. Whatever the preceding weather, tomorrow's indicator equals $1$ with probability $0.7$ after a wet day and $0.2$ after a dry day. The *state space* is the set of values the state can take, here $\{0,1\}$ for the rain indicator and $\R$ for the log variance.

**Definition 15.2 (Markov process).** Let $I\subseteq\Z$ be the time index set and let $(S_t)_{t\in I}$ be a family of random variables on a common probability space with values in a set $E$. Writing $\mathcal L(Y\mid Z)$ for the conditional distribution of $Y$ given $Z$, we call $(S_t)_{t\in I}$ *Markov* if, for every integer $k\geq1$ and all time points $t_0<\cdots<t_k<u$ in $I$,
$$
\mathcal L(S_u\mid S_{t_0},\ldots,S_{t_k})
=
\mathcal L(S_u\mid S_{t_k}).
$$
Once we condition on the most recent state $S_{t_k}$, learning the earlier states leaves the conditional distribution of $S_u$ unchanged.

In the return model, the Markov property applies to $S_t=h_t$. What does the model say about the returns $X_t$? Given the log variances $h_1,\ldots,h_T$, the only randomness left in each return $X_t=\delta+e^{h_t/2}Z_t$ comes from its own shock $Z_t$. Those shocks are independent of each other and of the log variances. Given the log variances, the returns are therefore independent, each with distribution $N(\delta,e^{h_t})$.

**Definition 15.3 (State-space representation).** Let $(S_t,X_t)_{t\in I}$ be a joint process with time index set $I\subseteq\Z$. The pair gives a *state-space representation* of $(X_t)$ with state $(S_t)$ if $(S_t)$ is Markov and, for every $k\geq1$, every $t_1<\cdots<t_k$ in $I$, and all events $\{X_{t_j}\in B_j\}$,
$$
\P(X_{t_1}\in B_1,\ldots,X_{t_k}\in B_k\mid S_{t_1},\ldots,S_{t_k})
=\prod_{j=1}^k \P(X_{t_j}\in B_j\mid S_{t_j}).
$$

We observe one record $x_{1:T}=(x_1,\ldots,x_T)$. We model it by introducing a state $S_t$ for each time $t$. We choose what $S_t$ represents and how it relates to $X_t$. Each candidate parameter value $\theta$ then specifies the initial distribution of $S_1$, the transition law of $S_t$ given $S_{t-1}$, and the observation law of $X_t$ given $S_t$.

Fix $\theta$ and draw $S_1$ from its initial law. Given $S_1$, draw $X_1$ from its observation law. Draw $S_2$ from its transition law given $S_1$, then $X_2$ from its observation law given $S_2$. Repeating these draws determines a unique joint law on the window $1,\ldots,T$. When we also require a stationary solution on $\Z$, we must establish that such a process exists, as we did for the log variance in Example 15.1.

The horizontal arrows show the transitions. Given the states, the observations are independent, and each downward arrow gives the observation law for $X_t$ conditional on its own state $S_t$.
$$
\begin{array}{ccccc}
 S_{t-1} & \longrightarrow & S_t & \longrightarrow & S_{t+1} \\
 \downarrow && \downarrow && \downarrow \\
 X_{t-1} && X_t && X_{t+1}
\end{array}
$$

### Likelihood of the observed record

What we fit is the family of distributions that the model assigns to the observed record. As $\theta$ varies, the marginal laws for the observation process $(X_t)$ form its *state-space model* (Workbook 2, Definition 2.5). Unknown parameters can occur in any of the three laws. For example, we never see a transition between log variances, but changing the persistence can change how often large returns occur close together. The realized state $S_t$ is a different unknown. Even if we knew $\theta$, we would still have to infer which states were plausible on day $t$.

To get the likelihood, first fix $\theta$. Suppose for a moment that we knew the entire state path $s_1,\ldots,s_T$. Write $g_{\theta,t}(x_t\mid s_t)$ for the observation mass or density at $x_t$ given $S_t=s_t$. Conditional on the path, the observations are independent, so the mass or density of the recorded observations would be
$$
\prod_{t=1}^T g_{\theta,t}(x_t\mid s_t).
$$

We do not know that path, but the initial and transition laws tell us how to weight the possible paths. Averaging over them gives the *observed-data likelihood*,
$$
L_T(\theta;x_{1:T})=p_\theta(x_{1:T})
=\E_\theta\!\left[\prod_{t=1}^T g_{\theta,t}(x_t\mid S_t)\right].
$$

A maximum likelihood estimate maximizes $L_T(\theta;x_{1:T})$ over $\theta$ for the record we observed. If we give $\theta$ a prior mass or density $\pi(\theta)$, the posterior mass or density is proportional to $\pi(\theta)L_T(\theta;x_{1:T})$, as in Workbook 4.

## Updating a hidden state {#how-a-forecast-changes-with-a-new-observation}

Suppose a lamp has low, medium, and high settings controlling how often it flashes. A sensor records the number of flashes each second, but the settings are unrecorded. The first count is seven. We want to infer the setting during that second.

Write $S_1$ for the first setting and $X_1$ for the first count. For this simulation, we give each setting probability $1/3$ and specify mean rates of $2$, $6$, and $12$ flashes per second for low, medium, and high. With $\lambda_s$ denoting the rate at setting $s$, the model is
$$
\P(S_1=s)=\frac13,
\qquad X_1\mid S_1=s\sim\operatorname{Poisson}(\lambda_s).
$$
The rate is an average. Even with the setting held fixed, the count varies, so seven flashes are possible at all three settings.

For a setting with rate $\lambda_s$, the probability of seven flashes is $e^{-\lambda_s}\lambda_s^7/7!$. This gives the probabilities in the second column below.

```{r}
#| label: first-count-state-update
#| echo: false
setting_names <- c("Low", "Medium", "High")
first_count_prob <- dpois(7, c(2, 6, 12))
first_state_prob <- first_count_prob / sum(first_count_prob)
knitr::kable(data.frame(setting_names, first_count_prob, first_state_prob),
             digits = c(0, 4, 3), row.names = FALSE,
             col.names = c("Setting", "Probability of seven flashes",
                           "Setting probability after observing seven"))
```

Since the settings were initially equally likely, Bayes' rule makes their updated probabilities proportional to their probabilities of producing seven flashes. For medium, the common factor $1/3$ cancels from the numerator and denominator,
$$
\P(S_1=\text{medium}\mid X_1=7)
\approx\frac{0.1377}{0.0034+0.1377+0.0437}
\approx0.745.
$$
Observing seven flashes raises the probability of medium from one third to about three quarters. The last column gives the resulting distribution of the setting. More generally, the conditional distribution of the current state given the observations so far is called its *filtering distribution*.

Before the first count, the probability of seven flashes was the average of the three entries in the second column, $(0.0034+0.1377+0.0437)/3\approx0.0616$. For a later observation $X_t$, the preceding observations $x_{1:t-1}=(x_1,\ldots,x_{t-1})$ tell us how to weight the states. Fix $\theta$. To obtain the conditional mass or density of $X_t$ at $x$, average the state-specific masses or densities $g_{\theta,t}(x\mid s)$ over the distribution of $S_t$ given those observations,
$$
p_\theta(x\mid x_{1:t-1})
=\E_\theta[g_{\theta,t}(x\mid S_t)\mid x_{1:t-1}].
$$

To find the probability that $S_t$ lies in a set $A$ after observing $x_t$, keep only the contributions from states in $A$ in the numerator. When $0<p_\theta(x_t\mid x_{1:t-1})<\infty$, Bayes' rule gives
$$
\P_\theta(S_t\in A\mid x_{1:t})=
\frac{\E_\theta[\mathbf{1}\{S_t\in A\}\,g_{\theta,t}(x_t\mid S_t)\mid x_{1:t-1}]}{p_\theta(x_t\mid x_{1:t-1})}.
$$
To predict $S_{t+1}$, average its transition law over this updated distribution of $S_t$. The denominator in Bayes' rule is the conditional mass or density at the observed $x_t$. Multiplying these factors for $t=2,\ldots,T$ by the mass or density of the first observation gives the likelihood of the record,
$$
L_T(\theta;x_{1:T})=p_\theta(x_1)\prod_{t=2}^T p_\theta(x_t\mid x_{1:t-1}).
$$

## Fitting persistence and forecasting a large return {#fitting-persistence-and-forecasting-a-large-return}

In Example 15.1 a large return was evidence of high volatility that day. Does it also change the forecast for tomorrow? The answer depends on $\phi$ in the transition mean $\mu+\phi(h_t-\mu)$. We compare $\phi=0$ with $\phi=0.95$ using the returns.

Here we fit only $\phi$. We take $\delta=0.02$, $\mu=-1$, and the stationary variance $v=\Var(h_t)=0.3^2/(1-0.95^2)$ of the log variance process as given. From the solution in Example 15.1, $\Var(h_t)=\sigma_\eta^2/(1-\phi^2)$. Holding $v$ fixed therefore sets $\sigma_\eta^2=v(1-\phi^2)$ for each candidate. The Gaussian solution gives $h_1\sim N(-1,v)$ as the starting law in both cases. Both candidates then give the same distribution for a single return. With $\phi=0$, even after a large return, tomorrow's log variance still has distribution $N(-1,v)$. With $\phi=0.95$, evidence of high volatility today changes the distribution of tomorrow's log variance. We can therefore compare the candidates by the forecasts they make along the series.

Fix one candidate $\phi$. Let $f_{t,\phi}(s\mid x_{1:t-1})$ be the conditional density of the log variance $h_t$ at $s$ given the preceding returns. On day 1 we start with the $N(-1,v)$ density. If $h_t=s$, the return has the $N(0.02,e^s)$ density, which we denote by $g(x\mid s)$. For the lamp we found the probability of seven flashes by averaging the three Poisson probabilities, one for each possible setting. Here the log variance can take any real value, so we average the return densities by integration,
$$
p_\phi(x\mid x_{1:t-1})=\int_{\R} g(x\mid s)f_{t,\phi}(s\mid x_{1:t-1})\,ds.
$$

Suppose we observe $X_t=x_t$. Bayes' rule updates the density of today's log variance using the return we just saw,
$$
f_{t,\phi}(s\mid x_{1:t})
=\frac{g(x_t\mid s)f_{t,\phi}(s\mid x_{1:t-1})}{p_\phi(x_t\mid x_{1:t-1})}.
$$
The update gives more weight to log variances $s$ for which $g(x_t\mid s)$ exceeds the forecast density $p_\phi(x_t\mid x_{1:t-1})$. At $h_t=s$, tomorrow's log variance has distribution $N(-1+\phi(s+1),v(1-\phi^2))$. To predict tomorrow's log variance we average these transition densities over the updated distribution of $h_t$. Repeating this prediction and update as returns arrive is called *filtering*.

Which persistence value better explains the first 100 returns? Under each candidate $\phi$, we start with the density of the first return. For each later return we use its conditional density given the earlier returns. The log likelihood adds the logarithms of these densities at the values we observed,
$$
\ell_{100}(\phi)=\log p_\phi(x_1)+\sum_{t=2}^{100}\log p_\phi(x_t\mid x_{1:t-1}).
$$
The [volatility code](R/volatility.R) approximates these integrals with weighted sums on a finite grid of log variances.

```{r}
#| label: sv-fit-persistence
#| code-fold: true
#| code-summary: "Compare the two candidates using days 1–100"
source("R/volatility.R")
sv_x <- sv_record()
sv_v <- 0.3^2 / (1 - 0.95^2)
sv_phi <- c(0, 0.95)
sv_loglik <- vapply(sv_phi, function(phi) {
  fit <- sv_filter(sv_x[1:100], phi = phi, sigma_eta = sqrt(sv_v * (1 - phi^2)))
  sum(fit$log_score)
}, numeric(1))
sv_selected <- sv_phi[which.max(sv_loglik)]
knitr::kable(data.frame(Persistence = sv_phi, `Log likelihood` = sv_loglik,
                        check.names = FALSE), digits = 2)
```

The 100 recorded returns have the larger joint density under $\phi=0.95$ (log likelihood $-94.38$), so they give the estimate $\hat\phi=0.95$. A different simulated record at the same true $\phi$ could select the other candidate. For this record we keep $\hat\phi$ fixed in the forecast calculation. Substituting the estimate in this way gives a *plug-in forecast*.

A large move in Example 15.1 is a return more than one unit from $\delta=0.02$. Suppose tomorrow's return standard deviation were $0.5$. A large move would require a standard normal shock of magnitude greater than $2$, which has probability $0.046$. With standard deviation $2$, a shock of magnitude greater than $0.5$ suffices, with probability $0.617$. For any log variance $s$, the return standard deviation is $e^{s/2}$. At that volatility, a move of more than one unit has probability
$$
\P(|X_{t+1}-\delta|>1\mid h_{t+1}=s)
=\P(|Z_{t+1}|>e^{-s/2})=2\Phi(-e^{-s/2}),
$$
where $\Phi$ is the standard normal distribution function.

For the day-141 forecast, we use the first 140 returns to find the predicted weights for $h_{141}$. We multiply the tail probability at each grid value by its predicted weight and sum. For day 156 we repeat the calculation using returns through day 155. Every $N(\delta,e^s)$ return distribution has mean $\delta$, so the forecast mean is always $0.02$. Its variance is the weighted average of the state-specific variances, $\E[e^{h_{t+1}}\mid x_{1:t}]$.

```{r}
#| label: sv-forecast-origins
#| code-fold: true
#| code-summary: "Calculate the forecasts using each history"
sv_at <- function(t, phi = sv_selected) {
  sv_filter(sv_x[seq_len(t)], phi = phi, sigma_eta = sqrt(sv_v * (1 - phi^2)))
}
sv_origins <- c(140L, 155L)
sv_forecasts <- lapply(sv_origins, sv_at)
sv_risks <- vapply(sv_forecasts, function(fit) {
  tail(sv_tail(fit, threshold = 1), 1)
}, numeric(1))
```

```{r}
#| label: fig-sv-forecast-densities
#| echo: false
#| fig-height: 3.3
#| fig-alt: "The forecast for day 141 is concentrated near 0.02. The forecast for day 156 spreads farther into both tails."
#| fig-cap: "Forecast densities for the simulated returns. The shaded regions mark returns below −0.98 or above 1.02. Both panels use the same scales."
sv_forecast_plot(sv_forecasts)
```

After the quiet stretch ending on day 140, the probability of a move into either shaded tail on day 141 is `r sprintf("%.1f%%", 100 * sv_risks[1])`. The history through day 155 includes the $-3.64$ return on day 153. The fitted model gives higher log variances more weight on day 156, and the forecast probability of a large move is `r sprintf("%.1f%%", 100 * sv_risks[2])`.

How much did the day-153 drop change the forecast? This time keep the target fixed at day 154. Using returns through day 152, calculate the probability of a large move on day 154. Then include the return of $-3.64$ on day 153 and forecast the same event again. The middle column is the conditional mean of day-153 volatility, $e^{h_{153}/2}$, given the returns listed in that row. The last column gives the forecast probability for day 154.

```{r}
#| label: sv-day153-update
#| code-fold: true
#| code-summary: "Update day 153 and calculate both forecasts for day 154"
fit152 <- sv_at(152)
prior153 <- fit152$next_weight[, 152]
density153 <- dnorm(sv_x[153], mean = fit152$delta, sd = exp(fit152$h / 2))
updated153 <- prior153 * density153
updated153 <- updated153 / sum(updated153)

weights153 <- cbind(prior153, updated153)
weights154 <- fit152$transition %*% weights153
large_move <- 2 * pnorm(-exp(-fit152$h / 2))
day154_prob <- drop(crossprod(large_move, weights154))
knitr::kable(data.frame(
  "Returns used" = c("Through day 152", "Through day 153"),
  "Estimated day-153 volatility" = round(drop(crossprod(exp(fit152$h / 2), weights153)), 2),
  "Chance of a large move on day 154" = sprintf("%.1f%%", 100 * day154_prob),
  check.names = FALSE), row.names = FALSE)
```

The day-153 return raises our estimate of that day’s volatility from $0.66$ to $1.30$. With $\hat\phi=0.95$, the updated distribution also puts more weight on high volatility the following day. The forecast probability of a large move on day 154 rises from $13.5\%$ to $41.3\%$.

This update takes us *half a step towards being Bayesian*. We used Bayes' rule to revise the distribution of $h_{153}$ and held $\phi$ at the value selected by likelihood. To account for uncertainty about persistence as in Workbook 4, we could give $\phi=0$ and $\phi=0.95$ prior probabilities before seeing the returns. Each prior probability times its model's likelihood through day 153 gives a weight. Dividing by the sum gives the posterior probabilities we would use to average the two day-154 forecasts.

## Do large moves occur as often as forecast? {#how-often-do-the-forecast-events-occur}

The model gave a $34.1\%$ chance that the return on day 156 would differ from $0.02$ by more than one unit. The return was $-1.26$, so the event occurred. Both persistence candidates gave valid stationary return processes, and the first 100 returns favored $\hat\phi=0.95$. How well do its forecast probabilities describe later outcomes? We check *calibration* by asking how often a large move occurred among days assigned similar probabilities.

Keep $\hat\phi=0.95$, fitted to the first 100 returns. As in Workbook 14, we advance the forecast origin one day at a time. At the end of each day $t=100,\ldots,399$ we forecast
$$
p_t=\P_{\hat\phi}\bigl(|X_{t+1}-0.02|>1\mid x_{1:t}\bigr).
$$
We pair $p_t$ with a $1$ if the next day's return falls outside $[-0.98,1.02]$ and a $0$ otherwise. Within each range of forecast probabilities, the table compares the mean of the forecasts with the fraction of returns that were large moves.

```{r}
#| label: sv-event-frequencies
#| code-fold: true
#| code-summary: "Pair each forecast with the next return and group the results"
origins <- 100:399
sv_rolling <- sv_at(max(origins))
# Entry t uses only the returns through day t.
prob <- sv_tail(sv_rolling, threshold = 1)[origins]
event <- abs(sv_x[origins + 1L] - 0.02) > 1
group <- cut(prob, c(0, 0.1, 0.2, 1), right = FALSE, include.lowest = TRUE,
             labels = c("Below 10%", "10% to below 20%", "20% or more"))
checks <- do.call(rbind, lapply(split(seq_along(prob), group), function(i) {
  c(length(i), sum(event[i]), 100 * mean(prob[i]), 100 * mean(event[i]))
}))
knitr::kable(data.frame(group = rownames(checks), checks, row.names = NULL),
             digits = 1, row.names = FALSE,
             col.names = c("Forecast range", "Forecasts", "Large moves",
                           "Mean forecast (%)", "Observed frequency (%)"))
```

Large moves were more frequent in each successive forecast group: $3.4\%$, $17.3\%$, and $23.0\%$ of days. Within the middle group, the observed frequency of $17.3\%$ exceeded the average forecast of $14.7\%$. In the lowest and highest groups the average forecasts, $5.6\%$ and $26.7\%$, exceeded the observed frequencies.

## Prediction intervals for JPMorgan Chase returns {#jpm-return-intervals}

In the simulation we fixed a threshold and forecast the probability that tomorrow's return would fall more than one unit from its mean. For JPMorgan Chase stock, let us reverse the question: at today's close, how wide an interval should we give for tomorrow's return if we want a $5\%$ chance of a miss on each side? Write $C_t$ for the closing price on trading day $t$, adjusted for splits and dividends. The return we forecast is $X_{t+1}=100\log(C_{t+1}/C_t)$. A return of $1$ is a rise of approximately $1\%$. We use [JPM daily adjusted closes from Yahoo Finance](https://finance.yahoo.com/quote/JPM/history/), with 2010–2019 for fitting and 2020–2023 for checking the forecasts.

First fit a constant-variance model: daily returns are independent $N(\delta,\sigma_X^2)$. Estimate $\delta$ by the training mean $\hat\delta$ and $\sigma_X^2$ by the average squared departure from that mean, $\hat\sigma_X^2$. Its central $90\%$ prediction interval is the same every day: $\hat\delta\pm 1.645\hat\sigma_X$. Now use the stochastic volatility model of Example 15.1. Given tomorrow's log variance $h_{t+1}$, the return has mean $\hat\delta$ and variance $e^{h_{t+1}}$. Recent large returns can make high values of $h_{t+1}$ more likely, even though the forecast mean stays the same.

We give both models the same variance for a single day's return before conditioning on the recent returns. If the stationary log variance has distribution $N(\mu,v)$, then $\Var(X_t)=\E[e^{h_t}]=e^{\mu+v/2}$. For each candidate pair $(\phi,v)$, we therefore set $\mu=\log(\hat\sigma_X^2)-v/2$ and $\sigma_\eta^2=v(1-\phi^2)$, starting the log variance in its stationary law. We compare five choices of persistence $\phi$ and three choices of the stationary variance $v$ of $h_t$. The training likelihood is largest at $\hat\phi=0.95$ and $\hat v=1$. We keep these estimates fixed and update the distribution of the log variance after each new return.

```{r}
#| label: jpm-training-fit
#| code-fold: true
#| code-summary: "Read the adjusted closes and compare the 15 training likelihoods"
source("R/volatility.R")
source("R/return-intervals.R")
jpm <- read.csv("data/jpm-adjusted-close.csv")
day <- as.Date(jpm$date)[-1L]
x_jpm <- 100 * diff(log(jpm$adjusted_close))
train <- day < as.Date("2020-01-01")
jpm_fit <- sv_fit_return_scale(x_jpm[train],
  phi = c(0.80, 0.90, 0.95, 0.975, 0.99), v = c(0.25, 1, 2.25))
```

For each trading day in 2020–2023, we use returns only through the preceding close. At each possible log variance $s$, tomorrow's return has distribution $N(\hat\delta,e^s)$. We average these distributions with the filter's predicted weights, then take the $5$th and $95$th percentiles of the resulting return distribution as our endpoints.

```{r}
#| label: jpm-next-intervals
#| code-fold: true
#| code-summary: "Find the endpoints for each trading day in 2020–2023"
best <- jpm_fit$best
jpm_filter <- sv_filter(x_jpm, jpm_fit$delta, best$mu, best$phi,
                        best$sigma_eta, nodes = 161L, span = 6)
target <- which(!train)
# For target t+1 use the weights from the close of trading day t.
half <- sv_interval_halfwidth(jpm_filter, origin = target - 1L)
fixed <- qnorm(.95) * sqrt(jpm_fit$variance)
bands <- data.frame(day = day[target], return = x_jpm[target],
  lower = jpm_fit$delta - half, upper = jpm_fit$delta + half,
  fixed_lower = jpm_fit$delta - fixed, fixed_upper = jpm_fit$delta + fixed)
```

```{r}
#| label: fig-jpm-return-intervals
#| echo: false
#| fig-height: 3.8
#| fig-alt: "Daily JPM returns with two 90% prediction intervals. The constant limits stay near minus and plus 2.6 percent. The volatility interval is narrower in January, widens after large returns in March and contracts in April."
#| fig-cap: "JPMorgan Chase daily log returns, January–May 2020. Each dot is a realized return; the shaded interval and dashed limits were calculated before that return was observed."
jpm_interval_plot(subset(bands, day < as.Date("2020-06-01")))
```

For February 10 the volatility model gave an interval from $-2.16$ to $2.27$, narrower than the constant limits of $-2.54$ and $2.65$. The first sharp declines fell outside the shaded band. After several large moves in both directions, the interval for March 19 extended from $-12.09$ to $12.20$. Both models still forecast the same mean return of about $0.06$.

A wider interval can cover more returns simply by making a less precise forecast. We compare *coverage*, the fraction of observed returns inside their intervals, together with width and the *interval score*. For a central $90\%$ interval, the score is its width on a day the return falls inside. On a miss we add $2/0.10=20$ times the distance from the return to the nearer endpoint. For example, an interval from $-2$ to $2$ scores $4$ when the return is $1$, and $4+20(3-2)=24$ when it is $3$. Lower scores are better.

```{r}
#| label: jpm-interval-scores
#| code-fold: true
#| code-summary: "Check coverage, width, and interval score over all 1,006 forecasts"
interval_check <- function(lower, upper, observed) {
  width <- upper - lower
  miss <- pmax(lower - observed, observed - upper, 0)
  c(100 * mean(miss == 0), mean(width), mean(width + 20 * miss))
}
jpm_checks <- with(bands, rbind(
  "Constant variance" = interval_check(fixed_lower, fixed_upper, return),
  "Stochastic volatility" = interval_check(lower, upper, return)))
knitr::kable(data.frame(Model = rownames(jpm_checks), unname(jpm_checks)),
  col.names = c("Model", "Coverage (%)", "Mean width (percentage points)", "Mean interval score"),
  row.names = FALSE, digits = c(0, 1, 2, 2))
```

Across all 1,006 forecasts, both intervals fell short of $90\%$ coverage. The changing interval covered $88.2\%$ of returns, compared with $85.6\%$ for the fixed interval. It was wider on average, but the reduction in the penalty for misses more than paid for the extra width: the mean score fell from $10.19$ to $8.21$.

## Gaussian dynamic linear models {#forecasting-with-a-hidden-level-and-slope}

Suppose a gauge measures water height once a day, in centimeters relative to a fixed reference mark. We want to forecast its next eight readings. The straight-line regression from Workbook 8 gives one possible model,
$$
X_t=a+bt+Y_t,
$$
where $X_t$ is the gauge reading and $Y_t$ is a mean-zero, weakly stationary deviation from the line. The mean reading changes by the same amount $b$ each day throughout the record. We estimate the fixed coefficients $a$ and $b$ from the readings.

To allow the rate of rise or fall to vary, let $L_t$ be the actual water height on day $t$ and let $B_t=L_{t+1}-L_t$ be its change over the next day. For example, $B_t=2$ means the height rises by $2$ cm from day $t$ to day $t+1$. We call $B_t$ the *slope*, measured in centimeters per day. The gauge measures the height with error,
$$
X_t=L_t+\varepsilon_t.
$$
Only $X_t$ is recorded. To forecast the next reading, we need to infer both the current height $L_t$ and the daily change $B_t$. We take $S_t=(L_t,B_t)^\top$ as the state, with state space $\R^2$.

**Definition 15.4 (Gaussian dynamic linear model).** A *Gaussian dynamic linear model* (DLM) for $(X_t)_{t\in I}$, with $I$ an interval of integers, is specified by a family of state-space representations with $S_t\in\R^m$ and $X_t\in\R^d$. In each representation, on every finite window $r,\ldots,u$ in $I$, the state $S_r$ is Gaussian and
$$
\begin{aligned}
S_t &= G_t S_{t-1}+Z_t, & Z_t &\sim N_m(0,Q_t), && r<t\leq u,\\
X_t &= F_t S_t+\varepsilon_t, & \varepsilon_t &\sim N_d(0,H_t), && r\leq t\leq u.
\end{aligned}
$$
The vectors $S_r,Z_{r+1},\ldots,Z_u,\varepsilon_r,\ldots,\varepsilon_u$ are mutually independent. The $m\times m$ transition matrices $G_t$ and $d\times m$ observation matrices $F_t$ are deterministic. The state covariance $\Var(S_r)$ and the disturbance covariances $Q_t$ may be singular. We require $H_t$ to be positive definite.

For a record beginning at $t=1$, write $S_1\sim N_m(m_1,P_1)$. The vector $m_1$ and matrix $P_1$ are the state mean and covariance before the first reading.

Given $S_{t-1}=s$, the state $S_t$ has mean $G_ts$ and covariance $Q_t$. Given $S_t=s$, the observation $X_t$ has mean $F_ts$ and covariance $H_t$.

**Example 15.5 (Water height with a changing slope).** We let the daily change $B_t$ follow a random walk. Its increments $U_t=B_t-B_{t-1}$ are independent $N(0,q)$ variables, with one unknown variance $q\geq0$ shared by all days. We fix the gauge's measurement-error standard deviation at $1$ cm. The model equations are
$$
\begin{aligned}
L_t &= L_{t-1}+B_{t-1},\\
B_t &= B_{t-1}+U_t, & U_t &\sim N(0,q),\\
X_t &= L_t+\varepsilon_t, & \varepsilon_t &\sim N(0,1).
\end{aligned}
$$
We assume the disturbances and initial state are mutually independent, as in Definition 15.4.

Suppose yesterday's actual height was $10$ cm and it was rising at $2$ cm per day. Today's height is then $12$ cm. If today's disturbance is $U_t=0.5$, the new slope is $2.5$ cm per day, giving a height of $14.5$ cm tomorrow. The gauge readings fluctuate around these actual heights because of measurement error.

For $S_t=(L_t,B_t)^\top$, the equations have the DLM matrices
$$
\begin{aligned}
G&=\begin{pmatrix}1&1\\0&1\end{pmatrix},
& Q&=\begin{pmatrix}0&0\\0&q\end{pmatrix},\\[.5em]
F&=\begin{pmatrix}1&0\end{pmatrix},
& H&=1.
\end{aligned}
$$
The first row of $G$ adds the previous height and slope. The second carries the previous slope forward before adding $U_t$. The observation row $F$ selects the height, and $H=1$ gives the gauge-error variance. The zeros in the first row and column of $Q$ express that the only new state disturbance is the increment to the slope.

We start the simulation with independent initial height and slope,
$$
L_1\sim N(0,4),\qquad B_1\sim N(0,0.25).
$$
Thus $m_1=(0,0)^\top$ and $P_1=\operatorname{diag}(4,0.25)$. We use this same initial law for every candidate $q$. For a fixed $q$, drawing the initial state and the independent disturbances determines all successive heights, slopes, and readings. This gives a unique joint law on the observed window. The heights and slopes are random quantities we infer from the gauge readings; $q$ is one fixed parameter we estimate from the whole fitting period.

If $q=0$, then $U_t=0$ and $B_t=B_1$ throughout the record. Given the initial height and slope, the readings follow a straight line with measurement error,
$$
X_t=L_1+(t-1)B_1+\varepsilon_t=(L_1-B_1)+tB_1+\varepsilon_t.
$$
This is the regression form $a+bt+\varepsilon_t$ with $a=L_1-B_1$ and $b=B_1$. Here the line varies across simulated records because we draw a new initial height and slope for each record. Positive $q$ allows the slope to change within a record. Exercise 15.3 works out how $q$ affects differences between successive gauge readings. $\diamond$

[**Example 15.6 (The Kalman filter for the gauge).**]{#kalman-filter-for-the-instrument} Fix a candidate value of $q$ in Example 15.5. The *Kalman filter* uses each new reading to update the Gaussian distribution of the water height and slope and to forecast the next reading. Consider the first simulated reading, $x_1\approx0.626$. Before it arrived, the height had mean zero and variance four. Since $X_1=L_1+\varepsilon_1$ adds independent error of variance one, the forecast for $X_1$ had mean zero and variance five.

The forecast error $x_1-0$ is the first *innovation*. How much of that error should we attribute to the water height? The covariance of the height $L_1$ with $X_1$ is four. The initial slope $B_1$ is independent of $X_1$, so its covariance with $X_1$ is zero. By Gaussian conditioning (Workbook 3, Proposition 3.5), the regression coefficients for $S_1=(L_1,B_1)^\top$ on $X_1$ give the *Kalman gain*,
$$
K_1=\frac{\Cov(S_1,X_1)}{\Var(X_1)}
   =\frac{1}{5}\begin{pmatrix}4\\0\end{pmatrix}
   =\begin{pmatrix}0.8\\0\end{pmatrix}.
$$
To update the state mean, add $K_1$ times the innovation to the initial mean $m_1=(0,0)^\top$. The estimated height rises to $0.8x_1\approx0.501$. The estimated slope remains zero.

The conditional variance of the height falls to $4-4^2/5=0.8$. The initial slope retains variance $0.25$ and remains independent of the height given $X_1=x_1$. Since the second reading is $X_2=L_1+B_1+\varepsilon_2$, its Gaussian forecast has mean $0.8x_1\approx0.501$ and variance $0.8+0.25+1=2.05$. The second reading was $x_2\approx0.943$, giving an innovation of about $0.442$. This reading carries information about the slope through the change in water height from day 1 to day 2.

For a later day $t$, suppose we have the Gaussian distribution of $S_t$ given the readings $x_{1:t}$. The state equation gives the mean vector $u$ and covariance matrix $V$ of tomorrow's state given those readings,
$$
\begin{aligned}
u &= G_{t+1}\E[S_t\mid x_{1:t}],\\
V &= G_{t+1}\Var(S_t\mid x_{1:t})G_{t+1}^{\top}+Q_{t+1}.
\end{aligned}
$$
The observation equation then gives the forecast for tomorrow's reading,
$$
X_{t+1}\mid x_{1:t}\sim N_d\!\left(F_{t+1}u,\ F_{t+1}VF_{t+1}^{\top}+H_{t+1}\right).
$$
The term $F_{t+1}VF_{t+1}^{\top}$ is the covariance of $F_{t+1}S_{t+1}$ given the readings through $t$. Measurement error adds $H_{t+1}$. Given those same readings, the next state and the next reading are jointly Gaussian. When $x_{t+1}$ arrives, we condition this joint distribution on $X_{t+1}=x_{t+1}$ to update the mean and covariance of $S_{t+1}$. This updated distribution supplies the state mean and covariance for the next prediction.

The first slope increment $U_2$ enters the height on day 3 through $L_3=L_2+B_2$. Its variance $q$ therefore first affects the forecast for $X_3$. We fit $q$ from the first 80 readings by running the filter separately at $q=0,0.01,\ldots,0.30$. Within each run, $q$ stays fixed and the filter revises the distribution of the height and slope as each reading arrives. Comparing runs, we select the $q$ with the largest log likelihood, the sum of the 80 log forecast densities evaluated at the recorded readings. The accompanying [DLM code](R/dlm.R) generates the readings and carries out the filter.

```{r}
#| label: dlm-fit
source("R/dlm.R")
trend <- trend_sample()
x_train <- trend$x[1:80]
q_grid <- seq(0, 0.30, by = 0.01)
loglik <- vapply(q_grid, function(q) trend_filter(x_train, q)$loglik, numeric(1))
q_hat <- q_grid[which.max(loglik)]
trend_fit <- trend_filter(x_train, q_hat)
round(c(q_hat = q_hat, trend_fit$m), 3)
```

The straight trend ($q=0$) gives a log likelihood of $-1448.81$. Allowing the slope to change raises the log likelihood to $-135.75$ at $\hat q=0.05$, the maximum on the grid. The readings give nearly the same log likelihood, $-135.84$, at the value $q=0.04$ used in the simulation.

With $\hat q=0.05$, the estimated height on day 80 is $-55.72$ cm and the estimated slope is $0.37$ cm per day. Adding them gives a forecast mean of $-55.34$ cm for the day-81 gauge reading, using unrounded values. The predicted height has variance $0.964$ square centimeters given the first 80 readings. The reading adds independent measurement error of variance one, so its forecast standard deviation is $\sqrt{0.964+1}\approx1.40$. With $q$ fixed at its fitted value, a 95% prediction interval extends about $1.96$ forecast standard deviations on either side of the mean.

```{r}
#| label: dlm-forecast
dlm_future <- trend_forecast(trend_fit, horizons = 1:8)
knitr::kable(dlm_future[c(1, 8), c("h", "mean", "lower", "upper")],
             row.names = FALSE, digits = 2,
             col.names = c("Days ahead", "Mean", "95% lower", "95% upper"))
```

All eight forecasts in the figure use the same first 80 readings. The forecast mean for day $80+h$ adds $h$ times the day-80 mean slope to the day-80 mean height because each future slope increment has mean zero.

```{r}
#| label: fig-dlm-look-ahead
#| echo: false
#| fig-alt: "Daily readings approach a low point before day 80. The fitted forecast rises gradually over the next eight days and its shaded interval widens."
#| fig-height: 3.5
#| fig-cap: "Simulated daily gauge readings, in centimeters relative to the reference mark. The forecasts use the observations through day 80 (vertical line). Shading gives a 95% prediction interval for each later reading with the parameter $q$ fixed at its fitted value. Open circles show the readings on days 81–88."
trend_plot(trend, origin = 80, forecasts = dlm_future)
```

The day-80 prediction interval for the reading eight days later is $[-61.70,-43.78]$. Uncertainty about the slope on day 80 affects all eight daily changes in water height. A new slope increment on day 81 first affects the height on day 82 and enters every subsequent change in height through day 88.

The first open circle in the figure is the day-81 reading, $x_{81}=-56.35$. It falls below the forecast of $-55.34$. How does it change our estimate of the day-81 height and slope? With $\hat q=0.05$ held fixed, the two conditional means of $S_{81}=(L_{81},B_{81})^\top$ are
$$
\begin{aligned}
m_{81\mid80} &= \E[S_{81}\mid x_{1:80}]\approx(-55.34,\;0.37)^\top,\\
m_{81\mid81} &= \E[S_{81}\mid x_{1:81}]\approx(-55.84,\;0.21)^\top.
\end{aligned}
$$
In both rows the first index is 81 because we are estimating the height and slope on day 81. The index after the bar gives the last day whose reading enters the estimate, day 80 in the first row and day 81 in the second.

Why does a reading of the height revise the slope? Both $L_{81}=L_{80}+B_{80}$ and $B_{81}=B_{80}+U_{81}$ depend on the previous slope. Given the first 80 readings, their covariance is $0.313$. The gauge adds independent error to the height, so $X_{81}$ has the same covariance with $B_{81}$ and has variance $0.964+1=1.964$. The slope component of the day-81 Kalman gain is therefore $0.313/1.964$. Multiplying by the innovation gives the change in the estimated slope,
$$
\frac{\Cov(B_{81},X_{81}\mid x_{1:80})}{\Var(X_{81}\mid x_{1:80})}
\bigl(x_{81}-\E[X_{81}\mid x_{1:80}]\bigr)
\approx \frac{0.313}{1.964}(-1.01)\approx -0.16.
$$
This accounts for the drop from $0.37$ to $0.21$ cm per day.

How much uncertainty about the day-81 height remains? Write $P_{81\mid80}$ and $P_{81\mid81}$ for the conditional covariance matrices of $S_{81}$ before and after the reading. Their first diagonal entries give the conditional variance of $L_{81}$. The reading reduces that variance from $0.96$ to $0.49$ square centimeters. $\diamond$

## Updating a forecast of Nile flow {#nile-flow}

In 1966 the Nile carried considerably less water past Aswan than in 1965. How much should that low reading lower our forecast for 1967? The observed variable $X_t$ is the volume of water that passes Aswan in year $t$, measured in hundreds of millions of cubic metres. We model $L_t$, the underlying mean flow for that year, as a hidden state. Part of a low reading can reflect a fall in $L_t$ that persists into the next year. Part can be a fluctuation specific to that year.

We use a DLM with one state coordinate. With the independence assumptions of Definition 15.4, its equations are
$$
\begin{aligned}
L_t &= L_{t-1}+Z_t, & Z_t &\sim N(0,q),\\
X_t &= L_t+\varepsilon_t, & \varepsilon_t &\sim N(0,r).
\end{aligned}
$$
The two parameters to fit are $q$, the variance of a yearly change in mean flow, and $r$, the variance of annual flow around that mean. We fit $q$ and $r$ by likelihood using the 95 flows recorded in 1871–1965. R starts the filter at the first recorded flow with a large initial variance for the level. The fitted standard deviations are $\sqrt{\hat q}\approx35$ and $\sqrt{\hat r}\approx124$ flow units. In this fit, the river's flow in a single year can vary considerably more than its underlying mean changes from one year to the next.

```{r}
#| label: nile-fit
#| code-fold: true
#| code-summary: "Fit the two variances and forecast 1966"
nile_train <- window(datasets::Nile, end = 1965)
nile_fit <- StructTS(nile_train, type = "level")
nile_before <- predict(nile_fit, n.ahead = 1)
nile_mean <- as.numeric(nile_before$pred)
nile_flow_var <- as.numeric(nile_before$se)^2
nile_level_var <- nile_flow_var - unname(nile_fit$coef["epsilon"])
nile_gain <- nile_level_var / nile_flow_var
```

Before seeing the 1966 reading, we forecast a flow of $960.2$. This is also our estimate of the 1966 level, since the annual fluctuation has mean zero. The recorded flow was $746$. The fitted Kalman gain is $0.243$, the predicted variance of the hidden level divided by the predicted variance of the flow. We therefore attribute about a quarter of the shortfall to a lower level. Its updated estimate is
$$
m_{1966\mid1966}
\approx 960.2+0.243(746-960.2)\approx 908.
$$

The state equation carries the updated level forward with a mean-zero increment. The observation equation adds another mean-zero fluctuation. The forecast for the **observed 1967 flow** is therefore also $908$, down from $960.2$ before we saw the low reading. The 1967 flow was $919$.

```{r}
#| label: fig-nile-next-flow
#| code-fold: true
#| code-summary: "Use the 1966 flow to update the forecast for 1967"
#| fig-height: 3.3
#| fig-alt: "Annual Nile flows with forecasts for the last two years. The 1966 forecast is 960 and the recorded flow is 746. The 1967 forecast is 908 and the recorded flow is 919."
#| fig-cap: "Annual flow at Aswan from R’s Nile series. Each triangle gives a forecast made using all recorded flows through the preceding year, with both variances fitted through 1965. The example follows [Raquel Prado’s local-level analysis](https://www2.stat.duke.edu/~mwest/TSFCourseSoftware/)."
nile_1966 <- as.numeric(window(datasets::Nile, 1966, 1966))
nile_after <- nile_mean + nile_gain * (nile_1966 - nile_mean)
nile_flow_plot(window(datasets::Nile, 1957, 1967), c(nile_mean, nile_after))
```

## Finite-state hidden Markov models {#when-the-hidden-state-is-a-setting}

A sensor counted seven flashes from the lamp in the first second, raising the probability of the medium setting from $1/3$ to about $0.75$. To forecast the next count, we need to know how often the lamp stays on the same setting. Can we estimate that probability from the sequence of flash counts?

**Definition 15.7 (Finite-state hidden Markov model).** A *finite-state hidden Markov model* for $(X_t)_{t\in I}$ is specified by a family of state-space representations (Definition 15.3) whose states $S_t$ take values in a common finite set $E$. With $K$ states we can label its elements $1,\ldots,K$.

**Example 15.8 (Will the lamp flash at least ten times in the next second?).** Let $S_t$ be the random setting in second $t$, labeled $1,2,3$ for low, medium, and high. The sensor records the number of flashes $X_t$ in that second. We fix the initial setting probabilities at $1/3$ each and the Poisson means at $\lambda=(2,6,12)$ flashes per second. The unknown parameter is the probability $p\in[0,1]$ of staying in the current setting. Otherwise the lamp advances one place in the cycle low, medium, high, low. The transition matrix $A(p)$ and the probability $g_j(x)$ of count $x$ in setting $j$ are
$$
A(p)=\begin{pmatrix}
p&1-p&0\\
0&p&1-p\\
1-p&0&p
\end{pmatrix},
\qquad g_j(x)=e^{-\lambda_j}\frac{\lambda_j^x}{x!},\quad x=0,1,\ldots.
$$
For every $p\in[0,1]$, the entries in each row of $A(p)$ are nonnegative and sum to one. Starting from $S_1$, row $i$ gives the probabilities $a_{ij}$ for drawing the next setting $j$ from setting $i$. Drawing the counts independently given the settings, with probabilities $g_j(x)$, completes the joint law of the setting and count sequences on the record.

Can we distinguish $p=0.50,0.80,0.95$ from the counts? Suppose each setting has probability $1/3$ before a transition. After the transition, each setting again has probability $(1/3)p+(1/3)(1-p)=1/3$ because the lamp either stayed there or advanced from its predecessor. Thus a count from any one second has the same mixture of three Poisson laws for every candidate $p$. Consecutive counts can distinguish the candidates. A count of at least ten favors the high setting. The larger $p$ is, the more likely the lamp is to stay high and produce another count of at least ten.

Fix a candidate $p$. We repeat the update we used for the first count of seven, now allowing the setting to change between counts. Write $w_{t\mid r}(j)=\P_p(S_t=j\mid x_{1:r})$ for the probability of setting $j$ in second $t$ given counts through second $r$. We start with $w_{1\mid0}(j)=1/3$. The probability of $x$ flashes in second $t$ averages the count probabilities over the three possible settings,
$$
\P_p(X_t=x\mid x_{1:t-1})=\sum_{j=1}^3 w_{t\mid t-1}(j)g_j(x).
$$
After observing $x_t$, Bayes' rule gives the share of that probability contributed by setting $j$,
$$
w_{t\mid t}(j)=\frac{w_{t\mid t-1}(j)g_j(x_t)}{\P_p(X_t=x_t\mid x_{1:t-1})}.
$$
To predict setting $j$ for the next second, add the contributions from the possible current settings $i$. This gives $w_{t+1\mid t}(j)=\sum_i w_{t\mid t}(i)a_{ij}$. These probabilities weight the Poisson laws in the forecast for the next count.

We simulate 81 seconds at $p=0.80$ with seed 1507 and use the first 80 counts to fit $p$. For each candidate value of $p$, the [count model code](R/hmm.R) uses the initial setting probabilities for the first count and conditions on the preceding counts thereafter. Summing the logarithms of the 80 count probabilities gives the log likelihood for that candidate.

```{r}
#| label: hmm-persistence-fit
source("R/hmm.R")
flashes <- wb15_count_record()
rate <- c(2, 6, 12)
stay <- c(0.50, 0.80, 0.95)
lamp_fits <- lapply(stay, function(p) wb15_count_fit(flashes[1:80], p, rate))
lamp_ll <- vapply(lamp_fits, function(fit) fit$loglik, numeric(1))
winner <- which.max(lamp_ll)
knitr::kable(data.frame(stay, lamp_ll), digits = 2, row.names = FALSE,
             col.names = c("Probability of staying", "Log likelihood"))
```

The largest log likelihood gives $\hat p=0.80$ among the three candidates. We observed 16 flashes in second 79 and 11 in second 80. Given all 80 counts, the fitted probability that the lamp was in the high setting in second 80 is 0.986. If it was high, the lamp has a 0.20 chance of switching to low for second 81. Accounting for transitions from all three possible settings, the predicted probabilities of low, medium and high in second 81 are $(0.197,0.011,0.792)$.

```{r}
#| label: hmm-count-forecast
lamp_fit <- lamp_fits[[winner]]
next_w <- lamp_fit$next_state
setting_tail <- ppois(9, rate, lower.tail = FALSE)
forecast_tail <- vapply(lamp_fits, function(fit) {
  sum(fit$next_state * setting_tail)
}, numeric(1))
```

```{r}
#| label: fig-hmm-count-forecast
#| fig-alt: "The flash counts fluctuate between stretches of low and high activity. The next-count forecast puts mass near both two and twelve flashes."
#| echo: false
#| fig-width: 9
#| fig-height: 3.35
#| fig-cap: "Simulated flash counts through second 80 (left) and the fitted probabilities for the count in second 81 (right). Bars at ten flashes and above are shaded in brown."
wb15_plot_count_forecast(flashes[1:80], next_w, rate)
```

The peaks near two and twelve flashes come from the low and high settings. If the lamp is in its high setting, the probability of at least ten flashes is $0.758$ under the Poisson law with mean 12. To forecast the same event from the first 80 counts, average the three Poisson tail probabilities, weighting each by the predicted probability of its setting in second 81,
$$
\P_{\hat p}(X_{81}\ge10\mid x_{1:80})
=\sum_{j=1}^3 w_{81\mid80}(j)\,\P(X_{81}\ge10\mid S_{81}=j)
\approx0.601.
$$

Using the same 80 counts with $p=0.50$ gives a probability of `r sprintf("%.3f", forecast_tail[1])` for at least ten flashes in second 81. With $p=0.95$ it is `r sprintf("%.3f", forecast_tail[3])`. The higher probability of staying puts more weight on the high setting for the next second.

The sensor counted five flashes in second 81. In Exercise 15.5 we use those five flashes to update the setting probabilities and forecast the probability of at least ten flashes in second 82. ◊

## Waiting for the next eruption of Old Faithful {#old-faithful-waits}

An eruption of Old Faithful has just begun. How long until the next one starts? Azzalini and Bowman (1990) report 299 consecutive waits at Yellowstone in August 1985. Write $X_t$ for the time in minutes between the starts of eruptions $t-1$ and $t$. At the start of eruption $t$ we know $x_t$ and want to forecast $X_{t+1}$. The waiting times are available in [`MASS::geyser`](https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/geyser.html).

We model each wait as coming from one of two hidden states. In state $S_t=j$, the waiting time $X_t$ has a normal distribution with mean $\mu_j$ and standard deviation $\sigma_j$. State 1 has the shorter mean wait and state 2 the longer mean wait. The states follow a Markov chain, and given the states the waits are independent.

We use the first 200 waits to fit the two means, the two standard deviations, and the two probabilities of moving to the long state, one from each current state. For each candidate set of parameters we start the chain in its stationary distribution and compute the likelihood from the successive forecast densities. We maximize this likelihood numerically, keeping each standard deviation at or above two minutes so a state cannot concentrate its density arbitrarily tightly around a single recorded wait.

```{r}
#| label: geyser-training-fit
#| code-fold: true
#| code-summary: "Fit the two states using only the first 200 waits"
source("R/hmm.R")
source("R/geyser.R")
wait <- MASS::geyser$waiting
geyser_fit <- wb15_fit_waits(wait[1:200])
knitr::kable(data.frame(
  state = c("Short", "Long"),
  mean = round(geyser_fit$means, 1),
  sd = round(geyser_fit$sds, 1),
  next_long = ifelse(geyser_fit$transition[, 2] > .999, "> 0.999",
                     sprintf("%.3f", geyser_fit$transition[, 2]))),
  col.names = c("Current state", "Mean wait (min)", "SD (min)",
                "Probability next state is long"), row.names = FALSE)
```

With the parameters fixed at their estimates, consider a forecast made at the start of eruption 208. The probability of the long state for wait 209 was $0.472$. Applying the transition matrix once more gave a probability of $0.618$ for the long state for wait 210. The mean forecast for wait 210 was therefore $0.382(59.1)+0.618(82.2)\approx73.4$ minutes.

Eruption 209 began after 69 minutes. Both states can produce that wait, but the fitted density at 69 is $0.0241$ per minute in the short state and only $0.00658$ in the long state. Bayes' rule gives the updated probability for the wait that just ended,
$$
w_{209\mid209}(2)
\approx \frac{0.472(0.00658)}{0.528(0.0241)+0.472(0.00658)}
\approx 0.196.
$$

The 69-minute wait makes the short state more likely for wait 209. Should we raise or lower our forecast of $73.4$ minutes for wait 210? In the fitted chain a short state almost always leads to a long state, while a long state leads to another long state with probability $0.192$. Weighting these transitions by the updated state probabilities gives $w_{210\mid209}(2)\approx0.8037(1.000)+0.1963(0.192)\approx0.841$. The mean forecast for wait 210 rises to $78.6$ minutes. The next eruption began 92 minutes later.

```{r}
#| label: geyser-next-wait
#| code-fold: true
#| results: hide
#| code-summary: "Forecast wait 210 before and after observing wait 209"
prior <- wb15_wait_filter(wait[1:208], geyser_fit)$next_state
before210 <- drop(prior %*% geyser_fit$transition)
weights <- prior * dnorm(wait[209], geyser_fit$means, geyser_fit$sds)
updated <- weights / sum(weights)
after210 <- drop(updated %*% geyser_fit$transition)
forecast210 <- rbind(through_208 = before210, through_209 = after210)
round(cbind(long_probability = forecast210[, 2],
            mean_minutes = drop(forecast210 %*% geyser_fit$means)), 3)
```

## Exercises

::: {.exercise}
**Exercise 15.1 (Can one day reveal persistence?).** In Example 15.1 the log variance follows $h_t-\mu=\phi(h_{t-1}-\mu)+\eta_t$ and the return is $X_t=\delta+e^{h_t/2}Z_t$. Compare $\phi=0$ with $\phi=0.95$, taking $\delta=0.02$ and $\mu=-1$ in both models. Take $v>0$ as given and start each process with $h_1\sim N(\mu,v)$. The return shocks $(Z_t)_{t\geq1}$ are standard normal and the state shocks $(\eta_t)_{t\geq2}$ are $N(0,v(1-\phi^2))$. These shocks and $h_1$ are mutually independent.

(a) Explain why $h_t$ has the same marginal distribution in both candidates and why the marginal distribution of $X_t$ therefore agrees. What happens to the joint law of the returns when $\phi=0$?

(b) Set $Y_t=\log((X_t-\delta)^2)$. Show that $\Cov(Y_t,Y_{t-1})=\phi v$. Which feature of the observed record can distinguish the models in part (a)? *Hint.* Square the return's deviation from $\delta$ and take logs to express $Y_t$ as the sum of $h_t$ and an error depending only on $Z_t$.

(c) At the end of days $t=140$ and $t=155$, forecast the probability that the next return differs from $0.02$ by more than one unit. Use the simulated returns and $v=0.3^2/(1-0.95^2)$ from [Fitting persistence and forecasting a large return](#fitting-persistence-and-forecasting-a-large-return). Run that section's R code to define `sv_at()`. Use `sv_at(t, phi = 0)` to filter the returns through day $t$ and take the last probability returned by `sv_tail()`. Compare these probabilities with the forecasts for days 141 and 156 under $\phi=0.95$. Under which model does the return history affect the forecast probability, and why?
:::


::: {.exercise}
**Exercise 15.2 (Does the latest position suffice?).** A point starts at
$L_0=0$ or $2$, each with probability $1/2$. Independently of $L_0$, draw a
step $V=-1$ or $+1$, each with probability $1/2$. Use this same step throughout the path. Thus $L_{t+1}=L_t+V$. Only the position
$X_t=L_t$ is recorded.

(a) Take $S_t=(L_t,V)$ with state space $\Z\times\{-1,+1\}$. The recorded position gives the first coordinate of this pair. Give the transition from a state $(l,v)$
and verify the Markov property.

(b) An observer sees only $X_1=1$. Calculate the probability that
$X_2=2$. How does that probability change if the observer also knows
that $X_0=0$? What if $X_0=2$?

(c) Is $(X_t)$ Markov? Express the next position in terms of the two
most recent positions, and explain what the earlier position tells us.
:::


::: {.exercise}
**Exercise 15.3 (What can differences of the readings reveal?).** In the water-height model of Example 15.5, the actual height, slope, and gauge reading satisfy
$$
L_t=L_{t-1}+B_{t-1},\qquad B_t=B_{t-1}+U_t,\qquad X_t=L_t+\varepsilon_t.
$$
The slope increments $U_t$ have variance $q$. Allow the measurement errors $\varepsilon_t$ to have an unknown variance $r>0$. All increments and errors are mutually independent. For $t\geq3$, write $D_t=X_t-2X_{t-1}+X_{t-2}$ for the second difference of the readings.

(a) Express $D_t$ in terms of one slope increment and three measurement errors. Which day's slope increment enters?

(b) Find $\Var(D_t)$ and $\Cov(D_t,D_{t-2})$ for $t\geq5$. *Hint: list the measurement errors shared by the two differences.*

(c) Suppose the population variance is $\Var(D_t)=6.4$ and the population lag-two covariance is $\Cov(D_t,D_{t-2})=1$. Recover $q$ and $r$. What information does the lag-two covariance add to the variance alone?
:::


::: {.exercise}
**Exercise 15.4 (Two forecasts for day 82).** Use `trend` and `q_hat` from the [model fit in Example 15.6](#kalman-filter-for-the-instrument), where $\hat q=0.05$. Forecast the reading on day 82 from the end of day 80 and again from the end of day 81. Keep $\hat q$ fixed for both forecasts.

(a) Use `trend_filter()` with the readings through day 80 and `trend_forecast()` to obtain the mean forecasts for days 81 and 82 and a 95% prediction interval for the day-82 reading.

(b) Observe the day-81 reading, `trend$x[81]`. Subtract the day-81 forecast mean in (a) from this reading to obtain the forecast error. Use the readings through day 81 to update the distribution of the height and slope. Give the new mean forecast and 95% prediction interval for day 82.

(c) Use $m_{81\mid80}$ and $m_{81\mid81}$ from Example 15.6. How much of the revision to the day-82 mean forecast comes from the change in the estimated height, and how much from the change in the estimated slope? Which source of uncertainty in $X_{82}$ remains even if the day-82 state were known?
:::


::: {.exercise}
**Exercise 15.5 (One more count).** In Example 15.8 the number of flashes in one second has a Poisson distribution with mean 2, 6, or 12 according to whether the lamp is in the low, medium, or high setting. Keep the fitted probability of staying in the same setting fixed at $0.80$. Given the first 80 counts, the probabilities of the low, medium, and high settings in second 81 are approximately $(0.197,0.011,0.792)$. The [forecast code in Example 15.8](#when-the-hidden-state-is-a-setting) saves their unrounded values in `next_w`. The observed count in second 81 is five.

(a) Calculate the probability that the count would be exactly five under each state. Which setting gives that count the largest probability?

(b) Multiply each probability in (a) by the corresponding setting probability in `next_w`. What is the probability of five flashes given the first 80 counts? Use the three products and their sum to find the setting probabilities given all 81 counts. Explain which setting now has the largest probability.

(c) Apply the transition matrix $A(0.80)$ to the updated setting probabilities to predict the setting in second 82. What is the probability of at least ten flashes in second 82? Compare it with the probability of at least ten flashes in second 81, as forecast from the first 80 counts.
:::


::: {.exercise}
**Exercise 15.6 (What do successive counts tell us about persistence?).**
In Example 15.8 the lamp starts in the low, medium, or high setting with equal probabilities. Each second it stays in its setting with probability $p$ or advances to the next setting in the cycle low, medium, high, low. Write $B_t=\{X_t\ge10\}$ for a count of at least ten flashes in second $t$, and $q_j=\P(B_t\mid S_t=j)$.

(a) Compute the three $q_j$ for Poisson counts with mean rates $(2,6,12)$. Show that $\P(B_1)$ has the same value for every $p$.

(b) Give a formula for $\P(B_1\cap B_2)$ by conditioning on the settings in seconds 1 and 2. Suppose all we observe in second 1 is that the lamp flashed at least ten times. Compute $\P(B_2\mid B_1)$ for $p=0.50$ and $p=0.95$. Explain how the probability of staying in the same setting affects this forecast.

(c) Suppose all three rates equal six. Show that the likelihood of any count record $x_1,\ldots,x_T$ is the same for every $p$. Give the conditional distribution of the next count $X_{T+1}$ given the counts through second $T$.
:::


## Further reading

Prado, Ferreira and West, *Time Series: Modeling, Computation, and Inference*
(2nd ed., 2021), Chapter 4 develops forecasting and filtering for dynamic
linear models. Sections 7.1 and 7.5 treat Markov switching and stochastic
volatility. The
[Stan User's Guide](https://mc-stan.org/docs/stan-users-guide/time-series.html#stochastic-volatility-models)
gives the Gaussian log-volatility specification, and the
[finite-state HMM reference](https://mc-stan.org/docs/functions-reference/hidden_markov_models.html)
describes the likelihood obtained by summing over states.