suppReferences

# Warped Dynamic Linear Models for Time Series of Counts

Brian King and Daniel R. Kowal
Department of Statistics, Rice University Brian King is PhD Graduate,
Department of Statistics, Rice University, Houston, TX (bking@rice.edu).
Daniel R. Kowal is Dobelman Family Assistant Professor, Department of
Statistics, Rice University, Houston, TX (Daniel.Kowal@rice.edu). This
material is based upon work supported by the National Science Foundation: a
Graduate Research Fellowship under Grant No. 1842494 (King) and SES-2214726
(Kowal).

###### Abstract

Dynamic Linear Models (DLMs) are commonly employed for time series analysis
due to their versatile structure, simple recursive updating, ability to handle
missing data, and probabilistic forecasting. However, the options for count
time series are limited: Gaussian DLMs require continuous data, while Poisson-
based alternatives often lack sufficient modeling flexibility. We introduce a
novel semiparametric methodology for count time series by _warping_ a Gaussian
DLM. The warping function has two components: a (nonparametric) transformation
operator that provides distributional flexibility and a rounding operator that
ensures the correct support for the discrete data-generating process. We
develop conjugate inference for the warped DLM, which enables analytic and
recursive updates for the state space filtering and smoothing distributions.
We leverage these results to produce customized and efficient algorithms for
inference and forecasting, including Monte Carlo simulation for offline
analysis and an optimal particle filter for online inference. This framework
unifies and extends a variety of discrete time series models and is valid for
natural counts, rounded values, and multivariate observations. Simulation
studies illustrate the excellent forecasting capabilities of the warped DLM.
The proposed approach is applied to a multivariate time series of daily
overdose counts and demonstrates both modeling and computational successes.

Keywords: Bayesian statistics; state-space model; particle filter; selection
normal

##  1 Introduction

Count time series data inherit all the complexities of continuous time series
data: the time-ordered observations may be multivariate, seasonal, dependent
on exogenous variables, and exhibit a wide variety of autocorrelation
structures. At the same time, count data often present uniquely challenging
distributional features, including zero-inflation, over-/underdispersion,
boundedness or censoring, and heaping. Additionally, discrete data require
distinct strategies for probabilistic forecasting, uncertainty quantification,
and evaluation. As modern datasets commonly feature higher resolutions and
lengthier time series, computational tools for both online and offline
inference and forecasting are in demand. Fundamentally, the goals in count
time series modeling are similar to those in the continuous setting:
forecasting, trend filtering/smoothing, seasonal decomposition, and
characterization of inter- and intra-series dependence, among others.

In this paper, we develop methods, theory, and computing tools for a broad
class of multivariate state space models that address each of these challenges
and objectives. The core model is defined by warping a Gaussian Dynamic Linear
Model (DLM; West and Harrison,, 1997), which we refer to as a _warped DLM_
(warpDLM):

$\displaystyle\bm{y}_{t}$ $\displaystyle=h\circ g^{-1}(\bm{z}_{t})$ (warping)
(1) $\displaystyle\\{\bm{z}_{t}\\}_{t=1}^{T}$ $\displaystyle\sim\text{DLM}$
(see (3) and (4)) (2)

where $\bm{y}_{t}\in\mathbb{N}^{n}$ is the observed count data and
$\bm{z}_{t}\in\mathbb{R}^{n}$ is continuous latent data. The warping has two
components: a rounding operator $h:\mathcal{T}\rightarrow\mathbb{N}^{n}$,
which ensures the correct support for the discrete data-generating process,
and a monotone transformation function
$g:\mathcal{T}\rightarrow\mathbb{R}^{n}$, which endows (possibly
nonparametric) flexibility in the marginal distributions. The latent DLM
enables straightforward embedding of familiar dynamic modeling structures such
as local levels, seasonality, or covariates, along with a natural way of
dealing with missing data. Note that this article focuses on modeling counts,
but the framework is easily adaptable for general integer-valued or rounded
data.

To isolate the rounding and transformation operations, the warping operation
may be decomposed into $\bm{y}_{t}=h(\bm{y}_{t}^{*})$ and
$\bm{z}_{t}=g(\bm{y}_{t}^{*})$ for a latent continuous variable
$\bm{y}_{t}^{*}$. We emphasize that the warpDLM is fundamentally distinct from
the “transformed DLM” strategy of (i) fitting a Gaussian DLM to (possibly
transformed, e.g., logarithmically) count data and (ii) rounding the resulting
(continuous) forecasts. The “transformed DLM” fails to account for the
discreteness of the data in the model-fitting process and introduces a
critical mismatch between the _fitted model_ and the _data-generating
process_. If the terminal rounding step is omitted, then the data-generating
process is not discrete; yet unless the rounding step is included within the
model-fitting process—as in the warpDLM—then the fitted model fails to account
for crucial features in the data. In particular, the rounding operation is a
nontrivial component of the model: within the warpDLM, $h$ provides the
capability to model challenging discrete distributional features, such as
zero-inflation, boundedness, or censoring (see Section 2). Transformation-only
models are known to be ineffective in many settings, such as low counts
(O’Hara and Kotze,, 2010), yet the warpDLM excels precisely in this case (see
Section 5.2).

A key contribution of this article is to develop conjugate inference for the
warpDLM. In particular, we show that the warpDLM likelihood is conjugate to
the _selection normal distribution_ (e.g., Arellano-Valle et al.,, 2006).
Based on this result, we derive analytic and recursive updates for the warpDLM
filtering and smoothing distributions. Crucially, we provide direct Monte
Carlo simulators for these distributions—as well as the count-valued
forecasting distribution—and construct an optimal particle filter for online
inference. These models, derivations, and algorithms remain valid in the
multivariate setting, and provide significant advancements over existing
latent models for count time series (e.g., Jia et al.,, 2021). To the best of
our knowledge, these results are unique for multivariate count time series
models.

The warpDLM offers a unified framework for several discrete data models
(dynamic as well as static) and incorporates strategies which have proven
successful in other related methods. In the non-dynamic realm, Siegfried and
Hothorn, (2020) demonstrated the benefits of learned transformations for
discrete data linear regression, while Kowal and Wu, (2021) adopted a related
transformation and rounding strategy to model heaped count data. Both took a
frequentist approach to estimation. Among Bayesian methods, Canale and Dunson,
(2011) and Canale and Dunson, (2013) similarly applied rounded Gaussian and
Dirichlet processes, respectively, without the transformation considerations.
Kowal and Canale, (2020) showcased the advantage of both the rounding _and_
transformation components within a static regression setting, and Kowal and
Wu, (2022) extended this framework to incorporate Bayesian nonparametric
estimation of the transformation. Looking to dynamic models, when $g$ is
viewed as a copula, the warpDLM resembles the count time series model proposed
by Jia et al., (2021). Unlike Jia et al., (2021), we do not focus on
stationary latent Gaussian processes, but rather incorporate DLMs to enable
nonstationary modeling, dynamic covariates, and Bayesian inference and
forecasting within a familiar setting. The warpDLM also generalizes recent
work done in the binary data space, in particular the dynamic probit model of
Fasano et al., (2021). In doing so, we construct novel theory utilizing a
broader class of distributions and requiring distinct computations of relevant
posterior quantities.

The warpDLM framework falls in the category of generalized state space models,
which are traditionally separated into two classes, stemming from Cox et al.,
(1981): parameter-driven and observation-driven. In observation-driven models,
the state process is treated as an explicit function of past data values, and
methods often utilize likelihood-based inference born out of GLM theory (see
Fokianos, (2015) for a review). Bayesian methods, including the warpDLM
framework, are more commonly parameter-driven, which means the latent state
parameter is treated as stochastic; model learning typically proceeds via
recursive updating, e.g., using the Kalman filter (Kalman,, 1960). The DLM is
the most well-known model in this category, but relies on Gaussian assumptions
unmet by count data. Dynamic Generalized Linear Models (DGLMs) were developed
to adapt the state space framework for non-Gaussian data within the
exponential family (West et al.,, 1985). For most count data, Poisson is the
only available observational density that belongs to the exponential family.
Binomial DGLMs are applicable for bounded data, yet perform quite poorly in
the case of zero-inflation or heaping on the boundary (see Section 5.2). The
negative binomial distribution with fixed dispersion parameter is also
exponential family, but there is often little guidance for determining the
dispersion parameter. Berry and West, (2019) recently extended the DGLM family
by mixing Bernoulli and Poisson DGLMs to better model zero-inflated and
overdispersed count data.

The Poisson DGLM and extensions provide “data coherent” (Freeland and McCabe,,
2004) inferences and forecasts, in the sense that predictions are appropriate
to the type of data (count time series), but the closed-form Kalman filtering
results are generally unavailable. In most DGLM specifications, the evolution
equation is assumed Gaussian (see (4)), which necessitates linearized
approximations or MCMC algorithms for smoothing, filtering, and forecasting.
West et al., (1985) prioritized closed-form and conjugate recursions, yet
require approximate and moment-based (rather than distributional) updates for
state parameters. Fahrmeir, (1992) designed an extended Kalman Filter to
estimate posterior modes in a multivariate setting. Durbin and Koopman, (2000)
used importance sampling based on a linear approximation. Frühwirth-Schnatter
and Wagner, (2006) constructed an approximate Gibbs sampler using data
augmentation and mixture sampling. Another strategy is to replace the Gaussian
evolution equation. Gamerman et al., (2013) used a Poisson model with a
multiplicative state update to preserve analytic and recursive inference.
However, this updating structure has limited dynamic flexibility, for example
to include seasonality or covariates. Aktekin et al., (2018) proposed a
multivariate extension of this multiplicative model, but the analytic updating
results were not preserved and the multivariate structure only accommodated
positive correlations among the series. Another exception is Bradley et al.,
(2018), who proposed log-Gamma processes that are conditionally conjugate to
the Poisson distribution. This model still requires Gibbs sampling for all
state smoothing, filtering, and forecasting distributions.

The common limitations among existing state space models for count data are
(i) a lack of exact, coherent, and recursive updates for filtering, smoothing,
and forecasting distributions and (ii) restricted options for count-valued
distributions. The warpDLM framework directly addresses and overcomes both
limitations.

The paper is organized as follows. Section 2 introduces DLMs, the proposed
model, and examples for the rounding and transformation. In Section 3, we
derive the smoothing, filtering, and forecasting distributions for the
warpDLM. We discuss computing strategies for online and offline inference in
Section 4. Finally, we present forecasting results on simulated data as well
as a real-data application in Section 5 before concluding. Supplementary
material includes proofs of presented theorems, additional simulations, and
further details on the application data and model specification. Code to
reproduce all findings is also available on GitHub. warpDLM functionality is
included in the `R` package `countSTAR` (King and Kowal,, 2023), which is
available on CRAN and documented in an online vignette.

##  2 Dynamic Linear Models and Time Series of Counts

The broad success of Bayesian time series analysis has largely been driven by
Dynamic Linear Models (DLMs), also known as linear state space models. The DLM
framework subsumes ARIMA models and provides decomposition of time series,
dynamic regression analysis, and multivariate modeling capabilities.
Additionally, the sequential updating structure provides a simple way to deal
with missing observations (cf. Durbin and Koopman,, 2012, Section 2.7). DLMs
are widely popular not only because of their versatility, but also because the
filtering, smoothing, and forecasting distributions are available in closed-
form via the recursive Kalman filter.

A DLM is defined by two equations: (i) the observation equation, which
specifies how the observations are related to the latent state vector and (ii)
the state evolution equation, which describes how the states are updated in a
Markovian fashion. We present the Gaussian DLM for an _observable_ continuous
$n$-dimensional time series $\\{\bm{z}_{t}\\}_{t=1}^{T}$, but note that these
continuous variables are latent (i.e., unobservable) within the warpDLM:

$\displaystyle\bm{z}_{t}$
$\displaystyle=\bm{F}_{t}\bm{\theta}_{t}+\bm{v}_{t},\quad\bm{v}_{t}\sim
N_{n}(\bm{0},\bm{V}_{t})$ (3) $\displaystyle\bm{\theta}_{t}$
$\displaystyle=\bm{G}_{t}\bm{\theta}_{t-1}+\bm{w}_{t},\quad\bm{w}_{t}\sim
N_{p}(\bm{0},\bm{W}_{t})$ (4)

for $t=1,\ldots,T$, where $\\{\bm{v}_{t},\bm{w}_{t}\\}_{t=1}^{T}$ are mutually
independent and $\bm{\theta}_{0}\sim N_{p}(\bm{a}_{0},\bm{R}_{0})$. Depending
on the DLM specification, the $p$-dimensional state vector $\bm{\theta}_{t}$
could describe a local level, regression coefficients, or seasonal components,
among other features. Common choices of the $n\times p$ observation matrix
$F_{t}$ include the identity $(n=p)$, a matrix of indicators that selects
certain elements of $\bm{\theta}_{t}$ (e.g., for structural time series
models), and covariate values (e.g., for dynamic regression analysis). The
$p\times p$ state evolution matrix $\bm{G}_{t}$ is frequently the identity but
can be more complex, for example to capture seasonality. The observation and
evolution covariance matrices are $\bm{V}_{t}$ ($n\times n$) and $\bm{W}_{t}$
($p\times p$), respectively. Taken together, the quadruple
$\\{\bm{F}_{t},\bm{G}_{t},\bm{V}_{t},\bm{W}_{t}\\}_{t=1}^{T}$ defines the DLM.
Often, these matrices will be time-invariant.

