Approximate filtering via discrete dual processes

Hidden Markov models are widely used statistical models for time series that assume an unobserved Markov process $(X_{t})_{t\geq 0}$, or hidden signal, driving the process that generates the observations $(Y_{t_{i}})_{i=0,\dots,n}$, e.g., by specifying the dynamics of one or more parameters of the observation density $f_{X_{t}}(y)$, called emission distribution. See [11] for a general treatment of hidden Markov models. In this framework, the first task is to estimate the trajectory of the signal given observations collected at discrete times $0=t_{0}<t_{1}<\cdots<t_{n}=T$, which amounts to performing sequential Bayesian inference by computing the so-called filtering distributions $p(x_{t_{i}}|y_{t_{0}},\ldots,y_{t_{i-1}})$, i.e., the marginal distributions of the signal at time $t_{i}$ conditional on observations collected up to time $t_{i-1}$. Originally motivated by real-time tracking and navigation systems and pioneered by [46, 47], classical and widely known explicit results for this problem include: the Kalman–Bucy filter, when both the signal and the observation process are formulated in a gaussian linear system; the Baum–Welch filter, when the signal has a finite state-space as the observations are categorical; the Wonham filter, when the signal has a finite state-space and the observations are Gaussian. These scenarios allow the derivation of so-called finite-dimensional filters, i.e., a sequence of filtering distributions whose explicit identification is obtained through a parameter update based on the collected observations and on the time intervals between the collection times. In such case, the resulting computational cost of the algorithm increases linearly with the number of observations. Other explicit results include [56, 26, 27, 55, 30, 31, 16]. Outside these classes, explicit solutions are difficult to obtain, and their derivation typically relies on ad hoc computations. This is especially true when the map $x\mapsto f_{x}$ is non-linear and when the signal transition kernel is known up to a series expansion, often intractable, as is the case for many widely used stochastic models. When exact solutions are not available, one must typically make use of approximate strategies, whose state of the art is most prominently based on extensions of the Kalman and particle filters. See, for example, [6, 15].

A somewhat weaker but useful notion with respect to that of a finite-dimensional filter was formulated in [13], who introduced the concept of computable filter. This extends the former class to a larger class of filters whose marginal distributions are finite mixtures of elementary kernels, rather than single kernels. Unlike the former case, such a scenario entails a higher computational cost, usually polynomial in the number of observations, but avoids the infinite-dimensionality typically implied by series expansion of the signal transition kernel. See [13, 14] for two examples.

Recently, [53] derived sufficient conditions for computable filtering based on duality. A dual process is a process $D_{t}$ which enjoys the identity

Here the expectation on the left-hand side is taken with respect to the transition law of the signal $X_{t}$, and that on the right hand side with respect to that of $D_{t}$, while the class of functions $h(x,d)$ which satisfy the above identity are called duality functions. See [44] for a review and for the technical details we have overlooked here. The use of duality is largely established in probability and statistical physics, see for example [21, 51, 7, 33, 23, 34, 35, 24, 52, 12, 25, 38, 1, 8, 18, 42, 28, 43]. The use of duality for inference was initiated by [53], who showed that for a reversible signal whose marginal distributions are conjugate to the emission distribution (i.e., the Bayesian update at a fixed $t$ can be computed in closed-form), computable filtering is guaranteed if the stochastic part of the dual process evolves on a finite state space. Examples of such scenario were given for non-linear hidden Markov models involving signals driven by Cox–Ingersoll–Ross (CIR) and $K$-dimensional Wright–Fisher (WF) diffusions, for which recursive formulae for the filtering distributions were derived. Along similar lines, duality was exploited for computable smoothing in [50], whereby the signal is also conditioned on data collected at later times, and for nonparametric hidden Markov models driven by Fleming–Viot and Dawson–Watanabe diffusions in [54, 2, 3].

In this paper, we investigate filtering problems for hidden Markov models when the dual process takes the more general form of a continuous-time Markov chain on a discrete state space. This is of interest for example in some population genetic models with selection [7] or interaction [4, 25] whose known dual processes are of birth-and-death (B&D) type, whose specific filtering problems are currently under investigation by some of the authors. When the dual process evolves in a countable state space, the filtering distributions can in general be expected to be countably infinite mixtures. This leads to inferential procedures which are not computable in the sense specified above, since the computation of the filtering distribution can no longer be exact. It is thus natural to wonder how the inferential procedures obtained in such duality-based scenario, possibly aided by some suitable approximation strategies, would perform.

