Unbiased inference for discretely observed hidden Markov model diffusions
Neil K. Chada Jordan Franks Ajay Jasra Kody J. H. Law Matti Vihola

 







  




























hypA
hypAhypA
  
  (AhypA)
hypD
hypDhypD
  
  (D)


 
 









































t  
ℝ
𝒩

t  
ℝ
ℙ
𝔼
𝒩

plain
theoremTheorem
corollary[theorem]Corollary
proposition[theorem]Proposition
lemma[theorem]Lemma

definition
definition[theorem]Definition
assumption[theorem]Assumption
condition[theorem]Condition
example[theorem]Example
recommendationRecommendation


assumptionxCondition 
  

remark
remark[theorem]Remark





 
 
 

 
 
 

 
 
 
 

 
 

 

 
 

 
 
 
 

 
 

 
 
 
 
 
 




 
 
 
 


 
 
 
 


 



 

 

 




 
 








Unbiased inference for diffusions]Unbiased inference for discretely observed hidden Markov model diffusions

N. K. Chada, J. Franks, A. Jasra, K. J. H. Law & M. Vihola]Neil K. Chada, Jordan Franks, Ajay Jasra, Kody J. H. Law &  Matti Vihola
N. K. Chada & A. Jasra, Applied Mathematics and Computational Science,
King Abdullah University of Science and Technology, Thuwal, 23955, KSA. Email:  &   
J. Franks, School of Mathematics, Statistics and Physics,
  Herschel Building, NE1 7RU Newcastle University, UK.  Email: 
M. Vihola, Department of Mathematics and Statistics
P.O.Box 35, FI-40014 University of Jyväskylä, FI. Email: 
K. J. H. Law, School of Mathematics,
University of Manchester, Manchester, M13 9PL, UK. Email: 

[2010]65C05 (primary);
60H35, 65C35, 65C40 (secondary)









 
  We develop a Bayesian inference method for diffusions observed discretely and with noise, which is free of discretisation bias. Unlike existing unbiased inference methods, our method does not rely on exact simulation techniques.  Instead, our method uses standard time-discretised approximations of diffusions, such as the Euler-Maruyama scheme.  Our approach is based on particle marginal Metropolis-Hastings, a particle filter, randomised multilevel Monte Carlo, and importance sampling type correction of approximate Markov chain Monte Carlo.  The resulting estimator leads to inference without a bias from the time-discretisation as the number of Markov chain iterations increases.  We give convergence results and recommend allocations for algorithm inputs.  Our method admits a straightforward parallelisation, and can be computationally efficient.  The user-friendly approach is illustrated on three examples, where the underlying diffusion is an Ornstein-Uhlenbeck process, a geometric Brownian motion, and a 2d non-reversible Langevin equation.
   