The Gaussian DLM (3)–(4) is data-incoherent for discrete time series: the
forecasting distribution does not match the support of the data. DGLMs offer
one resolution by replacing the Gaussian observation equation with an
exponential family distribution. However, the primary option for count data is
the Poisson distribution, which is often inadequate and requires additional
modeling layers for common discrete data features such as zero-inflation,
over/underdispersion, boundedness or censoring, and heaping. These additional
layers introduce significant modeling and computational complexity. By
comparisons, the warpDLM is capable of modeling _each_ of these distributional
features under default specifications via the warping (rounding and
transformation) operation, yet maintains the useful and familiar state space
formulation through the latent DLM. By leveraging Gaussian state space models,
the warpDLM builds on the long history of theoretical and computational tools
(West and Harrison,, 1997; Petris et al.,, 2009; Prado et al.,, 2021), and
operates within a familiar setting for practitioners.

The warpDLM framework links count data $\bm{y}_{t}$ with a (latent) Gaussian
DLM for $\bm{z}_{t}$ in (3)–(4) via the warping operation (1). The rounding
operation $h:\mathcal{T}\rightarrow\mathbb{N}^{n}$ serves as the connection
mechanism between the real-valued latent space and the non-negative integers,
setting $h(\bm{y}^{*})=\bm{j}$ for any $\bm{y}^{*}\in\mathcal{A}_{j}$. These
pre-image sets $\mathcal{A}_{j}$ form a disjoint partition of the space
$\mathcal{T}$. In the univariate setting, $\mathcal{A}_{j}=[a_{j},a_{j+1})$ is
simply an interval. The warpDLM likelihood can thus be written as

$\mathbb{P}(\bm{y}_{t}=j|\bm{\theta})=\mathbb{P}\\{\bm{z}_{t}\in
g(\mathcal{A}_{j})|\bm{\theta}\\},\quad t=1,\ldots,T$ (5)

for $j\in\mathbb{N}^{n}$, where $\bm{z}_{t}\in g(\mathcal{A}_{j})$ is defined
elementwise when $\bm{y}_{t}$ is multivariate.

Both components of the warping operation serve important purposes that lead
directly to desirable model properties. The rounding function ensures that the
warpDLM has the correct support for the (possibly bounded or censored) count
data. For simplicity, suppose $n=1$; generalizations occur by applying these
specifications elementwise. By default, we take the rounding function to be
the floor function, so $\mathcal{A}_{j}=[j,j+1)$. In addition, we include the
zero modification $g(\mathcal{A}_{0})=(-\infty,0)$ so that $y_{t}=0$ whenever
$z_{t}<0$. This specification maps much of the latent space to zero, with
persistence of zeros determined by the DLM (3)–(4), for example,
$\mathbb{P}(y_{t}=0|y_{t-1}=0)=\mathbb{P}(z_{t}<0|z_{t-1}<0)$. Similarly, when
there is a known upper bound $y_{max}$ due to natural bounds or censoring, we
may simply set $\mathcal{A}_{y_{max}}=[y_{\max},\infty)$ so that the warpDLM
has the correct support, $\mathbb{P}(y_{t}\leq y_{max}|\bm{\theta})=1$. These
useful rounding operation properties are formalized in Kowal and Canale,
(2020) and Kowal and Wu, (2021). Importantly, the rounding operator does not
require any modification to the computing algorithms: once it is specified,
inference proceeds the exact same way for all choices of $h$.

While the rounding operation matches the discreteness and support of the data,
the transformation enables (nonparametric) distributional flexibility. We
apply the transformation elementwise,
$g(\bm{y}_{t}^{*})=(g_{1}(y_{1,t}^{*}),\ldots,g_{n}(y_{n,t}^{*}))^{\prime}$,
and again present the $n=1$ case for simplicity. The only requirement of the
transformation $g$ is that it be strictly monotonic, which preserves ordering
in the latent data space and ensures an inverse exists. We present both
_parametric_ and _nonparametric_ modeling strategies. Parametric examples
include classical transformations of count data, such as logarithmic, square-
root, or identity transformations, and introduce no additional parameters into
the model.

The nonparametric strategy uses a flexible and data-driven approach to infer
the transformation based on the marginal distribution of each component of
$\bm{y}$. Specifically, let $\mathcal{A}_{j}=[a_{j},a_{j+1})$ as above and
consider $n=1$. The cumulative distribution function (CDF) of $y$ and $z$ are
linked via $F_{y}(j)=F_{z}\\{g(a_{j+1})\\},$ which suggests the transformation

$\hat{g}_{0}(a_{j+1})=\bar{y}+\hat{s}_{y}\Phi^{-1}\\{\tilde{F}_{y}(j)\\},$ (6)

where $\tilde{F}_{y}$ is an estimate or model for $F_{y}$ and $\bar{y}$ and
$\hat{s}_{y}$ are the sample mean and sample standard deviation, respectively,
of $\\{y_{t}\\}_{t=1}^{T}$, to match the marginal moments of $y$ and $z$
(Kowal and Wu,, 2021). We smoothly interpolate
$(y_{t},\hat{g}_{0}(a_{y_{t}+1}))$ using a monotonic spline, which ensures
that the warpDLM is supported on $\mathbb{N}$ instead of only the observed
data values. Hence, a (nonparametric) model for $g$ may be equivalently
specified by a (nonparametric) model for $F_{y}$. We adopt the (rescaled)
empirical CDF $\tilde{F}_{y}(j)=(T+1)^{-1}\sum_{t=1}^{T}\mathbb{I}(y_{t}\leq
j)$, which implies a _semiparametric_ model for the warpDLM. Many other models
for $F_{y}$ are compatible within the warpDLM, including Bayesian
nonparametric models and parametric distributions (e.g., Poisson or Negative
Binomial marginals). By design, the model for the transformation $g$ is
decoupled from both the rounding operator $h$ and the DLM (3)–(4), so the
subsequent derivations and algorithms require only trivial modifications for
distinct choices of $g$.

The nonparametric transformation serves as a reasonable default across a
variety of scenarios, as suggested by our results in Section 5 and further
confirmed in other (non-time series) settings (e.g. Kowal and Wu,, 2021). To
compare models with different transformations, there are a variety of choices.
In the simulation study of Section 5.2, we select the best model based on out-
of-sample forecasting performance using leave-future-out cross-validation.
However, in practice, any Bayesian model comparison metric could be used, such
as the widely applicable or Watanabe-Akaike information criterion (WAIC,
Watanabe,, 2010) or approaches more specifically tailored to time series
analysis (Bürkner et al.,, 2020).

##  3 Exact Filtering and Smoothing

In this section, we derive the recursive updates for the warpDLM filtering,
smoothing, and forecasting distributions. Currently, these _exact_ results
stemming from a _coherent_ joint distribution (1)–(2) are unique among state
space models for multivariate count (or rounded) data, and crucially enable
MCMC-free inference and forecasting.

###  3.1 Selection Distributions and the warpDLM

Consider the first time step $t=1$ of the warpDLM. Here, we omit the time
subscripts for simplicity. The latent data (2) are described by the two DLM
equations (3)–(4), which can be rewritten as a single equation

$\bm{z}=\bm{F}\bm{\theta}+\bm{v},\quad\bm{v}\sim N_{n}(\bm{0},\bm{V})$ (7)

with $\bm{\theta}\coloneqq\bm{\theta}_{1}$ and the prior $\bm{\theta}\sim
N_{p}(\bm{\mu}_{\theta}=\bm{G}a_{0},\bm{\Sigma}_{\theta}=\bm{G}\bm{R}_{0}\bm{G}^{\prime}+\bm{W})$.
Using the likelihood (5), the posterior distribution is

$p(\bm{\theta}|\bm{y})=p(\bm{\theta}|\bm{z}\in\mathcal{C})=\frac{p(\bm{\theta})p\\{\bm{z}\in\mathcal{C}|\bm{\theta}\\}}{p\\{\bm{z}\in\mathcal{C}\\}}$
(8)

which is also the first-step filtering distribution. Within the warpDLM,
conditioning on $\bm{y}$ is equivalent to conditioning on $\bm{z}$ belonging
to a set $\mathcal{C}=g(\mathcal{A}_{\bm{y}})$. Although the representation is
general, this set is typically simple: the default floor operator for $h$
implies that
$\mathcal{A}_{\bm{y}}=[y_{1},y_{1}+1)\times\cdots\times[y_{n},y_{n}+1)$, so
that $\bm{z}\in\mathcal{C}$ implies that each element of $\bm{z}$ belongs to a
(transformed) interval.

A distribution of the form (8) is known as a selection distribution (Arellano-
Valle et al.,, 2006). When the two variables $\bm{\theta}$ and $\bm{z}$ are
jointly normal—as in the warpDLM—the resulting distribution is a _selection
normal_ (SLCT-N). More formally, given the joint distribution

$\begin{pmatrix}\bm{z}\\\ \bm{\theta}\end{pmatrix}\sim
N_{n+p}\left\\{\begin{pmatrix}\bm{\mu}_{z}\\\
\bm{\mu}_{\theta}\end{pmatrix},\begin{pmatrix}\bm{\Sigma}_{z}&\bm{\Sigma}_{z\theta}\\\
\bm{\Sigma}_{z\theta}^{\prime}&\bm{\Sigma}_{\theta}\end{pmatrix}\right\\}$

we denote the conditional random variable
$[\bm{\theta}|\bm{z}\in\mathcal{C}]\sim\mbox{SLCT-N}_{n,p}(\bm{\mu}_{z},\bm{\mu}_{\theta},\bm{\Sigma}_{z},\bm{\Sigma}_{\theta},\bm{\Sigma}_{z\theta},\mathcal{C})$
for constraint region $\mathcal{C}$. This random variable has density

$p(\bm{\theta}|\bm{z}\in\mathcal{C})=\phi_{p}(\bm{\theta};\bm{\mu}_{\theta},\bm{\Sigma}_{\theta})\frac{\bar{\Phi}_{n}(\mathcal{C};\bm{\Sigma}_{z\theta}\bm{\Sigma}_{\theta}^{-1}(\bm{\theta}-\bm{\mu}_{\theta})+\bm{\mu}_{z},\bm{\Sigma}_{z}-\bm{\Sigma}_{z\theta}\bm{\Sigma}_{\theta}^{-1}\bm{\Sigma}_{z\theta}^{\prime})}{\bar{\Phi}_{n}(\mathcal{C};\bm{\mu}_{z},\bm{\Sigma}_{z})}$
(9)

where $\phi_{p}(\cdot;\bm{\mu},\bm{\Sigma})$ denotes the Gaussian density
function of a Gaussian random variable with mean $\bm{\mu}$ and covariance
$\bm{\Sigma}$ and
$\bar{\Phi}_{n}(\mathcal{C};\bm{\mu},\bm{\Sigma})=\int_{\mathcal{C}}\phi_{n}(\bm{x};\bm{\mu},\bm{\Sigma})d\bm{x}$.
The density (9) is somewhat unwieldy in practice, but there is also a
constructive representation which allows for direct Monte Carlo simulation
from the posterior density $p(\bm{\theta}|\bm{y})$ (see Section 4).

For the first-step model (7) with prior $\bm{\theta}\sim
N_{p}(\bm{\mu}_{\theta},\bm{\Sigma}_{\theta})$, we report the exact posterior
distribution:

######  Theorem 1.

Under (7), the posterior distribution is
$[\bm{\theta}|\bm{y}]\sim\mbox{SLCT-N}_{n,p}(\bm{\mu}_{z}=\bm{F}\bm{\mu}_{\theta},\bm{\mu}_{\theta},\bm{\Sigma}_{z}=\bm{F}\bm{\Sigma}_{\theta}\bm{F}^{\prime}+\bm{V},\bm{\Sigma}_{\theta},\bm{\Sigma}_{z\theta}=\bm{F}\bm{\Sigma}_{\theta},\mathcal{C}=g(\mathcal{A}_{\bm{y}}))$.

The Gaussian prior in Theorem 1 is conjugate: Gaussian distributions can be
viewed as a limiting case of a $\mbox{SLCT-N}_{n,p}$ where $n=0$. However, the
conjugacy suggests a more general prior for $\bm{\theta}$, namely, the SLCT-N
distribution:

######  Lemma 1.

Consider the prior
$\bm{\theta}\sim\mbox{SLCT-N}_{n_{0},p}(\bm{\mu}_{z_{0}},\bm{\mu}_{\theta},\bm{\Sigma}_{z_{0}},\bm{\Sigma}_{\theta},\bm{\Sigma}_{z_{0}\theta},\mathcal{C}_{0})$
with the latent-data observation equation (7). The posterior is