The paper is organized as follows. In Section 2 we identify sufficient conditions for filtering based on discrete duals and provide a general description of the filtering operations in this setting. In Section 3, we apply these results to devise practical algorithms which allow to evaluate in recursive form filtering and smoothing distributions under this formulation. Section 4 and 5 investigate hidden Markov models driven by a Cox–Ingersoll–Ross diffusion, which admits a dual given by a one-dimensional B&D process, and by a $K$-dimensional Wright–Fisher diffusions, which is shown to admit a dual given by a $K$-dimensional Moran model. The latter can be seen as a multidimensional B&D process with constant total population size. Section 6 discusses several approximation strategies used to implement the above algorithms with these dual processes, and compares their performance with exact filtering based on the results in [53] and with a bootstrap particle filter as a general benchmark. Finally, we conclude with some brief remarks.

Consider a hidden signal $(X_{t})_{t\geq 0}$ given by a diffusion process on $\mathcal{X}\subset\mathbb{R}^{K}$, for $K\geq 1$. This takes here the role of a temporally evolving parameter which is the main target of estimation. Observations $Y_{t_{i}}\in\mathcal{Y}\subset\mathbb{R}^{D}$, $D\geq 1$, are collected at discrete times $0=t_{0}<t_{1}<\cdots<t_{n}=T$ with distribution $Y_{t_{i}}\overset{\text{ind}}{\sim}f_{x}(\cdot)$, given $X_{t_{i}}=x$. Knowledge of the signal state $x$ thus makes the observations conditionally independent. Given an observation $Y=y$, define the update operator $\phi_{y}$ acting on densities $\xi$ on $\mathcal{X}$ by

Here and later we assume all densities of interest exist with respect to an appropriate dominating measure. In (2), $\xi$ acts as a prior distribution on the signal state, which encodes the current knowledge (or lack thereof) on $X_{t}$, whereas $\mu_{\xi}(y)$ is interpreted as the marginal likelihood of a data point $y$ when $X_{t}$ has distribution $\xi$. The update operator thus amounts to an application of Bayes’ theorem for conditioning $\xi$ on a new observation $y$, leading to the updated density $\phi_{y}(\xi)$. For example, if $\xi$ is a $\text{Beta}(a,b)$ density on $[0,1]$, and $f_{x}$ is $\text{Bern}(x)$, then $\phi_{y}(\xi)$ is $\text{Beta}(a+y,b+1-y)$ as in a classical Beta-Bernoulli update.

Define also the propagation operator $\psi_{t}$ by

where $P_{t}$ is the transition density of the signal. Here $\psi_{t}(\xi)$ is the probability density at time $t$ obtained by propagating forward the density $\xi$ of the signal at time 0 by means of the signal semigroup.

We will make three assumptions, the first two of which are the same as in [53].

Assumption 1 (Reversibility). The signal $X_{t}$ is reversible with respect to the density $\pi$, i.e., $\pi(x)P_{t}(x^{\prime}|x)=\pi(x^{\prime})P_{t}(x|x^{\prime})$.

See the discussion in [53] on the possibility of relaxing the above assumption. For $K\in\mathbb{Z}_{+}$, define now the space of multi-indices $\mathcal{M}=\mathbb{Z}_{+}^{K}$ to be

whose origin is denoted $\mathbf{0}=(0,\ldots,0)$.

Assumption 2 (Conjugacy). For $\Theta\subset\mathbb{R}^{l}$, $l\in\mathbb{Z}_{+}$, let $h:\mathcal{X}\times\mathcal{M}\times\Theta\rightarrow\mathbb{R}_{+}$ be such that $\sup_{x\in\mathcal{X}}h(x,\mathbf{m},\theta)<\infty$ for all $\mathbf{m}\in\mathcal{M},\theta\in\Theta$ and $h(x,\mathbf{0},\theta^{\prime})=1$ for some $\theta^{\prime}\in\Theta$. Then $f_{x}(\cdot)$ is conjugate to distributions in the family