[
    [
    Received: date / Accepted: date
===================================





§ INTRODUCTION
 

Hidden Markov models (HMMs) are widely used in real applications, for example, for financial and physical systems modeling; see <cit.>. We focus on the case where the hidden Markov chain arises from a diffusion process that is observed with noise at some number of discrete points in time; see e.g. <cit.>.
The parameters associated to the model are static and assigned a prior density. 
Bayesian inference involves expectations with respect to (w.r.t.) the joint posterior distribution of parameters and states, and is important in problems of model calibration and uncertainty quantification. A difficult aspect of Bayesian inference for these models is simulation or evaluation of the diffusion dynamics. Unless the transition probability is explicitly known (see Section 4.4 of <cit.>),  one often resorts to time discretisation, leading to biased inference.  
If one is to seek inference where such discretisation is avoided, there are broadly two schemes (an exception is <cit.>, which we shall discuss below) which one can
follow, the first are the elegant exact simulation methods in <cit.> or the debiasing methods of  
<cit.>, which is the direction followed in this article. We call inference when there is no time-discretisation error, unbiased inference and this is the objective of this article.

Approaches to unbiased inference using the afore-mentioned exact simulation methodology apply for a certain class of diffusions. First, the existence of the Lamperti transformation, after which the process has unit diffusion matrix, and second, the drift of the transformed process has to be of gradient form. Although this includes some important diffusion processes, often, these conditions do not hold for multivariate diffusion processes; this limits the scope of the application of these novel schemes. A more recent and general approach can be found in <cit.>,
which focuses on a type of continuous-time importance sampling method for continuous-time Markov processes, including diffusion processes. The method is, in essence, a type of continuous-time sequential importance sampling algorithm that produces a signed approximation of laws of diffusion processes. Whilst the methodology applies to a reasonably wide class of diffusions, the signed approximation can introduce several algorithmic issues. This includes a large cost associated to the final time of the diffusion process, which at the very least can lead to linear-in-time errors. This latter issue is problematic in our context, as if the methodology is used for HMMs, the time will relate to the number of data and the errors reported in <cit.> will be prohibitive for our application.
It should also be noted that this methodology requires a good proposal based upon a process which is analytically soluble. This limits the applicability of the method, which will not be the case in our context. As the utility of the method has not been explored for the problem of interest in this article, we have considered an alternative method. We proceed with an Euler-Maruyama time-discretisation (see <cit.>), referred henceforth as Euler which is generally applicable and is combined with the debiasing schemes of <cit.>.

Traditional inference approaches based on time-discretisations face a trade-off between bias and computational cost. 
Once the user has decided on a suitably fine discretisation size, one can run, for example, the particle marginal Metropolis-Hastings (PMMH) <cit.>.  This algorithm uses a particle filter (PF) (see <cit.>), where proposals between time points are generated by the approximation scheme, and ultimately accepted or rejected according to a Metropolis-Hastings type acceptance ratio; see <cit.>.  As the discretisation size adopted must be quite fine, a PMMH algorithm can be computationally intensive.

To deal with the computational cost of PMMH, <cit.> develop a PMMH based method
which uses (deterministic) multilevel Monte Carlo (dMLMC) <cit.>. The basic premise of MLMC is to introduce a telescoping sum representation of the posterior expectation associated to the most precise time
discretisation. Then, given an appropriate coupling of posteriors with `consecutive' time discretisations, the cost associated to a target mean square error is reduced,
relative to exact sampling from the most precise (time-discretised) posterior. In the HMM diffusion context, the standard MLMC method (for diffusions without observations) is not applicable, so based upon a PF coupling approach and PMMH, an MLMC method is devised in <cit.>, which achieves fine-level, though biased, inference. 



 §.§ Method

The unbiased and computationally efficient inference method suggested in this paper is built firstly on PMMH, using Euler type discretisations, but using a PMMH targeting a coarse-level model, which is less computationally expensive. This does not yield unbiased inference yet, but it can be achieved by an importance sampling (IS) type correction; see <cit.>.

We suggest an IS type correction that is based on a single-term (randomised) MLMC type estimator <cit.> and the PF coupling approach of <cit.>.
The rMLMC correction is based on randomising the running level in the multilevel context of a certain PF, which we refer to as the `delta PF (ΔPF)' (Algorithm <ref>).  In short, the ΔPF uses the PF coupling introduced in <cit.>, but here an estimator is used for unbiased estimation of the difference of unnormalised integrals corresponding to two consecutive discretisation levels, over the latent states with parameter held fixed (see Section <ref>). 
In <cit.>,
each term in the difference of PMMH averages 
is individually self-normalised (at each level) because of the unknown normalising constants.

The resulting IS type estimator leads to unbiased inference over the joint posterior distribution, and is highly parallelisable, as the more costly (randomised) ΔPF corrections may be performed independently en masse given the PMMH base chain output.  We are also able to suggest optimal choices for algorithm inputs in a straightforward manner (Recommendation <ref> and Figure <ref>).  This is because there is no bias, and therefore the difficult cost-variance-bias trade-off triangle associated with dMLMC is not present.  Besides being unbiased and efficient, our method is user-friendly, as it is a combination of well-known and relatively straightforward components: PMMH, Euler approximations, PF, rMLMC, and an IS type estimator.
For more about the strengths of the method, see Remark <ref> later, as well as <cit.> for more discussion about IS (type) estimators based on approximate Markov chain Monte Carlo (MCMC).

Key to verifying consistency of the method is a finite variance assumption for the rΔPF estimator.  We verify a parameter-uniform bound for the variance under a simple set of HMM diffusion conditions in Section <ref>.
Note, however, that consistency of our method is likely to hold more generally.  This is in contradistinction to methods based on exact simulation, which require analytically tractable transformations to unit covariance diffusion term and computable bounds in the rejection sampler, in order to even apply the method (see for example the review in the recent preprint <cit.>).

We consider a non-reversible  Langevin equation  in Section <ref>, where, to the authors' best knowledge, exact simulation is not applicable.
If an exact simulation method is applicable, the obvious question arises whether our method or the exact simulation method should be applied.  The efficiency of exact simulation type methods is dependent upon several and different factors than our method.  These factors for exact simulation include proper tuning and tight computable bounds for the rejection sampler.  In an ideal scenario for exact simulation, a method based on exact simulation is likely to perform better than our method.  However, in the reverse case, our method can perform better, if the efficiency of exact simulation is poor.  For instance, the efficiency of exact simulation decreases to zero as the analytically computed upper bound of the IS weight used in the rejection sampler  increases to infinity. 

We remark that in principle our algorithm may require simulations with

arbitrarily fine discretisation sizes and corresponding arbitrarily large cost;

this can be infeasible. 
However, it should be noted that the user specifies the chance that this might occur, 
so one can ensure that the probability of 'very expensive' simulations is 

arbitrarily small.
In addition, our method typically has finite expected cost, or cost that is finite with high probability: we give conditions (see Section <ref>), typical elsewhere in the context of HMM diffusions <cit.>, which ensure this (see Section <ref>).








Although we have mostly in mind the case of Euler approximation schemes for the diffusion dynamics approximation, which are generally implementable, other schemes could be possibly be used as well; see <cit.>.  However, suitable couplings for these schemes in dimensions greater than one may not be trivial.  For the sake of theory and proof of consistency, ideally these would have also known weak and strong order convergence rates; see <cit.>.
Indeed, assuming a coupling exists, such higher-order schemes can improve convergence of our method (see Sections <ref> and <ref>).  More generally, our approach based on PMMH or other approximate MCMC, increasingly fine families of approximations, MLMC, and IS correction, could be applied beyond the HMM diffusion context, for example, to HMM jump-diffusions; see <cit.>.    


 §.§ Outline

Section <ref> introduces the aforementioned ΔPF (Algorithm <ref>) and subsequently discusses some applications of randomisation techniques.  The
theoretical properties of the ΔPF in the HMM diffusion context are summarised in Section <ref>. Section <ref> presents the suggested IS type estimator (Algorithm <ref>), based on PMMH with rMLMC (i.e. rΔPF) correction, and details its consistency and a corresponding central limit theorem (CLT). Section <ref> suggests suitable allocations in the ΔPF based on rMLMC efficiency considerations. The numerical experiments in Section <ref> illustrate our method in practice in the setting of an Ornstein-Uhlenbeck process, geometric Brownian motion, and a non-reversible Langevin equation.
Proofs for the technical results of Sections <ref>, <ref> and <ref> are given in Appendix <ref>, <ref> and <ref>, respectively.





 §.§ Notation


Let (, ) be a measurable space.  Functions φ:→ will be assumed measurable.
We denote by 𝒫() the collection of probability measures on (,), and by ℬ_b() the set of φ:→ℝ with 
φsup_x∈ E|φ(x)|<∞.
For a measure μ on (,), we set μ(φ)∫_φ(x)μ( x) whenever well-defined.
For K:×→[0,1] a Markov kernel and μ∈𝒫(), we set
μ K( y) ∫_μ( x) K(x, y),
and 
K(φ)(x) ∫_φ(y) K(x, y),
whenever well-defined.
We use the convention ∏_∅ 1, and p:q{r∈ℤ: p≤ r ≤ q}.  






§ DELTA PARTICLE FILTER FOR UNBIASED ESTIMATION OF LEVEL DIFFERENCES

 


Consider the (Itô) diffusion process

    Z_t     =     a_(Z_t) t + b_(Z_t) W_t,      t≥ 0,

with Z_t∈𝖷ℝ^d, model parameter θ∈, a:𝖷×𝖳→𝖷, b:𝖷×𝖳→𝖷×𝖷, {W_t}_t≥ 0 a d-dimensional Brownian motion, and the initial value Z_0=z_0∈𝖷 a fixed value.
We suppose that there are data 
{Y_p=y_p}_p=0^, y_p ∈^m,
which are observed at equally spaced discrete times, 1:(+1) for simplicity. 
We shall consider the discrete time skeleton of the diffusion (<ref>) Z_1,…,Z_n+1,
where we shall set X_p=Z_p+1, for p=0:.
The Markov transition between X_p-1 and X_p, p=1:n, 
is given by the transition kernel M^(θ,∞)(x_p-1,  x_p) of the diffusion process over unit time, 
with initial
distribution η_0^(θ,∞)(dx_0)=M^(θ,∞)(z_0,dx_0).
It is assumed that conditional on X_p, 
Y_p is independent of random variables {X_i, Y_i}_i≠ p and has density g_(y_p|x_p) G_p^()(x_p).  
Setting M_0^(θ,∞)=η_0^(θ,∞) and M_p^(θ,∞)=M^(θ,∞), p=1:n,
the resulting pair (M_p^(θ,∞), G_p^(θ)) defines the HMM diffusion, and is an example of a so-called Feynman-Kac model (see <cit.>) described below.  As the results of this section can just as easily be stated in terms of Feynman-Kac models, we do so in the following, which shows the generality of our approach.


  In many situations of practical interest, M^(θ,∞)(x_p-1,  x_p), p=1:n exists, but one cannot even simulate from it and/or evaluate a non-negative unbiased estimator of it.
One can often consider the HMM diffusion where one works with a time discretisation of M^(θ,∞)(x_p-1,  x_p). For instance for ℓ∈{0}×ℕ, h_ℓ=2^-ℓ, one
could consider the Euler approximation for t=p,p+h_ℓ,…,p+1-h_ℓ, with Z_p=x_p-1 given

    Z_t+h_ℓ^ℓ = Z_t^ℓ + a_θ(Z_t^ℓ) h_ℓ + b_θ(Z_t^ℓ) _t+h_ℓ^ℓ,

where _t+h_ℓ^ℓi.i.d.∼𝒩(0,h_ℓ) (Gaussian distribution 0 mean, covariance h_ℓ I) and one will set X_p=Z_p+1^ℓ.
The induced transition kernel over unit time is written M^(θ,ℓ)(x_p-1, dx_p). A similar remark can be made for η_0^(θ,∞) with discretisation η_0^(θ,ℓ).





 §.§ Particle filters

A Feynman-Kac model (M_n,G_n) on spaces (_n,_n) arises when


  * _n(_0:n-1, _n) are (regular) probability `transition' kernels from
  _0:-1 to _ for n≥ 1, and M_0( x_0) η_0(_0)∈𝒫(_0), and 
    
  * G_n(_0:) are [0,∞)-valued (measurable) `potential' functions for n≥ 0.

The Particle filter (Algorithm <ref>) (see <cit.>) generates sets of samples
and weights corresponding to the Feynman-Kac model, which for :E_0:→ lead to an unbiased estimator for the (unnormalised) smoother _(G_), defined here in terms of the (unnormalised) predictor 

    _() ∫(x_0:)
        (∏_t=0^-1_t( _0:t)
        )
        η_0(_0)
        ∏_t=1^_t(_0:t-1, _t).

We remark that Step <ref> of Algorithm <ref> refers to the resampling step, which can be multinomial, residual,
stratified or systematic; see e.g. <cit.>.


 

 
     
Suppose that 
:E_0:→ is such that _(G_n )<∞.  Then, the output of Algorithm <ref>
satisfies

    [ ∑_i=1^ V^(i)(^(i)) ]
        = _(_).

 
Proposition <ref> is a restatement of Theorem
7.4.2 of <cit.> in case A_t-1^(i) are sampled independently
(`multinomial resampling').
The extension to the general unbiased case
is straightforward; see
<cit.>.





 §.§ Level difference estimation
 
Suppose that we have two Feynman-Kac models
(M_n^F, G_n^F) and (M_n^C, G_n^C) defined on common spaces (_n,_n). The models correspond to `finer' and `coarser'
Euler type discretised HMM diffusions.
We are interested in estimating (unbiasedly) the difference

    ^F_(_^F)
    -^C_(_^C).

If the models are close to each other, as they will be in the
multilevel (diffusion) context, we would like the estimator
also to be typically small.  In many contexts, if one can estimate the difference using a coupling, it is possible to obtain a variance reduction.  The particular coupling approach we use here is based on using a combined Feynman-Kac model as in <cit.>, which provides a simple, general and effective coupling of PFs, and which we will use to estimate the level difference of unnormalised smoother <ref>. 

Hereafter, we denote _ = (^F_,
^C_)∈ E_× E_, and for _0:n=(_0,…, _n)∈_0^2×…_n^2, we set _0:n^s ( _0^s,…, _n^s)∈_0:n for s∈{F,C}.  

 
     
Suppose that (_t, _t) is a Feynman-Kac
model on the product spaces (E_t× E_t, _t⊗_t), such that:


    
  *  _t is a 
      coupling of the probability measures _t^F and _t^C, i.e. for all A∈_t, we have 
      
    ∫_A×_t _t(_0:t-1, _t) 
    = _t^F( _0:t-1^F, A)
    ,    ∫__t× A_t(_0:t-1, _t) 
    = _t^C( _0:t-1^C, A),

and for A∈_0, we have
η̌_0(A×_0) = η_0^F(A)
and
η̌_0(_0 × A) = η_0^C(A).

  * 
  _t(_0:t) 1/2[_t^F(^F_0:t)
+ G_t^C(^C_0:t)]. 


 

Algorithm <ref> presents a methodology to unbiasedly estimate the level differences (<ref>). In the context of hidden Markov model diffusions, we explain in Remark <ref> how to satisfy Assumption <ref>.


    



 
     
Under Assumption <ref>, the output of Algorithm 
<ref> satisfies

    [ ∑_i=1^2 V^(i)(^(i))]
        =
        _^F(G_^F ) - _^C(G_^C ),

whenever both integrals on the right are well-defined and finite.
 
 
Applying Assumption <ref> one can use the unbiasedness property of PF in Algorithm <ref>, to yield 

    [ ∑_i=1^ V^(i)(^(i))]
           = ∫
        w^F(_0:) (_0:^F)
        (∏_t=0^_t(_0:t))
        η̌( x_0)
        ∏_t=1^_t(_0:t-1, _t)
         
       = ∫(_0:^F)
         (∏_t=0^_t^F(_0:t^F))
         η^F_0( x_0)
         _t^F(_0:t-1^F, _t^F)
        =
        _^F(G_^F),

where, specifically, Assumption <ref><ref>
guarantees _t>0
whenever _t^F>0, and Assumption <ref> implies the
marginal law of ∏_t=0^_t is ∏_t=0^_t^F.  Similarly,
[ ∑_i=+1^2 V^(i)(^(i))] =-_n^C(_^C).
 


Regarding Algorithm <ref>:


    
  * 
   In the hidden Markov model diffusion context, one could consider F to correspond to an Euler discretisation
at level ℓ step-size h_ℓ=2^-ℓ and C  to correspond to an Euler discretisation
at level ℓ-1 step-size h_ℓ-1=2^-(ℓ-1). The couplings 
_p, p=1:n, (of the two Euler transitions over unit time - see Remark <ref>) could be based on using the same underlying Brownian motion; see <cit.>.  That is, for t=p,p+2h^ℓ,…,p+1-2h^ℓ,
(Z_p^ℓ,Z_p^ℓ-1)=(x_p-1^ℓ,x_p-1^ℓ-1) given,

    Z_t+h_ℓ^ℓ   = Z_t^ℓ + a_θ(Z_t^ℓ) h_ℓ + b_θ(Z_t^ℓ) _t+h_ℓ^ℓ
    
      Z_t+2h_ℓ^ℓ   = Z_t+h_ℓ^ℓ + a_θ(Z_t+h_ℓ^ℓ) h_ℓ + b_θ(Z_t+h_ℓ^ℓ) _t+2h_ℓ^ℓ,

with _t+k h_ℓ^ℓi.i.d.∼𝒩(0, h_ℓ), k=1:2, then we can use

    Z_t+h_ℓ-1^ℓ-1 = Z_t^ℓ-1 + a_θ(Z_t^ℓ-1) h_ℓ-1 + b_θ(Z_t^ℓ-1) (_t+h_ℓ^ℓ + _t+2h_ℓ^ℓ),

for the coarser Euler discretisation.  We set (X_p^ℓ,X_p^ℓ-1)=(z_p+1^ℓ,z_p+1^ℓ-1). A similar remark can be made for the initialisation.
The potentials _t^ℓ and _t^ℓ-1 are then simply the conditional likelihood functions _t^ℓ=_t^ℓ-1=G_t^(θ).

  * The choice of _t in Assumption <ref><ref> 
  provides a safe `balance' in between the
approximations, as w^F and w^C are upper bounded by 2^ +1.
Indeed, the coupled Feynman-Kac model can be thought
as an `average' of the two extreme cases–with the choice
_t(x_0:t)=_t^F(^F_0:t) the coupled PF
would coincide marginally with the Feynman-Kac model with dynamics
M_t^F. What is the optimal choice for _t is an interesting question.

  * 
  Clearly, the choice of _0:t can be made also in other
ways.  It is sufficient
for unbiasedness to choose _t(_0:t) such that it 
is strictly positive whenever either the _t^F(_0:t^F) or _t^C(_0:t^C) product is positive, but choices which make w^F and w^C bounded are safer, for instance _0:t(_0:t) =
max{_t^F(_0:t^F), _t^C(_0:t^C)}.  This was the original choice made in <cit.> for approximation of normalised smoother differences. This PF coupling approach based on change of measure and weight corrections w^F and w^C, has been further used also, for example, in <cit.>.

  * Later, in the HMM diffusion context, we set _t^F = _t^C, corresponding to common observational densities, but the method is also of interest with differing potentials.






 §.§ Unbiased latent inference

We show here how the randomisation techniques of <cit.> can be used with the output of Algorithm  <ref> and <ref> to provide an unbiased estimator according to the true model, even though the PFs are only run according to approximate models.

Let us index the transitions ^(ℓ)_p and potentials ^(ℓ)_p by ℓ≥ 0.  They are assumed throughout to be increasingly refined approximations, in the (weak) sense that

    _^(ℓ)(G_^(ℓ))
    ⟶_^(∞)(G_^(∞)),
        as  ℓ→∞,

for all ∈ℬ_b(E_0:), where

    _^(ℓ)() ∫(x_0:)
        (∏_t=0^-1_t^(ℓ)( _0:t)
        )
        η_0^(ℓ)(_0)
        ∏_t=1^_t^(ℓ)(_0:t-1, _t).

In Assumption <ref> we set symbols (F,C) to be (ℓ,ℓ-1) for ℓ≥ 1. We will write the potentials and kernels of the coupled Feynman-Kac model (in the sense
of Assumption <ref>) as (M̌_0:n^(ℓ),Ǧ_0:n^(ℓ)). As a result of <ref>, Algorithm <ref> can provide unbiased estimation of
_^(∞)(_^(∞)), leading to unbiased inference w.r.t. the normalised smoother

    ↦_^(∞)(_) /_^(∞)(_)_n^(∞)(),

which is stated as Proposition <ref> below.

We remark that in step (<ref>) of Algorithm <ref>, in principle, one may have to run Algorithm <ref>
for arbitrarily large L. However, it should be noted that the user specifies the probability 𝐩=(p_ℓ)_ℓ∈, so one can ensure the probability of simulating `very large' values of L is arbitrarily small.

 
  Assumption <ref> holds, 𝐩=(p_ℓ)_ℓ∈ is a probability on ℕℤ_≥ 1 with p_ℓ>0 for all ℓ≥ 1,
  :_0:→ is a function, and 
  
    s_∑_ℓ≥ 0Δ_ℓ^2()/p_ℓ <∞,

where

    Δ_ℓ() ∑_i=1^2 V^(i)(^(i)),

  is formed from the output (V^(1:2), ^(1:2)) of Algorithm <ref> with (^(ℓ)_0:, ^(ℓ)_0:,).
  


  Under Assumption <ref>, the estimator

    ζ()∑_i=1^ V^(i)'(^(i)')
    +
    1/p_LΔ_L(g),

formed from the 
  output of Algorithm <ref> satisfies

    [ζ()]
    =   _^(∞)(G_^(∞)),

    whenever _^(0)(G_) and _^(∞)(G_) are both finite.
  

 Under Assumption <ref>, 
  
  we have (see <cit.>)
  
    [ p_L^-1Δ_L(g)]= _^(∞)(_^(∞)) - _^(0)(_n^(0)),

    so the result follows by Proposition <ref> and linearity of the expectation.
  
The following suggests a fully parallelisable algorithm for unbiased inference over the normalised smoother, and is an unbiased alternative to the particle independent Metropolis-Hastings (PIMH) <cit.> run at some fine level of discretisation.

  Suppose 𝐩 on  satisfies Assumption <ref> for functions ∈{1,}, with
  _^(0)(_n^(0)) and
  _^(∞)(_n^(∞)) finite,
  and
  ^(∞)(_^(∞))>0.
  For each k∈{1:}, if one runs independently Algorithm <ref>, forming ζ_k() from the output as in <ref> for each k, then
  
    E_,,𝐩()
      ∑_k=1^ζ_k()/∑_k=1^ζ_k(1)_n^(∞)()
              almost surely.

Moreover, with - _n^(∞)(),

    √()[E_,,𝐩() - _n^(∞)()]
        𝒩(0,σ^2)
            in distribution,

    where
    
    σ^2=  s_ - ( ^(∞)(_^(∞)) - ^(0)(_^(0)))^2/ [^(∞)(_^(∞))]^2.


The above result follows directly from the results of Section <ref>.  It can also be seen as a multilevel version of Proposition 23 of <cit.>, with straightforward estimators for σ^2. See Section <ref> for suggested choices for 𝐩 and number of particles run at each level.  





§ A VARIANCE BOUND FOR THE DELTA PARTICLE FILTER


In this section we give theoretical results for the ΔPF (Algorithm <ref>) in the setting of HMM diffusions, which can be used to verify finite variance and therefore consistency of related estimators. In particular, Corollary <ref> below can be used to verify Assumption <ref>.


 §.§ Hidden Markov model diffusions

We consider an HMM diffusion and corresponding Feynman-Kac model as in Section <ref>.
We omit θ from the notation in the following, which is allowed as the remaining conditions and results in this Section <ref> will hold uniformly in  (i.e. any constants are independent of ).  The following will be assumed throughout.


The coefficients a^j, b^j,k are twice differentiable for j,k= 1,…, d, and





  (i) uniform ellipticity: b(x)b(x)^T is uniformly positive definite;

  (ii) globally Lipschitz:

there is a C>0 such that 
|a(x)-a(y)|+|b(x)-b(y)| ≤ C |x-y| 
for all x,y ∈ℝ^d; 





Let (x,  y) _p(x, y) for p=0: denote the Markov transition of the unobserved diffusion <ref>, i.e. the distribution of the solution X_1 of <ref> started at X_0=x.  With similar setup from Section <ref>, with _0:n𝖷^n+1, we have that <ref> takes the form

    _n()
    =
    ∫(x_0:)
    (∏_p=0^-1 G_p(x_p))
    η_0( x_0)
    ∏_p=1^(x_p-1,  x_p).

In practice one usually must approximate the true dynamics (x, y) of the underlying diffusion with a simpler transition M^(ℓ)(x, y), based on some Euler type scheme using a discretisation parameter h_ℓ= 2^-ℓ for ℓ≥ 0; see <cit.>.
The scheme allows for a coupling of the diffusions (X_t^(ℓ), X_t^(ℓ-1))_t≥ 0 running at discretisation levels ℓ and ℓ-1 (based on using the same Brownian path W_t),
such that for some β∈{1,2}, we have

    _(x,y)[|X_1^(ℓ)-X_1^(ℓ-1)|^2] ≤ M(|x-y|^2 + h_ℓ^β),

where M<∞ does not depend on ℓ≥ 1.
In particular, if the diffusion coefficient b(X_t) in <ref> is constant or if a Milstein scheme can be applied otherwise, then β=2; otherwise β=1; see Proposition D.1 of <cit.>.



 §.§ Variance bound

Assume we are in the above HMM diffusion setting, and that the coupling of Assumption <ref> holds,
with symbols (F,C) equal to (ℓ,ℓ-1) for ℓ≥ 1, and G_p^(ℓ)=G_p^(ℓ-1) G_p for p=0:.
Running Algorithm <ref>, we recall that Δ_ℓ(), defined in <ref>, satisfies, by Proposition <ref>,

    [Δ_ℓ()]
    =
    _n^(ℓ)(G_n ) - _n^(ℓ-1)(G_n ),

regardless of the number ≥ 1 of particles.

A (measurable) function :𝖷→ is Lipschitz, denoted ∈Lip(𝖷), if for some C'<∞, |(x) - (y)| ≤ C' |x-y| for all x,y∈𝖷.    
 

  The following conditions hold for the model (_n,_n):
    
  


  (i)G_n<∞ for each n≥ 0.

  (ii)G_n∈Lip(𝖷) for each n≥ 0.

  (iii)inf_x∈𝖷G_n(x)>0 for each n≥ 0.


  
    For every n≥ 1, φ∈Lip(𝖷)∩ℬ_b(𝖷)
    there exist a C'<∞ such that for s∈{F,C}, we have for every (x,y)∈𝖷×𝖷 that
    |M_n^s(φ)(x)-M_n^s(φ)(y)|≤ C'|x-y|.



In the following results for Δ_ℓ(), the constant M<∞ may change
from line-to-line. It will not depend upon  or ℓ (or θ), but may depend on the time-horizon  or the function . 
𝔼 denotes expectation w.r.t. the law associated to the ΔPF started at (,), with ∈𝖷.
Below we only consider multinomial resampling in the ΔPF for simplicity, though Theorem <ref> and Corollary <ref> can be proved also assuming other resampling schemes. 


Assume (A<ref>-<ref>).  Then for any φ∈ℬ_b(𝖷^+1)∩Lip(𝖷^+1), 
there exists a M<∞ such that

    𝔼[(Δ_ℓ(φ) - 𝔼[Δ_ℓ(φ)])^2]
    ≤Mh_ℓ^2∧β/,
        with β as in <ref>.




  Assume (A<ref>-<ref>).  Then for any φ∈ℬ_b(𝖷^+1)∩Lip(𝖷^+1),
there exists a M<∞ such that

    𝔼[(Δ_ℓ(φ))^2]
    ≤ M(h_ℓ^2∧β/ +  h_ℓ^2 ),
            with β as in <ref>.


The proofs are given in Appendix <ref>.

Based on Corollary <ref>, Recommendation <ref> of Section <ref> suggests allocations for 𝐩 and _ℓ in the ΔPF (Algorithm <ref>) to optimally use resources and minimise variance <ref>.

  



§ UNBIASED JOINT INFERENCE FOR HIDDEN MARKOV MODEL DIFFUSIONS

 

We are interested in unbiased inference for the Bayesian model posterior

    π^(∞)(θ,  x_0:)
        ∝pr(θ) G_^(θ)(x_) γ_^(θ,∞) ( x_0:),

    where pr(θ) = pr(θ) θ is the prior on the model parameters, and 
    
    γ_^(θ,∞)( x_0:) =
        (∏_t=0^-1 G_t^(θ)(x_t)
        )
        η_0^(θ)( x_0)
        ∏_t=1^
        M_t^(θ,∞)(x_t-1, x_t).
    
Here, M^(θ,∞)_t
corresponds to the transition density of the
diffusion model of interest.
The dependence of the HMM on θ is made explicit in this section.  As in Section <ref>, we assume the transition densities 
M_t^(θ,∞) cannot be simulated, 
but that there are increasingly refined discretisations M_t^(θ,ℓ) approximating M_t^(θ,∞) in the sense of <ref> (with E_0:𝖷^+1).  



 §.§ Randomised MLMC IS type estimator based on coarse-model PMMH

We now consider Algorithm <ref> for joint inference w.r.t. the above Bayesian posterior.  Algorithm <ref> uses the following ingredients:

    
  * ^(θ,ℓ)_0: satisfying Assumption
<ref><ref> with M_0:^F = M_0:^(θ,ℓ), 
and M_0:^C = M_0:^(θ,ℓ-1).

  * Ǧ^(θ)_0: defined as in Assumption 
<ref><ref>, with G_0:^F = G_0:^C =
G_0:^(θ).

  * Metropolis-Hastings proposal distribution q(|θ) 
  for parameters.

  * Algorithm constant ϵ≥ 0 (e.g. ϵ=10^-10; see Remark <ref><ref> below).

  * Number of MCMC iterations ∈ and number of particles ∈.

  * Probability mass 𝐩 = (p_ℓ)_ℓ∈ on  with p_ℓ>0 for all
  ℓ∈.




 




Before stating consistency and central limit theorems, we briefly
discuss various aspects of this approach, which are appealing from
a practical perspective, and we also mention certain algorithmic
modifications which could be further considered.
  


  * 
  The first phase (P<ref>) of Algorithm 
<ref> implements a PMMH type algorithm
<cit.>. If ϵ=0, this is exactly PMMH
targeting the model
π^(0)(θ,  x_0:)∝pr(θ) G_^(θ)(x_) γ_^(θ,0)( x_0:). It is generally safer to choose
ϵ>0 <cit.>, which ensures that the IS type
correction 
in phase (P<ref>) will yield consistent inference for the ideal model 

    π^(∞)(θ,  x_0:)∝pr(θ) G_^(θ)(x_)
    γ_^(θ,∞)(x_0:)
 (Theorem <ref>). Setting ϵ>0 may be helpful otherwise in terms of improved mixing, as the PMMH will target marginally an averaged probability between a `flat' prior and a `multimodal' ℓ=0 marginal posterior.

  * It is only necessary to implement PMMH for the coarsest level.
  This is typically relatively cheap, and therefore allows for 
  a relatively long MCMC run. Consequently, relative cost of burn-in is
  small, and if the proposal q is adapted (see <cit.>), it has time to converge.

  * The (potentially costly) rΔPFs 
  are applied independently for each Θ_k, which allows for 
  efficient parallelisation.

  * We suggest that the number of particles N, here referred to as `_0', used in the PMMH be chosen based on <cit.>, while the number of particles `N_ℓ' (and p_ℓ) can be optimised for each level ℓ based on Recommendation <ref> of Section <ref>, or kept constant.  One can also afford to increase the number of particles when a `jump chain' representation is used (see the following remark).

  * The rΔPF corrections may be calculated only once for each
  accepted state <cit.>. That is, 
  suppose (Θ̃_k, Ṽ_k^(1:), 
  𝐗̃_k^(1:))_k=1^ are the accepted
  states, (D_k)_k=1^ are the corresponding holding
  times, and (Ṽ_k,L_k^(1:2_L_k), 
  𝐗̃_k,L_k^(1:2_L_k))_k=1^ are corresponding ΔPF outputs, then the estimator is formed as in Algorithm <ref> using (Θ̃_k, Ṽ_k^(1:), 
  𝐗̃_k^(1:)), and accounting for the holding times
  in the weights defined as
  W_k^(i)
   D_k Ṽ_k^(i) / (∑_j=1^Ṽ_k^(j)+ϵ)
   and W_k,L_k^(i)Ṽ_k,L_k^(i)/[p_L_k(∑_j=1^Ṽ_k^(j)+ϵ)].

  * In case the Markov chain in (P<ref>) phase is slow
  mixing, (further) thinning may be applied (to the jump chain) before the (P<ref>) phase.

  * In practice, Algorithm <ref> may be implemented in
  an on-line fashion w.r.t. the number of
  iterations , and by progressively refining the estimator
  E_,,𝐩(). The rΔPF corrections
  may be calculated in parallel with the Markov chain.

  * In Algorithm <ref>, the rΔPFs
  depend only on Θ_k. They could depend also on 
  V_k^(i) and 𝐗_k^(i), but it is not clear how such
  dependence could be used in practice to achieve better performance.
  Likewise, the `zeroth level' estimate in Algorithm <ref>
  is based solely on particles in (P<ref>), but it
  could also be based on (additional) new particle filter output.

  * In order to save memory, it is possible also to `subsample' only one 
  trajectory 𝐗_k^* from 𝐗_k^(1:), such that
  [𝐗_k^* = 𝐗_k^(i)]∝ V_k^(i), and set W_k,0^* ∑_i=1^ W_k,0^(i), and similarly in Algorithm
  <ref> find 𝐗̌^* such that
  [𝐗̌^* = 𝐗̌^(i)]∝V̌^(i), setting 𝐗_k,L_k^*(1:2)𝐗̌^*, and defining from the usual output of Algorithm <ref>, W^*(1)_k,L_k∑_i=1^ W_k,L_k^(i) and W^*(2)_k,L_k∑_i=+1^2 W_k,L_k^(i).  The subsampling output estimator then takes the form,
  
    E_,,𝐩^subsample()
          ∑_k=1^[W_k,0^* (Θ_k, 𝐗_k^*) +
            ∑_i=1^2 W_k,L_k^*(i)(Θ_k,
            𝐗_k,L_k^*(i)) ]
            /∑_k=1^[W_k,0^* + 
             ∑_i=1^2 W_k,L_k^*(i)] 
            .

  Note, however, that the asymptotic variance of this estimator is higher,  because 

  E_,,𝐩() may be viewed as a Rao-Blackwellised version
  of E_,,𝐩^subsample().








 §.§ Consistency and central limit theorem




 
     
Assume that the algorithm constant ϵ≥ 0 is chosen positive, and that the Markov chain 
(Θ_k,
X_k^(1:), V_k^(1:)))_k≥ 1 is ψ-irreducible, and that
π^(0)() and π^(∞)() are finite.
For each θ∈, suppose Assumption
<ref> holds for  
≡ 1 and = ^(θ)(θ,), with
_0:^(ℓ)_0:^(θ,ℓ) and _0:^(ℓ)_0:^(θ).
Assume

    ∫pr(θ) (√(s_1(θ)) +
        √(s_^(θ)(θ)))θ  <∞.

Then, the estimator of Algorithm
<ref> is strongly consistent:

    E_,,𝐩() ∫π^(∞)(θ,  x_0:) (θ, x_0:)
            (a.s.)

 



 
Regarding Theorem <ref>, whose proof is given in Appendix <ref>:
  
    


  * If all potentials G_t are strictly positive, the
  algorithm constant ϵ may be taken to be zero. However, if
  ϵ=0 and Algorithm
  <ref> with 
  (M_0:^(Θ̂_k,0), G_0:^(Θ̂_k),)  
  can produce an estimate with ∑_i=1^ V^(i)=0 with
  positive probability, the consistency may be lost
  <cit.>.


 





 
  Suppose that the conditions of Theorem <ref> hold.  Suppose additionally that π^(∞)(^2) <∞ and that the base chain (Θ_k, V^(1:)_k,𝐗^(1:)_k)_k≥ 1 is aperiodic, with transition probability denoted by P.  Then,

    √()[ E_,,𝐩() - π^(∞)() ]
        𝒩(0,σ^2),
            in distribution,

whenever the asymptotic variance 

    σ^2 = (P, μ_) +
          (σ^2_ξ)/c^2,

is finite.  Here, - π^(∞)(), c>0 is a constant (equal to (μ_1)), and 

    σ^2_ξ(θ, v^(1:), 𝐱^(1:))    ( ξ_k()Θ_k=θ,V_k^(1:)=v^(1:),
        𝐗_k^(1:)=𝐱^(1:)) 
        = s_^(θ)(θ) -
          (_^(θ,∞)(G_n ^(θ) ) - 
    _^(θ,0)(G_n ^(θ))  )^2
          /(∑_i=1^ v^(i) + ϵ)^2.

 


Proposition <ref> follows from Theorem 7 <cit.>.
  We suggest that =_0 for (P<ref>) be chosen based on
<cit.> to minimise (P, μ_), and that (p_ℓ) and =_ℓ in (P<ref>) for the rΔPF be chosen as in Recommendation <ref> of Section <ref>, to minimise σ_ξ^2, subject to cost constraints, in order to jointly minimise σ^2.  However, the question of the optimal choice for N_0 in the IS context is not yet settled.

  
 Regarding ψ-irreducibility and aperiodicity in Theorem <ref> and Proposition <ref>, these are inherited by the coarse PMMH chain <cit.> if the corresponding idealised marginal Metropolis-Hastings for the coarse model has these properties.  





§ ASYMPTOTIC EFFICIENCY AND RANDOMISED MULTILEVEL CONSIDERATIONS
 

We summarise the results of this section by suggesting the following safe allocations for probability 𝐩=(p_ℓ)_ℓ∈ and number =_ℓ of particles at level ℓ≥ 1 in the ΔPF (Algorithm <ref>) used in Algorithm <ref> and Algorithm <ref>, and Proposition <ref>, with β given in <ref> in the HMM diffusion context of Section <ref>, or, indeed, with β given in the abstract framework under Assumption <ref> given later. See also Figure <ref> for the recommended allocations.

With strong error convergence rate  β given in <ref>, we suggest the following for 𝐩=(p_ℓ)_ℓ∈ and _ℓ∈ in ΔPF (Algorithm <ref>):
  
  
  
  * (e.g. Euler scheme).
Choose p_ℓ = (1/2)^ℓ and N_ℓ∝ 1 constant.
  
  * (e.g. Milstein scheme).
    Choose p_ℓ∝ 2^-1.5 ℓ≈ (1/3)^ℓ and _ℓ∝ 1 constant. 
  
  

The suggestions are based on Corollary <ref> of Section <ref>, and 
Proposition <ref> (β=2) and <ref> (β=1) given below (with weak convergence rate α=1; see Figure <ref> for general α).  In the Euler case, although the theory below gives the same computational complexity order by choosing any ρ∈ [0,1] and setting p_ℓ∝ 2^-(1+ρ)ℓ and N_ℓ∝ 2^ρℓ, the experiment in Section <ref> gave a better result using simply ρ=0, corresponding to no scaling.  However, this may depend on the implementation.





 §.§ Efficiency framework

The asymptotic efficiency of simulation-based estimators has been considered theoretically in <cit.>; see <cit.> in the dMLMC context.  The developments of this section follow <cit.> for rMLMC (originally in the i.i.d. setting without observations), while also giving some extensions (also applicable to that setting).  We will see that the basic rMLMC results carry over to our setting involving MCMC and randomised estimators based on PF outputs, but also discover a novelty in the common Euler case (β=1 in Figure <ref>).  Proofs are given in Appendix <ref>.

We are interested in modeling the computational costs involved in running Algorithm <ref>; the algorithm of Proposition <ref> is recovered with ={θ}.
Let τ_Θ_k,L_k represent the combined cost at iteration k of the base Markov chain and weight calculation in Algorithm <ref>, so that the
total cost () of Algorithm <ref> with  iterations is
  
    ()∑_k=1^τ_Θ_k,L_k.

The following assumption seems natural in our setting.


  For Θ_k∈𝖳, a family {τ_Θ_k,ℓ}_k, ℓ≥ 1 consists of positive-valued random variables that are independent of {L_k}_k≥ 1, where L_k∼𝐩 i.i.d., and that are conditionally independent given {Θ_k}_k≥ 1, such that τ_Θ_k,ℓ depends only on Θ_k∈𝖳 and ℓ∈ℕ.  

  Under a budget constraint κ>0, the realised length of the chain is (κ) iterations, where
  
    (κ)
      max{≥ 1: () ≤κ}.

Under a budget constraint, the CLT of Proposition <ref> takes the following altered form, where here _m(θ) denotes the θ-marginal of the invariant probability measure (given as <ref> in Appendix <ref>)  of the base Markov chain (equal to the θ-marginal posterior of the ℓ=0 model).

  If the assumptions of Proposition <ref> hold with σ^2<∞,
 and if __m⊗𝐩[]<∞ with (θ,ℓ)[τ_Θ_k,L_k| Θ_k =θ, L_k =ℓ],
  then
  
    √(κ)[ E_(κ),,𝐩() - π^(∞)()]
        𝒩(0,[]σ^2),
            in distribution.

  
 The quantity []σ^2 is called the `inverse relative efficiency' by <cit.>, and is considered a more accurate quantity than the asymptotic variance (σ^2 here) for comparison of Monte Carlo algorithms run on the same computer, as it takes into account also the average computational time.  
  
In the following we consider (possibly) variance reduced (if >0) versions of Δ_ℓ() of Assumption <ref>, denoted Δ_ℓ, where =^(θ), based on running the ΔPF (Algorithm <ref>) with parameters θ, ℓ fixed.  The constant C<∞ may change line-to-line, but does not depend on , ℓ, or θ, but may depend on the time-horizon  and function .          
  
Assumption <ref> holds, and constants 2 α≥β>0, >0, and ≥ 0 are such that the following hold:
  
    
    
  * (Mean cost)
[ τ_θ,ℓ] ≤ C 2^ℓ(1+) 
    
  * (Strong order)
      [Δ_ℓ^2] ≤ C 2^-ℓ(β + )  + C 2^-2 αℓ
    
  * (Weak order)
|Δ_ℓ| ≤ C 2^-αℓ
    
    
  
   Regarding Assumption <ref>: 
    
    
    
  * We only assume bounded mean cost in Assumption <ref>, rather than the almost sure cost bound commonly used.  This generalisation allows for the setting where occasional algorithmic runs may take a long time.
    
  * In the original MLMC setting, the cost scaling  in Assumption <ref> is taken to be =1 <cit.>.  However, in settings involving uncertainty quantification, and where the forward solver may involve non-sparse matrix inversions, often ≥ 1 <cit.>.
      
  * We assume in Assumption <ref> that the mean cost to form Δ_ℓ is bounded by the -scaled product of the number of samples or particles _ℓ times the number of Euler time steps 2^ℓ + 2^ℓ-1 together with the O(_ℓ)-resampling cost, where there are _ℓ∝ 2^ℓ particles at level ℓ.  Here, we recall that the stratified, systematic, and residual resampling algorithms have O(_ℓ) cost, as does an improved implementation of multinomial resampling; see <cit.>.
      
  * With =0, by Jensen's inequality one sees why α≥β/2 can be assumed, and that Assumption <ref> becomes Δ_ℓ^2 ≤ C 2^-ℓβ.  
    
  * ≥ 0 in Assumption <ref> and <ref> corresponds to using an average of
      _ℓ⌈ 2^ℓ⌉
      i.i.d samples of Δ_ℓ^(1), i.e. Δ_ℓ = 1/_ℓ∑_i=1^_ℓΔ_ℓ^(i), or, of more present interest to us, to increasing the number of particles used in a PF by a factor of _ℓ instead of the default lower number.  The former leads to 
      Δ_ℓ^2 = 1/_ℓ(Δ_ℓ^(1)) + [Δ_ℓ^(1)]^2,
justifying Assumption <ref>, as does Corollary <ref>, with β∈{1,2} and α=1, for the ΔPF (Algorithm <ref>) in the HMM diffusion context (Section <ref>).    
      
      
    
  

  Suppose Assumption <ref> and the assumptions of Proposition <ref> hold, with (P, μ_)<∞.  If p_ℓ∝ 2^-r ℓ for some r∈((1+), min(β + , 2α)), then <ref> holds, i.e.
    
    √(κ)[ E_(κ),,𝐩() - π^(∞)()]
        𝒩(0,σ^2),
            in distribution.


 Regarding Proposition <ref>, in the common case =1 for simplicity:
  
  
    
  * If β>1 (`canonical convergence regime') and =0, then a choice for r∈ (1,β) exists.  See also Theorem 4 of <cit.> for a discussion of the theoretically optimal 𝐩. 

  * If β≤ 1 (`subcanonical convergence regime'), then β + ≤ 1 + and so no choice for r exists.
  
  
  



 §.§ Subcanonical convergence

When β>1, within the framework above we have seen that a canonical convergence rate holds (Proposition <ref>) because <∞ and σ^2<∞.  When β≤ 1, this is no longer the case, and one must choose between a finite asymptotic variance and infinite expected cost, or vice versa.  Assuming the former, and that a CLT holds (Proposition <ref>), for ϵ>0 and 0<δ <1 the Chebyshev inequality implies that the number of iterations of Algorithm <ref> so that
  
    [  E_,,𝐩() -π^(∞)()≤ϵ]≥ 1-δ,

  holds implies that  must be of the order O(ϵ^-2).  The question is then how to minimise the total cost 𝒞(), or computational complexity, involved in obtaining the  samples.  This will involve optimising for (p_ℓ) and _ℓ to minimise 𝒞(), while keeping the asymptotic variance finite.

  
    Suppose that the assumptions of Proposition <ref> hold with σ^2<∞, and Assumption <ref> holds with
    [τ_Θ_k_0,L_k_0]
    =∞
    for some k_0≥ 1.  If
    
    ∑_k≥ 1sup_j≥ 1[τ_Θ_j,L_j > a_k]
        <
        ∞,

    with
    a_k = O( k^c_1 (log_2 k)^c_0)
    for some constants c_0 >0 and c_1≥ 1, then <ref> can be obtained with computational complexity
    
    O( ϵ^-2 c_1log_2 ϵ^c_0)
              as  ϵ→ 0.

  
  The above result shows that even for costs with unbounded tails, reasonable confidence intervals and complexity order may be possible.  This may be the case for example when a rejection sampler or adaptive resampling mechanism is used within Algorithm <ref> or <ref>, which may lead to large costs for some Θ_k, for example a cost with a geometric tail.
  
  
  The next results are as in Proposition 4 and 5 of <cit.> in the standard rMLMC setting, and shows how one can choose 𝐩, assuming an additional almost sure cost bound,  so that σ^2<∞, with reasonable complexity.  
  
    
      Suppose that the assumptions of Proposition <ref> hold with (P, μ_f̅)<∞, and that Assumption <ref> holds with β≤ 1, where moreover τ_θ,ℓ≤ C 2^γℓ(1+) almost surely, uniformly in Θ_k=θ∈𝖳.  For all q>2 and η>1, the choice of probability
      
    p_ℓ∝
          2^-2 b ℓℓ [ log_2 (ℓ +1)]^η,

      where b min((β + )/2,α), leads to σ^2 <∞, and <ref> can be obtained with computational complexity
      
    O(
          ϵ^-(1+)/blog_2 ϵ^q (1+)/2b)
              as  ϵ→ 0.

      
     Regarding Proposition <ref>, with =1:
      
      

  * Under Assumption <ref> with =0, the usual setup in MLMC before variance reduced estimators are used, the above proposition shows that finite variance and <ref> can be obtained without increasing the number of particles at the higher levels, even in the subcanonical regime.  We have in this case b=β/2≤α and complexity
O(
      ϵ^-2/βlog_2 ϵ^q/β).
      When β=1 (borderline case), dMLMC gives complexity O(ϵ^-2log_2 ϵ^2) <cit.>, which is negligibly better (recall q>2), but is biased inference.
     
  * When α>β/2, which is the usual case in the subcanonical regime (β≤ 1); see <cit.>, a more efficient use of resources can be obtained by increasing the number of particles (see Proposition <ref> below).


      
    
Suppose the assumptions of Proposition <ref> hold, where moreover ≥ 0 may vary as a free parameter without changing the constant C>0.  Then,
      for all q>2, η>1 constants, the choice = 2α-β and probability
      
    p_ℓ∝ 2^-2αℓℓ [log_2(ℓ +1)]^η,

leads to σ^2<∞, and <ref> can be obtained with computational complexity

    O( ϵ^[-2 - (1-β)/α ]log_2 ϵ^[q + (1-β)/2α ])
              as  ϵ→ 0.

    




    


    

§ NUMERICAL SIMULATIONS


Now the theoretical results relating to the method herein introduced 
will be demonstrated on three examples.  
We will consider one example in the canonical regime, and two in the sub-canonical.  
In the first two experiments, the likelihoods can be computed exactly, 
so that the ground truth π^(∞)() can be easily calculated to arbitrary precision.  We run each example with 100 independent replications, and calculate the mean squared error (MSE) when the chain is at length  as

    MSE()
    =
    1/100∑_i=1^100| E_,,𝐩^(i)() - π^(∞)()|^2,

which is depicted as the thick red line, average of the thin lines, in Figure <ref> below.
The error decays with the optimal rate of cost^-1 and log(cost)cost^-1 in the canonical and sub-canonical cases, respectively, 
where cost is the realised cost of the run, () from Section <ref>, measured in seconds, with  iterations of the Markov chain.

Three examples will be considered. 
First we consider two models where the exact marginal likelihood can be computed. 
This way a reliable ground truth can be computed with a long MCMC chain, 
providing a strong verification of the theoretical results.
In Section <ref> the Ornstein–Uhlenbeck (OU) process is considered,
where Euler-Maruyama provides the canonical convergence regime.
In Section <ref> the Geometric Brownian motion is considered,
where Euler-Maruyama provides sub-canonical convergence regime.
In Section <ref> we consider a more complicated 2d model
which does not allow exact computation of the marginal likelihood.
  
It is of interest to compare our methodology to existing unbiased methods. 
 The method we consider for comparison is PMMH using 
 the exact method introduced by Fearnhead et al. <cit.>. 
 In their work they provide a way to construct unbiased estimates without approximating the transition density.
 
 The key idea is to assign to each particle a random positive weight which is an unbiased estimator of the true weight. 
 This method is referred to as the random weight particle filter.
 
 It has been later extended to continuous-time observations
 
 in <cit.>.
The method of <cit.> is implemented when it is applicable 
(in particular, for the models of Sections <ref> and <ref>)
and the corresponding MSE is plotted in comparison to 
our method. 
That method is not amenable to the example of Section <ref>.
In all the examples, for the sake of comparison with the standard approach, we also implement a finely-discretised PMMH.  This shows the benefit of our alternative approach based on a coarsely-discretised PMMH with  multilevel IS correction.



 §.§ Ornstein–Uhlenbeck process

Consider the OU  process

    Z_t = -a Z_t  t + b  W_t   ,
         t≥ 0,

with initial condition Z_0=0, model parameter θ = (θ_1,θ_2) ∼ N(0,σ^2 I), and a a_θ=exp(θ_1) and b b_θ=exp(θ_2).
The process is discretely observed for k=1,…,,

    Y_p = X_p + ξ_p   ,

where ξ_p ∼ N(0,γ^2) i.i.d. and recall that X_p=Z_p+1.  
Therefore, 

    G_p(x) = exp(-1/2γ^2 |x - y_p|^2 )   .


The marginal likelihood is given by 

    [y_1:|θ] = ∏_p=1^[y_p|y_1:p-1,θ]   ,

and each factor can be computed as the marginal of the joint on the prediction and current observation, i.e.

    [y_p|y_1:p-1,θ] = ∫_[y_p|x_p,θ] [x_p | y_1:p-1,θ]  x_p   .



In this example the ground truth can be computed exactly via the Kalman filter.  
In particular, the solution of <ref> is given by 

    Z_1 = e^-aX_0 + W_1   ,    W_1 ∼𝒩 (0, b^2/2a(1 - e^-2a)  )   .

The filter at time p is given by the following simple recursion

    m_p = c_p  (y_p/γ^2 + m̂_p/ĉ_p ), c_p = (γ^-2 + ĉ_p^-1)^-1  ,   m̂_p = e^-a m_p-1  ,   ĉ_p = e^-2a c_p-1 + b^2/2a(1 - e^-2a)   .
 
Additionally, the incremental marginal likelihoods <ref> can be computed exactly

    [y_p|y_1:p-1,θ] = √(c_p/2πĉ_pγ^2)exp{-1/2 [ y_p^2/γ^2 + m̂_p^2/ĉ_p - c_p  ( y_p/γ^2 + m̂_p/ĉ_p)^2 ] }  .


The parameters are chosen as γ=1, σ^2=0.1, =5, and 
the data is generated with θ= (0,0)^T.
Our aim is to compute (θ| y_1:) (or [ (a,b)^T | y_1: ], etc., but we will content ourselves with the former).
This is done via a brute force random walk MCMC for =10^8 steps using the exact likelihood [y_1:|θ] as above.  
The IACT is around 10, so this gives a healthy limit for MSE computations.

For the numerical experiment, we use Euler-Maruyama method at resolution h_ℓ=2^-ℓ
to solve <ref> as follows

    Z_p+1 = (1-a h_ℓ) Z_p + b B_p+1  ,    B_p+1∼𝒩(0,h_ℓ)    i.i.d.

for p=1,…, K_ℓ=h^-1_ℓ. Levels ℓ and ℓ-1 are coupled in the simulation of Δ_ℓ
by defining B^C_1:K_ℓ/2 = B^F_1:2:K_ℓ-1 + B^F_2:2:K_ℓ
Algorithm <ref> is then run using the standard bootstrap particle filter (Algorithm <ref>) with =20 particles and O()-complexity multinomial resampling; see <cit.>.  Theorem <ref> provides a rate of β=2 for Algorithm <ref>, because the diffusion coefficient is constant, which implies we are essentially running a Milstein scheme (see <ref> and <cit.>).  Recommendation <ref> (or Proposition <ref>) of Section <ref> suggests arbitrary precision can be obtained by Algorithm <ref> with p_ℓ∝ 2^-3ℓ/2 and no scaling of particle numbers based on ℓ in this canonical β=2 regime (with weak rate α=1).  We choose a positive PMMH algorithm constant ϵ=10^-6 (see Remark <ref><ref>).  We run Algorithm <ref> for 10^4 steps, 
with 100 replications. For the finely-discretised PMMH experiment we run 10^4 steps, with 100 replications with a discretisation of h_ℓ = 2^-5.  The results are presented in Figure <ref>, where it is clear that the theory holds and the MSE decays according to 1/cost.
The variance of the run-times is very small over replications.
The method of <cit.>, within PMMH,
also converges with the theoretically-predicted canonical rate, 
but with a slightly smaller constant. 
This is not unexpected. 
The important point is that both methods achieve the same canonical rate, 
while our method is quite generally applicable, 
in particular, to a wide range of models inaccessible to methods of the type of <cit.>.
Also with the finely discretised PMMH from Figure <ref>, i.e. the curve titled as MSE3, we see the bias kicking in at the end. We also see that for the same level of cost the MSE is higher than that of the other methodologies. 








 §.§ Geometric Brownian motion

We next consider the following stochastic differential equation 

    Z_t = 
    
    
     a Z_t W_t,

with initial condition Z_0 = 1, and a a_θ=exp(θ) with θ∼𝒩(0,σ^2).  

This equation is analytically tractable as well, and the solution
of the transformed equation U = log Z is given via Itô's formula by

    U_t = -a^2/2 t + a  W_t.

Defining W_p ∼𝒩(0,1) i.i.d., one has that 

    U_p+1 = U_p +  -a^2/2 + a  W_p   ,    with    U_0 = log z_0 = 0,

and the solution of <ref> can be obtained via exponentiation: Z_p = e^U_p.
Moreover, noisy observations 
are introduced on the form , with X_p=Z_p+1,







    Y_p = log(X_p) + ξ_p,

where ξ_p ∼𝒩(0,γ^2) i.i.d. as above.
Therefore we have

    G_p(x) = exp(-1/2γ^2 |log(x) - y_p|^2 ).

Again [y_1:|θ] can be computed analytically. 
The parameters γ=1, σ^2=0.1, =5 are chosen 
the same as in the previous example and the true observations are generated again with θ=0.


In order to investigate the theoretical sub-canonical rate, we return to <ref>
and approximate this directly using Euler-Maruyama method <ref>, which introduces
artificial approximation error.  This problem suffers from stability problems when X<0, 
so we take h_ℓ=2^-6-ℓ.
 Algorithm <ref> is then used along with the selection functions (<ref>).



Here the diffusion coefficient is not constant, and Euler-Maruyama method yields a rate of β=1=α,
the borderline case, which is expected to give a logarithmic penalty.
Based on Recommendation <ref> (or Proposition <ref>) of Section <ref>, 

we consider scaling the particles as 2^ℓ with =2α-β=1 and =0, 
with p_ℓ∝ 2^-2ℓℓlog(ℓ)^2 in both cases.  Again we let ϵ=10^-6, and 
 the standard bootstrap particle filter is used, with =20× 2^ℓ particles.  
Algorithm <ref> is run for 10^4 steps, 
with 100 replications. Again, for the finely-discretised PMMH, we run 10^4 steps, with 100 replications with a discretisation of h_ℓ = 2^-4 
For this sub-canonical case we impose an artificial upper bound ℓ≤ 10, 
corresponding to an induced bias of ≈ 10^-5.
The results are presented in Figure <ref>, and they 
show good agreement with the theory, in terms of rate. 
On the other hand, the cost for =0 is apparently smaller than that of 
=1 by a factor of approximately 100. 
The method of <cit.> is not expected to suffer from 
a logarithmic penalty on the MSE convergence, i.e. it achieves canonical rate also in this example.
This can be seen in Figure <ref>, in addition to a slightly better constant, as before.
For the finely-discretised PMMH from Figure <ref>, with geometric Brownian motion, we again notice the effect of the bias arising from the discretisation, and the overall higher MSE.

    

 §.§ 2d Non-reversible  Langevin equation

   
    We now consider a 2d example which is not amenable to approaches of the type <cit.>.
    Consider a target distribution of the type ρ(z) ∝exp(-Φ(z)/a_2), a_2>0, 
    and the following non-reversible Langevin equation
    
    Z_t =
        (A - I_2) ∇Φ(Z_t)  t + √(2 a_2) W_t,
             t≥ 0   ,

and noisy observations Y_p ∼𝒩(X_p, γ^2 I_2), with γ=1, n=10 observations where A ∈^2× 2
is anti-symmetric and parameterised by a_1 ≥ 0, I_2 is the 2-dimensional identity matrix,
and Φ(z) = a_3/2(z_1^2+z_2^2-1)^2 is the ring potential, parameterised by a_3 >0.
The initial condition is specified as Z_0=[1,1]^T. 
It is easy to see that the right-hand side of the Fokker-Planck equation vanishes for the invariant distribution
ρ given above, so the dynamics are well-behaved <cit.>. We let a_i=exp(θ_i), for i=1,…, 3, and the prior is given by 
θ∼ N(0,ε^2 I_3), ε^2=0.1. The resolution of the Euler–Maruyama scheme for this experiment is set as 2^-ℓ+1.
This problem is no longer analytically soluble, so a high-resolution simulation is used with a large sample size 
as ground truth.





For our setup it follows similarly to that of the previous experiment for the OU process, where we are working in a canonical regime. Again we run Algorithm <ref> with 10^4 steps, and 100 replications, 
and Algorithm <ref> is run using the standard bootstrap particle filter (Algorithm <ref>) with =20 particles. We compare our results to the single level PMMH, with a similar setup where we specify its resolution as 2^-4. As before we choose a positive PMMH algorithm constant ϵ=10^-6. As we can see from Figure <ref> the MSE of the proposed methodology in the paper decays at the rate of 1/cost, which is as expected, which outperforms that of the single-level PMHH.

For multidimensional diffusions, it is well known that the exact methodology works only on specific diffusions, which require strong assumptions <cit.>.  It is not so clear how the exact methodology can be applied to our non-reversible diffusion (due to the difficult drift term) <cit.>.  An alternative to this is the work of Blanchet el al. <cit.>, which does not rely on the assumption of drift term equal to the gradient of a suitable potential function, or on Lamperti's transformation for that matter. However, despite this, the major drawback is that the running time of their methodology, although finite with probability one, has infinite mean. As a result, the comparison would not be practical due to the cost of the experiment.




§ ACKNOWLEDGMENTS
 
JF, AJ, KL and MV have received support from the Academy of Finland (grants 274740, 312605 and 315619) and from the Institute for Mathematical Sciences, Singapore (2018 programme `Bayesian Computation for High-Dimensional Statistical Models'). NC and AJ have received support from KAUST baseline funding, JF and KL from The Alan Turing Institute, AJ from the Singapore Ministry of Education (R-155-000-161-112), and KL from the University of Manchester (School of Mathematics). This research made use of the Rocket High Performance Computing service at Newcastle University. We thank Paul Fearnhead, Santeri Karppinen, Anthony Lee and the anonymous referees for their many insightful remarks.






§ ANALYSIS OF THE DELTA PARTICLE FILTER


We now give our analysis that is required for the proofs of Theorem <ref> and Corollary <ref> of Section <ref> regarding the ΔPF (Algorithm <ref>) for HMM diffusions.  
The structure of the appendix is as follows. 
In Section <ref> we introduce some more Feynman–Kac notations, following <cit.>, emphasising that here we consider standard HMMs that can be coupled.
In Section <ref> we recall the ΔPF stated earlier.  A general variance bound for quantities such as Δ_ℓ(φ) is given in Section
<ref>.  This is particularised to the HMM diffusion case in Section <ref>, where we supply the proofs for the results of Section <ref>.




 §.§ Models


Let (𝖷,𝒳) be a measurable space and {G_n}_n≥ 0 a sequence of non-negative, bounded and measurable functions such that G_n:𝖷→ℝ_+.
Let η_0^F,η_0^C∈𝒫(𝖷) and {M_n^F}_n≥ 1,   {M_n^C}_n≥ 1 be two sequences of Markov kernels, i.e. M_n^F:𝖷→𝒫(𝖷),
M_n^C:𝖷→𝒫(𝖷). 
Set 𝖤_n 𝖷^n+1 for n≥ 0, and for x_0:n∈𝖤_n,

    G_n(x_0:n) = G_n(x_n)

and for n≥ 1, s∈{F,C}, x_0:n-1∈𝖤_n-1

    M_n^s(x_0:n-1, x_0:n') = δ_{x_0:n-1}( x_0:n-1') M_n^s(x_n-1', x_n').

Define for s∈{F,C}, φ∈ℬ_b(𝖤_n), u_n∈𝖤_n

    γ_n^s(φ) = ∫_𝖤_0×⋯×𝖤_nφ(u_n) (∏_p=0^n-1G_p^s(u_p)) η_0^s( u_0)∏_p=1^n M_p^s(u_p-1, u_p)

and

    η_n^s(φ) = γ_n^s(φ)/γ_n^s(1).


Throughout this appendix, we assume Assumption <ref>, and that
Assumption <ref>(<ref>) holds, i.e. there exists η̌_0∈𝒫(𝖷×𝖷) such that
for any A∈𝒳

    η̌_0(A×𝖷) = η_0^F(A)     η̌_0(𝖷× A) = η_0^C(A)

and moreover for any n≥ 1 there exists Markov kernels {M̌_n}, M̌_n:𝖷×𝖷→𝒫(𝖷×𝖷)
such that for any A∈𝒳, (x,x')∈𝖷×𝖷:

    M̌_n(A×𝖷)(x,x') = M_n^F(A)(x)     M̌_n(𝖷× A)(x,x') = M_n^C(A)(x').






 §.§ Delta particle filter

Define x_p=(x_p^F,x_p^C)∈𝖷×𝖷 and

    Ǧ_p(x_p) = 1/2(G_p(x_p^F)+G_p(x_p^C)),

as in Assumption <ref>(<ref>).
Set, for n≥ 0, x_0:n∈𝖷^2(n+1)

    Ǧ_n(x_0:n) = Ǧ_n(x_n)

and for n≥ 1, x_0:n-1∈𝖷^2n

    M̌_n(x_0:n-1, x_0:n') = δ_{x_0:n-1}( x_0:n-1') M̌_n(x_n-1', x_n'),.

Note that coupling assumption (<ref>) for M̌_n can be equivalently formulated for M̌_n.

For n≥ 0, φ∈ℬ_b(𝖤_n×𝖤_n), u_n∈𝖤_n×𝖤_n, we have

    γ̌_n(φ) = ∫_𝖤_0^2×⋯×𝖤_n^2φ(u_n) (∏_p=0^n-1Ǧ_p(u_p)) η̌_0( u_0)∏_p=1^n M̌_p(u_p-1, u_p)

and

    η̌_n(φ) = γ̌_n(φ)/γ̌_n(1).


As noted in <cit.> it is simple to establish that for φ∈ℬ_b(𝖤_n), if

    ψ(x_0:n) =
    Ǧ_n(x_0:n)(φ(x_0:n^F)∏_p=0^n G_p(x_0:p^F)/Ǧ_p(x_0:p)
    -
    φ(x_0:n^C)∏_p=0^n G_p(x_0:p^C)/Ǧ_p(x_0:p))

then

    γ̌_n(ψ) = γ̌_n(1)η̌_n(ψ)
    =
    γ_n^F(G_nφ) - γ_n^C(G_nφ).

Note 

    γ̌_n(1) = ∏_p=0^n-1η̌_p(Ǧ_p).


 In order to approximate γ̌_n(ψ) one can run the following abstract version of Algorithm <ref> (recall from Section <ref> that we will only consider multinomial resampling). Define for n≥ 1, μ∈𝒫(𝖤_n-1×𝖤_n-1), φ∈ℬ_b(𝖤_n×𝖤_n)

    ϕ̌_n(μ)(φ) = μ(Ǧ_n-1M̌_n(φ))/μ(Ǧ_n-1).

The algorithm begins by generating u_0^i∈𝖤_0×𝖤_0, i∈{1,…,N} with joint law

    ∏_i=1^N η̌_0( u_0^i) = ∏_i=1^N η̌_0( u_0^i).

Defining 

    η̌_0^N( u_0) = 1/N∑_i=1^N δ_u_0^i( u_0)

we then generate  u_1^i∈𝖤_1×𝖤_1, i∈{1,…,N} with joint law

    ∏_i=1^N ϕ̌_1(η̌_0^N)( u_1^i).

This proceeds recursively, so the joint law of the particles up to time n is

    (∏_i=1^N η̌_0( u_0^i))(∏_p=1^n∏_i=1^N ϕ̌_p(η̌_p-1^N)( u_p^i)).

Hence we have the estimate

    γ̌^N_n(ψ) = (∏_p=0^n-1η̌_p^N(Ǧ_p))η̌^N_n(ψ).


Note that γ̌^N_n(ψ) corresponds to the quantity Δ_ℓ(φ) in (<ref>) from the ΔPF output (Algorithm <ref>).




 §.§ General hidden Markov model case


Define for p≥ 1 the semigroup

    Q̌_p(x_0:p-1, x_0:p') = Ǧ_p-1(x_0:p-1) M̌_p(x_0:p-1, x_0:p')

with the definition for 0≤ p ≤ n, φ∈ℬ_b(𝖤_n×𝖤_n)

    Q̌_p,n(φ)(u_p) = ∫φ(u_n) ∏_j=p+1^n  Q̌_j(u_j-1, u_j)

if p=n clearly Q̌_n,n is the identity operator. For any 0≤ n, φ∈ℬ_b(𝖤_n×𝖤_n)
we set Q̌_-1,n(φ)(u_-1)=0.

Now following <cit.> we have the following martingale (w.r.t. the natural filtration of the particle system), φ∈ℬ_b(𝖤_n×𝖤_n):

    γ̌^N_n(φ) - γ̌_n(φ) = 
    ∑_p=0^n  γ̌_p^N(1)[η̌_p^N-ϕ̌_p(η̌_p-1^N)](Q̌_p,n(φ))

with the convention that ϕ̌_p(η̌_p-1^N)=η̌_0 if p=0. The representation immediately establishes
that 

    𝔼[γ̌^N_n(φ)]=  γ̌_n(φ)

where the expectation is w.r.t. the law associated to the particle system. 
We will use the following convention that C' is a finite positive constant that does not depend upon n,N or any of the G_n, M_n^s (s∈{F,C}, M̌_n.
The value of C' may change from line-to-line.
Define for 0≤ p ≤ n<∞

    G_p,n = ∏_q=p^nG_q

with the convention that if p=0 we write G_n.
We have the following result.


  Suppose that G_n<∞ for each n≥ 0. Then there exist a C'<∞ such that for any n≥ 0, φ∈ℬ_b(𝖤_n×𝖤_n)  

    𝔼[(γ̌^N_n(φ) - γ̌_n(φ))^2]
    ≤C'/N∑_p=0^n G_p-1^2 𝔼[Q̌_p,n(φ)(u_p^1)^2].




Set 

    Š_p,n^N(φ)
    =
    γ̌_p^N(1)[η̌_p^N-ϕ̌_p(η̌_p-1^N)](Q̌_p,n(φ))

By (<ref>), one can apply the Burkholder-Gundy-Davis inequality to obtain

    𝔼[(γ̌^N_n(φ) - γ̌_n(φ))^2] ≤ 
    C'
    ∑_p=0^n 
    𝔼[Š_p,n^N(φ)^2].

Now, we have that

    𝔼[Š_p,n^N(φ)^2]
    ≤G_p-1^2𝔼[[η̌_p^N-ϕ̌_p(η̌_p-1^N)](Q̌_p,n(φ))^2].

Application of the (conditional) Marcinkiewicz-Zygmund inequality yields

    𝔼[Š_p,n^N(φ)^2]
    ≤C' G_p-1^2/N𝔼[
    (
      Q̌_p,n(φ)(u_p^1) -
    ϕ̌_p(η̌_p-1^N)( Q̌_p,n(φ) )
      )^2].

After applying C_2 and Jensen inequalities, we then conclude by (<ref>).




 §.§ Diffusion case


We now consider the model of Section <ref>, where we recall that θ is omitted from the notation. 
A series of technical results are given and the proofs for Theorem <ref> and Corollary <ref> are given at the end of this section.

We recall that the joint probability density of the observations and the unobserved diffusion at the observation times 
is given by

    ∏_p=0^n G_p(x_p) Q^(∞)(x_p-1,x_p).

As the true dynamics can not be simulated, in practice we work with

    ∏_p=0^n G_p(x_p) Q^(ℓ)(x_p-1,x_p).

Recall an (Euler) approximation scheme
with discretisation h_ℓ=2^-ℓ, ℓ≥ 0.
In our context then, M_n^F corresponds Q^(ℓ) (ℓ≥ 1) and
M_n^C corresponds Q^(ℓ-1). 
The initial distribution η_0 is simply the (Euler) kernel started at some given x_0.
As noted earlier in Remark <ref>(<ref>), a natural coupling of M_n^F and M_n^C (and hence of η_0) exists. 
As established in <cit.> one has (uniformly in  as Assumption <ref> holds with  independent constants)
for C'<∞

    sup_𝒜sup_x∈𝖷|M_n^F(φ)(x)-M_n^C(φ)(x)| ≤ C' h_ℓ

where 𝒜={φ∈ℬ_b(𝖷)∩Lip(𝖷): φ≤ 1|}.
We also recall that (<ref>) holds (recall Assumption <ref> is assumed).

We will use M<∞ to denote a constant that may change
from line-to-line. It will not depend upon θ nor N, ℓ, but may depend on the time parameter or a function.
The following result will be needed later on. The proof is given after the proof of Lemma <ref> below.



Assume (A<ref> (i)-(ii),<ref>). Then for any n≥ 0 and φ∈ℬ_b(𝖷^n+1)∩Lip(𝖷^n+1)
there exists a M<∞ such that

    |γ_n^F(G_nφ) - γ_n^C(G_nφ)| ≤ M h_ℓ




We write expectations w.r.t. the time-inhomogeneous Markov chain associated to the sequence of kernels (M_p^F)_p≥ 1 (resp. (M_p^C)_p≥ 1) 
as 𝔼^F, (resp. 𝔼^C). 


Assume (A<ref>(i)-(ii),<ref>). Let s∈{F,C} and φ∈ℬ_b(𝖷^n+1)∩Lip(𝖷^n+1),
then, define the function for 0≤ p ≤ n

    φ_p,n^s(x_0:p) := 𝔼^s[φ(x_0:p,X_p+1:n)∏_q=p+1^n G_q(X_q)|x_p].

Then we have that φ_p,n^s∈ℬ_b(𝖷^p+1)∩Lip(𝖷^p+1). 



The case p=n follows immediately from φ∈ℬ_b(𝖷^n+1)∩Lip(𝖷^n+1).
We will use a backward inductive argument on p. Suppose p=n-1 then we have for any (x_0:n-1,x_0:n-1')∈𝖷^n×𝖷^n

    |φ_n-1,n^s(x_0:n-1)-φ_n-1,n^s(x_0:n-1')| =


    |𝔼^s[φ(x_0:n-1,X_n)G_n(X_n)|x_n-1]-𝔼^s[φ(x_0:n-1',X_n)G_n(X_n)|x_n-1']| ≤


    |𝔼^s[φ(x_0:n-1,X_n)G_n(X_n)|x_n-1]-𝔼^s[φ(x_0:n-1',X_n)G_n(X_n)|x_n-1]| +


    |𝔼^s[φ(x_0:n-1',X_n)G_n(X_n)|x_n-1]-𝔼^s[φ(x_0:n-1',X_n)G_n(X_n)|x_n-1']|

By φ∈Lip(𝖷^n+1) it easily follows via (A<ref>(i)) that

    |𝔼^s[φ(x_0:n-1,X_n)G_n(X_n)|x_n-1]-𝔼^s[φ(x_0:n-1',X_n)G_n(X_n)|x_n-1]| ≤ M
    ∑_j=0^n-1 |x_j-x_j'|.

By (A<ref>(ii)) and φ∈Lip(𝖷^n+1), φ(x_0:n)G_n(x_n) is Lipschitz in x_n and hence by (A<ref>)

    |𝔼^s[φ(x_0:n-1',X_n)G_n(X_n)|x_n-1]-𝔼^s[φ(x_0:n-1',X_n)G_n(X_n)|x_n-1']|
    ≤ M|x_n-1-x_n-1'|.

Hence it follows

    |φ_n-1,n^s(x_0:n-1)-φ_n-1,n^s(x_0:n-1')| ≤ M ∑_j=0^n-1 |x_j-x_j'|.

The induction step follows by almost the same argument as above and is hence omitted.



We have the following standard collapsing sum representation:

    γ_n^F(G_nφ) - γ_n^C(G_nφ)      =    ∑_p=0^n(
    𝔼^F[∏_q=0^p G_q(X_q)𝔼^C[φ(X_0:n)∏_q=p+1^n G_q(X_q)|X_p]]
     - 
          𝔼^F[∏_q=0^p-1 G_q(X_q)𝔼^C[φ(X_0:n)∏_q=p^n G_q(X_q)|X_p-1]])


The summand is

    T_p:=𝔼^F[(∏_q=0^p-1 G_q(X_q)) (𝔼^F-𝔼^C)(𝔼^C[φ(X_0:n)∏_q=p+1^n G_q(X_q)|X_p]G_p(X_p)|X_p-1)].

By Lemma <ref>, 𝔼^C[φ(x_0:p,X_p+1:n)∏_q=p+1^n G_q(X_q)|x_p]∈ℬ_b(𝖷^p+1)∩Lip(𝖷^p+1)
and by (A<ref>) (i) and (ii) G_p∈ℬ_b(𝖷)∩Lip(𝖷). So by (<ref>)

    |(𝔼^F-𝔼^C)(𝔼^C[φ(X_0:n)∏_q=p+1^n G_q(X_q)|X_p]G_p(X_p)|X_p-1)|
    ≤


    M h_ℓsup_x_0:p∈𝖷^p+1|𝔼^C[φ(x_0:p,X_p+1:n)∏_q=p+1^n G_q(X_q)|

and hence

    |T_p| ≤ M h_ℓ𝔼^F[∏_q=0^p-1 G_q(X_q)] sup_x_0:p∈𝖷^p+1|𝔼^C[φ(x_0:p,X_p+1:n)∏_q=p+1^n G_q(X_q)|x_p]G_p(x_p)|.

Application of (A<ref>) (i) gives |T_p| ≤ M h_ℓ and the proof is hence concluded. 



Assume (A<ref>). Then for any n≥ 0 there exists a M<∞ such that for any x_0:n∈𝖷^2(n+1)

    |∏_p=0^n G_p(x_p^F)/Ǧ_p(x_p) - ∏_p=0^n G_p(x_p^C)/Ǧ_p(x_p)| ≤ M ∑_p=0^n |x_p^F-x_p^C|.




The is proof by induction. The case n=0:

    |G_0(x_0^F)/Ǧ_0(x_0) - G_0(x_0^C)/Ǧ_0(x_0)| = 1/Ǧ_0(x_0)|G_0(x_0^F) - G_0(x_0^C)|.

Application of (A<ref>) (ii) and (iii) yield that

    |G_0(x_0^F)/Ǧ_0(x_0) - G_0(x_0^C)/Ǧ_0(x_0)| ≤ M |x_0^F-x_0^C|.

The result is assumed to hold at rank n-1, then 

    |∏_p=0^n G_p(x_p^F)/Ǧ_p(x_p) - ∏_p=0^n G_p(x_p^C)/Ǧ_p(x_p)| ≤


    |G_n(x_n^F)/Ǧ_n(x_n) - G_n(x_n^C)/Ǧ_n(x_n)|·∏_p=0^n-1G_p(x_p^F)/Ǧ_p(x_p)
    + |∏_p=0^n-1G_p(x_p^F)/Ǧ_p(x_p) - ∏_p=0^n-1G_p(x_p^C)/Ǧ_p(x_p)|
    ·G_n(x_n^C)/Ǧ_n(x_n).

For the first term of the R.H.S. one can follow the argument at the initialisation and apply (A<ref>) (i) and (iii). For the second term of the R.H.S., the induction
hypothesis and (A<ref>) (i) and (iii) can be used. That is one can deduce that 

    |∏_p=0^n G_p(x_p^F)/Ǧ_p(x_p) - ∏_p=0^n G_p(x_p^C)/Ǧ_p(x_p)| ≤ M ∑_p=0^n |x_p^F-x_p^C|.



Recall (<ref>) for the definition of ψ and that
x_p=(x_p^F,x_p^C)∈𝖷×𝖷.


Assume (A<ref>-<ref>). Then 
for any 0≤ p < n, φ∈ℬ_b(𝖷^n+1)∩Lip(𝖷^n+1) there exists a M<∞ such that for any
x_0:p∈𝖤_p×𝖤_p

    |Q̌_p,n(ψ)(x_0:p)| ≤  M (∑_j=0^p |x_j^F-x_j^C| + h_ℓ)




We have 

    Q̌_p,n(ψ)(x_0:p) = Ǧ_p(x_p) ×(   ∏_q=0^p G_q(x_q^F)/Ǧ_q(x_q)𝔼^F[φ(x_0:p^F,Y_p+1:n) ∏_s=p+1^n G_s(X_s^F)|x_p^F] 
       -∏_q=0^p G_q(x_q^C)/Ǧ_q(x_q)𝔼^C[φ(x_0:p^C,Y_p+1:n)∏_s=p+1^n G_s(X_s^C)|x_p^C]).

It then follows that Q̌_p,n(ψ)(x_0:p) = Ǧ_p(x_p) ( T_1+T_2) where

    T_1     =    (∏_q=0^p G_q(x_q^F)/Ǧ_q(x_q) - ∏_q=0^p G_q(x_q^C)/Ǧ_q(x_q)) 𝔼^F[φ(x_0:p^F,Y_p+1:n)∏_s=p+1^n G_s(X_s^F)|x_p^F] 
    
    T_2     =    ∏_q=0^p G_q(x_q^C)/Ǧ_q(x_q)(𝔼^F[φ(x_0:p^F,Y_p+1:n)∏_s=p+1^n G_s(X_s^F)|x_p^F]-𝔼^C[φ(x_0:p^C,Y_p+1:n)∏_s=p+1^n G_s(X_s^C)|x_p^C]).

By Lemma <ref>, φ∈ℬ_b(𝖷^n+1)∩Lip(𝖷^n+1) and (A<ref>) (i)

    |T_1| ≤ M ∑_j=0^p |x_j^F-x_j^C|.


Now T_2 = T_3 + T_4 where

    T_3     =    ∏_q=0^p G_q(x_q^C)/Ǧ_q(x_q)(𝔼^F[φ(x_0:p^F,Y_p+1:n)∏_q=p+1^n G_s(X_s^F)|x_p^F] - 
    𝔼^F[φ(x_0:p^F,Y_p+1:n)∏_s=p+1^n G_s(X_s^F)|x_p^C]) 
    
    T_4     =    ∏_q=0^p G_q(x_q^C)/Ǧ_q(x_q)(𝔼^F[φ(x_0:p^F,Y_p+1:n)∏_s=p+1^n G_s(X_s^F)|x_p^C]-𝔼^C[φ(x_0:p^C,Y_p+1:n)∏_s=p+1^n G_s(X_s^C)|x_p^C]).

For T_3 one can use Lemma <ref> (along with (A<ref>) (i) and (iii)) to get that

    |T_3| ≤ M ∑_j=0^p |x_j^F-x_j^C|.

For T_4 a similar collapsing sum argument that is used in the proof of Proposition <ref> can be used to deduce that

    |T_4| ≤ M h_ℓ.

One can then conclude the proof via the above bounds (along with (A<ref>) (i)).


Below 𝔼 denotes expectation w.r.t. the particle system described in Section <ref> started at position (x,x) at time n=0 with x∈𝖷, in the diffusion case of Section <ref>.  Recall the particle U_n^i∈𝖤_n×𝖤_n at time n≥ 0 in path space. We denote by U_n^i,s(j)∈𝖷 as the j∈{0,…,n} component of
particle i∈{1,…,N} at time n≥ 0 of s∈{F,C} component. Recall (U_n^i,F(n),U_n^i,C(n)) for n≥ 1 is sampled from the kernel M̌_n((u̅_n-1^i,F(n-1),u̅_n-1^i,C(n-1)),)  where the u̅ denotes post-resampling and the component (U_n^i,F(j),U_n^i,C(j))=(u̅_n-1^i,F(j),u̅_n-1^i,C(j)) for j∈{0,…,n-1} is kept the same for the earlier components of the particle.


Assume (A<ref> (i) (iii), <ref>). Then for any n≥ 0 there exists a M<∞ such that

    𝔼[∑_j=0^n |U_n^1,F(j)-U_n^1,C(j)|^2] ≤ M h_ℓ^β.

where β is as in (<ref>).



Our proof is by induction, the case n=0 following by (<ref>). Assuming the result at n-1 we have

    𝔼[∑_j=0^n |U_n^1,F(j)-U_n^1,C(j)|^2] = 𝔼[∑_j=0^n-1 |U̅_n-1^1,F(j)-U̅_n-1^1,C(j)|^2 + |U_n^1,F(n)-U_n^1,C(n)|^2 ].


Now

    𝔼[∑_j=0^n-1 |U̅_n-1^1,F(j)-U̅_n-1^1,C(j)|^2]     =     N∑_j=0^n-1𝔼[
    Ǧ_n-1(U_n-1^1,F(n-1),U_n-1^1,C(n-1))/∑_j=1^N Ǧ_n-1(U_n-1^j,F(n-1),U_n-1^j,C(n-1))×
           
    |U_n-1^1,F(j)-U_n-1^1,C(j)|^2] 
       ≤   
    M𝔼[∑_j=0^n-1|U_n-1^1,F(j)-U_n-1^1,C(j)|^2]

where we have used (A<ref>) (i) and (iii). 
Applying the induction hypothesis along with (<ref>) yields

    𝔼[∑_j=0^n |U_n^1,F(j)-U_n^1,C(j)|^2] ≤ M(h_ℓ^β + 𝔼[|U̅_n-1^1,F(n-1)-U̅_n-1^1,C(n-1)|^2])

Now

    𝔼[|U̅_n-1^1,F(n-1)-U̅_n-1^1,C(n-1)|^2] =


    N 𝔼[Ǧ_n-1(U_n-1^1,F(n-1),U_n-1^1,C(n-1))/∑_j=1^N Ǧ_n-1(U_n-1^j,F(n-1),U_n-1^j,C(n-1))|U_n-1^1,F(n-1)-U_n-1^1,C(n-1)|^2]

Then by (A<ref>) (i) and (iii)

    𝔼[Ǧ_n-1(U_n-1^1,F(n-1),U_n-1^1,C(n-1))/∑_j=1^N Ǧ_n-1(U_n-1^j,F(n-1),U_n-1^j,C(n-1))|U_n-1^1,F(n-1)-U_n-1^1,C(n-1)|^2] ≤


    M/N𝔼[|U_n-1^1,F(n-1)-U_n-1^1,C(n-1)|^2] ≤M/N𝔼[∑_j=0^n-1 |U_n-1^1,F(j)-U_n-1^1,C(j)|^2].

Hence via the induction hypothesis, one has

    𝔼[|U̅_n-1^1,F(n-1)-U̅_n-1^1,C(n-1)|^2]  ≤ M h_ℓ^β

and the proof is concluded.


Recall Remark <ref>.


This follows first by applying Proposition <ref>, followed by Lemma <ref> and then some standard calculations followed by Lemma <ref>.



Easily follows by adding and subtracting γ̌_n(ψ) the C_2 inequality along with Theorem <ref>, and
then using (<ref>) combined with Proposition <ref>.




§ PROOF OF CONSISTENCY OF THE MARKOV CHAIN MONTE CARLO


 
  Denote
  
    ξ_k() (∑_i=1^ V_k^(i) +
          ϵ)^-1[∑_i=1^
          V_k^(i)(Θ_k, X_k^(i)) +
          _k(^(Θ_k))],

       where
      ^(θ)(x) (θ, x) and
       _k(^(θ)) p_L_k^-1∑_i=1^2  V_k,L_k^(i)^(θ)( 𝐗_k,L_k^(i)).
      Then
      E_,,𝐩() = ∑_k=1^ξ_k()/∑_j=1^ξ_k(1).
Furthermore, by Assumption <ref> <cit.>,
we have

    [_k^2() |Θ_k=θ]     = s_(θ), 
    [_k() |Θ_k=θ]    =
    _^(θ,∞)(_)
    -_^(θ,0)(_)

for = 1 and = ^(θ).
This implies for = and = 1, 

    μ_(θ, v^(1:),  𝐱^(1:))    [ξ_k()| (Θ_k, V_k^(1:), 𝐗_k^(1:))
        =(θ,v^(1:),     𝐱^(1:))] 
       = 1/∑_j=1^ v^(j) + ϵ[∑_i=1^
          v^(i) g(θ, x^(i)) 
          - _^(θ,0)(_)
          + _^(θ,∞)(_)], 
    
        m_^(1)(θ, v^(1:),  𝐱^(1:))
           [|ξ_k(g)|| (Θ_k, V_k^(1:), 𝐗_k^(1:))
        =(θ,v^(1:),     𝐱^(1:))]  
       ≤1/∑_j=1^ v^(j) + ϵ[∑_i=1^
          v^(i) |(θ, x^(i))|
          + √(s_^(θ)(θ))].

It is direct to check that the PMMH type chain (Θ_k,
X_k^(1:), V_k^(1:))) is reversible with respect to the
probability 

    (θ,  x^(1:),  v^(1:))
        = c_0 pr(θ) θ_θ^(0)( x^(1:),
         v^(1:)) (∑_i=1^ v^(i) + ϵ),

where c_0>0 is a normalisation constant and _θ^(0)()
stands for the law of the output of Algorithm <ref> with
(M_0:^(θ, 0), G_0:^(θ,0),), and
therefore is Harris recurrent as a full-dimensional
Metropolis–Hastings that is ψ-irreducible <cit.>. It is direct to check that 
(m_^(1))<∞, 
(m_1^(1))<∞, (μ_) = c
π^(∞)() and (μ_1) = c, where c>0
is a constant, so the result follows from <cit.>.
 





§ PROOFS ABOUT ASYMPTOTIC EFFICIENCY AND ALLOCATIONS



  By Harris ergodicity,
  ^-1()
  →
almost surely. 
Dividing the inequality
  
    ((κ)) ≤κ < ((κ) +1)

by (κ) and taking the limit κ→∞, which implies (κ)→∞, we get that κ/(κ)→ almost surely.  Also, by Proposition <ref>,

    √((κ))[ E_(κ),,𝐩() - π^(∞)() ]
        𝒩(0,σ^2),    in distribution,

so the result follows by Slutsky's theorem.    


  We have that
  
    [()]
      =
      ∑_k=1^[τ_Θ_k,L_k]
      =
          ∑_k=1^∑_ℓ=1^∞[τ_Θ_k,ℓ] p_ℓ.
          ≤
    C ∑_k=1^∞ p_ℓ 2^ℓ(1+),

by Assumption <ref>(<ref>), 
which is finite if r>(1+).  Also,

    s_(θ)
    =
    [_k^2()| Θ_k=θ]
    =
    ∑_ℓ≥ 1Δ_ℓ^2/p_ℓ≤
    C ∑(
    2^-ℓ(β +  - r) + 2^-ℓ(2α -r)),
      
which is finite if r< min(β+, 2α).  Therefore, σ^2<∞, and the CLT follows by Proposition <ref>.
  



Let {X_k}_k≥ 1 be a sequence of independent random variables with [X_k_0]=∞ for at least one k_0, and let {a_k}_k≥ 1 be a sequence of monotonically increasing real numbers with a_k/k ⟶∞.  Suppose one of the following assumptions holds:


  * 
  ∑_k≥ 1[ X_k > a_k]<∞,
and {X_k}_k≥ 1 are also identically distributed, or

  * 
  ∑_k≥ 1sup_≥ 1[ X_ > a_k]
<
∞.

Then

    [ ∑_k=1^ X_k > a_ infinitely many ∈ℕ]=0.



(<ref>) is <cit.> since [X_k_0]=∞ implies [X_k]=∞ for all k≥ 1 as {X_k}_k≥ 1 are i.i.d.  For (<ref>), note that if X_k has c.d.f. denoted F_k, then it is straightforward to check that
  
    F^*(x)
      inf_k≥ 1 F_k(x)

  is a c.d.f. also.  With X^*_k∼ F^* i.i.d. for k≥ 1, we have
  
    [X_k^* > a_k]
      =
      1- F^*(a_k)
      =
      sup_≥ 1 1-F_(a_k)
      =
      sup_≥ 1 P[X_ > a_k].

  Summing over k≥ 1, we obtain ∑_k≥ 1[X_k^* > a_k] < ∞.
  In addition,
  
    [X_k^*]
      =
      ∫[X_k^* > x]  x
      ≥∫[X_k_0 > x]  x
      =
      ∞,

for all k≥ 1.  Hence, we can apply (<ref>) for i.i.d. random variables, obtaining

    0
    =
    [∑_k=1^ X^*_k > a_ infinitely many ]
    ≥[∑_k=1^ X_k > a_ infinitely many ],

where the first equality comes from (<ref>), and so we conclude.



  Conditional on output {Θ_k}_k ≥ 1 of Algorithm <ref>, {τ_Θ_k,L_k}_k≥ 1 are independent random variables.  Our assumptions imply Lemma <ref>(<ref>) holds, so
  
    [() > a_ infinitely many ] =0,

which means that () is asymptotically bounded by a_.  Setting =O(ϵ^-2) allows us to conclude.
  

The proofs below of Proposition <ref> and <ref> are similar to that of <cit.>.   

  With the prescribed choice of p_ℓ we have finite variance, as 
  
    s_(θ)
      =
      ∑_ℓ≥ 1Δ_ℓ^2 /p_ℓ≤ C
      ∑_ℓ≥11/ℓ [ log_2 (ℓ +1)]^η
      <∞,

  uniformly in θ∈𝖳.  To determine the order of complexity, we would like to apply Lemma <ref>(<ref>) to the i.i.d sequence {τ_L_k^*}_k≥ 1, where τ_ℓ^* C2^ℓ(1+). For any k≥ 1, where a_k>0 is some positive real number, we have,
  
    [τ_L_k^* > a_k]
    =
      ∑_ℓ≥ 1[τ_ℓ^* > a_k] p_ℓ
      =
      ∑_ℓ≥ 11{ℓ > 1/(1+)log_2a_k/C} p_ℓ.

  Because ∑_ℓ≥ 1 p_ℓ = 1 and p_ℓ is monotonically decreasing, we have ∑_ℓ≥ℓ_* p_ℓ is O(p_ℓ_*).  Setting ℓ_*=⌊1/(1+)log_2 a_k/C⌋, we therefore obtain that (<ref>) is of order  
  
    a_k^-2b/(1+)(log_2 a_k)( log_2log_2 a_k)^η.

Setting

    a_k [k (log_2 k)^q]^(1+)/2 b

then ensures that ∑_k ≥ 1[τ_L_k^* > a_k] <∞.  As β≤ 1, it is easy to check that [τ_L_k^*]=∞.  We then apply Lemma <ref>(<ref>), obtaining

    0
    =
    [ ∑_k=1^τ_L_k^* > a_ infinitely many ]
    ≥[ ∑_k=1^τ_Θ_k, L_k > a_ infinitely many ].

and conclude as before, by using that () is asymptotically bounded by a_ and setting =O(ϵ^-2).  
  


We are in the basic setting of Proposition <ref> as before, but additionally may choose ≥ 0 as we please.  
The growth of a_k given in (<ref>) is essentially determined by (1+)/2b, which can be made small when = 2α-β, implying b=α.
  


abbrvnat


34


urlstyle

[Andrieu and Thoms(2008)]andrieu-thoms
C. Andrieu and J. Thoms.
A tutorial on adaptive MCMC.
Statist. Comput., 180 (4):0 343–373, Dec.
  2008.

[Andrieu et al.(2010)Andrieu, Doucet, and
  Holenstein]andrieu-doucet-holenstein
C. Andrieu, A. Doucet, and R. Holenstein.
Particle Markov chain Monte Carlo methods.
J. R. Stat. Soc. Ser. B Stat. Methodol., 720
  (3):0 269–342, 2010.
(with discussion).

[Beskos and Roberts(2005)]beskos-roberts
A. Beskos and G. Roberts.
Exact simulation of diffusions.
Ann. Appl. Probab., 150 (4):0 2422–2444, 11
  2005.
10.1214/105051605000000485.

[Beskos et al.(2006)Beskos, Papaspiliopoulos, Roberts, and
  Fearnhead]beskos-papaspiliopoulos-roberts-fearnhead
A. Beskos, O. Papaspiliopoulos, G. O. Roberts, and P. Fearnhead.
Exact and computationally efficient likelihood-based estimation for
  discretely observed diffusion processes.
J. R. Stat. Soc. Ser. B Stat. Methodol., 680
  (3):0 333–382, 2006.
(with discussion).

[Blanchet and Zhang(2020)]blanchet-zhang
J. Blanchet and F. Zhang.
Exact simulation for multivariate Itô diffusions.
Adv. Appl. Probab., 520 (4), 2020.

[Cappé et al.(2005)Cappé, Moulines, and Ryden]cappe
O. Cappé, E. Moulines, and T. Ryden.
Inference in Hidden Markov Models.
Springer, New York, 2005.

[Cliffe et al.(2011)Cliffe, Giles, Scheichl, and
  Teckentrup]cliffe2011multilevel
A. Cliffe, M. Giles, R. Scheichl, and A. Teckentrup.
Multilevel Monte Carlo methods and applications to elliptic
  PDEs with random coefficients.
Comput. Vis. Sci., 140 (1):0 3, 2011.

[Del Moral(2004)]del-moral
P. Del Moral.
Feynman-Kac Formulae.
Springer, New York, 2004.

[Douc et al.(2005)Douc, Cappé, and Moulines]douc-cappe-moulines
R. Douc, O. Cappé, and E. Moulines.
Comparison of resampling schemes for particle filtering.
In Proc. Image and Signal Processing and Analysis, 2005, pages
  64–69, 2005.

[Doucet et al.(2015)Doucet, Pitt, Deligiannidis, and
  Kohn]doucet-pitt-deligiannidis-kohn
A. Doucet, M. Pitt, G. Deligiannidis, and R. Kohn.
Efficient implementation of Markov chain Monte Carlo when using
  an unbiased likelihood estimator.
Biometrika, 1020 (2):0 295–313, 2015.

[Fearnhead et al.(2008)Fearnhead, Papaspiliopoulos, and
  Roberts]fearnheadPR
P. Fearnhead, O. Papaspiliopoulos, and G. O. Roberts.
Particle filters for partially observed diffusions.
J. R. Stat. Soc. Ser. B Stat. Methodol., 700
  (4):0 755–777, 2008.

[Fearnhead et al.(2010)Fearnhead, Papaspiliopoulos, Roberts, and
  Stuart]fearnheadPRS
P. Fearnhead, O. Papaspiliopoulos, G. O. Roberts, and A. Stuart.
Random-weight particle filtering of continuous time processes.
J. R. Stat. Soc. Ser. B Stat. Methodol., 720
  (4):0 497–512, 2010.

[Fearnhead et al.(2017)Fearnhead, Latuszynski, Roberts, and
  Sermaidis]fearnheadLRS
P. Fearnhead, K. Latuszynski, G. Roberts, and G. Sermaidis.
Continuous-time importance sampling: Monte Carlo methods which
  avoid time-discretisation error.
Preprint arXiv:1712.06201, 2017.

[Feller(1946)]feller-infinite
W. Feller.
A limit theorem for random variables with infinite moments.
Amer. J. Math., 680 (2):0 257–262, 1946.

[Franks and Vihola(2020)]franks-vihola
J. Franks and M. Vihola.
Importance sampling correction versus standard averages of reversible
  MCMCs in terms of the asymptotic variance.
Stochastic Process. Appl., 1300 (10), 2020.

[Giles and Szpruch(2014)]giles-szpruch
M. Giles and L. Szpruch.
Antithetic multilevel Monte Carlo estimation for
  multi-dimensional SDEs without Lévy area simulation.
Ann. Appl. Probab., 240 (4):0 1585–1620,
  2014.

[Giles(2008)]giles-or
M. B. Giles.
Multilevel Monte Carlo path simulation.
Oper. Res., 560 (3):0 607–617, 2008.

[Glynn and Whitt(1992)]glynn-whitt
P. Glynn and W. Whitt.
The asymptotic efficiency of simulation estimators.
Oper. Res., 400 (3):0 505–520, 1992.

[Golightly and Wilkinson(2011)]golightly-wilkinson-bayesian-2011
A. Golightly and D. Wilkinson.
Bayesian parameter inference for stochastic biochemical network
  models using particle Markov chain Monte Carlo.
Interface focus, 10 (6):0 807–820, 2011.

[Heinrich(2001)]heinrich
S. Heinrich.
Multilevel Monte Carlo methods.
In Large-scale scientific computing, pages 58–67. Springer,
  2001.

[Jasra et al.(2017)Jasra, Kamatani, Law, and Zhou]mlpf
A. Jasra, K. Kamatani, K. J. H. Law, and Y. Zhou.
Multilevel particle filters.
SIAM J. Numer. Anal., 55:0 3068–3096, 2017.

[Jasra et al.(2018a)Jasra, Kamatani, Law, and
  Zhou]jasra-kamatani-law-zhou
A. Jasra, K. Kamatani, K. J. H. Law, and Y. Zhou.
Bayesian static parameter estimation for partially observed
  diffusions via multilevel Monte Carlo.
SIAM J. Sci. Comp., 40:0 A887–A902,
  2018a.

[Jasra et al.(2018b)Jasra, Kamatani, Law, and
  Zhou]jasra2018multi
A. Jasra, K. Kamatani, K. J. H. Law, and Y. Zhou.
A multi-index Markov chain Monte Carlo method.
Intern. J. Uncertainty Quantif., 80 (1),
  2018b.

[Jasra et al.(2019)Jasra, Law, and Osei]jasraLO
A. Jasra, K. J. Law, and P. P. Osei.
Multilevel particle filters for Lévy-driven stochastic
  differential equations.
Stat. Comp., 29:0 775–789, 2019.

[Kloeden and Platen(1999)]kloeden-platen
P. Kloeden and E. Platen.
Numerical solution of stochastic differential equations.
Springer, Berlin Heidelberg, 3rd edition, 1999.

[McLeish(2011)]mcleish
D. McLeish.
A general method for debiasing a Monte Carlo estimator.
Monte Carlo Methods Appl., 170 (4):0 301–315,
  2011.

[Pavliotis(2016)]pavliotis
G. Pavliotis.
Stochastic Processes and Applications.
Springer, 2016.

[Rhee and Glynn(2015)]rhee-glynn
C.-H. Rhee and P. W. Glynn.
Unbiased estimation with square root convergence for SDE models.
Oper. Res., 630 (5):0 1026–1043, 2015.

[Roberts and Rosenthal(2006)]roberts-rosenthal-harris
G. Roberts and J. Rosenthal.
Harris recurrence of Metropolis-within-Gibbs and
  trans-dimensional Markov chains.
Ann. Appl. Probab., 160 (4):0 2123–2139,
  2006.

[Sherlock et al.(2015)Sherlock, Thiery, Roberts, and
  Rosenthal]sherlock-thiery-roberts-rosenthal
C. Sherlock, A. H. Thiery, G. O. Roberts, and J. S. Rosenthal.
On the efficiency of pseudo-marginal random walk Metropolis
  algorithms.
Ann. Statist., 430 (1):0 238–275, 2015.

[Sørensen(2004)]sorensen
H. Sørensen.
Parametric inference for diffusion processes observed at discrete
  points in time: a survey.
Intern. Statist. Review, 720 (3):0 337–354,
  2004.

[Vihola(2018)]vihola-unbiased
M. Vihola.
Unbiased estimators and multilevel Monte Carlo.
Oper. Res., 660 (2):0 448–462, 2018.

[Vihola et al.(2020)Vihola, Helske, and Franks]vihola-helske-franks
M. Vihola, J. Helske, and J. Franks.
Importance sampling type estimators based on approximate marginal
  MCMC.
Scand. J. Statist., 470 (4), 2020.

[Wang et al.(2019)Wang, Rao, and Teh]wang-rao-teh
Q. Wang, V. Rao, and Y. W. Teh.
An exact auxiliary variable Gibbs sampler for a class of
  diffusions.
Preprint arXiv:1903.10659, 2019.










