Online Smoothing for Diffusion Processes Observed with Noise

Research in Hidden Markov Models (HMMs) has – thus far – provided effective online algorithms for the estimation of expectations of the smoothing distribution for the case of a class of additive functionals of the underlying signal. Such methods necessitate knowledge of the transition density of the Markovian part of the model between observation times. We carry out a related exploration for the (common in applications) case when the signal corresponds to a diffusion process, thus we are faced with the challenge that such transition densities are typically unavailable. Standard data augmentation schemes that work with the multivariate density of a large enough number of imputed points of the continuous-time signal will lead to ineffective algorithms. The latter will have the abnormal characteristic that – for given Monte-Carlo iterates – the variability of the produced estimates will increase rapidly as the resolution of the imputation becomes finer. One of the ideas underpinning the work in this paper is that development of effective algorithms instead requires respecting the structural properties of the diffusion process, thus we build up imputation schemes on the infinite-dimensional diffusion pathspace itself. As a consequence, the time-discretised algorithm used in practice on a PC will be stable under mesh-refinement.

We consider continuous-time jump-diffusion models observed at discrete-time instances. The $d_{x}$-dimensional process, $X=\{X_{t};t\geq 0\}$, $d_{x}\geq 1$, is defined via the following time-homogeneous stochastic differential equation (SDE), with $X_{t-}:=\lim_{s\uparrow t}X_{t}$,

Solution $X$ is driven by the $d_{w}$-dimensional Brownian motion, $W=\{W_{t};t\geq 0\}$, $d_{w}\geq 1$, and the compound Poisson process, $J=\{J_{t};t\geq 0\}$. The SDE involves a drift function $b=b_{\theta}:\mathbb{R}^{d_{x}}\mapsto\mathbb{R}^{d_{x}}$ and coefficient matrix $\sigma=\sigma_{\theta}:\mathbb{R}^{d_{x}}\mapsto\mathbb{R}^{d_{x}\times d_{w}}$, for parameter $\theta\in\mathbb{R}^{p}$, $p\geq 1$. Let $\{N_{t};t\geq 0\}$ be a Poisson process with intensity function $\lambda_{\theta}(\cdot)$, and $\{\xi_{k}\}_{k\geq 1}$ i.i.d. sequence of random variables with Lebesgue density $h_{\theta}(\cdot)$; the càdlàg process $J$ is determined as $J_{t}=J_{t,\theta}=\sum_{i=1}^{N_{t}}\xi_{i}$. We work under standard assumptions (e.g. linear growth, Lipschitz continuity for $b$, $\sigma$) that guarantee a unique global solution of (1.1), in a weak or strong sense, see e.g. Øksendal and Sulem, (2007).

SDE (1.1) is observed with noise at discrete-time instances $0=t_{0}<t_{1}<t_{2}<\cdots<t_{n}$, $n\geq 1$. Without loss of generality, we assume equidistant observation times, with $\Delta:=t_{1}-t_{0}$. We consider data $Y_{t_{0}},\ldots,Y_{t_{n}}$, and for simplicity we set,

Let $\mathcal{F}_{i}$ be a filtration generated by $X_{s}$ for $s\in[t_{i-1},t_{i}],0<i\leq n$. We assume,

for conditional distribution $g_{\theta}(\cdot|Y_{t_{i-1}},\mathcal{F}_{i})$ on $\mathbb{R}^{d_{y}}$, $d_{y}\geq 1$, under convention $Y_{t_{-1}}=y_{-1}=\emptyset$. We write,

where $f_{\theta}(dx_{i}|x_{i-1})$ is the transition distribution of the driving SDE process (1.1). We consider the density functions of $g_{\theta}(dy_{i}|y_{i-1},\mathcal{F}_{i})$ and $f_{\theta}(dx_{i}|x_{i-1})$, and – with some abuse of notation – we write $g_{\theta}(dy_{i}|y_{i-1},\mathcal{F}_{i})=g_{\theta}(y_{i}|y_{i-1},\mathcal{F}_{i})dy_{i}$, $f_{\theta}(dx_{i}|x_{i-1})=f_{\theta}(x_{i}|x_{i-1})dx_{i}$, where $dy_{i}$, $dx_{i}$ denote Lebesgue measures. Our work develops under the following regime.