$$\mathcal{F}=\{g(x,\mathbf{m},\theta)=h(x,\mathbf{m},\theta)\pi(x),\ \mathbf{m}\in\mathcal{M},\theta\in\Theta\},$$

i.e., there exist functions $t:\mathcal{Y}\times\mathcal{M}\rightarrow\mathcal{M}$ and $T:\mathcal{Y}\times\Theta\rightarrow\Theta$ such that if $X_{t}\sim g(x,\mathbf{m},\theta)$ and $Y_{t}|X_{t}=x\sim f_{x}$, we have $X_{t}|Y_{y}=y\sim g(x,t(y,\mathbf{m}),T(y,\theta))$.

Here $g(x,\mathbf{m},\theta)$ takes the role of “current” prior distribution, i.e., the prior information on the signal state, possibly based on past observations through previous conditioning and propagations, and $g(x,t(y,\mathbf{m}),T(y,\theta))$ takes the role of the posterior, i.e., $g(x,\mathbf{m},\theta)$ conditional on a data point $y$ observed from $f_{x}$ for $X_{t}=x$. The functions $t(y,\mathbf{m}),T(y,\theta)$ provide the transformations that update the parameters based on $y$. In absence of data, the condition $h(x,\mathbf{0},\theta^{\prime})=1$ reduces $g(x,\mathbf{m},\theta)$ to $\pi(x)$.

The third assumption weakens Assumption 3 in [53] by assuming the dual process has finite activity on a discrete state space, and possibly has a deterministic companion.

Assumption 3 (Duality). Given a deterministic process $\Theta_{t}\in\Theta$ and a regular jump continuous-time Markov chain $M_{t}$ on $\mathbb{Z}_{+}^{K}$ with transition probabilities

$$p_{\mathbf{m},\mathbf{n}}(t;\theta):=\mathbb{P}(M_{t}=\mathbf{n}|M_{0}=\mathbf{m},\Theta_{0}=\theta),$$

(5)

equation (1) holds with $D_{t}=(M_{t},\Theta_{t})$ and $h$ as in Assumption 2.

The following result provides a full description of the propagation and update steps which allow to compute the filtering distribution.

The expression (6), together with (7), provides a general recipe on how to compute the forward propagation of the current marginal distribution of the signal $g(x,\mathbf{m},\theta)$, based on the transition probabilities of the dual continuous-time Markov chain. Since the update operator (2) can be easily applied to the resulting distribution, Proposition 2.1 then shows that under these assumptions all filtering distributions are countable mixtures of elementary kernels indexed by the state space of the dual process, with mixture weights proportional to the dual process transition probabilities $p_{\mathbf{m},\mathbf{n}}(t;\theta)$. When these transition probabilities happen to give positive mass only to points $\{\mathbf{n}\in\mathcal{M}:\ \mathbf{n}\leq\mathbf{m}\}$, as is the case for a pure-death process, then the right-hand side of (6) reduces to a finite sum, and one can construct an exact filter with a computational cost that is polynomial in the number of observations, as shown in [53].

The above approach can be seen as an alternative to deriving the filtering distribution of the signal by leveraging on a spectral expansion of the transition function $P_{t}$ in (3), which typically requires ad hoc computations and does not lend itself easily to explicit update operations through (2). Note also that expressions like (6) can be used, by taking appropriate limits of $p_{\mathbf{m},\mathbf{n}}(t;\theta)$ as $t\rightarrow 0$, to identify the transition kernel of the signal $P_{t}$ itself, see, e.g., [7, 53, 38].

In order to translate Proposition 2.1 into practical recursive formulae for filtering and smoothing, we are going to assume for simplicity of exposition that the time intervals between successive data collections $t_{i}-t_{i-1}$ equal $\Delta$ for all $i$. For ease of reading, we will therefore use the symbol $P_{\Delta}$ instead of $P_{t_{i}-t_{i-1}}$ for the signal transition function over the interval $\Delta=t_{i}-t_{i-1}$. We will also use the established notation whereby $i|0:i-1$ indicates that the argument refers to time $t_{i}=i\Delta$, and we are conditioning on the data collected at times from $0$ to $t_{i-1}=(i-1)\Delta$.