$\displaystyle[\bm{\theta}|\bm{y}]\sim\mbox{SLCT-N}_{n_{0}+n,p}\Big{\\{}$
$\displaystyle\bm{\mu}_{z_{1}}=\begin{pmatrix}\bm{\mu}_{z_{0}}\\\
\bm{F}\bm{\mu}_{\theta}\end{pmatrix},\bm{\mu}_{\theta},\bm{\Sigma}_{z_{1}}=\begin{pmatrix}\bm{\Sigma}_{z_{0}}&\bm{\Sigma}_{z_{0}\theta}\bm{F}^{\prime}\\\
\bm{F}\bm{\Sigma}_{z_{0}\theta}^{\prime}&\bm{F}\bm{\Sigma}_{\theta}\bm{F}^{\prime}+\bm{V}\end{pmatrix},\bm{\Sigma}_{\theta},$
$\displaystyle\bm{\Sigma}_{z_{1}\theta}=\begin{pmatrix}\bm{\Sigma}_{z_{0}\theta}\\\
\bm{F}\bm{\Sigma}_{\theta}\end{pmatrix},\mathcal{C}_{1}=\mathcal{C}_{0}\times
g(\mathcal{A}_{\bm{y}})\Big{\\}}.$

Thus the selection normal distribution is conjugate with the warpDLM
likelihood. Moreover, since the Gaussian distribution is closed under linear
transformations, this property is preserved for the selection normal
distribution (Arellano-Valle et al.,, 2006). These results link the selection
normal distribution with the warpDLM prior-to-posterior updating mechanism,
and provide the key building blocks for deriving the warpDLM filtering and
smoothing distributions.

###  3.2 Filtering and Smoothing for Warped DLMs

We now proceed to the most general setting for warpDLM, using the conjugacy
and closure under linear transformation results to derive the exact forms of
the filtering, smoothing, and forecasting distributions. Throughout, we assume
the quadruple $\\{\bm{F}_{t},\bm{G}_{t},\bm{V}_{t},\bm{W}_{t}\\}$ for
$t=1,\ldots,T$ is known. Estimation of variance parameters is addressed in
Section 4.

####  3.2.1 Filtering

Adding back the appropriate subscripts to Theorem 1, the first-step filtering
distribution is
$(\bm{\theta}_{1}|\bm{y}_{1})\sim\mbox{SLCT-N}_{n,p}(\bm{\mu}_{z_{1}}=\bm{F}_{1}\bm{a}_{1},\bm{\mu}_{\theta}=\bm{a}_{1},\bm{\Sigma}_{z_{1}}=\bm{F}_{1}\bm{R}_{1}\bm{F}_{1}^{\prime}+\bm{V}_{1},\bm{\Sigma}_{\theta}=\bm{R}_{1},\bm{\Sigma}_{z_{1}\theta}=\bm{F}_{1}\bm{R}_{1},\bm{C}_{1}=g(\mathcal{A}_{\bm{y}_{1}})),$
where $\bm{a}_{1}=\bm{G}_{1}\bm{a}_{0}$ and
$\bm{R}_{1}=\bm{G}_{1}\bm{R}_{0}\bm{G}_{1}^{\prime}+\bm{W}_{1}$. Given the
first-step filtering distribution, we proceed inductively for time $t$:

######  Theorem 2.

Let
$(\bm{\theta}_{t-1}|\bm{y}_{1:t-1})\sim\mbox{SLCT-N}_{n(t-1),p}(\bm{\mu}_{z_{t-1|t-1}},\bm{\mu}_{\theta_{t-1|t-1}},\bm{\Sigma}_{z_{t-1|t-1}},\bm{\Sigma}_{\theta_{t-1|t-1}},\\\
\bm{\Sigma}_{(z\theta)_{t-1|t-1}},C_{t-1|t-1})$ be the filtering distribution
at time $t-1$ under the warpDLM. Then the one-step-ahead state predictive
distribution at $t$ is

$(\bm{\theta}_{t}|\bm{y}_{1:t-1})\sim\mbox{SLCT-N}_{n(t-1),p}(\bm{\mu}_{z_{t|t-1}},\bm{\mu}_{\theta_{t|t-1}},\bm{\Sigma}_{z_{t|t-1}},\bm{\Sigma}_{\theta_{t|t-1}},\bm{\Sigma}_{(z\theta)_{t|t-1}},C_{t|t-1})$
(10)

with $\bm{\mu}_{z_{t|t-1}}=\bm{\mu}_{z_{t-1|t-1}}$,
$\bm{\mu}_{\theta_{t|t-1}}=\bm{G}_{t}\bm{\mu}_{\theta_{t-1|t-1}}$,
$\bm{\Sigma}_{z_{t|t-1}}=\bm{\Sigma}_{z_{t-1|t-1}}$,
$\bm{\Sigma}_{\theta_{t|t-1}}=\bm{G}_{t}\bm{\Sigma}_{\theta_{t-1|t-1}}\bm{G}_{t}^{\prime}+\bm{W}_{t}$,
$\bm{\Sigma}_{(z\theta)_{t|t-1}}=\bm{\Sigma}_{(z\theta)_{t-1|t-1}}\bm{G}_{t}^{\prime}$,
and $C_{t|t-1}=C_{t-1|t-1}$. Furthermore, the filtering distribution at time
$t$ is

$(\bm{\theta}_{t}|\bm{y}_{1:t})\sim\mbox{SLCT-N}_{n(t),p}(\bm{\mu}_{z_{t|t}},\bm{\mu}_{\theta_{t|t}},\bm{\Sigma}_{z_{t|t}},\bm{\Sigma}_{\theta_{t|t}},\bm{\Sigma}_{(z\theta)_{t|t}},C_{t|t})$
(11)

with $\bm{\mu}_{z_{t|t}}=\begin{pmatrix}\bm{\mu}_{z_{t|t-1}}\\\
\bm{F}_{t}\bm{\mu}_{\theta_{t|t-1}}\end{pmatrix}$,
$\bm{\mu}_{\theta_{t|t}}=\bm{\mu}_{\theta_{t|t-1}}$,
$\bm{\Sigma}_{z_{t|t}}=\begin{pmatrix}\bm{\Sigma}_{z_{t|t-1}}&\bm{\Sigma}_{(z\theta)_{t|t-1}}\bm{F}_{t}^{\prime}\\\
\bm{F}_{t}\bm{\Sigma}_{(z\theta)_{t|t-1}}^{\prime}&\bm{F}_{t}\bm{\Sigma}_{\theta_{t|t-1}}\bm{F}_{t}^{\prime}+\bm{V}_{t}\end{pmatrix}$,
$\bm{\Sigma}_{\theta_{t|t}}=\bm{\Sigma}_{\theta_{t|t-1}}$,
$\bm{\Sigma}_{(z\theta)_{t|t}}=\begin{pmatrix}\bm{\Sigma}_{(z\theta)_{t|t-1}}\\\
\bm{F}_{t}\bm{\Sigma}_{\theta_{t|t-1}}\end{pmatrix}$ , and
$C_{t|t}=C_{t|t-1}\times g(\mathcal{A}_{\bm{y}_{t}})$.

The state predictive distribution updating in (10) is analogous to that of the
Kalman filter, but the filtering distribution update in (11) has a different
form. In the Kalman filter, the filtering distribution parameters are updated
by taking the previous best estimate for the states (from the state predictive
distribution) and “correcting” them using information gained from the new
observation. Under the warpDLM, the previous parameters are not corrected;
instead, the dimension of the new filtering distribution increases, with the
latest data point controlling the bounds over which the latent state can vary.
The growing dimension ensures we are retaining all the information about the
states, although it also means that sampling can become demanding for longer
time series (see Section 4 for discussion). The selection region $C$ also
grows in size with each update, but from a storage perspective this is not
burdensome: the region $C_{t|t}$ adds an interval for each component of
$\bm{y}_{t}$, which simply requires appending a row to two $(t-1)\times n$
matrices. Importantly, the parameter updates depend on the new observation
only through $C$, which implies that the system matrices do not need to be
updated recursively. This point is better clarified by considering the
smoothing distribution.

####  3.2.2 Smoothing

The joint smoothing distribution is $p(\bm{\theta}_{1:T}|\bm{y}_{1:T})$ where
$T$ is the terminal time point. Most commonly, the smoothing distribution in a
DLM is obtained by first filtering “forward” recursively and then smoothing
“‘backward” in time (Rauch et al.,, 1965). This strategy targets the marginal
smoothing distributions $p(\bm{\theta}_{t}|\bm{y}_{1:T})$ and requires
modifications to access the _joint_ smoothing distribution (e.g., Durbin and
Koopman,, 2002).

Within the warpDLM, we instead target the joint smoothing distribution
directly and analytically. Crucially, this process does not require any
preliminary passes through the data: all smoothing parameters can be
constructed _a priori_ through matrix multiplications. Let
$\bm{G}_{1:t}=\bm{G}_{t}\bm{G}_{t-1}\cdots\bm{G}_{1}$ and
$\bm{\mu}_{\theta}=(\bm{G}_{1:1}\bm{a}_{0},\dotsi,\bm{G}_{1:T}\bm{a}_{0})$, a
$pT$-dimensional vector. Now, let
$\bm{R}_{t}=\bm{G}_{t}\bm{R}_{t-1}\bm{G}_{t}^{\prime}+\bm{W}_{t}$ for
$t=1,\ldots,T$. Then the covariance matrix of $\bm{\theta}$ is a $pT\times pT$
matrix with $p\times p$ diagonal block entries of
$\bm{\Sigma}_{\theta}[t,t]=\bm{R}_{t}=\bm{G}_{1:t}\bm{R}_{0}\bm{G}_{1:t}^{\prime}+\sum_{q=2}^{t}\bm{G}_{q:t}\bm{W}_{q-1}\bm{G}_{q:t}^{\prime}+\bm{W}_{t}$
and cross covariance entries (for $t>q$)
$\bm{\Sigma}_{\theta}[t,q]={\bm{\Sigma}_{\theta}[q,t]}^{\top}=\bm{G}_{(q+1):t}\bm{\Sigma}_{\theta}[q,q].$
We also define two block diagonal matrices, $\bm{\mathfrak{F}}$ and
$\bm{\mathfrak{V}}$, with diagonal entries of $\bm{F}_{t}$ and $\bm{V}_{t}$,
respectively, for $t=1,\dots,n$, so $\bm{\mathfrak{F}}$ is a $nT\times pT$
matrix and $\bm{\mathfrak{V}}$ is a $nT\times nT$ matrix. Finally, let
$\mathcal{C}_{1:T}=(g(\mathcal{A}_{\bm{y}_{1}}),\dotsi,g(\mathcal{A}_{\bm{y}_{T}}))$
be the matrix of selection intervals through time $T$. The joint smoothing
distribution is given below.

######  Theorem 3.

Under the warpDLM, the joint smoothing distribution is

$(\bm{\theta}_{1:T}|\bm{y}_{1:T})\sim\mbox{SLCT-N}_{nT,pT}\\{\bm{\mu}_{z}{=}\bm{\mathfrak{F}}\bm{\mu}_{\theta},\bm{\mu}_{\theta},\bm{\Sigma}_{z}{=}\bm{\mathfrak{V}}+\bm{\mathfrak{F}}\bm{\Sigma}_{\theta}\bm{\mathfrak{F}}^{\prime},\bm{\Sigma}_{\theta},\bm{\Sigma}_{z\theta}{=}\bm{\mathfrak{F}}\bm{\Sigma}_{\theta},\mathcal{C}{=}\mathcal{C}_{1:T}\\}$
(12)

These results can be applied to the time $s$ filtering distribution
$p(\bm{\theta}_{1:s}|\bm{y}_{1:s})$ to provide a more concise representation.
Using the analogous smoothing parameters but applied up to time $s$ instead of
$T$, we can rewrite the parameters of the filtering distribution as
$\bm{\mu}_{z_{s|s}}=\bm{\mathfrak{F}}\bm{\mu}_{\theta}\;,\;\bm{\mu}_{\theta_{s|s}}=\bm{G}_{1}^{s}\bm{a}_{0}\;,\;\bm{\Sigma}_{z_{s|s}}=\bm{\mathfrak{V}}+\bm{\mathfrak{F}}\bm{\Sigma}_{\theta}\bm{\mathfrak{F}}^{\prime}\;,\;\bm{\Sigma}_{(z\theta)_{s|s}}=\bm{\mathfrak{F}}\bm{\Sigma}_{\theta}\;,\;\bm{\Sigma}_{\theta_{s|s}}=\bm{R}_{s}.$
Although these covariance terms are quite complex in terms of notation,
constructing such matrices is computationally straightforward via matrix
multiplications and additions.

The marginal smoothing distribution $p(\bm{\theta}_{t}|\bm{y}_{1:T})$ at each
time $t$ is readily available from Theorem 3. In particular, the SLCT-N family
is closed under marginalization, so parameters for the observations stay the
same as in the joint case, and we simply pick off the correct block of the
state parameters. This is formalized in the corollary below.

######  Corollary 1.

Under the warpDLM, the marginal smoothing distribution at time $t$ is

$\displaystyle(\bm{\theta}_{t}|\bm{y}_{1:T})\sim\mbox{SLCT-N}_{nT,p}$
$\displaystyle\\{\bm{\mu}_{z_{t|T}}=\bm{\mathfrak{F}}\bm{\mu}_{\theta},\bm{\mu}_{\theta_{t|T}}=\bm{\mu}_{\theta}[t],\bm{\Sigma}_{z_{t|T}}=\bm{\mathfrak{V}}+\bm{\mathfrak{F}}\bm{\Sigma}_{\theta}\bm{\mathfrak{F}}^{\prime},\bm{\Sigma}_{\theta_{t|T}}=\bm{\Sigma}_{\theta}[t,t],$
(13)
$\displaystyle\bm{\Sigma}_{(z\theta)_{t|T}}=\bm{\Sigma}_{z\theta}[,t],\mathcal{C}=\mathcal{C}_{1:T}\\}\;.$

