% ============================================================================
% Appendix C — Mathematical derivations
% Written by Agent L (Appendices B/C instance), 2026-11-03. Derivations and
% proofs referenced by Part IV (sections/part04_methodology.tex); notation
% follows Part IV. Each derivation is checked against the implementation in
% econometrics/structural_change.py, bibliometrics/03_citation_impact.py,
% bibliometrics/04_diversity_cagr.py, knowledge_graph/02_category_flows.py,
% taxonomy/classify.py, taxonomy/validate_sample.py and
% experiments/emerging_detection.py; where the code and the usual textbook
% convention differ, the code's convention is stated and the difference noted.
% ============================================================================

\section{Mathematical derivations}
\label{sec:appC}

Part~IV defines the statistics of this study and states their properties;
this appendix supplies the derivations behind those statements, in the order
in which Part~IV uses them. Each subsection fixes the object as implemented,
derives the property that the main text relies on, and records any point at
which the implementation departs from the convention of the cited source.
The derivations are elementary and self-contained; their purpose is not
novelty but the assurance that the quantities tabulated in Part~IX are what
the text says they are. Throughout, $n$ is the length of a monthly series,
$y_t = \log(N_t + 1)$ its logarithm, and $\mathbf{1}[\cdot]$ the indicator
function.

\subsection{The interrupted time series as a linear regression}
\label{sec:appC-its}

\paragraph{Design matrix.}
Equation~\eqref{eq:its} of Part~IV is the linear model $y = X\beta +
\varepsilon$ with $y \in \R^n$ and $X \in \R^{n \times 15}$ whose columns
are, for $t = 0, \dots, n-1$ and event index $t_e$,
\begin{equation}
  x_t \;=\; \bigl(1,\; t,\; \mathbf{1}[t \ge t_e],\; (t - t_e)\,\mathbf{1}[t \ge t_e],\;
  \mathbf{1}[\mathrm{month}(t) = 2],\; \dots,\; \mathbf{1}[\mathrm{month}(t) = 12]\bigr)^{\top},
  \label{eq:appC-design}
\end{equation}
exactly the data frame built by \texttt{\_design()} in
\texttt{econometrics/structural\_change.py} (columns \texttt{const},
\texttt{t}, \texttt{post}, \texttt{t\_post}, \texttt{m2}, \dots,
\texttt{m12}). The post-event trend is centred at $t_e$ rather than at zero
so that $\beta_2$ is the level change \emph{at} the event month and not the
extrapolated intercept difference at $t = 0$; with an uncentred post-trend
the two parametrisations are related by $\beta_2^{\mathrm{unc}} = \beta_2 -
\beta_3 t_e$, and only the centred one has the interpretation used in the
text. The eleven month dummies are linearly independent of the constant
because January is omitted; the full-column-rank condition requires at least
one observation in every month of the year, satisfied by any series of
twelve months or more, and at least one post-event month with $t > t_e$,
guaranteed by the driver's rule that twelve months follow every event.
The ordinary least-squares estimator is $\hat\beta = (X^{\top}X)^{-1}
X^{\top}y$ and the residuals $\hat\varepsilon = y - X\hat\beta$.

\paragraph{The HAC covariance estimator.}
Under the assumption that $\varepsilon_t$ is a stationary, serially
correlated and possibly heteroskedastic sequence with
$\E[x_t \varepsilon_t] = 0$, the asymptotic covariance of $\sqrt{n}(\hat\beta
- \beta)$ is $Q^{-1}\Omega Q^{-1}$ with $Q = \lim n^{-1}X^{\top}X$ and
$\Omega = \sum_{\ell=-\infty}^{\infty} \Gamma_\ell$, where
$\Gamma_\ell = \E[x_t \varepsilon_t \varepsilon_{t-\ell} x_{t-\ell}^{\top}]$
is the autocovariance of the score at lag $\ell$. The estimator of
\citet{newey1987hac} replaces the infinite sum by a kernel-weighted finite
one,
\begin{equation}
  \hat\Omega_L \;=\; \hat\Gamma_0 + \sum_{\ell=1}^{L} w_\ell\,\bigl(\hat\Gamma_\ell + \hat\Gamma_\ell^{\top}\bigr),
  \qquad
  \hat\Gamma_\ell = \frac{1}{n}\sum_{t=\ell}^{n-1} \hat\varepsilon_t \hat\varepsilon_{t-\ell}\, x_t x_{t-\ell}^{\top},
  \qquad
  w_\ell = 1 - \frac{\ell}{L+1},
  \label{eq:appC-nw}