Define the filtering density

i.e., the law of the signal at time $t_{i}$ given data up to time $t_{i}$, obtained by integrating out the past trajectory. Define also the predictive density

i.e, the marginal density of the signal at time $t_{i+1}$, given data up to time $t_{i}$. This can be expressed recursively as a function of the previous filtering density $p(x_{i}|y_{0:i})$, as displayed. Finally, define the marginal smoothing density

where the signal is evaluated at time $t_{i}$ conditional on all available data. The first two distributions above are natural objects of inferential interest, whereas the latter is typically used to improve previous estimates once additional data become available. Finally, for $\Theta_{\Delta}$ as in Assumption 3 and $t(\cdot,\cdot),T(\cdot,\cdot)$ as in Assumption 2, define for $i=0,\dots,n$ the quantities

Here, $\vartheta_{i|0:i-1}$ denotes the state of the deterministic component of the dual process at time $i$, after the propagation from time $i-1$ and before updating with the datum collected at time $i$, and $\vartheta_{i|0:i}$ the state after such update.

The following Corollary of Proposition 2.1 extends a result of [53] (see also Theorem 1 in [50] for an easier comparison in a similar notation).

Input: $Y_{0:n}$, $t_{0:n}$

Result: $\vartheta_{i|0:i}$, $\mathbf{M}_{i|0:i}$ and $W_{i}=\{w_{\mathbf{m}}^{(i)},\mathbf{m}\in\mathbf{M}_{i|0:i}\}$

Initialise

      

      Set $\vartheta_{0|0}=T(Y_{0},\theta_{0})$ with $T$ as in Assumption 2

      

      Set $\mathbf{M}_{0|0}=\{t(Y_{0},\mathbf{0})\}=\{\mathbf{m}^{*}\}$ and $W_{0}=\{1\}$ with $t$ as in Assumption 2

      

      Compute $\vartheta_{1|0}$ from $\vartheta_{0|0}$ as in (17)

      

      Set $\mathbf{M}^{*}=\mathcal{B}(\mathbf{M}_{0|0})$ and $W^{*}=\{p_{\mathbf{m}^{*},\mathbf{n}}(\Delta,\vartheta_{0|0}),\mathbf{n}\in\mathbf{M}^{*}\}$ with $p_{\mathbf{m},\mathbf{n}}$ as in (5)

for $i$ from $1$ to $n$ do

      

      

      Update

            

            Set $\vartheta_{i|0:i}=T(Y_{i},\vartheta_{i|0:i-1})$

            

            Set $W_{i}=\left\{\frac{w_{\mathbf{m}}^{*}\mu_{\mathbf{m},\vartheta_{i|0:i-1}}(Y_{i})}{\sum_{\mathbf{n}\in\mathbf{M}^{*}}w_{\mathbf{n}}^{*}\mu_{\mathbf{n},\vartheta_{i}}(Y_{i})},\mathbf{m}\in\mathbf{M}^{*}\right\}$ with $\mu_{\mathbf{m},\theta}$ defined as in (9)

            Set $\mathbf{M}_{i|0:i}=\{t(Y_{i},\mathbf{m}),\mathbf{m}\in\mathbf{M}^{*}\}$ and update the labels in $W_{i}$

             Copy $\vartheta_{i|0:i}$, $\mathbf{M}_{i|0:i}$ and $W_{i}$ to be reported as the output

      

      

      

      Propagation

            

            Compute $\vartheta_{i+1|0:i}$ from $\vartheta_{i|0:i}$

            

            Set $\mathbf{M}^{*}=\mathcal{B}(\mathbf{M}_{i|0:i})$ and $W^{*}=\bigg{\{}\displaystyle{\sum_{\mathbf{m}\in\mathbf{M}_{i|0:i}}}w_{\mathbf{m}}^{(i)}p_{\mathbf{m},\mathbf{n}}(\Delta,\vartheta_{i|0:i}),\mathbf{n}\in\mathbf{M}^{*}\bigg{\}}$

      

      

end for