where $\bm{\Sigma}_{z\theta}[,t]$ refers to the $t$-th block of $p$ columns in
$\bm{\Sigma}_{z\theta}$.

The joint smoothing distribution also enables direct computation of the
warpDLM marginal likelihood:

######  Corollary 2.

Under the warpDLM, the marginal likelihood is

$\displaystyle
p(\bm{y}_{1:T})=\bar{\Phi}_{nT}(\mathcal{C}_{1:T};\bm{\mu}_{z}=\bm{\mathfrak{F}}\bm{\mu}_{\theta},\bm{\Sigma}_{z}=\bm{\mathfrak{V}}+\bm{\mathfrak{F}}\bm{\Sigma}_{\theta}\bm{\mathfrak{F}}^{\prime})\;.$
(14)

The marginal likelihood is useful for marginal maximum likelihood estimation
of the variance parameters in (3)–(4) and other model comparison metrics.

####  3.2.3 Forecasting

A primary advantage of the warpDLM is the availability of discrete and
correctly-supported forecasting distributions. Let $\bm{\mu}_{z_{t+1|t+1}}$
and $\bm{\Sigma}_{z_{t+1|t+1}}$ denote parameters of $p(\bm{y}_{1:(t+1)})$
analogous to (14), and let $\mathcal{C}_{1:t}$ be the vector of selection
intervals up to the known time $t$. Then the one step forecasting distribution
is

$p(\bm{y}_{t+1}|\bm{y}_{1:t})=\frac{p(\bm{y}_{1:(t+1)})}{p(\bm{y}_{1:t})}=\frac{\bar{\Phi}_{n(t+1)}\\{\mathcal{C}_{1:t}\times
g(\mathcal{A}_{y_{t+1}});\bm{\mu}_{z_{t+1|t+1}},\bm{\Sigma}_{z_{t+1|t+1}}\\}}{\bar{\Phi}_{nt}(\mathcal{C}_{1:t};\bm{\mu}_{z_{t|t}},\bm{\Sigma}_{z_{t|t}})}\;.$
(15)

The $h$-step ahead forecasting distribution can be defined similarly. In most
applications, the forecasting distribution is defined over all non-negative
integers, and thus evaluating it at every point is not possible. However, in
practice most time series of counts take values in a small range (especially
locally), and thus the forecasting distribution only has significant mass on a
small range of values. For short to medium time series applications,
evaluating these probabilities is quite fast and does not pose a computational
burden.

More directly, it is straightforward to simulate from the forecasting
distribution given draws from the filtering distribution. Specifically, let
$\bm{\theta}_{t}^{*}\sim p(\bm{\theta}_{t}|\bm{y}_{1:t})$ denote a draw from
the filtering distribution. By passing this draw through the DLM, the
evolution equation (4) samples $\bm{\theta}_{t+1}^{**}\sim
N_{p}(\bm{G}_{t+1}\bm{\theta}_{t}^{*},\bm{W}_{t+1})$ and the observation
equation (3) samples $\bm{z}_{t+1}^{*}\sim
N_{n}(\bm{F}_{t}\bm{\theta}_{t+1}^{**},\bm{V}_{t+1})$, which equivalently
yields a draw from $(\bm{z}_{t+1}|\bm{y}_{1:t})$. By retaining
$\bm{\tilde{y}}_{t+1}=h\circ g^{-1}(\bm{z}_{t+1}^{*})$, we obtain a draw from
the one-step forecasting distribution of the warpDLM; multi-step forecasting
proceeds similarly. Hence, the primary computational cost is the draw from the
filtering distribution, which is detailed in Section 4.

####  3.2.4 Missing Data

The warpDLM seamlessly handles missing data in any component of $\bm{y}_{t}$
at any time $t$. Notably, the filtering, smoothing, and forecasting
distributions depend on $\bm{y}_{t}$ only through the implied constraint
region (i.e., the collection of intervals). Hence, the absence of an
observation $y_{jt}$ implies that the interval, such as
$[g(y_{jt}),g(y_{jt}+1))$, is simply replaced by the interval
$(-\infty,\infty)$: there is no information about $y_{jt}$ and therefore no
constraint for the latent $z_{jt}$. As a result, the modifications for the
filtering, smoothing, and forecasting distributions in the presence of
missingness are trivial both analytically and computationally. This result
also suggests that the warpDLM is readily applicable to discrete time series
data observed at mixed frequencies.

##  4 Computing

We develop computing strategies for offline and online inference and
forecasting. For offline inference, we propose direct Monte Carlo sampling
from the relevant filtering or smoothing distributions. We also introduce a
Gibbs sampler for offline inference, which decreases raw computing time for
large $T$, but sacrifices some Monte Carlo efficiency due to the
autocorrelated samples. Online inference is enabled by an optimal particle
filter. The _Monte Carlo_ sampler and the _optimal_ particle filter are
uniquely available as a consequence of the results in Section 3.

###  4.1 Direct Monte Carlo Sampling

The analytic filtering and smoothing distributions for the warpDLM unlock the
potential for direct Monte Carlo sampling. This task requires sampling from a
selection normal distribution. Crucially, the selection normal distribution
admits a constructive representation using a multivariate normal distribution
and a multivariate truncated normal distribution (Arellano-Valle et al.,,
2006), which is converted to a sampler in Algorithm 1. We then apply Algorithm
1 to obtain Monte Carlo draws from the warpDLM state predictive distribution
(10), filtering distribution (11), joint smoothing distribution (12), marginal
smoothing distribution (13), or forecasting distribution (Section 3.2.3). This
computing strategy avoids the need for MCMC, which often requires lengthy
simulation runs and various diagnostics.

Result: One sample $\bm{\theta}$ from SLCT-N distribution

Given
$[\bm{\theta}]\sim\mbox{SLCT-N}_{d_{1},d_{2}}(\bm{\mu}_{z},\bm{\mu}_{\theta},\bm{\Sigma}_{z},\bm{\Sigma}_{\theta},\bm{\Sigma}_{z\theta},\mathcal{C})$;

  1. 1.

Sample truncated multivariate normal: Sample $\bm{V}_{0}$ from
$N_{d_{1}}(\bm{0},\bm{\Sigma}_{Z})$ truncated to region
$\mathcal{C}-\bm{\mu}_{z}$

  2. 2.

Sample multivariate normal: Sample $\bm{V}_{1}$ from
$N_{d_{2}}(\bm{0},\bm{\Sigma}_{\theta}-\bm{\Sigma}_{z\theta}^{\prime}\bm{\Sigma}_{z}^{-1}\bm{\Sigma}_{z\theta})$

  3. 3.

Combine results: Compute
$\bm{\theta}=\bm{\mu}_{\theta}+\bm{V}_{1}+\bm{\Sigma}_{z\theta}^{\prime}\bm{\Sigma}_{z}^{-1}\bm{V}_{0}$

Algorithm 1 Sampling from Selection Normal Distribution

These distributions treat the variance components $\bm{V}_{t}$ and
$\bm{W}_{t}$ as fixed and known. Unknown variance components may be
accommodated using point estimates (e.g., marginal maximum likelihood
estimators using (14)) or within a blocked Gibbs sampler that iterates draws
of the state parameters (Algorithm 1) and the variance components.

The computational bottleneck of the Monte Carlo sampler is the draw (a) from
the multivariate truncated normal. Efficient sampling from multivariate
truncated normal distributions is an active area of research; we apply the
minimax tilting strategy from Botev, (2017), which is implemented in the R
package `TruncatedNormal`. This method is highly efficient and accurate for
dimensions up to about $d_{1}\approx 500$. For the warpDLM, the dimension
$d_{1}$ corresponds to $nt$ for the filtering distribution at time $t$ or $nT$
for the smoothing distributions. Hence, when the combined data dimension
($nT$) is very large, other computational approaches are needed; we provide
these in the subsequent sections.

###  4.2 Gibbs sampling

For lengthy time series, we develop a Gibbs sampler for offline inference that
circumvents the computational bottlenecks of direct Monte Carlo sampling. The
Gibbs sampler iterates between a _latent data augmentation_ draw and a _state
parameter_ draw (plus the variance components $\bm{\psi}$, if unknown). This
is represented in Algorithm 2. Although not shown, one can also easily sample
from the posterior predictive and forecast distributions for $\bm{z}_{t}$ in
each pass. Draws from the corresponding distributions for $\bm{y}_{t}$ are
given by computing $h\left(g^{-1}(\bm{z}_{t})\right)$.

  1. 1.

Sample the latent data: draw
$[\bm{z}_{t}|\\{\bm{y}_{t},\bm{\theta}_{t}\\}_{t=1}^{T},\psi]$ from
$N_{n}(\bm{F}_{t}\bm{\theta}_{t},\bm{V}_{t})$ truncated to $\bm{z}_{t}\in
g(\mathcal{A}_{y_{t}})$ for $t=1,\ldots,T$

  2. 2.

Sample the states: draw $[\bm{\theta}_{1:T}|\bm{z}_{1:T},\bm{\psi}]$

  3. 3.

Sample the parameters: draw
$[\bm{\psi}|\\{\bm{z}_{t},\bm{\theta}_{t}\\}_{t=1}^{T}]$ (from the appropriate
full conditionals depending on priors)

Algorithm 2 Gibbs Sampler for WarpDLM

Crucially, the first step of latent data augmentation can be decomposed into
$T$ _independent_ draws from an $n$-dimensional truncated normal distribution;
hence the Gibbs sampling approach doesn’t suffer from the same dimensionality
problems as the direct sampling method. Furthermore, the latent state
parameter full conditional in step 2 is exactly the smoothing distribution for
a Gaussian DLM, and we can thus rely on already-developed Kalman filter-based
simulation techniques such as forward filtering backward sampling (FFBS;
Frühwirth‐Schnatter,, 1994; Carter and Kohn,, 1994) or the simulation smoother
(Durbin and Koopman,, 2002). Importantly, these sampling steps are _linear_ in
the length of the time series $T$. Additionally, both methods have R
implementations which can be readily leveraged: the package `dlm` (Petris,,
2010) performs FFBS and the `KFAS` package (Helske,, 2017) applies the
simulation smoother.

In a fully Bayesian framework, we place priors on all variances, and sample
from the appropriate full conditionals in step 3. For univariate models, we
put a Uniform$(0,A)$ prior with $A$ large on all standard deviations, leading
to relatively simple Gibbs updates. In the multivariate context, we place
independent inverse-Wishart priors on variance matrices as in Petris et al.,
(2009).

The sampler for the warpDLM framework is similar to that developed in Kowal
and Canale, (2020) for the static count regression scenario. The Bayesian
dynamic probit sampler (Albert and Chib,, 1993; Fasano et al.,, 2021) can be
regarded as a special case. In particular, if our data $y_{t}$ were binary, we
could let our transformation $g$ be the identity while defining our rounding
operator such that $\mathcal{A}_{y_{t}=0}=(-\infty,0)$ and
$\mathcal{A}_{y_{t}=1}=(0,\infty)$, thus recovering the probit Gibbs sampling
algorithm.

Overall, the Gibbs sampler provides scalable offline inference for the warpDLM
in the case of large $T$, although at the expense of Monte Carlo efficiency;
naturally, the direct Monte Carlo sampler (Algorithm 1) is preferred when the
data dimension permits.

###  4.3 Optimal Particle Filtering

We design an _optimal_ particle filter for online inference with the warpDLM.
The key sampling steps do _not_ depend on $T$, which provides scalability for
lengthy time series. Online particle filtering algorithms for state space
models usually apply Sequential Importance Sampling with or without resampling
(Doucet et al.,, 2000). For each time point, sample trajectories or particles
are drawn from an importance function, and each particle has an associated
importance weight. For the resampling step, draws from the filtering
distribution are obtained via a weighted resampling from the previously drawn
particles. The purpose of this step is to slow the effect of particle
degeneration.

A well-known drawback of particle filtering methods is that the variance of
the importance weights grows over time, and thus the trajectories will
degenerate, eventually leaving only one particle with nonzero weight. To avoid
quick degeneration, it is important to choose a good importance function; the
_optimal_ function is that which minimizes the conditional variance of
weights. Applying Doucet et al., (2000) to the warpDLM, the optimal importance
function is $p(\bm{\theta}_{t}|\bm{\theta}_{t-1},\bm{y}_{t})$ with weights
given by $p(\bm{y}_{t}|\bm{\theta}_{t-1})$. For many state space models, this
function is not known analytically; however, our results from Section 3
produce an exact form for the warpDLM, formalized in the following corollary.

######  Corollary 3.

Under the warpDLM, we have

$\displaystyle(\bm{\theta}_{t}|\bm{\theta}_{t-1},\bm{y}_{t})\sim\mbox{SLCT-N}_{n,p}$
$\displaystyle(\bm{\mu}_{z}=\bm{F}_{t}\bm{G}_{t}\bm{\theta}_{t-1},\bm{\mu}_{\theta}=\bm{G}_{t}\bm{\theta}_{t-1},\bm{\Sigma}_{z}=\bm{V}_{t}+\bm{F}_{t}\bm{W}_{t}\bm{F}_{t}^{\prime},$
(16)
$\displaystyle\bm{\Sigma}_{\theta}=\bm{W}_{t},\bm{\Sigma}_{z\theta}=\bm{F}_{t}\bm{W}_{t},C=g(\mathcal{A}_{\bm{y}_{t}}))$

