SPB Git

spb/wp44_uqo

Public
23.8 KB · 420 lines latex
1% ============================================================================2% Appendix C — Mathematical derivations3% Written by Agent L (Appendices B/C instance), 2026-11-03. Derivations and4% proofs referenced by Part IV (sections/part04_methodology.tex); notation5% follows Part IV. Each derivation is checked against the implementation in6% econometrics/structural_change.py, bibliometrics/03_citation_impact.py,7% bibliometrics/04_diversity_cagr.py, knowledge_graph/02_category_flows.py,8% taxonomy/classify.py, taxonomy/validate_sample.py and9% experiments/emerging_detection.py; where the code and the usual textbook10% convention differ, the code's convention is stated and the difference noted.11% ============================================================================1213\section{Mathematical derivations}14\label{sec:appC}1516Part~IV defines the statistics of this study and states their properties;17this appendix supplies the derivations behind those statements, in the order18in which Part~IV uses them. Each subsection fixes the object as implemented,19derives the property that the main text relies on, and records any point at20which the implementation departs from the convention of the cited source.21The derivations are elementary and self-contained; their purpose is not22novelty but the assurance that the quantities tabulated in Part~IX are what23the text says they are. Throughout, $n$ is the length of a monthly series,24$y_t = \log(N_t + 1)$ its logarithm, and $\mathbf{1}[\cdot]$ the indicator25function.2627\subsection{The interrupted time series as a linear regression}28\label{sec:appC-its}2930\paragraph{Design matrix.}31Equation~\eqref{eq:its} of Part~IV is the linear model $y = X\beta +32\varepsilon$ with $y \in \R^n$ and $X \in \R^{n \times 15}$ whose columns33are, for $t = 0, \dots, n-1$ and event index $t_e$,34\begin{equation}35  x_t \;=\; \bigl(1,\; t,\; \mathbf{1}[t \ge t_e],\; (t - t_e)\,\mathbf{1}[t \ge t_e],\;36  \mathbf{1}[\mathrm{month}(t) = 2],\; \dots,\; \mathbf{1}[\mathrm{month}(t) = 12]\bigr)^{\top},37  \label{eq:appC-design}38\end{equation}39exactly the data frame built by \texttt{\_design()} in40\texttt{econometrics/structural\_change.py} (columns \texttt{const},41\texttt{t}, \texttt{post}, \texttt{t\_post}, \texttt{m2}, \dots,42\texttt{m12}). The post-event trend is centred at $t_e$ rather than at zero43so that $\beta_2$ is the level change \emph{at} the event month and not the44extrapolated intercept difference at $t = 0$; with an uncentred post-trend45the two parametrisations are related by $\beta_2^{\mathrm{unc}} = \beta_2 -46\beta_3 t_e$, and only the centred one has the interpretation used in the47text. The eleven month dummies are linearly independent of the constant48because January is omitted; the full-column-rank condition requires at least49one observation in every month of the year, satisfied by any series of50twelve months or more, and at least one post-event month with $t > t_e$,51guaranteed by the driver's rule that twelve months follow every event.52The ordinary least-squares estimator is $\hat\beta = (X^{\top}X)^{-1}53X^{\top}y$ and the residuals $\hat\varepsilon = y - X\hat\beta$.5455\paragraph{The HAC covariance estimator.}56Under the assumption that $\varepsilon_t$ is a stationary, serially57correlated and possibly heteroskedastic sequence with58$\E[x_t \varepsilon_t] = 0$, the asymptotic covariance of $\sqrt{n}(\hat\beta59- \beta)$ is $Q^{-1}\Omega Q^{-1}$ with $Q = \lim n^{-1}X^{\top}X$ and60$\Omega = \sum_{\ell=-\infty}^{\infty} \Gamma_\ell$, where61$\Gamma_\ell = \E[x_t \varepsilon_t \varepsilon_{t-\ell} x_{t-\ell}^{\top}]$62is the autocovariance of the score at lag $\ell$. The estimator of63\citet{newey1987hac} replaces the infinite sum by a kernel-weighted finite64one,65\begin{equation}66  \hat\Omega_L \;=\; \hat\Gamma_0 + \sum_{\ell=1}^{L} w_\ell\,\bigl(\hat\Gamma_\ell + \hat\Gamma_\ell^{\top}\bigr),67  \qquad68  \hat\Gamma_\ell = \frac{1}{n}\sum_{t=\ell}^{n-1} \hat\varepsilon_t \hat\varepsilon_{t-\ell}\, x_t x_{t-\ell}^{\top},69  \qquad70  w_\ell = 1 - \frac{\ell}{L+1},71  \label{eq:appC-nw}72\end{equation}73and the variance estimate is $\hat V = n\,(X^{\top}X)^{-1} \hat\Omega_L74(X^{\top}X)^{-1}$, which is the expression of Part~IV once the factors of75$n$ are absorbed into the un-normalised $\hat\Gamma_\ell$. The Bartlett76weights $w_\ell$ are what make $\hat\Omega_L$ positive semi-definite for every77sample: $\hat\Omega_L$ can be written as $\frac{1}{n(L+1)}\sum_{t}78\bigl(\sum_{s=t-L}^{t} \hat u_s\bigr)\bigl(\sum_{s=t-L}^{t} \hat u_s\bigr)^{\top}$79up to end effects, with $\hat u_s = x_s \hat\varepsilon_s$, a sum of outer80products, whereas the truncated sum with $w_\ell \equiv 1$ is not guaranteed81to be positive semi-definite. The implementation calls82\texttt{statsmodels} with \texttt{cov\_type="HAC"} and \texttt{maxlags=6}, so83$L = 6$ and $w_\ell = 1 - \ell/7$ as stated in Part~IV, without the84small-sample degrees-of-freedom correction and with $p$-values from the85normal approximation to the $t$-ratio, which is the library's default for86robust covariances. With $n$ of the order of one hundred and fourteen87regressors, the difference between the normal and the $t_{n-15}$ reference88distribution is immaterial at the magnitudes of the $t$-statistics reported;89what matters is the placebo rank of the next paragraph, which does not90depend on any reference distribution.9192\paragraph{The placebo $p$-value as a rank statistic.}93Let $\mathcal{T}' = \{t' : 24 \le t' \le n - 13\}$ be the candidate event94indices (the implementation's \texttt{min\_pre = 24} and95\texttt{min\_post = 12}) and, for each $t' \in \mathcal{T}'$, let96$t^{(2)}_{t'} = \hat\beta_2(t') / \widehat{\mathrm{se}}(\hat\beta_2(t'))$ be97the HAC $t$-ratio of the level change when the model98\eqref{eq:appC-design} is fitted with $t_e = t'$. The placebo value of the99true event $t_e \in \mathcal{T}'$ is100\begin{equation}101  p^{\mathrm{plac}} \;=\; \frac{1}{|\mathcal{T}'|}\sum_{t' \in \mathcal{T}'}102  \mathbf{1}\bigl[\,|t^{(2)}_{t'}| \ge |t^{(2)}_{t_e}|\,\bigr]103  \;=\; \frac{\mathrm{rank}\bigl(|t^{(2)}_{t_e}|\bigr)}{|\mathcal{T}'|},104\end{equation}105where $\mathrm{rank}$ is the descending rank of the true event's absolute106$t$-ratio among the $|\mathcal{T}'|$ values, ties counted as not better107(\texttt{>=}). Three properties follow directly. Because $t_e \in108\mathcal{T}'$, the sum contains the term $t' = t_e$, so $p^{\mathrm{plac}}109\ge 1/|\mathcal{T}'|$, which for the default window ($n = 101$, hence110$|\mathcal{T}'| = 65$) is about $0.015$: a placebo rank cannot be smaller,111however large the true statistic, and a reported value of $1/|\mathcal{T}'|$112means ``largest of all candidate dates'', not ``significant at that level''.113Because the statistic is a rank, it is invariant to any monotone114transformation of the $|t|$-ratios and in particular does not depend on the115HAC bandwidth through anything but the ordering it induces. And if the116candidate dates were exchangeable under a null hypothesis of no break117anywhere, $p^{\mathrm{plac}}$ would be a valid permutation-type $p$-value118with $P(p^{\mathrm{plac}} \le \alpha) \le \alpha$; they are not exchangeable,119since all $|\mathcal{T}'|$ regressions share the same $n$ observations and120adjacent candidates give nearly identical fits, so the value is reported as121a rank and read as such. The same construction gives the placebo rank of122the slope change from $t^{(3)}_{t'}$. The driver additionally computes the123rank restricted to the baseline candidates $\{t' < t_e\}$, which answers the124protocol's original question (is the reference date special relative to the125pre-period?) and is reported next to the full-range rank.126127\subsection{Bai--Perron multiple breaks by dynamic programming}128\label{sec:appC-bp}129130\paragraph{Segment costs and additivity.}131For a segment $[i, j)$ of the series with $j - i \ge h$, define132\begin{equation}133  S(i, j) \;=\; \min_{a, b \in \R}\ \sum_{t=i}^{j-1} (y_t - a - b\,t)^2,134  \label{eq:appC-segcost}135\end{equation}136the residual sum of squares of a linear trend fitted to the segment alone137(the implementation solves the two-parameter least-squares problem for each138segment with \texttt{numpy.linalg.lstsq}; with \texttt{trend=False} the139regressor $t$ is dropped and $S(i,j)$ is the centred sum of squares). The140model with $m$ breaks $0 < \tau_1 < \dots < \tau_m < n$ assigns its own141$(a, b)$ to each of the $m+1$ segments $[\tau_{k-1}, \tau_k)$ with142$\tau_0 = 0$ and $\tau_{m+1} = n$. Because the parameters of different143segments are distinct and each observation enters exactly one segment, the144residual sum of squares of the whole model separates,145\begin{equation}146  \mathrm{SSR}(\tau_1, \dots, \tau_m) \;=\; \sum_{k=1}^{m+1} S(\tau_{k-1}, \tau_k),147  \label{eq:appC-additive}148\end{equation}149and the global minimisation over $\tau$ decomposes: the optimal partition of150the first $j$ observations into $m+1$ segments contains, as its first $m$151segments, an optimal partition of the first $\tau_m$ observations into $m$152segments. This is the principle of optimality, and it is the only property of153the cost that the dynamic programme uses; it holds for any cost that is a154sum over segments of a function of the segment alone, which is why the same155recursion serves the piecewise-constant model of Section~\ref{sec:appC-pelt}.156It fails for a model with parameters shared across segments (a common slope,157or the seasonal dummies of the ITS), which is why the segmented regression158carries no seasonal terms and the series is seasonally adjusted beforehand159by subtracting the baseline-period month means160(\texttt{run\_structural\_change.py}).161162\paragraph{The recursion.}163Let $C_m(j)$ be the minimum of \eqref{eq:appC-additive} over partitions of164observations $0, \dots, j-1$ into $m+1$ admissible segments. Then165\begin{equation}166  C_0(j) = S(0, j), \qquad167  C_m(j) = \min_{\,mh \le i \le j-h} \bigl\{ C_{m-1}(i) + S(i, j) \bigr\},168  \qquad j \ge (m+1)h,169  \label{eq:appC-dp}170\end{equation}171the lower limit of $i$ ensuring that $m$ segments of length at least $h$172fit before the last break and the upper limit that the last segment has173length at least $h$. The argument minimiser $i^{\ast}_m(j)$ is stored, and174the breaks of the $m$-break solution for the full series are recovered175backward: $\tau_m = i^{\ast}_m(n)$, $\tau_{m-1} = i^{\ast}_{m-1}(\tau_m)$,176and so on, which is what the dictionary \texttt{B[m][j]} of the177implementation accumulates forward instead. This is the dynamic programme178that \citet{bai2003computation} propose for the global minimiser of179\citet{bai1998breaks}; the earlier literature used it for the180piecewise-constant case, and its correctness requires nothing beyond181\eqref{eq:appC-additive} and the induction on $m$: if $C_{m-1}(i)$ is the182true minimum for every $i$, then any partition of $[0, j)$ into $m+1$183segments has a last segment $[i, j)$ and a prefix partition of cost at least184$C_{m-1}(i)$, so its cost is at least the right-hand side of185\eqref{eq:appC-dp}, which is attained.186187\paragraph{Complexity.}188The table of segment costs has one entry per admissible pair $(i, j)$, that189is $\binom{n-h+1}{2} + O(n)$ entries, of order $n^2/2$. The recursion190\eqref{eq:appC-dp} evaluates, for each $m \le M$ and each $j$, a minimum over191at most $n$ candidates, so its cost is $O(n^2 M)$ operations; the192back-tracking costs $O(M)$. The implementation fills the cost table by193solving each segment's least-squares problem independently, at a cost194proportional to the segment length, so the table costs $O(n^3)$ rather than195the $O(n^2)$ obtainable by updating running sums as in196\citet{bai2003computation}; for $n$ of the order of one hundred, as here,197the difference is a fraction of a second and the simpler code was preferred198for its transparency. Trimming enters only through the admissibility199condition $j - i \ge h$; with $h = \max\{8, \lfloor 0.15\,n \rfloor\}$ and200$M = 3$, the default window of $n = 101$ months gives $h = 15$ and a201feasible region $(m+1)h \le n$ for every $m \le 3$, so all four candidate202models exist.203204\paragraph{Information criteria as implemented.}205For $m = 0, \dots, M$ let $\mathrm{SSR}_m = C_m(n)$ and let $p_m = k(m+1) +206m$ be the parameter count, with $k = 2$ for the trend model ($k = 1$ for the207mean model): one intercept and one slope per segment, plus the $m$ break208positions, counted as parameters by convention. The two criteria computed by209the code are210\begin{align}211  \mathrm{BIC}(m) &= n \log\frac{\mathrm{SSR}_m}{n} + p_m \log n, \label{eq:appC-bic}\\212  \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}213\end{align}214and the selected number of breaks is $\hat m = \arg\min_m$ of each. The BIC215form \eqref{eq:appC-bic} is the Gaussian log-likelihood, maximised over the216error variance, in which $-2\log\hat L = n\log(\mathrm{SSR}/n) + \text{const}$,217plus the Schwarz penalty $p \log n$ \citep{schwarz1978bic}; the two218differences from a textbook BIC are that the break positions are counted219among the parameters and that the constant and the segment variances are220treated as common, both as in \citet{bai2003computation}. The LWZ form221\eqref{eq:appC-lwz} is the modified Schwarz criterion of222\citet{liu1997segmented}, whose penalty grows as $(\log n)^{2+\delta_0}$ with223$\delta_0 = 0.1$ and whose variance estimate divides by $n - p_m$; the224constant $c_0 = 0.299$ is the value they recommend and the one used in the225simulations of \citet{bai2003computation}, which find that LWZ selects fewer226spurious breaks than BIC under serial correlation at the cost of some227power. The comparison across $m$ is legitimate because every $\mathrm{SSR}_m$228is a global minimum over the same observations and the same admissible229region; a sequential procedure that fixed earlier breaks before adding the230next would not have this property and would need the sup-$F$ machinery that231is not implemented. Neither criterion yields a confidence set for the break232dates; the text reports the agreement of the two selections and the233sensitivity to the truncation month instead.234235\subsection{PELT: cost, penalty and pruning}236\label{sec:appC-pelt}237238\paragraph{Objective.}239With the $\ell_2$ cost $\mathcal{C}(y_{i:j}) = \sum_{t=i}^{j-1} (y_t -240\bar y_{i:j})^2$ for a segment $[i,j)$ and a penalty $\lambda > 0$ per241changepoint, the problem of equation~\eqref{eq:pelt} of Part~IV is242\begin{equation}243  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\},244  \label{eq:appC-peltobj}245\end{equation}246and, by the same additivity argument as in Section~\ref{sec:appC-bp}, it247satisfies the optimal-partitioning recursion of \citet{jackson2005optimal},248\begin{equation}249  F(0) = -\lambda, \qquad F(j) \;=\; \min_{0 \le i < j}\ \bigl\{ F(i) + \mathcal{C}(y_{i:j}) + \lambda \bigr\},250  \label{eq:appC-op}251\end{equation}252where the number of breaks is not fixed in advance but chosen by the253penalty; the convention $F(0) = -\lambda$ cancels the penalty charged to the254first segment. Solving \eqref{eq:appC-op} for all $j$ costs $O(n^2)$ cost255evaluations.256257\paragraph{The pruning condition.}258\citet{killick2012pelt} observe that if there exists a constant $K$ such that259for all $i < j < l$260\begin{equation}261  \mathcal{C}(y_{i:j}) + \mathcal{C}(y_{j:l}) + K \;\le\; \mathcal{C}(y_{i:l}),262  \label{eq:appC-Kcond}263\end{equation}264then any candidate $i$ for which265\begin{equation}266  F(i) + \mathcal{C}(y_{i:j}) + K \;\ge\; F(j)267  \label{eq:appC-prune}268\end{equation}269can never be the last changepoint before any later $l > j$ and may be270discarded from the candidate set. The proof is two lines. For $l > j$, by271\eqref{eq:appC-Kcond} and \eqref{eq:appC-prune},272$F(i) + \mathcal{C}(y_{i:l}) + \lambda \ge F(i) + \mathcal{C}(y_{i:j}) +273\mathcal{C}(y_{j:l}) + K + \lambda \ge F(j) + \mathcal{C}(y_{j:l}) +274\lambda$, so the partition that places a changepoint at $j$ after the275optimal prefix of $[0,j)$ is at least as good as any that ends its last276segment at $i$; hence $i$ never needs to be examined again. For the277$\ell_2$ cost, \eqref{eq:appC-Kcond} holds with $K = 0$, because adding a278split point to a segment can only reduce its residual sum of squares (the279split model nests the unsplit one). The recursion with pruning examines, at280each $j$, only the candidates not yet discarded; under the regularity281conditions of \citeauthor{killick2012pelt} (changepoints whose number grows282linearly with $n$) the expected number of surviving candidates is bounded and283the total cost is linear in $n$, whence the name. For a short series with284few changes, as here, the pruning saves little and the method is simply an285exact solution of \eqref{eq:appC-peltobj}; the point of using it is286exactness at a given penalty, not speed. The implementation287(\texttt{ruptures.Pelt}, \citealp{truong2020ruptures}) is called with288\texttt{model="l2"}, \texttt{min\_size=6} and \texttt{jump=1}, the last two289restricting segments to at least six months and allowing a changepoint at290every index; the library returns the segment end points including $n$, which291the wrapper drops.292293\paragraph{The penalty $2\sigma^2 \log n$.}294Suppose the observations are independent Gaussian with piecewise-constant295mean and a common variance $\sigma^2$. Twice the negative log-likelihood of296a partition, up to a constant, is $\sigma^{-2}\sum_k \mathcal{C}(y_{\tau_{k-1}:\tau_k})$.297Each additional changepoint adds two parameters, a segment mean and a298position, so the Schwarz criterion in this scale adds $2\log n$ per299changepoint. Multiplying through by $\sigma^2$ to express the criterion in300the units of the cost gives301\begin{equation}302  \sum_k \mathcal{C}(y_{\tau_{k-1}:\tau_k}) + 2\sigma^2 \log n \cdot m,303\end{equation}304which is \eqref{eq:appC-peltobj} with $\lambda = 2\sigma^2\log n$. The305implementation replaces $\sigma^2$ by the sample variance of the whole306series, $\widehat{\mathrm{Var}}(y) = n^{-1}\sum_t (y_t - \bar y)^2$307(\texttt{numpy.var}, population form). This is the conservative choice308discussed in Part~IV: in a series with a trend or with genuine level309shifts, the overall variance exceeds the within-segment variance, so the310penalty is larger than the within-segment Schwarz penalty and fewer311changepoints are declared. \citet{yao1988schwarz} proves that the Schwarz312criterion, with a known or consistently estimated $\sigma^2$, selects the313true number of changes with probability tending to one in the314piecewise-constant Gaussian model; the result does not extend to a series315whose mean is a trend, for which a piecewise-constant fit must spend316changepoints to track the slope, and the PELT dates of this study are read317with that in mind.318319\subsection{Bayesian online changepoint detection}320\label{sec:appC-bocpd}321322\paragraph{Conjugate updates.}323Within a run the observations are modelled as324$y \mid \mu, \sigma^2 \sim \mathcal{N}(\mu, \sigma^2)$ with the325Normal--Inverse-Gamma prior326$\mu \mid \sigma^2 \sim \mathcal{N}(\mu_0, \sigma^2/\kappa_0)$,327$\sigma^2 \sim \mathrm{IG}(\alpha_0, \beta_0)$, whose joint density is328proportional to329\begin{equation}330  (\sigma^2)^{-\alpha_0 - 3/2} \exp\Bigl\{ -\frac{\beta_0}{\sigma^2} - \frac{\kappa_0 (\mu - \mu_0)^2}{2\sigma^2} \Bigr\}.331\end{equation}332Multiplying by the likelihood of one observation $y$, proportional to333$(\sigma^2)^{-1/2}\exp\{-(y-\mu)^2/(2\sigma^2)\}$, and completing the square334in $\mu$,335\begin{equation}336  \kappa_0 (\mu - \mu_0)^2 + (y - \mu)^2337  \;=\; (\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},338\end{equation}339shows that the posterior is again Normal--Inverse-Gamma with340\begin{equation}341  \mu_1 = \frac{\kappa_0 \mu_0 + y}{\kappa_0 + 1}, \quad342  \kappa_1 = \kappa_0 + 1, \quad343  \alpha_1 = \alpha_0 + \tfrac12, \quad344  \beta_1 = \beta_0 + \frac{\kappa_0 (y - \mu_0)^2}{2(\kappa_0 + 1)},345  \label{eq:appC-nigupdate}346\end{equation}347which are the four one-step updates of the implementation348(\texttt{mu\_new}, \texttt{kappa\_new}, \texttt{alpha\_new},349\texttt{beta\_new}). Applying \eqref{eq:appC-nigupdate} $r$ times and350collecting terms gives the batch form quoted in Part~IV, with351$\kappa_r = \kappa_0 + r$, $\mu_r = (\kappa_0\mu_0 + r\bar y)/\kappa_r$,352$\alpha_r = \alpha_0 + r/2$ and353$\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)$;354the equality of the two forms is a matter of algebra (induction on $r$) and355is the reason the online recursion can be trusted to reproduce the batch356posterior.357358\paragraph{Student-$t$ predictive.}359The predictive density of the next observation integrates the Gaussian360likelihood against the current Normal--Inverse-Gamma posterior. Integrating361$\mu$ first gives $y \mid \sigma^2 \sim \mathcal{N}(\mu_r, \sigma^2 (1 +3621/\kappa_r))$; integrating $\sigma^2 \sim \mathrm{IG}(\alpha_r, \beta_r)$363against this scale mixture of normals gives364\begin{equation}365  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^2366  \;=\; t_{2\alpha_r}\Bigl( y;\ \mu_r,\ \frac{\beta_r(\kappa_r + 1)}{\alpha_r \kappa_r} \Bigr),367  \label{eq:appC-student}368\end{equation}369a Student-$t$ density with $2\alpha_r$ degrees of freedom, location $\mu_r$370and squared scale $\beta_r(\kappa_r+1)/(\alpha_r\kappa_r)$: the standard371result that a normal with inverse-gamma mixing on the variance is a372$t$ with degrees of freedom twice the shape and scale equal to the373rate over the shape times the normal's variance factor. The code evaluates374\eqref{eq:appC-student} through \texttt{scipy.stats.t.pdf} with375\texttt{df = 2*alpha} and \texttt{scale = sqrt(beta*(kappa+1)/(alpha*kappa))},376vectorised over the run lengths.377378\paragraph{Run-length recursion and the mass on $r = 0$.}379Let $R_t(r) = P(r_t = r \mid y_{1:t})$ be the filtered run-length posterior380after $t$ observations, with $R_0(0) = 1$, and let $\pi_t(r) = p(y_{t+1} \mid381r_t = r, y_{1:t})$ be the predictive \eqref{eq:appC-student} evaluated with382the hyperparameters of the run of length $r$ ending at $t$. With a constant383hazard $H$ (a geometric prior on run lengths), the joint recursion of384\citet{adams2007bocpd} is, for the growth and the change moves,385\begin{equation}386  \tilde R_{t+1}(r+1) = R_t(r)\,\pi_t(r)\,(1 - H), \qquad387  \tilde R_{t+1}(0) = \sum_{r=0}^{t} R_t(r)\,\pi_t(r)\, H,388  \label{eq:appC-rl}389\end{equation}390followed by normalisation, $R_{t+1}(r) = \tilde R_{t+1}(r) / \sum_{r'}391\tilde R_{t+1}(r')$. Write $Z_t = \sum_{r} R_t(r)\pi_t(r)$ for the marginal392predictive likelihood of $y_{t+1}$. Then $\sum_{r \ge 1}\tilde R_{t+1}(r) =393(1-H) Z_t$ and $\tilde R_{t+1}(0) = H Z_t$, so the normaliser is $Z_t$ and394\begin{equation}395  R_{t+1}(0) \;=\; \frac{H Z_t}{H Z_t + (1-H) Z_t} \;=\; H \qquad \text{for every } t .396  \label{eq:appC-r0}397\end{equation}398The filtered probability that a changepoint occurs \emph{at} the current399step therefore equals the prior hazard identically, whatever the data: the400run-length posterior carries information about where the last change was401only through the relative masses at $r \ge 1$, which are reweighted by the402predictive likelihoods, and not through the mass at $r = 0$, which the403constant hazard fixes. This is why the detection statistics of Part~IV are404not $R_t(0)$. The first, $P(r_t \le 3 \mid y_{1:t}) = \sum_{r=0}^{3} R_t(r)$,405is $H$ plus the mass on runs of length one to three, that is, on changes406that occurred one to three months earlier and have since accumulated407evidence; its value exceeds $H$ only when the recent observations are408better explained by a fresh run than by the continuation of the old one. The409second is the maximum-a-posteriori run length $\hat r_t = \arg\max_r R_t(r)$,410whose resets from a value above twelve to a value at most two are the411online detector's changepoints. The implementation keeps the full412$(n+1)\times(n+1)$ matrix $R$, adds $10^{-300}$ to the normaliser against413underflow, and initialises $\mu_0$ and $\beta_0$ from the first six414observations with $\kappa_0 = \alpha_0 = 1$; the prior mean run length415$1/H = 60$ months is the only tuning constant. The recursion is the filtered416one; the smoothed posterior of \citet{fearnhead2006exact} would require a417backward pass and answers the retrospective question that the offline418methods of Sections~\ref{sec:appC-bp} and~\ref{sec:appC-pelt} already419address.420
Files · 3951030904d1ec09140fb0be444dba70975f7bdf
Menu

Theme