\end{equation}
and the variance estimate is $\hat V = n\,(X^{\top}X)^{-1} \hat\Omega_L
(X^{\top}X)^{-1}$, which is the expression of Part~IV once the factors of
$n$ are absorbed into the un-normalised $\hat\Gamma_\ell$. The Bartlett
weights $w_\ell$ are what make $\hat\Omega_L$ positive semi-definite for every
sample: $\hat\Omega_L$ can be written as $\frac{1}{n(L+1)}\sum_{t}
\bigl(\sum_{s=t-L}^{t} \hat u_s\bigr)\bigl(\sum_{s=t-L}^{t} \hat u_s\bigr)^{\top}$
up to end effects, with $\hat u_s = x_s \hat\varepsilon_s$, a sum of outer
products, whereas the truncated sum with $w_\ell \equiv 1$ is not guaranteed
to be positive semi-definite. The implementation calls
\texttt{statsmodels} with \texttt{cov\_type="HAC"} and \texttt{maxlags=6}, so
$L = 6$ and $w_\ell = 1 - \ell/7$ as stated in Part~IV, without the
small-sample degrees-of-freedom correction and with $p$-values from the
normal approximation to the $t$-ratio, which is the library's default for
robust covariances. With $n$ of the order of one hundred and fourteen
regressors, the difference between the normal and the $t_{n-15}$ reference
distribution is immaterial at the magnitudes of the $t$-statistics reported;
what matters is the placebo rank of the next paragraph, which does not
depend on any reference distribution.