and

$\displaystyle
p(\bm{y}_{t}|\bm{\theta}_{t-1})=\bar{\Phi}_{n}(g(\mathcal{A}_{y_{t}});\bm{\mu}_{z}=\bm{F}_{t}\bm{G}_{t}\bm{\theta}_{t-1},\bm{\Sigma}_{z}=\bm{V}_{t}+\bm{F}_{t}\bm{W}_{t}\bm{F}_{t}^{\prime})$
(17)

We apply Corollary 3 to design an _optimal_ particle filter for the warpDLM in
Algorithm 3. Crucially, the draw from the selection normal in (16) is not
dependent on the length of the time series, which provides scalability in $T$.

Result: Draws ($\bm{\theta}_{t}^{(1)},\ldots,\bm{\theta}_{t}^{(S)}$) from
filtering distribution $p(\bm{\theta}_{t}|\bm{y}_{1:t})$ for $t=1,\ldots,T$

Given $S$ particle trajectories at time $t=0$ denoted by
$\bm{\theta}_{0}^{(s)}$ ;

for  _$t\leftarrow 1$ to $T$_ do

for  _$s\leftarrow 1$ to $S$_ do

Sample $\tilde{\bm{\theta}}_{t}^{(s)}$ from (16) conditional on
$\bm{\theta}_{t-1}=\bm{\theta}_{t-1}^{(s)}$ ;

Calculate weight up to normalizing constant using (17): $w_{t}^{(s)}\propto
p(\bm{y}_{t}|\bm{\theta}_{t-1}^{(s)})$

end for

Normalize importance weights:
$\tilde{w_{t}}^{(s)}=\frac{w_{t}^{(s)}}{\sum_{s=1}^{S}w_{t}^{(s)}}$ ;

Obtain final draws $\bm{\theta}_{t}^{(1)},\ldots,\bm{\theta}_{t}^{(S)}$ by
resampling
$\tilde{\bm{\theta}}_{t}^{(1)},\ldots,\tilde{\bm{\theta}_{t}}^{(S)}$ with
weights $\tilde{w_{t}}^{(s)}$

end for

Algorithm 3 Optimal Particle Filter for warpDLM

The (unnormalized) weights can be used to estimate the marginal likelihood
$p(\bm{y}_{1:T})$:

$\hat{p}(\bm{y}_{1:T})\coloneqq\hat{p}(\bm{y}_{1})\prod_{t=2}^{T}\hat{p}(\bm{y}_{t}|\bm{y}_{1:t-1})$
(18)

where $\hat{p}(\bm{y}_{1})=\frac{1}{S}\sum_{s=1}^{S}w_{1}^{(s)}$ and
$\hat{p}(\bm{y}_{t}|\bm{y}_{1:t-1})=\frac{1}{S}\sum_{s=1}^{S}w_{t}^{(s)}$. As
with equation (14), this can be used for model comparisons or estimation of
unknown parameters (Kantas et al.,, 2015). Such an approximation is also key
when designing more complex sequential Monte Carlo schemes or sampling based
on particle methods (Andrieu et al.,, 2010).

As both the cross-section dimension $n$ and the length of the time series $T$
increase, we expect that the quality of the particles will decline, despite
the optimality of our transition density. This effect is typically more
pronounced for the smoothing distribution than the filtering distribution. In
that case, we may seek to improve the particle filter (Algorithm 3) by
applying relevant techniques from sequential Monte Carlo sampling, such as a
mutation step, or resort to particle smoothing (Doucet et al.,, 2009) when $T$
is especially large. However, we show empirically that the proposed particle
filter is capable of accurate inference for $n=5$ and $T=4000$ (see the
supplementary material).

##  5 Warped DLMs in Action

In this section, we showcase several facets of the warpDLM framework. We first
perform a forecasting analysis on simulated data, showing that the warpDLM
offers better distributional forecasting than competing Bayesian methods. We
then apply the warpDLM to model a bivariate time series of drug overdose
counts in Cincinnati, and highlight both the offline and online computing
capabilities.

###  5.1 Evaluating Probabilistic Forecasts

A prominent advantage of the warpDLM is that it provides a forecasting
_distribution_ that matches the (discrete) support of the data. Distributional
forecasts are especially important for discrete data, whereas point forecasts
are less informative. Evaluation of probabilistic forecasts should first and
foremost verify that forecast distributions are calibrated in the sense that
observed values are consistent with samples from the distribution (Czado et
al.,, 2009). To assess calibration, we apply the randomized Probability
Integral Transform (rPIT). Let $H_{t}$ denote the true (discrete) cumulative
forecasting distribution, from which a count $y_{t}$ is observed. Each model
produces an estimate $\hat{H}_{t}$ of the forecasting distribution. The rPIT
is drawn from $\tilde{p}_{t}\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Uniform}[\hat{H}_{t}(y_{t}-1),\hat{H}_{t}(y_{t})].$ If the
model is well-calibrated, then $\hat{H}_{t}\approx H_{t}$ and consequently
$\tilde{p}_{t}\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Uniform}(0,1)$. Thus,
forecasts are evaluated by (i) drawing $\tilde{p}_{t}$ for each $t$ and (ii)
comparing the draws to a standard uniform distribution.

Among calibrated forecasts, we evaluate the sharpness of the forecasting
distribution, such as the tightness of the prediction intervals. Czado et al.,
(2009) propose to evaluate the sharpness via scoring rules, in particular the
logarithmic scoring rule $-\log(p_{x})$ where $p_{x}$ is the probability mass
of the predictive distribution at the observed count $x$. We will utilize
these diagnostics to evaluate various warpDLMs and DGLM competitors.

###  5.2 Forecasting on Simulated Data

Existing models for count time series are often unable to capture multiple
challenging discrete distributional features, such as zero-inflation,
boundedness, and heaping. By comparison, the warpDLM is particularly well-
suited for these features. At the same time, the warpDLM is sufficiently
flexible to model simpler count data settings for which existing models have
been designed, such as Poisson data. To showcase this breadth of modeling, we
designed simulation experiments for both scenarios and compared the
forecasting performance of various warpDLMs against Bayesian competitors.

In the first experiment, we simulated zero-inflated and bounded count time
series data. One real-world example of such data is analyzed in Ensor et al.,
(2014), where the quantity of interest is the number of hours per day which
exceed a certain pollutant threshold. Many days show no pollution, leading to
zero-inflation, but occasionally every hour in a day exceeds the threshold,
resulting in heaped values at the upper bound of 24. Specifically, we sampled
from a zero-inflated Poisson with time varying mean $\lambda_{t}$, where
$\lambda_{1}\sim\mbox{Uniform}(5,15)$ and $\lambda_{t+1}\sim
N(\lambda_{t},0.2)$ for $t=2,\ldots,T$. The zero-inflation component was drawn
uniformly from $[0.1,0.3]$. Any sampled value above 24 was rounded down to the
bound.

In the second scenario, we simulated data from the INGARCH class using the R
package `tscount` (Liboschik et al.,, 2017). An INGARCH($p,q$) model assumes
$y_{t}\sim\mbox{Poisson}(\lambda_{t})$ with
$\lambda_{t}=\beta_{0}+\sum_{k=1}^{p}\beta_{k}Y_{t-k}+\sum_{\ell=1}^{q}\alpha_{\ell}\lambda_{t-\ell}.$
In our simulation, we chose relatively simple dynamics of $p=1$ and $q=1$ and
set $\beta_{0}=0.3$, $\beta_{1}=0.6$, and $\alpha_{1}=0.2$, with coefficients
selected such that simulated series were low-count time series.

For both simulation scenarios, 30 time series of length $T=200$ were
simulated. Among all competing state space models, we applied the same
evolution equation: a local level DLM with the univariate state $\theta_{t}$
and the time-invariant system $\\{1,1,V,W\\}$ with initial state
$\theta_{0}\sim N(0,3)$. The warpDLMs included _parametric_ transformations
(identity and square-root) and the _nonparametric_ transformation (6). In each
case, the variance components $V$ and $W$ were estimated via maximization of
the marginal likelihood (14) and then Monte Carlo samples from the joint
smoothing distribution (12) were drawn using Algorithm 1. These draws were
used to simulate from the forecasting distribution, as outlined in Section
3.2.3. Competing methods included Poisson, negative binomial, and (for the
bounded data scenario) binomial DGLMs, as implemented in the R package `KFAS`
(Helske,, 2017). We estimated variance components with maximum likelihood and
used the `KFAS` default value of the dispersion parameter for the negative
binomial DGLM.

For the DGLM models, it took several attempts to identify an initial value
that resulted in convergence for the maximum likelihood estimation. In
general, the likelihoods of state space models often have complex landscapes
which can make parameter estimation difficult and sensitive to starting inputs
(Helske,, 2017). For the warpDLM methods, a suitable starting point was found
by running a a Gibbs sampler for 2000 iterations (treating the first 500 as
burn-in) and then using the posterior means as the input to the optimization
algorithm.

In order to evaluate the forecasting performances, we used time series cross-
validation for one-step-ahead predictions. Starting at $t=100$, each of the
five methodologies (three warpDLM methods and two DGLM methods) were fitted
and one-step-ahead forecasts were drawn. For 50 equally dispersed time points
from $t=100$ to $t=200$, the process of training and drawing one-step ahead
forecasts was repeated. Thus, we have 50 observations to compare to 50 one-
step-ahead forecasts for each of the 30 time series.

Since we evaluate a multitude of series, we follow Kolassa, (2016) in testing
the uniformity of the rPIT values using the data-driven smooth test of
Ledwina, (1994). Hence, the output is a p-value, where small values indicate
poor calibration. We then consider the sharpness by evaluating the logarithmic
scoring rule. Our output for these models is 5000 draws from the predictive
distribution, so we approximate the probability mass by calculating the
proportion of draws which predicted that value. If a value was not sampled in
any posterior draw, but was observed, then the logarithmic scoring rule would
be infinite. To avoid this, we instead set the mass in that case to 0.0001.
For each series, we find the mean log score across the 50 time points and
compute the percent difference over a baseline method, chosen to be the
Poisson DGLM. The logarithmic scoring rule is negatively oriented, so a
negative percent difference equates to improved forecasting.

Figure 1 presents the calibration and sharpness across simulations for the
zero-inflated and bounded Poisson example. The binomial DGLM exhibited
extremely poor forecasting performance and is omitted for clarity of
presentation; results including the binomial DGLM are in the supplement. A
clear takeaway is that the warpDLM with nonparametric transformation is well-
calibrated. The identity warpDLM is also reasonably well-calibrated for many
simulations, highlighting the importance of the rounding operator. The
competing models produce extremely low p-values for almost every simulation,
which decisively shows that these models are not well-calibrated; Czado et
al., (2009) argue that such poor calibration should automatically disqualify
these models. The results for sharpness are consistent: the warpDLM with
nonparametric transformation produces the best distributional forecasts (on
average 30% better than the Poisson DGLM), while the warpDLM with identity
transformation is the second best. Modifications for the Poisson and negative
binomial DGLMs to enforce the upper bound of 24 on forecasts did not
substantively change the results. Clearly, the warpDLM offers significant
forecasting gains over competing methods, and advertises both good calibration
and maximal sharpness.

Figure 1: Zero-inflated Poisson simulations; Left: p-values measuring
calibration (larger is better) across simulations, with a line at p=0.05;
Right: percent difference in log score compared to baseline (negative values
imply improved forecasting relative to Poisson DGLM)

The second simulation scenario is presented in Figure 2. This simpler scenario
should be more favorable for the Poisson DGLM. Nonetheless, the warpDLM has
the flexibility to match or improve upon the DGLM forecasts. All five models
show proper calibration for nearly all simulations. Each of warpDLMs exhibits
sharpness similar to the Poisson DGLM, and the warpDLM with square-root
transformation actually appears to improve on the forecasting for the majority
of the series. The negative binomial DGLM performs uniformly worse than all
other methods.

Figure 2: INGARCH simulations; Left: p-values measuring calibration (larger is
better) across simulations, with a line at p=0.05; Right: percent difference
in log score compared to baseline (negative values imply improved forecasting
relative to Poisson DGLM)

Additional simulation results for higher-dimensional data ($n=5$, $T=4000$)
are included in the supplementary material.

###  5.3 Real Data Analysis

We apply the warpDLM to a multivariate time series of drug overdose (OD)
counts in Cincinnati. The OD reports are contained within a dataset of all
fire and EMS incidents in the city’s open access data repository, and include
daily counts of EMS calls for heroin and non-heroin ODs from January 1st, 2018
to January 31st, 2021; see the supplement for data wrangling and cleaning
details. Notably, Cincinnati has struggled with a severe problem of drug ODs
in the city, and in particular heroin (Li et al.,, 2019). From a statistical
modeling perspective, these data are discrete, time-ordered, multivariate
($n=2$), and moderately lengthy ($T=1127$); see Figure 3.

Figure 3: Bivariate time series of daily counts of EMS calls for heroin
(larger values) and non-heroin ODs. The dotted vertical line represents the
transition between the online and offline datasets. The solid lines are the
median smoothed predictions from the warpDLM.