Note: $\mathbf{M}_{i|0:i}=\{\mathbf{m}\in\mathcal{M}:\ w_{\mathbf{m}}^{(i)}>0\}\subset\mathcal{M}$ is the support of the weights of $\nu_{i|0:i}$; $\mathcal{B}(\mathbf{m})$ denotes the states reached by the dual process from $\mathbf{m}$, and $\mathcal{B}(\mathbf{M})$ those reached from all $\mathbf{m}\in\mathbf{M}$.

Algorithm 1 outlines the pseudo-code for implementing the update and propagation steps of Corollary 3.1. How to use these results efficiently can depend on the model at hand. When the transition probabilities $p_{\mathbf{m},\mathbf{n}}(t;\theta)$ are available in closed form, their use could lead to the best performance, but can also at times face numerical instability issues (as is the case pointed out in Section 4 below). When the transition probabilities $p_{\mathbf{m},\mathbf{n}}(t;\theta)$ are not available in closed form, one can approximate them by simulating $N$ replicates of the dual component $M_{t}$, and then regroup probability masses according to the arrival states as done in (18). In our framework, the dual process is typically easier to simulate than the original process, given its discrete state space. For instance, pure-death or B&D processes are easily simulated using a Gillespie algorithm [36], whereby one alternates sampling waiting times and jump transitions for the embedded chain. Depending on the dual process, there might also be more efficient simulation strategies.

A different type of approximation of the propagation step (18) in Corollary 3.1 can be based on pruning the transition probabilities or the arrival weights under a given threshold, followed by a renormalisation of the weights. Both this approximation strategy and that outlined above assign positive weights only to a finite subset of $\mathcal{M}$, hence they overcome the infinite dimensionality of the dual process state space. [50] showed that the latter strategy allows to control the approximation error while retaining almost entirely the distributional information, thus affecting the inference negligibly. In the next sections we will investigate such strategies for two hidden Markov models driven by Cox–Ingersoll–Ross and $K$-dimensional Wright–Fisher diffusions.

Now, in order to describe the marginal smoothing densities (16), we need an additional assumption and some further notation.

Assumption 4 For $h$ as in Assumption 3, there exist functions $d:\mathcal{M}^{2}\to\mathcal{M}$ and $e:\Theta^{2}\to\Theta$ such that for all $x\in\mathcal{X},\ \mathbf{m},\mathbf{m}^{\prime}\in\mathcal{M},\ \theta,\theta^{\prime}\in\Theta$

$$h(x,\mathbf{m},\theta)h(x,\mathbf{m}^{\prime},\theta^{\prime})=C_{\mathbf{m},\mathbf{m}^{\prime},\theta,\theta^{\prime}}h(x,d(\mathbf{m},\mathbf{m}^{\prime}),e(\theta,\theta^{\prime})),$$

(20)

where $C_{\mathbf{m},\mathbf{m}^{\prime},\theta,\theta^{\prime}}$ is constant in $x$.

Denote by $\overset{\leftarrow}{\vartheta}_{i},\overset{\leftarrow}{\vartheta}_{i}^{\prime}$ the quantities defined in (17) computed backwards. Equivalently, these are computed as in (17) with data in reverse order, i.e. using $y_{n:0}$ in place of $y_{0:n}$, namely

The main difference between the above result and Theorem 4 in [50] lies in the fact that the support of the weights $\{\overset{\leftarrow}{\omega}_{\mathbf{m}}^{(i+1)},\mathbf{m}\in\mathcal{M}\}$ (which in [50] is denoted by $\overset{\leftarrow}{M}_{i|i+1:n}$) can be countably infinite and coincide with the whole of $\mathcal{M}$. Indeed, which points of $\mathcal{M}$ have positive weight are determined by the transition probabilities of the dual process, which in the present framework is no longer assumed to make only downward moves in $\mathcal{M}$. Section 6 will deal with this possibly infinite support for a concrete implementation of the inferential strategy.

The Cox–Ingersoll–Ross diffusion, also known as the square-root process, is widely used in financial mathematics for modelling short-term interest rates and stochastic volatility, see [10, 9, 29, 40, 37]. It also belongs to the class of continuous-state branching processes with immigration, arising as the large-population scaling limit of certain branching Markov chains [49] and as the time-evolving total mass of a Dawson–Watanabe branching measure-valued diffusion [20].

Let $X_{t}$ be a CIR diffusion on $\mathbb{R}_{+}$ that solves the one-dimensional SDE