Models defined via (1.1)-(1.3) are extensively used, e.g., in finance and econometrics, for instance for capturing the market microstructure noise, see Aït-Sahalia et al., (2005); Aït-Sahalia and Yu, (2008); Hansen and Lunde, (2006). The above setting belongs to the general class of HMMs, with a signal defined in continuous-time. See Cappé et al., (2005); Douc et al., (2014) for a general treatment of HMMs fully specified in discrete-time. A number of methods have been suggested in the literature for approximating the unavailable transition density – mainly in the case of processes without jumps – including: asymptotic expansion techniques (Aït-Sahalia et al.,, 2005; Aït-Sahalia,, 2002, 2008; Kessler,, 1997); martingale estimating functions (Kessler and Sørensen,, 1999); generalized method of moments (Hansen and Scheinkman,, 1993); Monte-Carlo approaches (Wagner,, 1989; Durham and Gallant,, 2002; Beskos et al.,, 2006). See, e.g., Kessler et al., (2012) for a detailed review.

For a given sequence $\{a_{m}\}_{m}$, we use the notation $a_{i:j}:=(a_{i},\ldots,a_{j})$, for integers $i\leq j$. Let $p_{\theta}(y_{0:n})$ denote the joint density of $y_{0:n}$. Throughout the paper, $p_{\theta}(\cdot)$ is used generically to represent probability distributions or densities of random variables appearing as its arguments.

Consider the maximum likelihood estimator (MLE),

Except for limited cases, one cannot obtain the MLE analytically for HMMs (even for discrete-time signal) due to the intractability of $p_{\theta}(y_{0:n})$.

We have set up the modelling context for this work. The main contributions of the paper in this setting – several of which relate with overcoming the intractability of the transition density of the SDE, and developing a methodology that is well-posed on the infinite-dimensional pathspace – will be as follows:

• (i)

  We present an online algorithm that delivers Monte-Carlo estimators of smoothing expectations,

  $$\mathcal{S}_{\theta,n}=\int S_{\theta}(\mathbf{x}_{0:n})p_{\theta}(d\mathbf{x}_{0:n}|y_{0:n}),\quad n\geq 1,$$

  (1.4)

  for the class of additive functionals $S_{\theta}(\cdot)$ of the structure,

  $$S_{\theta}(\mathbf{x}_{0:n})=\sum_{k=0}^{n}s_{\theta,k}(x_{k-1},\mathbf{x}_{k}),$$

  (1.5)

  under the conventions $x_{-1}=\emptyset$, $\mathbf{x}_{0}=x_{0}$. The bold type notation $\mathbf{x}_{k}$, $k\geq 0$, is reserved for carefully defined variables involving elements of the infinite-dimensional pathspace. The specific construction of $\mathbf{x}_{k}$ will depend on the model at hand, and will be explained in the main part of the paper. We note that the class of additive functionals considered supposes additivity over time, and therefore includes the general case of evaluating the score function. The online solution of the smoothing problem is often used as the means to solve some concrete inferential problems - see, e.g., (ii) below.

• (ii)

  We take advantage of the new approach to show numerical applications, with emphasis on carrying out online parameter inference for the designated class of models via a gradient-ascent approach (in a Robbins-Monro stochastic gradient framework). A critical aspect of this particular online algorithm (partly likelihood based, when concerned with parameter estimation; partly Bayesian, with regards to identification of filtering/smoothing expectations) is that it delivers estimates of the evolving score function, of the model parameters, together with particle representations of the filtering distributions, through a single passage of the data. This is a unique favourable algorithmic characteristic, when constrasted with alternative algorithms with similar objectives, such as, e.g., Particle MCMC (Andrieu et al.,, 2010), or SMC² (Chopin et al.,, 2013).

• (iii)

  In this work, we will not characterise analytically the size of the time-discretisation bias relevant to the SDE models at hand, and are content that: (I) the bias can be decreased by increasing the resolution of the numerical scheme (typically an Eyler-Maruyama one, or some other Taylor scheme, see e.g. Kloeden and Platen, (2013)); (II) critically, the Monte-Carlo algorithms are developed in a manner that the variance does not increase (in the limit, up to infinity) when increasing the resolution of the time-discretisation method; to achieve such an effect, the algorithms are (purposely) defined on the infinite-dimensional pathspace, and SDE paths are only discretised when implementing the algorithm on a PC (to allow, necessarily, for finite computations).

• (iv)

  Our method draws inspiration from earlier works, in the context of online filtering for discrete-time HMMs and infinite-dimensional pathspace MCMC methods. The complete construct is novel; one consequence of this is that it is applicable, in principle, for a wide class of SDEs, under the following weak assumption (relatively to restrictive conditions often imposed in the literature of infinite-dimensional MCMC methods).

  Assumption 2.

  The diffusion covariance matrix function,

  $$\Sigma_{\theta}(v):=\sigma_{\theta}(v)\sigma_{\theta}(v)^{\top}\in\mathbb{R}^{d_{x}\times d_{x}},$$

  is invertible, for all relevant $v$, $\theta$.

Thus, the methodology does not apply as defined here only for the class of hypoelliptic SDEs.