We specify a warpDLM with the nonparametric transformation from (6), the
default rounding operator, and a linear growth model for each series
$i\in\\{1,2\\}$:

$\displaystyle z_{i,t}$ $\displaystyle=\mu_{i,t}+v_{i,t}$ (19)
$\displaystyle\mu_{i,t}$
$\displaystyle=\mu_{i,t-1}+\beta_{i,t-1}+w_{i,t-1}^{\mu}$ (20)
$\displaystyle\beta_{i,t}$ $\displaystyle=\beta_{i,t-1}+w_{i,t-1}^{\beta}$
(21)

where $\mu_{i,t}$ is a local level, $\beta_{i,t}$ is a slope, and (20)–(21)
comprise the DLM evolution equation (4). The linear growth model allows for
more flexibility than the local level specification previously introduced, but
also can effectively reduce to the simpler local level if the data supports it
(i.e. if the slope variance parameters are estimated to be very small).
Dependence between the $n=2$ series is induced by the observation error
distribution $\bm{v}_{t}\sim N_{2}\left(\bm{0},\bm{V}\right)$ and the state
error distributions $\bm{w}_{t}^{\mu}\sim
N_{2}\left(\bm{0},\bm{W}_{\mu}\right)$ and $\bm{w}_{t}^{\beta}\sim
N_{2}\left(\bm{0},\bm{W}_{\beta}\right)$ for
$\bm{v}_{t}=(v_{1,t},v_{2,t})^{\prime}$,
$\bm{w}_{t}^{\mu}=(w_{1,t}^{\mu},w_{2,t}^{\mu})^{\prime}$, and
$\bm{w}_{t}^{\beta}=(w_{1,t}^{\beta},w_{2,t}^{\beta})^{\prime}$, with
$\bm{V},\bm{W}_{\mu},\bm{W}_{\beta}$ assigned inverse-Wishart priors. Model
(19)–(21) is expressed in the usual DLM form of (3)–(4) in the supplement. As
an alternative approach, we also constructed a seasonal warpDLM, but its
performance was indistinguishable from the linear growth model (see the
supplementary material).

We highlight the capabilities of the warpDLM for both offline and online
analysis of multivariate count time series data. Specifically, we perform an
offline analysis for the years of 2018 and 2019 ($T=730$) and switch to an
online particle filter for the 2020 and January 2021 data ($T=397$). Given the
lengthy time series, offline inference was conducted using the Gibbs sampler
(Algorithm 2), which returned 10,000 draws from the state filtering
distribution of $\bm{\theta}_{t}$ for $t=730$ and the posterior means of the
variance components. These quantities were input into the particle filter of
Algorithm 3 and run sequentially through each of the $T=397$ time points in
the online dataset, resulting in 10,000 sequential paths of the smoothed
states. These smoothed states were iterated one step ahead to obtain draws
from the one-step forecasting distribution as in Section 3.2.3. We also
considered longer forecast horizons (see the supplementary material).
Unsurprisingly, the forecast distribution accuracy decreased for longer
horizons.

To summarize the trends in Figure 3, we compute the pointwise medians of the
smoothing predictive distributions $p(\tilde{\bm{y}}_{t}|\bm{y}_{1:T})$ under
the warpDLM. For each sampled path of states $\bm{\theta}_{1:T}^{*}\sim
p(\bm{\theta}_{1:T}|\bm{y}_{1:T})$, we draw a corresponding path
$\tilde{\bm{z}}_{1:T}$ via the DLM observation equation (3), and set
$\bm{\tilde{y}}_{1:T}=h\circ g^{-1}(\tilde{\bm{z}}_{1:T})$, similar to the
forecasting procedure in Section 3.2.3. Notably, this quantity provides a
count-valued point estimate that “smooths” each time series conditional on the
complete data $\bm{y}_{1:T}$. The heroin ODs exhibit both higher values and
greater variability over time, although the two time series appear to converge
to similar levels in early 2021.

To evaluate the fitness of the warpDLM, we report calibration of the one-step
forecasts in Figure 4. The sorted rPIT values of both series are plotted
against standard uniform quantiles; the 45 degree line indicates calibration.
Since the rPIT values are random quantities, some variation is natural. We
illustrate this inherent variability by sampling 100 draws of size $396$ from
a standard uniform distribution, and plot their sorted values against the
uniform quantiles.

Figure 4: rPIT plots from one-step-ahead particle filter forecasts for both OD
count series. The dashed 45 degree line indicates perfect calibration. Gray
lines are computed from standard uniform draws, to illustrate inherent
variability in the randomized metric.

The warpDLM one-step forecasts are well-calibrated for non-heroin ODs. For the
heroin ODs, the calibration shows minor deviations from uniformity. The
histogram of the rPIT values (not shown here) shows a slightly inverse-U
shape, which implies forecasts are overdispersed relative to the actual
values. A closer investigation of the data shows that the online dataset
showed considerably less variation than the offline dataset. Since the
nonparametric transformation (6) was inferred using only the _offline_
dataset, this result suggests the need to occasionally update the
transformation during online inference, for example by updating the estimate
of the marginal CDF $F_{y}$.

Lastly, we evaluate the computational performance of the optimal particle
filter, both for Monte Carlo efficiency and raw computing time. Figure 5
presents the effective sample size,
$ESS=1/\sum_{s=1}^{S}(\tilde{w_{t}}^{(s)})^{2}$ using the normalized weights
$\tilde{w_{t}}$, along with the number of seconds needed to update the model
at each time point (on a laptop with an Intel Core i5-6200U CPU with 8 GB
RAM). The _optimal_ particle filter—uniquely obtained via the warpDLM
distributional results from Section 3—produces consistently large ESS values.
Crucially, the computation time does _not_ increase with the time index and
hovers around 13-14 seconds for most updates. Given this computation time, the
warpDLM may be applied to more demanding streaming data, such as minute-by-
minute data. In all, the optimal particle filter for the warpDLM enabled (i)
modeling of multivariate count time series data, (ii) well-calibrated
forecasts, and (iii) sufficient scalability for moderate- to high-frequency
online analysis.

Figure 5: Effective sample sizes (left; maximum is 10,000) and raw computing
time in seconds (right) for each time point in the online dataset.

##  6 Conclusion

We introduced a state space modeling framework for a multivariate time series
of counts by _warping_ a dynamic linear model (warpDLM). The warpDLM
advertises (i) a data-coherent model (and forecasting distribution) that
matches the discrete support of the data; (ii) semiparametric modeling
capabilities via an inferred transformation; (iii) the familiarity and dynamic
flexibility of DLMs, including the ability to handle missing data; (iv)
analytic filtering and smoothing recursions; and (v) customized algorithms for
forecasting and inference, including direct Monte Carlo sampling and a Gibbs
sampler for offline inference and an optimal particle filter for online
inference. These results and methods also apply for integer-valued and rounded
data.

Using simulated data, we showed that the warpDLM offers better distributional
forecasting than Bayesian competitors, particularly when the count data
exhibit multiple complexities such as zero-inflation and boundedness. Finally,
we showcased the online inference and forecasting capabilities of the warpDLM
for a multivariate time series of drug overdose counts.

The generality of the warpDLM enables several useful extensions. Successful
tools from (Gaussian) DLMs, such as scale-free variance modeling or discount
factors, may be incorporated into the warpDLM framework. These specifications
often result in $t$-distributed state updates and predictions, which may be
linked to our results via the selection-$t$ distribution. Similarly, other
increases in modeling complexity can be accompanied by appropriate algorithmic
advancements, such as sequential Monte Carlo methods that sample both the
unknown parameters and the states concurrently (Chopin et al.,, 2012). Despite
the inherent challenges in these more complex modeling and computing
environments, the analytic and recursive updates derived for the warpDLM offer
a promising pathway for efficient and perhaps optimal implementations.

The warpDLM can also be thought of as a (nonlinear) hierarchical state-space
model (HSSM; Gamerman and Migon,, 1993), with the latent data $z_{t}$
introducing a first-level hierarchy. In particular, connections may be made to
the very general HSSM structure presented in Katzfuss et al., (2020), with the
warping operation presented here acting as the transformation layer in their
setup. Exploring this relationship further might open up new possibilities for
applying the warpDLM in high-dimensional problems.

## References

  * Aktekin et al., (2018) Aktekin, T., Polson, N., and Soyer, R. (2018).  Sequential bayesian analysis of multivariate count data.  Bayesian Analysis, 13(2).
  * Albert and Chib, (1993) Albert, J. H. and Chib, S. (1993).  Bayes inference via gibbs sampling of autoregressive time series subject to markov mean and variance shifts.  Journal of Business and Economic Statistics, 11(1):1–15.
  * Andrieu et al., (2010) Andrieu, C., Doucet, A., and Holenstein, R. (2010).  Particle markov chain monte carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342.
  * Arellano-Valle et al., (2006) Arellano-Valle, R. B., Branco, M. D., and Genton, M. G. (2006).  A unified view on skewed distributions arising from selections.  Canadian Journal of Statistics, 34(4):581–601.
  * Berry and West, (2019) Berry, L. R. and West, M. (2019).  Bayesian forecasting of many count-valued time series.  Journal of Business & Economic Statistics, 38(4):872–887.
  * Botev, (2017) Botev, Z. I. (2017).  The normal law under linear restrictions: simulation and estimation via minimax tilting.  Journal of the Royal Statistical Society. Series B: Statistical Methodology, 79(1):125–148.
  * Bradley et al., (2018) Bradley, J. R., Holan, S. H., and Wikle, C. K. (2018).  Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion).  Bayesian Analysis, 13(1).
  * Bürkner et al., (2020) Bürkner, P.-C., Gabry, J., and Vehtari, A. (2020).  Approximate leave-future-out cross-validation for bayesian time series models.  Journal of Statistical Computation and Simulation, 90(14):2499–2523.
  * Canale and Dunson, (2011) Canale, A. and Dunson, D. B. (2011).  Bayesian kernel mixtures for counts.  Journal of the American Statistical Association, 106(496):1528–1539.
  * Canale and Dunson, (2013) Canale, A. and Dunson, D. B. (2013).  Nonparametric Bayes modelling of count processes.  Biometrika, 100(4):801–816.
  * Carter and Kohn, (1994) Carter, C. K. and Kohn, R. (1994).  On Gibbs sampling for state space models.  Biometrika, 81(3):541–553.
  * Chopin et al., (2012) Chopin, N., Jacob, P. E., and Papaspiliopoulos, O. (2012).  Smc2: an efficient algorithm for sequential analysis of state space models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(3):397–426.
  * Cox et al., (1981) Cox, D. R., Gudmundsson, G., Lindgren, G., Bondesson, L., Harsaae, E., Laake, P., Juselius, K., and Lauritzen, S. L. (1981).  Statistical analysis of time series: Some recent developments [with discussion and reply].  Scandinavian Journal of Statistics, 8(2):93–115.
  * Czado et al., (2009) Czado, C., Gneiting, T., and Held, L. (2009).  Predictive model assessment for count data.  Biometrics, 65(4):1254–1261.
  * Doucet et al., (2000) Doucet, A., Godsill, S., and Andrieu, C. (2000).  On sequential monte carlo sampling methods for bayesian filtering.  Statistics and computing, 10(3):197–208.
  * Doucet et al., (2009) Doucet, A., Johansen, A. M., et al. (2009).  A tutorial on particle filtering and smoothing: Fifteen years later.  Handbook of nonlinear filtering, 12(656-704):3.
  * Durbin and Koopman, (2000) Durbin, J. and Koopman, S. J. (2000).  Time series analysis of non-gaussian observations based on state space models from both classical and bayesian perspectives.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(1):3–56.
  * Durbin and Koopman, (2002) Durbin, J. and Koopman, S. J. (2002).  A simple and efficient simulation smoother for state space time series analysis.  Biometrika, 89(3):603–615.
  * Durbin and Koopman, (2012) Durbin, J. and Koopman, S. J. (2012).  Time Series Analysis by State Space Methods.  Oxford University Press.
  * Ensor et al., (2014) Ensor, K. B., Ray, B. K., and Charlton, S. J. (2014).  Point source influence on observed extreme pollution levels in a monitoring network.  Atmospheric Environment, 92:191–198.
  * Fahrmeir, (1992) Fahrmeir, L. (1992).  Posterior mode estimation by extended kalman filtering for multivariate dynamic generalized linear models.  Journal of the American Statistical Association, 87(418):501–509.
  * Fasano et al., (2021) Fasano, A., Rebaudo, G., Durante, D., and Petrone, S. (2021).  A closed-form filter for binary time series.  Statistics and Computing, 31(4).
  * Fokianos, (2015) Fokianos, K. (2015).  Statistical analysis of count time series models: A glm perspective.  Handbook of discrete-valued time series, Handbooks of Modern Statistical Methods, pages 3–28.
  * Freeland and McCabe, (2004) Freeland, R. and McCabe, B. (2004).  Forecasting discrete valued low count time series.  International Journal of Forecasting, 20(3):427–434.
  * Frühwirth-Schnatter and Wagner, (2006) Frühwirth-Schnatter, S. and Wagner, H. (2006).  Auxiliary mixture sampling for parameter-driven models of time series of counts with applications to state space modelling.  Biometrika, 93(4):827–841.
  * Frühwirth‐Schnatter, (1994) Frühwirth‐Schnatter, S. (1994).  Data Augmentation and Dynamic Linear Models.  Journal of Time Series Analysis, 15(2):183–202.
  * Gamerman et al., (2013) Gamerman, D., dos Santos, T. R., and Franco, G. C. (2013).  A non-gaussian family of state-space models with exact marginal likelihood.  Journal of Time Series Analysis, 34(6):625–645.
  * Gamerman and Migon, (1993) Gamerman, D. and Migon, H. S. (1993).  Dynamic Hierarchical Models.  Journal of the Royal Statistical Society: Series B (Methodological), 55(3):629–642.
  * Helske, (2017) Helske, J. (2017).  Kfas: Exponential family state space models in r.  Journal of Statistical Software, 78(10).
  * Jia et al., (2021) Jia, Y., Kechagias, S., Livsey, J., Lund, R., and Pipiras, V. (2021).  Latent gaussian count time series.  Journal of the American Statistical Association, page 1–11.
  * Kalman, (1960) Kalman, R. E. (1960).  A new approach to linear filtering and prediction problems.  Transactions of the ASME–Journal of Basic Engineering, 82(Series D):35–45.
  * Kantas et al., (2015) Kantas, N., Doucet, A., Singh, S. S., Maciejowski, J., and Chopin, N. (2015).  On particle methods for parameter estimation in state-space models.  Statistical Science, 30(3).
  * Katzfuss et al., (2020) Katzfuss, M., Stroud, J. R., and Wikle, C. K. (2020).  Ensemble kalman methods for high-dimensional hierarchical dynamic space-time models.  Journal of the American Statistical Association, 115(530):866–885.
  * King and Kowal, (2023) King, B. and Kowal, D. R. (2023).  countSTAR: Flexible Modeling of Count Data.  R package version 1.0.1.9000.
  * Kolassa, (2016) Kolassa, S. (2016).  Evaluating predictive count data distributions in retail sales forecasting.  International Journal of Forecasting, 32(3):788–803.
  * Kowal and Canale, (2020) Kowal, D. R. and Canale, A. (2020).  Simultaneous Transformation and Rounding (STAR) Models for Integer-Valued Data.  Electronic Journal of Statistics.
  * Kowal and Wu, (2021) Kowal, D. R. and Wu, B. (2021).  Semiparametric count data regression for self-reported mental health.  Biometrics.
  * Kowal and Wu, (2022) Kowal, D. R. and Wu, B. (2022).  Semiparametric discrete data regression with monte carlo inference and prediction.
  * Ledwina, (1994) Ledwina, T. (1994).  Data-driven version of neyman’s smooth test of fit.  Journal of the American Statistical Association, 89(427):1000–1005.
  * Li et al., (2019) Li, Z. R., Xie, E., Crawford, F. W., Warren, J. L., McConnell, K., Copple, J. T., Johnson, T., and Gonsalves, G. S. (2019).  Suspected heroin-related overdoses incidents in cincinnati, ohio: A spatiotemporal analysis.  PLoS medicine, 16(11):e1002956.
  * Liboschik et al., (2017) Liboschik, T., Fokianos, K., and Fried, R. (2017).  Tscount: An R package for analysis of count time series following generalized linear models.  Journal of Statistical Software, 82(5):1–51.
  * O’Hara and Kotze, (2010) O’Hara, R. and Kotze, J. (2010).  Do not log-transform count data.  Nature Precedings, pages 1–1.
  * Petris, (2010) Petris, G. (2010).  An r package for dynamic linear models.  Journal of Statistical Software, 36(12).
  * Petris et al., (2009) Petris, G., Petrone, S., and Campagnoli, P. (2009).  Dynamic linear models with R.  Springer.
  * Prado et al., (2021) Prado, R., Ferreira, M. A. R., and West, M. (2021).  Time series: Modeling, computation, and inference.
  * Rauch et al., (1965) Rauch, H. E., Tung, F., and Striebel, C. T. (1965).  Maximum likelihood estimates of linear dynamic systems.
  * Siegfried and Hothorn, (2020) Siegfried, S. and Hothorn, T. (2020).  Count transformation models.  Methods in Ecology and Evolution, 11(7):818–827.
  * Watanabe, (2010) Watanabe, S. (2010).  Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  Journal of Machine Learning Research, 11(Dec):3571–3594.
  * West and Harrison, (1997) West, M. and Harrison, J. (1997).  Bayesian Forecasting and Dynamic Models (2nd Ed.).  Springer-Verlag, Berlin, Heidelberg.
  * West et al., (1985) West, M., Harrison, P. J., and Migon, H. S. (1985).  Dynamic generalized linear models and bayesian forecasting.  Journal of the American Statistical Association, 80(389):73–83.