where $\delta,\gamma,\sigma>0$, which is reversible with respect to the Gamma density $\pi=\text{Ga}(\delta/2,\gamma/\sigma^{2})$. The following proposition identifies a B&D process as dual to the CIR diffusion.

Assign now prior $\nu_{0}=\text{Ga}(\delta/2,\gamma/\sigma^{2})$ to $X_{0}$, and assume Poisson observations are collected at equally spaced intervals of length $\Delta$. Specifically, $Y|X_{t}=x\overset{\text{iid}}{\sim}\text{Po}(\tau x)$, for some $\tau>0$. By the well-known conjugacy to Gamma priors, we have that $X_{t}|Y=y\sim\text{Ga}(\delta/2+y,\gamma/\sigma^{2}+\tau)$. Without loss of generality, we can set $\tau=1$, which allows to interpret the update of the gamma rate parameter as the size of the conditioning data set. The filtering algorithm starts by first updating the prior $\nu_{0}$ to $\nu_{0|0}:=\phi_{Y_{0}}(\nu_{0})$. If we observe $Y_{0}=(Y_{0,1},\ldots,Y_{0,k})$ at time $0$, then $\nu_{0|0}$ is the law of $X_{0}|\sum_{j=1}^{k}y_{0,j}=m\sim\text{Ga}(\delta/2+m,\gamma/\sigma^{2}+k)$. Then $\nu_{0|0}$ is propagated forward for a $\Delta$ time interval, yielding $\nu_{1|0}:=\psi_{\Delta}(\nu_{0|0})$. In light of Proposition 4.1, an application of (6) to $\nu_{0|0}$ yields the infinite mixture

where $p_{m,n}(t)$ are the transition probabilities of $M_{t}$ in Proposition 4.1. Hence, the law of the signal is indexed by $\mathbb{Z}_{+}$, the state space of the dual process. While after the update at time 0 mass one is assigned to the sum of the observations $\sum_{j=1}^{k}y_{0,j}=m$, after the propagation the mass is spread over the whole $\mathbb{Z}_{+}$ by the effect of the dual process. We then observe $Y_{1}\sim f_{X_{1}}$ at time 1, which is used to update $\nu_{1|0}$ to $\nu_{1|1}$ and has the effect of shifting the probability masses of the mixture weights. For example, the weight $p_{m,n}(\Delta)$ in (36) is assigned to $n\in\mathbb{Z}_{+}$, but after the update based on $Y_{1}=(Y_{1,1},\ldots,Y_{1,k^{\prime}})$ it will be assigned to $n+m^{\prime}$ if $\sum_{j=1}^{k^{\prime}}y_{1,j}=m^{\prime}$, on top of being transformed according to (11). We then propagate forward again and proceed analogously.

When the current distribution of the signal, after the update, is given by a mixture of type $\sum_{m\in\mathbb{Z}_{+}}w_{m}\text{Ga}(\delta/2+m,\gamma/\sigma^{2}+k)$, it is enough to rearrange the mixture weights after the propagation step as done in (7).

The main difference with qualitatively similar equations found in [50] is now given by the transition probabilities $p_{m,n}(t)$ in (36), which are those of the B&D process in Proposition 4.1. Before tackling the problem of how to use the above expressions for inference, we try to provide further intuition of the extent and implications of such differences. To this end, consider the simplified parameterization $\alpha=\delta/2,\beta=\gamma/\sigma^{2},\sigma^{2}=1/2,\tau=1$, whereby one can check that the embedded chain of the B&D process of Proposition 4.1 has jump probabilities

Here $m,k$ are the same as in the left-hand side of (36), so $m/k$ is the sample mean. It is easily verified that $p_{m,m+1}<p_{m,m-1}$ if $m/k>\alpha/\beta$ and viceversa. Therefore, the dual evolves on $\mathbb{Z}_{+}$ so that it reverts $m/k$ to the prior mean $\alpha/\beta$. Indeed, the dual has Negative Binomial ergodic distribution $\text{NBin}\left(\alpha,\beta/(\beta+k)\right)$, whose mean is $k\alpha/\beta$, i.e., such that $m/k$ on average coincides with $\alpha/\beta$.