\paragraph{The placebo $p$-value as a rank statistic.}
Let $\mathcal{T}' = \{t' : 24 \le t' \le n - 13\}$ be the candidate event
indices (the implementation's \texttt{min\_pre = 24} and
\texttt{min\_post = 12}) and, for each $t' \in \mathcal{T}'$, let
$t^{(2)}_{t'} = \hat\beta_2(t') / \widehat{\mathrm{se}}(\hat\beta_2(t'))$ be
the HAC $t$-ratio of the level change when the model
\eqref{eq:appC-design} is fitted with $t_e = t'$. The placebo value of the
true event $t_e \in \mathcal{T}'$ is
\begin{equation}
  p^{\mathrm{plac}} \;=\; \frac{1}{|\mathcal{T}'|}\sum_{t' \in \mathcal{T}'}
  \mathbf{1}\bigl[\,|t^{(2)}_{t'}| \ge |t^{(2)}_{t_e}|\,\bigr]
  \;=\; \frac{\mathrm{rank}\bigl(|t^{(2)}_{t_e}|\bigr)}{|\mathcal{T}'|},
\end{equation}
where $\mathrm{rank}$ is the descending rank of the true event's absolute
$t$-ratio among the $|\mathcal{T}'|$ values, ties counted as not better
(\texttt{>=}). Three properties follow directly. Because $t_e \in
\mathcal{T}'$, the sum contains the term $t' = t_e$, so $p^{\mathrm{plac}}
\ge 1/|\mathcal{T}'|$, which for the default window ($n = 101$, hence
$|\mathcal{T}'| = 65$) is about $0.015$: a placebo rank cannot be smaller,
however large the true statistic, and a reported value of $1/|\mathcal{T}'|$
means ``largest of all candidate dates'', not ``significant at that level''.
Because the statistic is a rank, it is invariant to any monotone
transformation of the $|t|$-ratios and in particular does not depend on the
HAC bandwidth through anything but the ordering it induces. And if the
candidate dates were exchangeable under a null hypothesis of no break
anywhere, $p^{\mathrm{plac}}$ would be a valid permutation-type $p$-value
with $P(p^{\mathrm{plac}} \le \alpha) \le \alpha$; they are not exchangeable,
since all $|\mathcal{T}'|$ regressions share the same $n$ observations and
adjacent candidates give nearly identical fits, so the value is reported as
a rank and read as such. The same construction gives the placebo rank of
the slope change from $t^{(3)}_{t'}$. The driver additionally computes the
rank restricted to the baseline candidates $\{t' < t_e\}$, which answers the
protocol's original question (is the reference date special relative to the
pre-period?) and is reported next to the full-range rank.

\subsection{Bai--Perron multiple breaks by dynamic programming}
\label{sec:appC-bp}

\paragraph{Segment costs and additivity.}
For a segment $[i, j)$ of the series with $j - i \ge h$, define
\begin{equation}
  S(i, j) \;=\; \min_{a, b \in \R}\ \sum_{t=i}^{j-1} (y_t - a - b\,t)^2,
  \label{eq:appC-segcost}
\end{equation}
the residual sum of squares of a linear trend fitted to the segment alone
(the implementation solves the two-parameter least-squares problem for each
segment with \texttt{numpy.linalg.lstsq}; with \texttt{trend=False} the
regressor $t$ is dropped and $S(i,j)$ is the centred sum of squares). The
model with $m$ breaks $0 < \tau_1 < \dots < \tau_m < n$ assigns its own
$(a, b)$ to each of the $m+1$ segments $[\tau_{k-1}, \tau_k)$ with
$\tau_0 = 0$ and $\tau_{m+1} = n$. Because the parameters of different
segments are distinct and each observation enters exactly one segment, the
residual sum of squares of the whole model separates,
\begin{equation}
  \mathrm{SSR}(\tau_1, \dots, \tau_m) \;=\; \sum_{k=1}^{m+1} S(\tau_{k-1}, \tau_k),
  \label{eq:appC-additive}
\end{equation}
and the global minimisation over $\tau$ decomposes: the optimal partition of
the first $j$ observations into $m+1$ segments contains, as its first $m$
segments, an optimal partition of the first $\tau_m$ observations into $m$
segments. This is the principle of optimality, and it is the only property of
the cost that the dynamic programme uses; it holds for any cost that is a
sum over segments of a function of the segment alone, which is why the same
recursion serves the piecewise-constant model of Section~\ref{sec:appC-pelt}.
It fails for a model with parameters shared across segments (a common slope,
or the seasonal dummies of the ITS), which is why the segmented regression
carries no seasonal terms and the series is seasonally adjusted beforehand
by subtracting the baseline-period month means
(\texttt{run\_structural\_change.py}).

\paragraph{The recursion.}
Let $C_m(j)$ be the minimum of \eqref{eq:appC-additive} over partitions of
observations $0, \dots, j-1$ into $m+1$ admissible segments. Then
\begin{equation}
  C_0(j) = S(0, j), \qquad
  C_m(j) = \min_{\,mh \le i \le j-h} \bigl\{ C_{m-1}(i) + S(i, j) \bigr\},
  \qquad j \ge (m+1)h,
  \label{eq:appC-dp}
\end{equation}
the lower limit of $i$ ensuring that $m$ segments of length at least $h$
fit before the last break and the upper limit that the last segment has
length at least $h$. The argument minimiser $i^{\ast}_m(j)$ is stored, and
the breaks of the $m$-break solution for the full series are recovered
backward: $\tau_m = i^{\ast}_m(n)$, $\tau_{m-1} = i^{\ast}_{m-1}(\tau_m)$,
and so on, which is what the dictionary \texttt{B[m][j]} of the
implementation accumulates forward instead. This is the dynamic programme
that \citet{bai2003computation} propose for the global minimiser of
\citet{bai1998breaks}; the earlier literature used it for the
piecewise-constant case, and its correctness requires nothing beyond
\eqref{eq:appC-additive} and the induction on $m$: if $C_{m-1}(i)$ is the
true minimum for every $i$, then any partition of $[0, j)$ into $m+1$
segments has a last segment $[i, j)$ and a prefix partition of cost at least
$C_{m-1}(i)$, so its cost is at least the right-hand side of
\eqref{eq:appC-dp}, which is attained.

\paragraph{Complexity.}
The table of segment costs has one entry per admissible pair $(i, j)$, that
is $\binom{n-h+1}{2} + O(n)$ entries, of order $n^2/2$. The recursion
\eqref{eq:appC-dp} evaluates, for each $m \le M$ and each $j$, a minimum over
at most $n$ candidates, so its cost is $O(n^2 M)$ operations; the
back-tracking costs $O(M)$. The implementation fills the cost table by
solving each segment's least-squares problem independently, at a cost
proportional to the segment length, so the table costs $O(n^3)$ rather than
the $O(n^2)$ obtainable by updating running sums as in
\citet{bai2003computation}; for $n$ of the order of one hundred, as here,
the difference is a fraction of a second and the simpler code was preferred
for its transparency. Trimming enters only through the admissibility
condition $j - i \ge h$; with $h = \max\{8, \lfloor 0.15\,n \rfloor\}$ and
$M = 3$, the default window of $n = 101$ months gives $h = 15$ and a
feasible region $(m+1)h \le n$ for every $m \le 3$, so all four candidate
models exist.

\paragraph{Information criteria as implemented.}
For $m = 0, \dots, M$ let $\mathrm{SSR}_m = C_m(n)$ and let $p_m = k(m+1) +
m$ be the parameter count, with $k = 2$ for the trend model ($k = 1$ for the
mean model): one intercept and one slope per segment, plus the $m$ break
positions, counted as parameters by convention. The two criteria computed by
the code are
\begin{align}
  \mathrm{BIC}(m) &= n \log\frac{\mathrm{SSR}_m}{n} + p_m \log n, \label{eq:appC-bic}\\
  \mathrm{LWZ}(m) &= n \log\frac{\mathrm{SSR}_m}{n - p_m} + p_m\, c_0\, (\log n)^{2.1}, \qquad c_0 = 0.299, \label{eq:appC-lwz}
\end{align}
and the selected number of breaks is $\hat m = \arg\min_m$ of each. The BIC
form \eqref{eq:appC-bic} is the Gaussian log-likelihood, maximised over the
error variance, in which $-2\log\hat L = n\log(\mathrm{SSR}/n) + \text{const}$,
plus the Schwarz penalty $p \log n$ \citep{schwarz1978bic}; the two
differences from a textbook BIC are that the break positions are counted
among the parameters and that the constant and the segment variances are
treated as common, both as in \citet{bai2003computation}. The LWZ form
\eqref{eq:appC-lwz} is the modified Schwarz criterion of
\citet{liu1997segmented}, whose penalty grows as $(\log n)^{2+\delta_0}$ with
$\delta_0 = 0.1$ and whose variance estimate divides by $n - p_m$; the
constant $c_0 = 0.299$ is the value they recommend and the one used in the
simulations of \citet{bai2003computation}, which find that LWZ selects fewer
spurious breaks than BIC under serial correlation at the cost of some
power. The comparison across $m$ is legitimate because every $\mathrm{SSR}_m$
is a global minimum over the same observations and the same admissible
region; a sequential procedure that fixed earlier breaks before adding the
next would not have this property and would need the sup-$F$ machinery that
is not implemented. Neither criterion yields a confidence set for the break
dates; the text reports the agreement of the two selections and the
sensitivity to the truncation month instead.

\subsection{PELT: cost, penalty and pruning}
\label{sec:appC-pelt}

\paragraph{Objective.}
With the $\ell_2$ cost $\mathcal{C}(y_{i:j}) = \sum_{t=i}^{j-1} (y_t -
\bar y_{i:j})^2$ for a segment $[i,j)$ and a penalty $\lambda > 0$ per
changepoint, the problem of equation~\eqref{eq:pelt} of Part~IV is
\begin{equation}
  F(n) \;=\; \min_{m,\ \tau_1 < \dots < \tau_m}\ \Bigl\{ \sum_{k=1}^{m+1} \mathcal{C}(y_{\tau_{k-1}:\tau_k}) + \lambda m \Bigr\},
  \label{eq:appC-peltobj}
\end{equation}
and, by the same additivity argument as in Section~\ref{sec:appC-bp}, it
satisfies the optimal-partitioning recursion of \citet{jackson2005optimal},
\begin{equation}
  F(0) = -\lambda, \qquad F(j) \;=\; \min_{0 \le i < j}\ \bigl\{ F(i) + \mathcal{C}(y_{i:j}) + \lambda \bigr\},
  \label{eq:appC-op}
\end{equation}
where the number of breaks is not fixed in advance but chosen by the
penalty; the convention $F(0) = -\lambda$ cancels the penalty charged to the
first segment. Solving \eqref{eq:appC-op} for all $j$ costs $O(n^2)$ cost
evaluations.

\paragraph{The pruning condition.}
\citet{killick2012pelt} observe that if there exists a constant $K$ such that
for all $i < j < l$
\begin{equation}
  \mathcal{C}(y_{i:j}) + \mathcal{C}(y_{j:l}) + K \;\le\; \mathcal{C}(y_{i:l}),
  \label{eq:appC-Kcond}
\end{equation}
then any candidate $i$ for which
\begin{equation}
  F(i) + \mathcal{C}(y_{i:j}) + K \;\ge\; F(j)
  \label{eq:appC-prune}
\end{equation}
can never be the last changepoint before any later $l > j$ and may be
discarded from the candidate set. The proof is two lines. For $l > j$, by
\eqref{eq:appC-Kcond} and \eqref{eq:appC-prune},
$F(i) + \mathcal{C}(y_{i:l}) + \lambda \ge F(i) + \mathcal{C}(y_{i:j}) +
\mathcal{C}(y_{j:l}) + K + \lambda \ge F(j) + \mathcal{C}(y_{j:l}) +
\lambda$, so the partition that places a changepoint at $j$ after the
optimal prefix of $[0,j)$ is at least as good as any that ends its last
segment at $i$; hence $i$ never needs to be examined again. For the
$\ell_2$ cost, \eqref{eq:appC-Kcond} holds with $K = 0$, because adding a
split point to a segment can only reduce its residual sum of squares (the
split model nests the unsplit one). The recursion with pruning examines, at
each $j$, only the candidates not yet discarded; under the regularity
conditions of \citeauthor{killick2012pelt} (changepoints whose number grows
linearly with $n$) the expected number of surviving candidates is bounded and
the total cost is linear in $n$, whence the name. For a short series with
few changes, as here, the pruning saves little and the method is simply an
exact solution of \eqref{eq:appC-peltobj}; the point of using it is
exactness at a given penalty, not speed. The implementation
(\texttt{ruptures.Pelt}, \citealp{truong2020ruptures}) is called with
\texttt{model="l2"}, \texttt{min\_size=6} and \texttt{jump=1}, the last two
restricting segments to at least six months and allowing a changepoint at
every index; the library returns the segment end points including $n$, which
the wrapper drops.

\paragraph{The penalty $2\sigma^2 \log n$.}
Suppose the observations are independent Gaussian with piecewise-constant
mean and a common variance $\sigma^2$. Twice the negative log-likelihood of
a partition, up to a constant, is $\sigma^{-2}\sum_k \mathcal{C}(y_{\tau_{k-1}:\tau_k})$.
Each additional changepoint adds two parameters, a segment mean and a
position, so the Schwarz criterion in this scale adds $2\log n$ per
changepoint. Multiplying through by $\sigma^2$ to express the criterion in
the units of the cost gives
\begin{equation}
  \sum_k \mathcal{C}(y_{\tau_{k-1}:\tau_k}) + 2\sigma^2 \log n \cdot m,
\end{equation}
which is \eqref{eq:appC-peltobj} with $\lambda = 2\sigma^2\log n$. The
implementation replaces $\sigma^2$ by the sample variance of the whole
series, $\widehat{\mathrm{Var}}(y) = n^{-1}\sum_t (y_t - \bar y)^2$
(\texttt{numpy.var}, population form). This is the conservative choice
discussed in Part~IV: in a series with a trend or with genuine level
shifts, the overall variance exceeds the within-segment variance, so the
penalty is larger than the within-segment Schwarz penalty and fewer
changepoints are declared. \citet{yao1988schwarz} proves that the Schwarz
criterion, with a known or consistently estimated $\sigma^2$, selects the
true number of changes with probability tending to one in the
piecewise-constant Gaussian model; the result does not extend to a series
whose mean is a trend, for which a piecewise-constant fit must spend
changepoints to track the slope, and the PELT dates of this study are read
with that in mind.

\subsection{Bayesian online changepoint detection}
\label{sec:appC-bocpd}

\paragraph{Conjugate updates.}
Within a run the observations are modelled as
$y \mid \mu, \sigma^2 \sim \mathcal{N}(\mu, \sigma^2)$ with the
Normal--Inverse-Gamma prior
$\mu \mid \sigma^2 \sim \mathcal{N}(\mu_0, \sigma^2/\kappa_0)$,
$\sigma^2 \sim \mathrm{IG}(\alpha_0, \beta_0)$, whose joint density is
proportional to
\begin{equation}
  (\sigma^2)^{-\alpha_0 - 3/2} \exp\Bigl\{ -\frac{\beta_0}{\sigma^2} - \frac{\kappa_0 (\mu - \mu_0)^2}{2\sigma^2} \Bigr\}.
\end{equation}
Multiplying by the likelihood of one observation $y$, proportional to
$(\sigma^2)^{-1/2}\exp\{-(y-\mu)^2/(2\sigma^2)\}$, and completing the square
in $\mu$,
\begin{equation}
  \kappa_0 (\mu - \mu_0)^2 + (y - \mu)^2
  \;=\; (\kappa_0 + 1)\Bigl(\mu - \frac{\kappa_0\mu_0 + y}{\kappa_0 + 1}\Bigr)^2 + \frac{\kappa_0 (y - \mu_0)^2}{\kappa_0 + 1},
\end{equation}
shows that the posterior is again Normal--Inverse-Gamma with
\begin{equation}
  \mu_1 = \frac{\kappa_0 \mu_0 + y}{\kappa_0 + 1}, \quad
  \kappa_1 = \kappa_0 + 1, \quad
  \alpha_1 = \alpha_0 + \tfrac12, \quad
  \beta_1 = \beta_0 + \frac{\kappa_0 (y - \mu_0)^2}{2(\kappa_0 + 1)},
  \label{eq:appC-nigupdate}
\end{equation}
which are the four one-step updates of the implementation
(\texttt{mu\_new}, \texttt{kappa\_new}, \texttt{alpha\_new},
\texttt{beta\_new}). Applying \eqref{eq:appC-nigupdate} $r$ times and
collecting terms gives the batch form quoted in Part~IV, with
$\kappa_r = \kappa_0 + r$, $\mu_r = (\kappa_0\mu_0 + r\bar y)/\kappa_r$,
$\alpha_r = \alpha_0 + r/2$ and
$\beta_r = \beta_0 + \tfrac12\sum_{i \le r}(y_i - \bar y)^2 + \kappa_0 r (\bar y - \mu_0)^2 / (2\kappa_r)$;
the equality of the two forms is a matter of algebra (induction on $r$) and
is the reason the online recursion can be trusted to reproduce the batch
posterior.

\paragraph{Student-$t$ predictive.}
The predictive density of the next observation integrates the Gaussian
likelihood against the current Normal--Inverse-Gamma posterior. Integrating
$\mu$ first gives $y \mid \sigma^2 \sim \mathcal{N}(\mu_r, \sigma^2 (1 +
1/\kappa_r))$; integrating $\sigma^2 \sim \mathrm{IG}(\alpha_r, \beta_r)$
against this scale mixture of normals gives
\begin{equation}
  p(y \mid r) \;=\; \int \mathcal{N}\Bigl(y;\ \mu_r,\ \sigma^2\frac{\kappa_r + 1}{\kappa_r}\Bigr)\,\mathrm{IG}(\sigma^2; \alpha_r, \beta_r)\, d\sigma^2
  \;=\; t_{2\alpha_r}\Bigl( y;\ \mu_r,\ \frac{\beta_r(\kappa_r + 1)}{\alpha_r \kappa_r} \Bigr),
  \label{eq:appC-student}
\end{equation}
a Student-$t$ density with $2\alpha_r$ degrees of freedom, location $\mu_r$
and squared scale $\beta_r(\kappa_r+1)/(\alpha_r\kappa_r)$: the standard
result that a normal with inverse-gamma mixing on the variance is a
$t$ with degrees of freedom twice the shape and scale equal to the
rate over the shape times the normal's variance factor. The code evaluates
\eqref{eq:appC-student} through \texttt{scipy.stats.t.pdf} with
\texttt{df = 2*alpha} and \texttt{scale = sqrt(beta*(kappa+1)/(alpha*kappa))},
vectorised over the run lengths.

\paragraph{Run-length recursion and the mass on $r = 0$.}
Let $R_t(r) = P(r_t = r \mid y_{1:t})$ be the filtered run-length posterior
after $t$ observations, with $R_0(0) = 1$, and let $\pi_t(r) = p(y_{t+1} \mid
r_t = r, y_{1:t})$ be the predictive \eqref{eq:appC-student} evaluated with
the hyperparameters of the run of length $r$ ending at $t$. With a constant
hazard $H$ (a geometric prior on run lengths), the joint recursion of
\citet{adams2007bocpd} is, for the growth and the change moves,
\begin{equation}
  \tilde R_{t+1}(r+1) = R_t(r)\,\pi_t(r)\,(1 - H), \qquad
  \tilde R_{t+1}(0) = \sum_{r=0}^{t} R_t(r)\,\pi_t(r)\, H,
  \label{eq:appC-rl}
\end{equation}
followed by normalisation, $R_{t+1}(r) = \tilde R_{t+1}(r) / \sum_{r'}
\tilde R_{t+1}(r')$. Write $Z_t = \sum_{r} R_t(r)\pi_t(r)$ for the marginal
predictive likelihood of $y_{t+1}$. Then $\sum_{r \ge 1}\tilde R_{t+1}(r) =
(1-H) Z_t$ and $\tilde R_{t+1}(0) = H Z_t$, so the normaliser is $Z_t$ and
\begin{equation}
  R_{t+1}(0) \;=\; \frac{H Z_t}{H Z_t + (1-H) Z_t} \;=\; H \qquad \text{for every } t .
  \label{eq:appC-r0}
\end{equation}
The filtered probability that a changepoint occurs \emph{at} the current
step therefore equals the prior hazard identically, whatever the data: the
run-length posterior carries information about where the last change was
only through the relative masses at $r \ge 1$, which are reweighted by the
predictive likelihoods, and not through the mass at $r = 0$, which the
constant hazard fixes. This is why the detection statistics of Part~IV are
not $R_t(0)$. The first, $P(r_t \le 3 \mid y_{1:t}) = \sum_{r=0}^{3} R_t(r)$,
is $H$ plus the mass on runs of length one to three, that is, on changes
that occurred one to three months earlier and have since accumulated
evidence; its value exceeds $H$ only when the recent observations are
better explained by a fresh run than by the continuation of the old one. The
second is the maximum-a-posteriori run length $\hat r_t = \arg\max_r R_t(r)$,
whose resets from a value above twelve to a value at most two are the
online detector's changepoints. The implementation keeps the full
$(n+1)\times(n+1)$ matrix $R$, adds $10^{-300}$ to the normaliser against
underflow, and initialises $\mu_0$ and $\beta_0$ from the first six
observations with $\kappa_0 = \alpha_0 = 1$; the prior mean run length
$1/H = 60$ months is the only tuning constant. The recursion is the filtered
one; the smoothed posterior of \citet{fearnhead2006exact} would require a
backward pass and answers the retrospective question that the offline
methods of Sections~\ref{sec:appC-bp} and~\ref{sec:appC-pelt} already
address.