This supplementary document contains the following:

  * •

Section A: Proposition and proof of SLCT-N closure under linear
transformations as well as proofs of main results in the paper

  * •

Section B: Results from simulations that include binomial DGLM

  * •

Section C: Further details on application dataset and additional model results

  * •

Section D: Results from particle filter applied to higher-dimensional
simulated dataset

##  Appendix A Further Theory and Proofs

As discussed in Section 3.1, the selection normal inherits closure under
linear transformations, a result which is used to show that the warpDLM state
predictive distribution stays within the selection normal family. We formalize
and prove this in the below proposition.

######  Proposition 1.

Suppose
$\bm{\theta}\sim\mbox{SLCT-N}_{n,p}(\bm{\mu}_{z},\bm{\mu}_{\theta},\bm{\Sigma}_{z},\bm{\Sigma}_{\theta},\bm{\Sigma}_{z\theta},\mathcal{C})$
and let $\bm{\theta}^{*}\stackrel{{\scriptstyle
d}}{{=}}\bm{A}\bm{\theta}+\bm{a}$ where $\bm{A}$ is a fixed $q\times p$ matrix
and $\bm{a}\sim N_{q}(\bm{\mu}_{a},\bm{\Sigma}_{a})$ is independent of
$\bm{\theta}$.

Then
$\bm{\theta}^{*}\sim\mbox{SLCT-N}_{n,q}(\bm{\mu}_{z},\bm{\mu}_{\theta^{*}},\bm{\Sigma}_{z},\bm{\Sigma}_{\theta^{*}},\bm{\Sigma}_{z\theta^{*}},\mathcal{C})$
where $\bm{\mu}_{\theta^{*}}=\bm{A}\bm{\mu}_{\theta}+\bm{\mu}_{a}$,
$\bm{\Sigma}_{\theta^{*}}=\bm{A}\bm{\Sigma}_{\theta}\bm{A}^{\prime}+\bm{\Sigma}_{a}$,
and $\bm{\Sigma}_{z\theta^{*}}=\bm{\Sigma}_{z\theta}\bm{A}^{\prime}$.

###### Proof.

First consider the distribution of $\bm{A}\bm{\theta}$. Since the moments of
the joint distribution $(\bm{z},\bm{\theta})$ are given by assumption, it
follows that $(\bm{z},\bm{A}\bm{\theta})$ is jointly Gaussian with moments
available by straightforward calculation, so
$\bm{A}\bm{\theta}\sim\mbox{SLCT-N}_{n,q}(\bm{\mu}_{z},\bm{A}\bm{\mu}_{\theta},\bm{\Sigma}_{z},\bm{A}\bm{\Sigma}_{\theta}\bm{A}^{\prime},\bm{\Sigma}_{z\theta}\bm{A}^{\prime},\mathcal{C})$.

As noted in ArellanoValle2006, selection normal distributions have a moment
generating function (MGF). For
$[\bm{\theta}|\bm{z}\in\mathcal{C}]\sim\mbox{SLCT-N}_{n,p}(\bm{\mu}_{z},\bm{\mu}_{\theta},\bm{\Sigma}_{z},\bm{\Sigma}_{\theta},\bm{\Sigma}_{z\theta},\mathcal{C})$,
the MGF is:

$M_{[\bm{\theta}|\bm{z}\in\mathcal{C}]}(\bm{s})=\exp\left(\bm{s}^{\prime}\bm{\mu}_{\theta}+\frac{1}{2}\bm{s}^{\prime}\bm{\Sigma}_{\theta}\bm{s}\right)\frac{\bar{\Phi}_{n}(\mathcal{C};\bm{\Sigma}_{z\theta}\bm{s}+\bm{\mu}_{z},\bm{\Sigma}_{z})}{\bar{\Phi}_{n}(\mathcal{C};\bm{\mu}_{z},\bm{\Sigma}_{z})}.$

Using the selection normal MGF and noting independence between
$\bm{A}\bm{\theta}$ and $\bm{a}$, it follows that
$M_{\bm{\theta}^{*}}(\bm{s})=M_{\bm{A}\bm{\theta}+\bm{a}}(\bm{s})=M_{\bm{A}\bm{\theta}}(\bm{s})M_{\bm{a}}(\bm{s})$
where $M_{\bm{A}\bm{\theta}}$ is given by inserting the appropriate parameters
into the equation above and
$M_{\bm{a}}(\bm{s})=\exp(\bm{s}^{\prime}\bm{\mu}_{a}+\frac{1}{2}\bm{s}^{\prime}\bm{\Sigma}_{a}\bm{s})$.
The result of this product is the MGF of the stated SLCT-N distribution. ∎

We now state the proofs of all results in the main article, with the exception
of Theorem 1, whose result follows directly from the definition of a selection
normal distribution and the warpDLM model setup. The concurrent work of
kowalSTARconjugate also explores Theorem 1 and Lemma 1, but with focus on non-
dynamic linear regression for discrete data.

###### Proof (Lemma 1).

The SLCT-N prior is equivalently defined by
$[\bm{\theta}|\bm{z}_{0}\in\mathcal{C}_{0}]$ for
$(\bm{z}_{0}^{\prime},\bm{\theta}^{\prime})^{\prime}$ jointly Gaussian with
moments given in the prior. The posterior is constructed similarly:
$[\bm{\theta}|\bm{y}]\stackrel{{\scriptstyle
d}}{{=}}[\bm{\theta}|\bm{z}_{0}\in\mathcal{C}_{0},\bm{z}\in
g(\mathcal{A}_{\bm{y}})]\stackrel{{\scriptstyle
d}}{{=}}[\bm{\theta}|\bm{z}\in\mathcal{C}_{0}\times g(\mathcal{A}_{\bm{y}})]$
where $\bm{z}_{1}=(\bm{z}_{0}^{\prime},\bm{z}^{\prime})^{\prime}$. It remains
to identify the moments of
$(\bm{z}_{1}^{\prime},\bm{\theta}^{\prime})^{\prime}=(\bm{z}_{0}^{\prime},\bm{z}^{\prime},\bm{\theta}^{\prime})^{\prime}$.
For each individual block of $\bm{z}_{0}$, $\bm{z}$, and $\bm{\theta}$ and the
pairs ($\bm{z}_{0},\bm{\theta}$) and ($\bm{z},\bm{\theta}$), the moments are
provided by either the prior or the posterior in Theorem 1. Finally, we have
cross-covariance
$\mbox{Cov}(\bm{z}_{0},\bm{z})=\mbox{Cov}(\bm{z}_{0},\bm{F}\bm{\theta}+\bm{v})=\bm{\Sigma}_{z_{0}\theta}\bm{F}^{\prime}$.
∎

###### Proof (Theorem 2).

From the model setup, we have
$\bm{\theta}_{t}=\bm{G}_{t}\bm{\theta}_{t-1}+\bm{w}_{t}$ and
$\bm{\theta}_{t-1}|\bm{y}_{1:t-1}$ follows a SLCT-N distribution as given.
Thus by Proposition 1, we have our result that
$\bm{\theta}_{t}|\bm{y}_{1:t-1}$ follows a SLCT-N with the given parameters.

Then, with the state predictive distribution $\bm{\theta}_{t}|\bm{y}_{1:t-1}$
as our prior combined with the warpDLM likelihood (5), we can directly apply
Lemma 1 to get the second result. ∎

###### Proof (Theorem 3).

We have that $p(\bm{\theta}_{1:T}|\bm{y}_{1:T})\propto
p(\bm{\theta}_{1:T})p(\bm{y}_{1:T}|\bm{\theta}_{1:T})$. By the model
specification, we know $\bm{\theta}_{1:T}\sim
N_{pT}(\bm{\mu}_{\theta},\bm{\Sigma}_{\theta})$ as defined above. Furthermore,
recall that given the latent states, the observations are conditionally
independent. Hence we can represent the likelihood as
$p(\bm{y}_{1:T}|\bm{\theta}_{1:T})=\prod_{t=1}^{T}\bar{\Phi}_{n}(g(\mathcal{A}_{y_{t}});\bm{F}_{t}\bm{\theta}_{t},\bm{V}_{t})$
which can be rewritten as
$\bar{\Phi}_{nT}(C;\bm{\mathfrak{F}}\bm{\theta}_{1:T},\bm{\mathfrak{V}})$.
Therefore, we have

$p(\bm{\theta}_{1:T}|\bm{y}_{1:T})\propto\phi_{pT}(\bm{\theta}_{1:T};\bm{\mu}_{\theta},\bm{\Sigma}_{\theta})\bar{\Phi}_{nT}(C;\bm{\mathfrak{F}}\bm{\theta}_{1:T},\bm{\mathfrak{V}})$

which is the kernel of a SLCT-N distribution with parameters given in (12). ∎

