Online sequential Monte Carlo smoother for partially observed diffusion processes

This paper introduces a new algorithm to approximate smoothed additive functionals of partially observed diffusion processes. This method relies on a new sequential Monte Carlo method which allows to compute such approximations online, i.e., as the observations are received, and with a computational complexity growing linearly with the number of Monte Carlo samples. The original algorithm cannot be used in the case of partially observed stochastic differential equations since the transition density of the latent data is usually unknown. We prove that it may be extended to partially observed continuous processes by replacing this unknown quantity by an unbiased estimator obtained for instance using general Poisson estimators. This estimator is proved to be consistent and its performance are illustrated using data from two models.


Introduction
This paper introduces a new algorithm to solve the smoothing problem for partially observed continuous time stochastic processes. In this setting, the hidden state process (X t ) t≥0 is assumed to be a solution to a stochastic differential equation (SDE) and the only information available is given by noisy observations (Y k ) 0≤k≤n of the states (X k ) 0≤k≤n at some discrete time points (t k ) 0≤k≤n . The bivariate stochastic process {(X k , Y k )} 0≤k≤n is a state space model such that conditional on the state sequence (X k ) 0≤k≤n the observations (Y k ) 0≤k≤n are independent and for all 0 ≤ ≤ n the conditional distribution of Y given {X k } 0≤k≤n depends on X only.
Statistical inference for partially observed state sequences often requires to solve bayesian filtering and smoothing problems, i.e. the computation of the posterior distributions of sequences of hidden states given observations. The filtering problem refers to the estimation, for each 0 ≤ k ≤ n, of the distributions of the hidden state X k given the observations (Y 0 , . . . , Y k ). Smoothing stands for the estimation of the distributions of the sequence of states (X k , . . . , X p ) given observations (Y 0 , . . . , Y ) with 0 ≤ k ≤ p ≤ ≤ n. These posterior distributions are crucial to compute maximum likelihood estimators of unknown parameters using the observations (Y 0 , . . . , Y n ) only. For instance, the E-step of the EM algorithm introduced in [7] boils down to the computation of a conditional expectation of an additive functional of the hidden states given all the observations up to time n. Similarly, by Fisher's identity, recursive maximum likelihood estimates may be computed using the gradient of the loglikelihood which can be written as the conditional expectation of an additive functional of the hidden states. See [5,Chapter 10 and 11], [15,19,20,27] for further references on the use of these smoothed expectations of additive functionals applied to maximum likelihood parameter inference in latent data models.
The exact computation of these expectations is usually not possible in the case of partially observed diffusions. In this paper, we propose to use Sequential Monte Carlo (SMC) methods to approximate smoothing distributions with random particles associated with importance weights. [13,18] introduced the first particle filters and smoothers for state space models by combining importance sampling steps to propagate particles with resampling steps to duplicate or discard particles according to their importance weights. Unfortunately, these methods cannot be applied directly to partially observed stochastic differential equations since some elementary quantities, such as transition densities of the hidden states, are not available explicitly. Discretization procedures may be used to approximate transition densities, for instance the Euler-Maruyama method, the Ozaki discretization which proposes a linear approximation of the drift coefficient between two observations [25,28], or Gaussian based approximations using Taylor expansions of the posterior mean and variance of an observation given the observation at the previous time step, [16,17,29]. Other approaches based on Hermite polynomials expansion were also introduced by [1, 2, 3] and extended in several directions recently, see [21] and all the references on the approximation of transition densities therein. However, even the most recent discretization based approximations of the transition densities induce a systematic bias of particle based approximations of posterior distributions, see for instance [6]. To overcome this difficulty, [11] proposed to solve the filtering problem by combining SMC methods with an unbiased estimate of the transition densities based on the generalized Poisson estimator (GPE). In this case, only the Monte Carlo error has to be controlled as there is no Taylor expansion to approximate unknown transition densities.
The only solution to solve the smoothing problem for partially observed SDE using SMC methods has been proposed in [23] and extends the fixed-lag smoother of [22]. Using forgetting properties of the hidden chain, the algorithm improves the performance of [11] to approximate smoothing distributions but at the cost of a bias that does not vanish as the number of particles grows to infinity. In the case of discrete time state space models, approximations of the smoothing distributions may also be obtained using the Forward Filtering Backward Smoothing algorithm (FFBS) and the Forward Filtering Backward Simulation algorithm (FFBSi) developed respectively in [18,14,9] and [12]. Both algorithms require first a forward pass which produces a set of particles and weights approximating the sequence of filtering distributions up to time n. Then, a backward pass is performed to compute new weights (FFBS) or sample trajectories (FFBSi) in order to approximate the smoothing distributions. Recently, [24] proposed a new SMC algorithm, the particle-based rapid incremental smoother (PaRIS), to approximate on-the-fly (i.e. using the observations as they are received) smoothed expectations of additive functionals. Unlike the FFBS algorithm, the complexity of this algorithm grows only linearly with the number of particles N and contrary to the FFBSi algorithm, no backward pass is required.
In this paper, we extend the use of PaRIS algorithm to partially observed SDE. The proposed algorithm allows to approximate smoothed expectations of additive functionals online and with a complexity growing only linearly with the number of particles. The crucial and simple result (Lemma 1) of the application of PaRIS algorithm to SDE is that the accept reject mechanism introduced in [8] ensuring the linear complexity of the procedure is still correct when the transition densities are replaced by unbiased estimates. The usual FFBS and FFBSi algorithms may not be extended this easily since they both require the computation of weights defined as ratios involving the transition densities, thus replacing these unknown quantities by unbiased estimates does not lead to unbiased estimators of the weights. The proposed Generalized Random version of PaRIS algorithm, hereafter named GRand PaRIS algorithm, may be applied to general hidden Markov models whose Markovian dynamics is ruled by a stochastic differential equation (one of the first two domains defined in [4]) but also to any general state space model where the transition density of the hidden chain may be estimated unbiasedly.
Section 2 describes the proposed algorithm to approximate smoothed additive functionals using unbiased estimates of the transition density of the hidden states and details the application of this algorithm when the transition density may be approximated using a GPE. In Section 3, classical convergence results for SMC smoothers are extended to the setting of this paper and illustrated with numerical experiments in Section 4. All proofs are postponed to Appendix A.
2 The Generalized Random PaRIS algorithm (X t ) t≥0 is defined as a weak solution to the following SDE in R d : where (W t ) t≥0 is a standard Brownian motion. It is assumed that α is of the form The solution to (1) is supposed to be partially observed at times t 0 = 0, . . . , t n through an observation process (Y k ) 0≤k≤n in (R m ) n+1 . For all 0 ≤ k ≤ n, the distribution of Y k given X k := X t k has a density with respect to a reference measure λ on R m given by g(X k , ·) = g k (X k ). The distribution of X 0 has a density with respect to a reference measure µ on R d given by χ. For all 0 ≤ k ≤ n − 1, the conditional distribution of X k+1 given X k has a density q k (X k , ·) with respect to µ. Let 0 ≤ k ≤ k ≤ n, the joint smoothing distributions of the hidden states are defined, for all measurable function h on (R d ) k −k+1 , by: For all 0 ≤ k ≤ n, φ k = φ k:k|k denote the filtering distributions. The aim of this section is to detail the extension of PaRIS algorithm to approximate expectations of the form when the transition density of the hidden states is not available explicitly and where {h k } n−1 k=0 are given functions on R d × R d . The algorithm is based on the following link between the filtering and smoothing distributions for additive functionals, see [24]: The approximation of (3) requires first to approximate the sequence of filtering distributions. Sequential Monte Carlo methods provide an efficient and simple solution to obtain these approximations using sets of particles are sampled independently according to ξ 0 ∼ η 0 , where η 0 is a probability density with respect to µ. Then, ξ 0 is associated with the importance weights ω 0 = χ(ξ 0 )g 0 (ξ 0 )/η 0 (ξ 0 ). For any bounded and measurable function h defined on R d , the expectation φ 0 [h] is approximated by Then, for 1 ≤ k ≤ n, using {(ξ k−1 , ω k−1 )} N =1 , the auxiliary particle filter of [26] samples pairs {(I k , ξ k )} N =1 of indices and particles using an instrumental transition density p k on R d × R d and an adjustment multiplier function ϑ k on R d . Each new particle ξ k and weight ω k at time k are computing following these steps: -sample ξ k using this chosen particle according to ξ k ∼ p k (ξ I k k−1 , ·) ; -associate the particle ξ k with the importance weight: . (4) PaRIS algorithm uses the same decomposition as the FFBS algorithm introduced in [10] and the FFBSi algorithm proposed by [12] to approximate smoothing distributions. It combines both the forward only version of the FFBS algorithm with the sampling mechanism of the FFBSi algorithm. It does not produce an approximation of the smoothing distributions but of the smoothed expectation of a fixed additive functional and thus may be used to approximate (2). Its crucial property is that it does not require a backward pass, the smoothed expectation is computed on-the-fly with the particle filter and no storage of the particles or weights is needed. PaRIS algorithm relies on the following fundamental property of T k [H k ] when H k is as in (2): Therefore, [24] introduces sufficient statistics τ i k (starting with τ i where Computing exactly these approximations would lead to a complexity growing quadratically with N because of the normalizing constant in (6). Therefore, PaRIS algorithm samples particles in the set {ξ j k−1 } N j=1 with probabilities Λ N k (i, ·) to approximate the expectation (5) and produce τ i k . ChoosingÑ ≥ 1, at each time step 0 ≤ k ≤ n − 1 these statistics are updated according to the following steps.
(i) Run one step of a particle filter to produce {(ξ k , ω k )} for 1 ≤ ≤ N .
Then, (2) is approximated by As proved in [24], the algorithm is asymptotically consistent (as N goes to infinity) for any precision parameterÑ . However, there is a significant qualitative difference between the casesÑ = 1 andÑ ≥ 2. As for the FFBSi algorithm, when there exists σ + such that 0 < q k < σ + , PaRIS algorithm may be implemented with O(N ) complexity using the accept-reject mechanism of [8].
In general situations, PaRIS algorithm cannot be used for stochastic differential equations as q k is unknown. Therefore, the computation of the importance weights ω k and of the acceptance ratio of [8] is not tractable. Following [11,23], filtering weights can be approximated by replacing q k (ξ k , ξ i k+1 ) by an unbiased estimator q k (ξ k , ξ i k+1 ; ζ k ), where ζ k is a random variable in R q such that: where, for all 0 ≤ k ≤ n, Practical choices for ζ k are discussed below, see for instance (9) which presents the choice made for the implementation of such estimators in our context. In the case where q k is unknown, the filtering weights in (4) then become: Therefore, to obtain a generalized random version of PaRIS algorithm, we only need to be able to sample from the discrete probability distribution Λ N k (i, ·) in the case of SDE based HMM. Consider the following assumption: for all 0 ≤ k ≤ n, there exists a random variableσ k + measurable with respect to G N k+1 such that, Note that Lemma 1 still holds if assumption (A1) is relaxed and replaced by one of the two following assumptions: It is worth noting that under assumptions (A1) or (A2), the linear complexity property of PaRIS algorithm still holds, whereas if only assumption (A3) holds, the algorithm has a quadratic complexity.
Bounded estimator of q k For x, y ∈ R d , by Girsanov and Ito's formulas, the transition density q k (x, y) of (1) satisfies, with ∆ k = t k+1 − t k , where W x,y,∆ k is the law of Brownian bridge starting at x at 0 and hitting y at ∆ k , (w t ) 0≤t≤∆ k is such a Brownian bridge, ϕ ∆ k (x, y) is the p.d.f. of a normal distribution with mean x and variance ∆ k , evaluated at y and φ : R d → R is defined as: with the Laplace operator. Assume that there exist random variables L w and U w such that for all 0 ≤ s ≤ ∆ k , L w ≤ φ(w s ) ≤ U w . Let κ be a random variable taking values in N with distribution µ and (U j ) 1≤j≤κ be independent uniform random variables on [0, ∆ k ], and ζ k = {κ, w, U 1 , . . . , U κ } . As shown in [11], a positive unbiased estimator is given by Interesting choices of µ are discussed in [11] and we focus here on the so called GPE-1, where µ is a Poisson distribution with intensity (U w − L w )∆ k . In that case, the estimator (8) becomes: On the r.h.s. of (9), the product over κ elements is bounded by 1, therefore, a sufficient condition to satisfy of the assumptions (A1)-(A3) is that the function: is upper bounded almost surely byσ k + . In particular, if L w is bounded almost surely, (10) always satisfies assumption (A3) and Algorithm 1 can be used. This condition is always satisfied for models in the domains D 1 and D 2 defined in [4], i.e. domains for which the exact algorithms EA1 and EA2 can be used. When (A1) or (A2) holds, it can be nonetheless of practical interest to choose the boundσ k + corresponding to (A3). Indeed, this might increase significantly the acceptance rate of the algorithm, and therefore reduce the number of drawings of the random variable ζ, which has a much higher cost than the computation of ρ, as it requires simulations of Brownian Bridges. Moreover, this latter option can also avoid numerical optimization if no analytical expression ofσ k + is available. In practice, we found this option more efficient in terms of computation time when N has moderate values.
Proof. See appendix A Proposition 1. Assume that H1 and H2 hold and that for all 1 ≤ k ≤ n, osc(h k ) < +∞. For all 0 ≤ k ≤ n and all N ≥ 1, there exist b k , c k > 0 such that for all N ≥ 1 and all ε ∈ R + , Proof. See appendix A

Numerical experiments
This section investigates the performance of the proposed algorithm with the sine and log-growth models. In both cases, the proposal distribution p k is chosen as the following approximation of the optimal filter (or the fully adapted particle filter in the terminology of [26]): the Euler approximation of equation (1). As the observation model is linear and Gaussian, the proposal distribution is therefore Gaussian with explicit mean and variance. In order to evaluate the performance of the proposed algorithm, the following strategy has been chosen. We compare the estimation of the EM intermediate quantity with the one obtained by the fixed lag method of [23], for different values of the lag (namely, 1,2,5,10,50). The particle approximation of Q(θ, θ) for each model is computed using each algorithm, see Figure 1 for the SINE model (and respectively Figure 3 for the log-growth model). This estimation is performed 200 times to obtain the estimates Q 1 , . . . , Q 200 , usingÑ = 2 particles for PaRIS algorithm, and M = 30 replications for the Monte Carlo approximation q k of each q k . Moreover, the E step requires the computation of a quantity such as (2) with h k = log g k +log q k . log q k is not available explicitly and is approximated using the unbiased estimator proposed in [23, Appendix B] based on 30 independent Monte Carlo simulations. The intermediate quantity of the EM algorithm is also estimated with our algorithm 30 times using N = 5000 particles, the reference value is then computed as the arithmetic mean of these 30 estimations, and denoted by Q . Figure 1 (resp. 3 ) shows this estimate for an example on one simulated data set. The GRand Paris algorithm is performed using N = 400 particles in both cases, the fixed lag technique using N = 1600 so that both estimations require similar computational times.

The SINE model
The performance of the GRand PaRIS algorithm are first highlighted using the SINE model, where (X t ) t≥0 is supposed to be the solution to: This simple model has no explicit transition density, however GPE estimators may be computed by simulating Brownian bridges. The process solution to (11) is observed regularly at times t 0 = 0, . . . , t 100 = 50 through the observation process (Y k ) 0≤k≤100 : where the (ε) 0≤k≤100 are i.i. d. N (0, 1). In the example displayed on Figure 1, we set θ = 0. In that case, the function ρ ∆ k defined in (10) can be upper bounded either on (x, y) or only on y, the GRand PaRIS algorithm has therefore a linear complexity. This same experiment was reproduced on 100 different simulated data sets. For each simulation s, the empirical absolute relative bias arb s and the empirical absolute coefficient of variation acv s are computed as where m( Q s ) and σ( Q s ) are the empirical mean and standard deviation of the sample Q s 1 , . . . , Q s 200 . For each estimation method, the resulting distributions of arb 1 , . . . , arb 100 and acv 1 , . . . , acv 100 are shown on Figure 2.
The GRand PaRIS algorithm outperforms the fixed lag methods for any value of the lag as the bias is the lowest (it is already negligible for N = 400) and with a lower variance than fixed lag estimates with negligible bias (i.e., in this case, lags larger than 10). Small lags lead to strongly biased estimates for the fixed lag method, and unbiased estimates are at the cost of a large variance. It is worth noting here that the lag for which the bias is small is model dependent.

Log-growth model
Following [4] and [24], the performance of the proposed algorithm are also illustrated with the log-growth model defined by:  In order to use the exact algorithms of [4] and the GPE of [11], we consider (15) after the Lamperti transform, i.e., the process defined by X t = η(Z t ), with η(z) := − log(z)/σ, which satisfies the following SDE: In this case, the conditions of the Exact Algorithm 2 defined in [4] are satisfied, as for any m ∈ R there exists U m such that for all x ≥ m, ψ(x) := α 2 (x)+α (x) ≤ U m . Moreover, ψ is lower bounded uniformly by L. Then, GPE estimators may be computed by simulating the minimum of a Brownian bridge, and simulating Bessel bridges conditionally to this minimum, as proposed by [4].
The process solution to (16) is observed regularly at times t 0 = 0, . . . , t 50 = 100 through the observation process (Y k ) 0≤k≤50 defined as: where the (ε k ) 0≤k≤50 are i.i.d. N (0, σ 2 obs ). The parameters are given by In that case, the ρ ∆ k function defined in (10)  The results for the fixed lag technique are similar to the ones presented in [23, Figure 1] on the same model. For small lags, the variance of the estimates is small, but the estimation is highly biased. The bias rapidly decreases as the lag increases, together with a great increase of variance. Again, the GRand PaRIS algorithm outperforms the fixed lag smoother as it shows a similar (vanishing) bias as the fixed lag for the largest lag and a smaller variance than the fixed lags estimates with negligible bias.

Conclusions
This paper presents a new online SMC smoother for partially observed differential equations. This algorithm relies on an accept-reject procedure inspired from the recent PaRIS algorithm. The main result of the article for practical applications is that the mechanism of this procedure remains valid when the transition density is approximated by a an unbiased positive estimator. The proposed procedure outperforms the existing fixed lag smoother for SDE of [23], as it does not introduce  an intrinsic and non vanishing bias. In addition, numerical simulations highlight a better variance using data from two different models. It can be implemented for the class of models D 1 and D 2 defined in [4] with a linear complexity in N .

A Proofs
Proof of Lemma 1. Let τ be the first time draws are accepted in the accept-reject mechanism. For all ≥ 1, write Let h be a function defined on {1, . . . , N }, which concludes the proof.
Proof of Lemma 2. The independence is ensured by the mechanism of SMC methods. By (7), Note that by Lemma 1, Since τ i k+1 and ζ k are independent conditionally to G N k+1 : .
Moreover, conditionally to F N k , the probability density function of (ξ i k+1 , I i k+1 ) is given by . Therefore, this yields: which concludes the proof.
Proof of Proposition 1. The results is proved by induction. At time k = 0, the result holds using that for all 1 ≤ i ≤ N , ρ i 0 = 0 and the convention T 0 [h 0 ] = 0. In addition, φ N 0 is a standard importance sampler estimator of φ 0 with ω i 0 ≤ | ω 0 | ∞ so that for any bounded function h on X, Assume the results holds for k ≥ 1 and that ϑ k+1 = 1 for simplicity. Write . By Lemma 2, the random variables { ω i k+1 τ i k+1 } N i=1 are independent conditionally on F N k and by H2, Therefore, by Hoeffding inequality, On the other hand, By [24,Lemma 11], φ k [Υ k ] = 0 which implies by the induction assumption that Then, P (|a N | ≥ ε) ≤ b k exp −c k N ε 2 . Similarly, as b N ≤ | ω k | ∞ , by Hoeffding inequality, Note that By the induction assumption, The proof is completed using Lemma 3. Then, Proof. See [8].