An elegant solution to the online smoothing problem posed above in (i), for the case of a standard HMM with discrete-time signal of known transition density $f_{\theta}(x^{\prime}|x)$, is given in Del Moral et al., (2010); Poyiadjis et al., (2011). Our own work overcomes the unavailability of the transition density in the continuous-time scenario by following the above literature but augmenting the hidden state with the complete continuous-time SDE path. Related augmentation approaches in this setting – though for different inferential objectives – have appeared in Fearnhead et al., (2008); Ströjby and Olsson, (2009); Gloaguen et al., (2018), where the auxiliary variables are derived via the Poisson estimator of transition densities for SDEs (under strict conditions on the class of SDEs; no jumps), introduced in Beskos et al., (2006), and in Särkkä and Sottinen, (2008) where the augmentation involves indeed the continuous-time path (the objective therein is to solve the filtering problem and the method is applicable for SDEs with additive Wiener noise; no jump processes are considered).

A Motivating Example: Fig. 1.1 (Top-Panel) shows estimates of the score function, evaluated at the true parameter value $\theta=\theta^{\dagger}$, for parameter $\theta_{3}$ of the Ornstein–Uhlenbeck (O–U) process, $dX_{t}=\theta_{1}(\theta_{2}-X_{t})dt+\theta_{3}dW_{t}$, $X_{0}=0.0$, for $n=10$ observations $y_{i}=x_{i}+\epsilon_{i}$, $\epsilon_{i}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mathcal{N}(0,0.1^{2})$, $1\leq i\leq n$. Data were simulated from $\theta^{\dagger}=(0.5,0.0,0.4)$ with an Euler-Maruyama scheme of $M^{\dagger}=1,000$ grid points per unit of time. Fig. 1.1 (Top-Panel) illustrates the ‘abnormal’ effect of a standard data-augmentation scheme, where for $N=100$ particles, the Monte-Carlo method (see later sections for details) produces estimates of increasing variability as algorithmic resolution increases with $M=10,50,100,200$ – i.e., as it approaches the ‘true’ resolution used for the data generation. In contrast, Fig. 1.1 (Bottom Panel) shows results for the same score function as estimated by our proposed method in this paper. Our method is well-defined on the infinite-dimensional pathspace, and, consequently, it manages to robustly estimate the score function under mesh-refinement.

The rest of the paper is organised as follows. Section 2 reviews the Forward-Only algorithm for the online approximation of expectations of a class of additive functionals described in Del Moral et al., (2010). Section 3 sets up the framework for the treatment of pathspace-valued SDEs, first for the conceptually simpler case of SDEs without jumps, and then proceeding with incorporating jumps. Section 4 provides the complete online approximation algorithm, constructed on the infinite-dimensional pathspace. Section 5 discusses the adaptation of the developed methodology for the purposes of online inference for unknown parameters of the given SDE model. Section 6 shows numerical applications of the developed methodology. Section 7 contains conclusions and directions for future research.

A bootstrap filter (Gordon et al.,, 1993) is applicable in the continuous-time setting, as it only requires forward sampling of the underlying signal $\{X_{t}\}_{t\geq 0}$; this is trivially possible – under numerous approaches – and is typically associated with the introduction of some time-discretisation bias (Kloeden and Platen,, 2013). However, the transition density is still required for the smoothing problem we have posed in the Introduction. In this section, we assume a standard discrete-time HMM, with initial distribution $p_{\theta}(dx_{0})$, transition density $f_{\theta}(x^{\prime}|x)$, and likelihood $g_{\theta}(y|x)$, for appropriate $x,x^{\prime}\in\mathbb{R}^{d_{x}}$, $y\in\mathbb{R}^{d_{y}}$, and review the online algorithm developed in Del Moral et al., (2010) for this setting. Implementation of the bootstrap filter provides an approximation of the smoothing distribution $p_{\theta}(dx_{0:n}|y_{0:n})$ by following the geneology of the particles. This method is studied, e.g., in Cappe, (2009); Dahlhaus and Neddermeyer, (2010). Let $\{x_{0:n}^{(i)},W_{n}^{(i)}\}_{i=1}^{N}$, $N\geq 1$, be a particle approximation of the smoothing distribution $p_{\theta}(dx_{0:n}|y_{0:n})$, in the sense that we have the estimate,

with $\delta_{x_{0:n}^{(i)}}(dx_{0:n})$ the Dirac measure with an atom at $x_{0:n}^{(i)}$. Then, replacing $p_{\theta}(x_{0:n}|y_{0:n})$ with its estimate in (2.1) provides consistent estimators of expectations of the HMM smoothing distributions. Though the method is online and the computational cost per time step is $\mathcal{O}(N)$, it typically suffers from the well-documented path-degeneracy problem – as illustrated via theoretical results or numerically (Del Moral et al.,, 2010; Kantas et al.,, 2015). That is, as $n$ increases, the particles representing $p_{\theta}(dx_{0:n}|y_{0:n})$ obtained by the above method will eventually all share the same ancestral particle due to the resampling steps, and the approximation collapses for big enough $n$. This is well-understood not to be a solution to the approximation of the smoothing distribution for practical applications.