###### Proof (Corollary 2).

The marginal likelihood is simply the normalizing constant of the joint
smoothing distribution, and its form is given by the denominator in (9). ∎

###### Proof (Corollary 3).

By Bayes’ Theorem, we have

$p(\bm{\theta}_{t}|\bm{\theta}_{t-1},\bm{y}_{t})=\frac{p(\bm{y}_{t}|\bm{\theta}_{t})p(\bm{\theta}_{t}|\bm{\theta}_{t-1})}{p(\bm{y}_{t}|\bm{\theta}_{t-1})}\;.$
(22)

Using the DLM and warpDLM properties, we know that

$p(\bm{y}_{t}|\bm{\theta}_{t})=\bar{\Phi}_{n}(g(\mathcal{A}_{y_{t}});\bm{F}_{t}\bm{\theta}_{t},\bm{V}_{t})\quad\mbox{and}\quad
p(\bm{\theta}_{t}|\bm{\theta}_{t-1})=\phi_{p}(\bm{G}_{t}\bm{\theta}_{t-1},\bm{W}_{t})\;.$

Putting this together, recognize that the numerator of (22) forms the kernel
of a selection normal distribution (9) with parameters given in (16).
Furthermore, the marginal distribution $p(\bm{y}_{t}|\bm{\theta}_{t-1})$ is
that kernel’s normalizing constant, leading us to the result of (17). ∎

##  Appendix B Binomial DGLM Simulation Results

In the main article, we did not show the binomial DGLM forecasting results for
the bounded, zero-inflated Poisson simulation because it performed much worse
than all other methods and therefore skewed the graphics. For completeness we
display here in Figures B.1 and B.2 the corresponding calibration and
sharpness plots with binomial DGLM included.

Figure B.1: Box plot of p-values measuring calibration (larger is better)
across simulations, with a line at p=0.05 Figure B.2: Violin plot of percent
difference in log score compared to baseline across simulations (negative
difference implies improved forecasting relative to Poisson DGLM)

##  Appendix C Application Details

###  C.1 Dataset Cleaning

The Cincinnati dataset contains all fire department data, but our focus in
this application is only on the overdose responses. To extract the appropriate
observations, we followed the data dictionary in selecting heroin overdose
responses as those with rows whose `INCIDENT_TYPE_ID` was one of ”HEROIF”,
”HEROIN-COMBINED”, ”HEROINF - FIRE ONLY”, ”HEROINF-FIRE ONLY”, ”HERON F”, or
”HERONF”. For the “other” overdose responses, we took all events which
correspond to EMS protocol 23 (Overdose/Poisoning), i.e. those rows whose
`INCIDENT_TYPE_ID` started with the number 23.

Once the appropriate responses were selected, the only cleaning necessary was
to count the number of occurrences per day and properly format dates. We did
so by binning the observations into days based on the `CREATE_TIME_INCIDENT`
and assigning zeroes on any days where no incidents were reported.

###  C.2 Linear Growth Model Details

The warpDLM latent specification for the application is an example of
Seemingly Unrelated Time Series Equations (SUTSE), where the inter-series
dependence is defined through the covariance matrix of the model errors. Given
the linear growth model for each series as defined in equations (19)–(21), the
full bivariate model can be rewritten in the general DLM formulation.
Specifically, with
$\bm{\theta}_{t}=(\mu_{1,t},\mu_{2,t},\beta_{1,t},\beta_{2,t})^{\prime}$, the
state evolution equation is

$\bm{\theta}_{t}=\begin{bmatrix}1&0&1&0\\\ 0&1&0&1\\\ 0&0&1&0\\\
0&0&0&1\end{bmatrix}\bm{\theta}_{t-1}+\bm{w}_{t}\quad\bm{w}_{t}\sim
N_{4}\left(\bm{0},\bm{W}=\left[\begin{array}[]{@{}c|c@{}}\Large\bm{W}_{\mu}&\mbox{\Large
0}\\\ \hline\cr\mbox{\Large 0}&\Large\bm{W}_{\beta}\end{array}\right]\right)$
(23)

and the observation equation is

$\bm{z}_{t}=\begin{bmatrix}1&0&0&0\\\
0&1&0&0\end{bmatrix}\bm{\theta}_{t}+\bm{v}_{t}\quad\bm{v}_{t}\sim
N_{2}\left(\bm{0},\bm{V}\right)$ (24)

for the latent $\bm{z}_{t}=(z_{1,t},z_{2,t})$.

###  C.3 Seasonal Model Specification

In addition to the linear growth model presented above, we also implemented a
SUTSE model with a seasonal term. There are two main ways to include
seasonality within DLMs: dummy variables or Fourier-form/trigonometric
seasonality. For a time series with period $s$, the dummy approach requires
$s-1$ state variables, whereas the trigonometric approach only requires
$\lfloor(s/2)\rfloor$, or even less, since the number of trigonometric terms
can be truncated to achieve a smoother and more parsimonious seasonal trend.
The overdose data appears to have a yearly seasonality and we have daily data,
so our period is $s=365$. For such a long period, Fourier-form seasonality is
most reasonable.

We used the same linear growth specification as described in Section 5.3, but
added the first 5 harmonics from a trigonometric seasonal component of period
365. A linear growth DLM with Fourier-form seasonality inherits the same state
evolution dynamics as (20)-(21) for the local level and slope, but the
observation equation adds a seasonal term $\gamma_{i,t}$ for each series
$i\in{1,2}$. Thus, equation (19) becomes

$z_{i,t}=\mu_{i,t}+\gamma_{i,t}+v_{i,t}\;.$ (25)

We write the seasonal term $\gamma_{i,t}$ dynamics as follows:

$\displaystyle\gamma_{i,t}$ $\displaystyle=\sum_{j=1}^{5}\gamma_{i,j,t}$ (26)
$\displaystyle\gamma_{i,j,t}$
$\displaystyle=\gamma_{i,j,t-1}\cos{\lambda_{j}}+\gamma^{*}_{i,j,t-1}\sin{\lambda_{j}}+w_{i,j,t-1}^{\gamma}$
(27) $\displaystyle\gamma^{*}_{i,j,t}$
$\displaystyle=-\gamma_{i,j,t-1}\sin{\lambda_{j}}+\gamma^{*}_{i,j,t-1}\cos{\lambda_{j}}+w_{i,j,t-1}^{\gamma^{*}}$
(28)

where $j=1,\ldots,5$ and $\lambda_{j}=2\pi j/365$. Inter-series dependence
between the seasonal components can come from the error terms:
$\bm{w}_{j,t-1}^{\gamma},\bm{w}_{j,t-1}^{\gamma^{*}}\overset{iid}{\sim}N_{2}\left(\bm{0},\bm{W}_{\gamma}\right)$.
For more details on Fourier-form seasonality and how it can be represented in
traditional DLM form, consult Section 3.2.3 of petrisDLM or Section 6.1 of
HelskeKFAS.

The above model results in a state space of dimension 24:
$\bm{\theta}_{t}=(\mu_{i,t},\beta_{i,t},\beta_{i,t},\gamma_{i,j,t},\gamma^{*}_{i,j,t})^{\prime}$
for $i\in\\{1,2\\}$ and $j=1,\ldots,5$. We set the variance of the seasonal
component $\bm{W}_{\gamma}$ to be 0, so that the periodic piece of the model
was static. Using this model specification, we performed the same analysis,
running first the offline Gibbs sampler followed by the particle filter. The
seasonal component does not appear to have much effect, as the forecasting
performance was much the same as the linear growth model, even at longer
horizons. This can be seen in Figure C.1, where we compare the performance of
the two model specifications on the heroin time series. To do so, we compute
the percent difference in log score between the seasonal model and the linear
growth model across different forecast horizons. Interestingly, the linear
growth model actually seems to produce slightly better forecasts on average at
longer horizons, although the difference is very minimal.

Figure C.1: Percent difference in log score of seasonal model forecasts for
heroin OD time series compared to linear growth model predictions, plotted
across forecast horizons. The positive values for longer horizons suggest that
the linear growth model outperforms the seasonal model for these longer-term
forecasting distributions.

###  C.4 Multi-Step Forecasting Results

In the main text, we primarily focus on the 1-step-ahead forecasts. Here, we
investigate longer forecast horizons, which may be more relevant for planners
and decision-makers in certain cases. We computed forecasts ranging from one
day ahead to two weeks ahead (14 steps ahead with our daily data). Calibration
for multi-step-ahead forecasts was usually less satisfactory, with the multi-
step-ahead forecast distributions often overdispersed relative to the observed
data. This is not surprising, given that the uncertainty of the forecast
distribution accumulates with each additional forecast step. Figure C.2 shows
the boxplot of log scores across forecast horizons from the linear growth
warpDLM model for both time series. It can be seen that the log scores
steadily increase with the forecast horizon.

Figure C.2: Boxplot of log scores from linear growth warpDLM model across
different forecast horizons

##  Appendix D Particle Filtering in Higher Dimensions

In Section 5.3, we combined offline inference using the Gibbs sampler with
online inference via the particle filter. We demonstrated that the particle
filter results in sensible smoothing and calibrated forecasts for a real-data
application with $n=2$ and $T=1127$. However, we also wanted to explore how
the particle filter performed in higher-dimensional settings, with larger $n$
and longer $T$. To do so, we simulated a dataset with $n=5$ and $T=4000$. For
the data-generating process, we used the multivariate Poisson-scaled Beta
(MPSB) class developed in Aktekin2018, as this provides a relatively simple
method of simulating multivariate count time series data. Briefly, this class
is described by the following hierarchical model:

$\displaystyle(Y_{jt}|\lambda_{j},\theta_{t})$
$\displaystyle\sim\mbox{Pois}(\lambda_{j}\theta_{t})\quad
j=1,\ldots,n\;;\;t=1,\ldots,T$ (29) $\displaystyle\theta_{t}$
$\displaystyle=\frac{\theta_{t-1}}{\gamma}\epsilon_{t}$ (30)
$\displaystyle(\epsilon_{t}|\bm{Y}_{1:t-1},\lambda_{1},\ldots,\lambda_{n})$
$\displaystyle\sim\mbox{Beta}(\gamma\alpha_{t-1},(1-\gamma)\alpha_{t-1})$ (31)
$\displaystyle\alpha_{t}$
$\displaystyle=\gamma\alpha_{t-1}+\sum_{j=1}^{n}Y_{jt}$ (32)

where $\alpha_{t-1}>0$ , $0<\gamma<1$ and
$D_{t-1}=\\{D_{t-2},Y_{1,t-1},\ldots,Y_{n,t-1}\\}$. The model is finalized by
setting the prior for $\theta_{0}$, which is of the form
$\theta_{0}\sim\mbox{Gamma}(\alpha_{0},\beta_{0})$. To simulate from this
model, it is also necessary to fix $\gamma$ and $\lambda_{j}$ for
$j=1,\ldots,n$. For our simulation, we generally followed the setup from
Section 4.1 of Aktekin2018, setting the series-specific factors $\lambda_{j}$
to $(2,2.5,3,3.5,4)$ and drawing the initial state from
$\mbox{Gamma}(\alpha_{0}=100,\beta_{0}=100)$. However, we modified the common
discount parameter $\gamma$, as the value of $0.3$ used in the Aktekin2018
simulation resulted in time series collapsing to zero for longer time scales.
To prevent this collapse, we increased its value considerably, in the end
fixing $\gamma=.93$. The resulting time series are plotted in Figure D.1.

Figure D.1: Time series plot of multivariate data simulated from the MPSB
class

We modeled this multivariate time series with the same linear growth model as
in the real-data application. Similarly, we split the dataset into offline and
online portions, running the Gibbs sampler for the first 1000 data points,
followed by the particle filter for the remaining 3000 time points. Our
primary objective was to ensure the particle filter is able to accurately
track the series in higher dimensions. We see that this goal is achieved by
visualizing the “push-forward” of our particles: for each particle
$\theta_{t}^{(s)}$ from our filter (not to be confused with the $\theta_{t}$
of the MPSB model), we draw from $N(F_{t}\theta_{t}^{(s)},V)$ and apply the
inverse transformation and rounding to arrive at a value on the scale of our
data. This is similar to the “smoothing predictive distribution” procedure
discussed in Section 5.3, only now we are looking at the “filtering predictive
distribution”. As with the smoothing, we can take pointwise medians of our
draws from this predictive filtering distribution; these are depicted in
Figure D.2.

Figure D.2: Multivariate online simulated dataset (light orange color)
overlaid with the median filtered predictive distribution (dark green color)

The particle filter is tracking each series well. However, it is worth noting
that the ESS in this setting is on average significantly lower than that of
our real-data application: the ratio of average ESS to number of particles
decreased from 58% to around 24%. This suggests that larger dimensions may
require adoption of further SMC techniques like mutation, as was discussed in
Section 4.3. As a final note, the computation time for each time point
increases with the number of series. This is due both to sampling from the
transition density (involving repeated simulation from an $n$-variate
multivariate truncated normal) and computation of the weights (evaluating an
$n$-variate multivariate normal CDF). These steps must be done for each
particle, although such computations could be parallelized for speed
improvements. For the $n=5$ example with $S=5000$ particles, each time point
is still executed in under 30 seconds, so our particle filter would remain
applicable for minute-by-minute streaming data of this dimension.

apalike refs.bib