An approach which overcomes path-degeneracy is the Forward Filtering Backward Smoothing (FFBS) algorithm of Doucet et al., (2000). We briefly review the method here, following closely the notation and development in Del Moral et al., (2010). In the forward direction, assume that a filtering algorithm (e.g. bootstrap) has provided a particle approximation of the filtering distribution $p_{\theta}(dx_{k-1}|y_{0:k-1})$ – assuming a relevant $k$ –,

for weighted particles $\{x_{k-1}^{(i)},W_{k-1}^{(i)}\}_{i=1}^{N}$. In the backward direction, assume that one is given the particle approximation of the marginal smoothing distribution $p_{\theta}(dx_{k}|y_{0:n})$,

One has that (Kitagawa,, 1987),

Using (2.2)-(2.3), and based on equation (2.4), we obtain the approximation,

Recalling the expectation of additive functionals in (1.4)-(1.5) – where now, in the discrete-time setting, we can ignore the bold $\mathbf{x}_{k}$ elements, and simply use $x_{k}$ instead – the above calculations give rise to the following estimator of the target quantity $\mathcal{S}_{\theta,n}$ in (1.5),

To be able to apply the above method, the marginal smoothing approximation in (2.3) is obtained via a backward recursive approach. In particular, starting from $k=n$ (where the approximation is provided by the standard forward particle filter), one proceeds as follows. Given $k$, the quantity for $k-1$ is directly obtained by integrating out $x_{k}$ in (2.5), thus we have,

for the normalised weights,

Notice that – in this version of FFBS – the same particles $\{x_{k}^{(i)}\}_{i=1}^{N}$ are used in both directions (the ones before resampling at the forward filter), but with different weights.

An important development made in Del Moral et al., (2010) is transforming the above offline algorithm into an online one. This is achieved by consideration of the sequence of instrumental functionals,

Notice that, first,

We also have that,

– see Proposition 2.1 of Del Moral et al., (2010) for a (simple) proof. This recursion provides an online – forward-only – advancement of FFBS for estimating the smoothing expectation of additive functionals. The complete method is summarised in Algorithm 1: one key ingredient is that, during the recursion, values of the functional $T_{\theta,n}(x_{n})$ are only required at the discrete positions $x_{n}^{(i)}$ determined by the forward particle filter.

In the SDE context, under Assumption 1, the transition density $f_{\theta}(\cdot|\cdot)$ is considered intractable, thus Algorithm 1 – apart from serving as a review of the method in Del Moral et al., (2010) – does not appear to be practical in the continuous-time case.

To overcome the intractability of the transition density $f_{\theta}(\cdot|\cdot)$ of the SDE, we will work with an algorithm that is defined in continuous-time and makes use of the complete SDE path-particles in its development. The new method has connections with earlier works in the literature. Särkkä and Sottinen, (2008) focus on the filtering problem for a class of models related to (1.1)-(1.3), and come up with an approach that requires the complete SDE path, for a limited class of diffusions with additive noise and no jumps. Fearnhead et al., (2008) also deal with the filtering problem, and – equipped with an unbiased estimator of the unknown transition density – recast the problem as one of filtering over an augmented space that incorporates the randomness for the unbiased estimate. The method is accompanied by strict conditions on the drift and diffusion coefficient (one should be able to transform the SDE – no jumps – into one of unit diffusion coefficient; the drift of the SDE must have a gradient form). Our contribution requires, in principle, solely the diffusion coefficient invertibility Assumption 2; arguably, the weakened condition we require stems from the fact that our approach appears as the relatively most ‘natural’ extension (compared to alternatives) of the standard discrete-time algorithm of Del Moral et al., (2010).

The latter discrete-time method requires the density $f_{\theta}(x^{\prime}|x)=f_{\theta}(dx^{\prime}|x)/dx^{\prime}$. In continuous-time, we obtain an analytically available Radon-Nikodym derivative of $p_{\theta}(d\mathbf{x}^{\prime}|x)$, for a properly defined variate $\mathbf{x}^{\prime}$ that involves information about the continuous-time path for moving from $x$ to $x^{\prime}$ within time $\Delta$. We will give the complete algorithm in Section 4. In this section, we prepare the ground via carefully determining $\mathbf{x}^{\prime}$ given $x$, and calculating the relevant densities to be later plugged in into our method.