2.1. Notation

Let $(\textsf{X},\mathcal{X})$ be a measurable space. For a given function $v\,:\,\textsf{X}\rightarrow[1,+\infty)$ and for measurable $\varphi\,:\,\textsf{X}\rightarrow\mathbb R$ ,

\[\|\varphi\|_v \,:\!= \sup_{x\in \textsf{X}}\frac{|\varphi(x)|}{v(x)}.\]

We define $\mathcal{L}_v(\textsf{X})=\{\varphi\,:\,\textsf{X}\rightarrow\mathbb R\,:\,\|\varphi\|_v<+\infty\}$ . When $v=1$ we write $\|\varphi\| \,:\!= \sup_{x\in \textsf{X}}|\varphi(x)|$ . We write $\mathcal{B}_b(\textsf{X})$ , $\mathcal{C}_b(\textsf{X})$ for the bounded measurable and continuous, bounded measurable real-valued functions respectively. $\mathcal{C}^2(\textsf{X})$ is the collection of twice continuously differentiable real-valued functions on $\textsf{X}$ . Let $\textsf{d}$ be a metric on $\textsf{X}$ ; then for $\varphi\in\mathcal{L}_v(\textsf{X})$ , we say that $\varphi\in\textrm{Lip}_{v,\textsf{d}}(\textsf{X})$ if there exists a $C<+\infty$ such that for every $(x,y)\in\textsf{X}\times\textsf{X}$ ,

\[|\varphi(x)-\varphi(y)| \leq C\textsf{d}(x,y)v(x)v(y).\]

If $v=1$ , we write $\varphi\in\textrm{Lip}_{\textsf{d}}(\textsf{X})$ and write $\|\varphi\|_{\textrm{Lip}}$ for the Lipschitz constant. $\mathscr{P}(\textsf{X})$ denotes the collection of probability measures on $(\textsf{X},\mathcal{X})$ . We also write $\|\mu\|_{v}\,:\!=\sup_{|\varphi|\leq v}|\mu(\varphi)|$ for $\mu\in\mathscr{P}(\textsf{X})$ . For a measure $\mu$ on $(\textsf{X},\mathcal{X})$ and a $\mu-$ integrable function $\varphi\,:\,\textsf{X}\rightarrow\mathbb{R}$ , the notation $\mu(\varphi)=\int_{\textsf{X}}\varphi(x)\mu(dx)$ is used. Let $K\,:\,\textsf{X}\times\mathcal{X}\rightarrow(0,+\infty)$ be a nonnegative kernel and $\mu$ a measure; then we use the notation $\mu K(dy) = \int_{\textsf{X}}\mu(dx) K(x,dy)$ , and for $K(x,\cdot)-$ integrable $\varphi\,:\,\textsf{X}\rightarrow\mathbb{R}$ , we write $K(\varphi)(x) = \int_{\textsf{X}} \varphi(y) K(x,dy).$ For $\mu,\nu\in\mathscr{P}(\textsf{X})$ , the total variation distance is denoted $\|\mu-\nu\|_{\textrm{tv}}=\sup_{B\in\mathcal{X}}|\mu(B)-\nu(B)|$ . For $B\in\mathcal{X}$ the indicator function is written $\mathbb{I}_B(x)$ and the Dirac measure $\delta_B(dx)$ . For two measures $\mu,\nu$ on $(\textsf{X},\mathcal{X})$ , the product measure is denoted $\mu\otimes \nu$ . For two measurable functions $\varphi$ , $\psi$ on $(\textsf{X},\mathcal{X})$ , the tensor product of functions is denoted $\varphi\otimes\psi$ . $\mathcal{U}_A$ denotes the uniform distribution on the set A. $\mathcal{N}_t(a,b)$ is the $t-$ dimensional Gaussian distribution of mean a and covariance b (if $t=1$ the subscript is dropped from $\mathcal{N}$ ). $\mathbb{P}$ and $\mathbb{E}$ are used to denote probability and expectation with respect to the law of the specified algorithm—the context will be clear in each instance. The symbols $\Rightarrow$ and $\rightarrow_{\mathbb{P}}$ are used to denote convergence in distribution and probability respectively. In the context of the article, this is as $N\rightarrow+\infty$ .

2.2. Models

Let $(\textsf{X},\mathcal{X})$ be a measurable space and $\{G_n\}_{n\geq 0}$ a sequence of nonnegative, bounded, and measurable functions such that $G_n\,:\,\textsf{X}\rightarrow\mathbb{R}_+$ . Let $\eta_0^{f},\eta_0^c\in\mathscr{P}(\textsf{X})$ and let $\{M_n^f\}_{n\geq 1}$ , $\{M_n^c\}_{n\geq 1}$ be two sequences of Markov kernels, i.e. $M_n^f\,:\,\textsf{X}\rightarrow\mathscr{P}(\textsf{X})$ , $M_n^c\,:\,\textsf{X}\rightarrow\mathscr{P}(\textsf{X})$ . For $s\in\{\,f,c\}$ , $B\in\mathcal{X}$ , define

\[\gamma_{n}^s(B) = \int_{\textsf{X}^{n+1}}\mathbb{I}_B(x_n) \Bigg(\prod_{p=0}^{n-1} G_p^s(x_p)\Bigg) \eta_0^s(dx_0)\prod_{p=1}^n M_p^s(x_{p-1},dx_p)\]

and

\[\eta_n^s(B) = \frac{\gamma_{n}^s(B)}{\gamma_{n}^s(1)}.\]

The objective is to consider Monte-Carlo–type algorithms which will approximate, for $\varphi\in\mathcal{B}_b(\textsf{X})$ and any $n\geq 0$ , quantities such as

(1) \begin{equation}\eta_n^f(\varphi) - \eta_n^c(\varphi)\end{equation}

or

(2) \begin{equation}\frac{\eta_n^f(G_n\varphi)}{\eta_n^f(G_n)} - \frac{\eta_n^c(G_n\varphi)}{\eta_n^c(G_n)}.\end{equation}

(1) corresponds to a predictor of a state-space model and (2) to the filter. We focus explicitly on the predictor from here on.

The major point is that one would like to approximate couplings of $(\eta_n^f,\eta_n^c)$ , say $\check{\eta}_n\in\mathscr{P}(\textsf{X}\times\textsf{X})$ , i.e. that for any $B\in\mathcal{X}$ and every $n\geq 0$ ,

\[\check{\eta}_n(B\times\textsf{X}) = \eta_n^f(B), \qquad \check{\eta}_n(\textsf{X}\times B) = \eta_n^c(B),\]

and we consider approximating

\[\check{\eta}_n(\varphi\otimes 1) - \check{\eta}_n(1\otimes \varphi).\]

An explanation of why coupling the pairs is of interest has been given in the introduction and will be further illuminated in Section 2.3.

Throughout the article, it is assumed that $\check{\eta}_0\in\mathscr{P}(\textsf{X}\times\textsf{X})$ is such that for any $B\in\mathcal{X}$ ,

\[\check{\eta}_0(B\times \textsf{X}) = \eta_0^f(B), \qquad \check{\eta}_0(\textsf{X}\times B) = \eta_0^c(B),\]

and moreover for any $n\geq 1$ there exist Markov kernels $\{\check{M}_n\}$ , $\check{M}_n\,:\,\textsf{X}\times\textsf{X}\rightarrow\mathscr{P}(\textsf{X}\times\textsf{X})$ , such that for any $B\in\mathcal{X}$ , $(x,x')\in\textsf{X}\times\textsf{X}$ ,

(3) \begin{equation}\check{M}_n(B\times \textsf{X})(x,x') = M_n^f(B)(x), \qquad \check{M}_n(\textsf{X}\times B)(x,x') = M_n^c(B)(x').\end{equation}

2.3. Example

The following example is from [Reference Jasra, Kamatani, Law and Zhou22] and there is some overlap with the presentation in that article. We start with a diffusion process

(4) \begin{eqnarray}dZ_t & \,=\, & a(Z_t)dt + b(Z_t)dW_t\end{eqnarray}

with $Z_t\in\mathbb{R}^d=\textsf{X}$ , $a\,:\,\mathbb{R}^d\rightarrow\mathbb{R}^d$ (jth element denoted $a^j$ ), $b\,:\,\mathbb{R}^d\rightarrow\mathbb{R}^{d\times d}$ ((j, k)th element denoted $b^{j,k}$ ), $t\geq 0$ , and $\{W_t\}_{t\geq 0}$ a $d-$ dimensional Brownian motion. The following assumption is made and is referred to as (D). We set $Z_0=x^*\in\textsf{X}$ .

The coefficients $a^j, b^{j,k}$ belong to $\mathcal{C}^2(\textsf{X})$ , for $j,k= 1,\ldots, d$ . Also, a and b satisfy the following properties:

1. (i) the uniform ellipticity property: $b(z)b(z)^T$ is uniformly positive definite;

2. (ii) the globally Lipschitz property: there is a $C>0$ such that $|a^j(z)-a^j(y)|+|b^{j,k}(z)-b^{j,k}(y)| \leq C |z-y|$ for all $(z,y)\in \textsf{X}\times\textsf{X}$ , $(j,k)\in\{1,\dots,d\}^2$ .

The data are observed at regular unit time-intervals (i.e. in discrete time) $y_1,y_2,\dots$ , $y_k \in \textsf{Y}$ . It is assumed that conditional on $Z_{k}$ , $Y_k$ is independent of all other random variables with density $G(z_{k},y_k)$ . Let M(z, dy) be the transition of the diffusion process (over unit time) and consider a discrete-time Markov chain $X_0,X_1,\dots$ with initial distribution $M(x^*,\cdot)$ and transition M(x, dy). Here we are creating a discrete-time Markov chain that corresponds to the discrete-time skeleton of the diffusion process at a time lag of 1. Now we write $G_{k}(x_{k})$ instead of $G(x_{k},y_{k+1})$ . Then we define, for $B\in\mathcal{X}$ ,

\[\gamma_n(B) \,:\!= \int_{\textsf{X}^{n+1}}\mathbb{I}_B(x_n)\Bigg(\prod_{p=0}^{n-1} G_{p}(x_{p})\Bigg)M(x^*,dx_{0})\prod_{p=1}^n M(x_{p-1},dx_{p}).\]

The predictor is $\eta_n(B)=\gamma_n(B)/\gamma_n(1)$ , which corresponds to the distribution associated to $Z_{n+1}|y_1,\dots,y_n$ .

In many applications, one must time-discretize the diffusion to use the model in practice. We suppose an Euler discretization with discretization $\Delta_l=2^{-l}$ , $l\geq 0$ , and write the associated transition kernel over unit time as $M^l(x,dy)$ . Note that, in practice, one may not be able to compute the density of the kernel, as it is a composition of $\Delta_l^{-1}-1$ Gaussians; however, one can certainly sample from $M^l$ in most cases. Hence we are interested in the Feynman–Kac model for $B\in\mathcal{X}$ ,

\[\gamma_n^l(B) \,:\!= \int_{\textsf{X}^{n+1}}\mathbb{I}_B(x_n)\Bigg(\prod_{p=0}^{n-1} G_{p}(x_{p})\Bigg)M^l(x^*,dx_{0})\prod_{p=1}^n M^l(x_{p-1},dx_{p}),\]

with associated predictor $\eta_n^l(B)=\gamma_n^l(B)/\gamma_n^l(1)$ .

Below, we will explain why one may wish to compute, for $l\geq 1$ , $\eta_n^l(\varphi)-\eta_n^{l-1}(\varphi)$ , $\varphi\in\mathcal{B}_b(\textsf{X})$ . That is, f as used above relates to the predictor associated to the discretization $\Delta_l$ and c to the predictor with discretization $\Delta_{l-1}$ . A natural coupling of $M_n^f=M^l$ and $M_n^c=M^{l-1}$ exists (e.g. [Reference Giles13]) so one also has a given $\check{\eta}_0$ and $\check{M}_n$ .

Before continuing, we note some results which will help in our discussion below. As established in [Reference Jasra, Kamatani, Law and Zhou22, Eq. (32)], for $C<+\infty$ one has

(5) \begin{equation}\sup_{\mathcal{A}}\sup_{x\in\textsf{X}}|M_n^f(\varphi)(x)-M_n^c(\varphi)(x)| \leq C \Delta_l,\end{equation}

where $\mathcal{A}=\{\varphi\in\mathcal{B}_b(\textsf{X})\cap\textrm{Lip}_{\textsf{d}_1}(\textsf{X})\,:\, \|\varphi\|\leq 1\}$ , with $\textsf{d}_1$ the $L_1-$ norm. In addition, [Reference Jasra, Kamatani, Law and Zhou22, Proposition D.1] states that for $p>0$ , $C<+\infty$ , and any $(x,y)\in\textsf{X}\times\textsf{X}$ ,

(6) \begin{equation}\int_{\textsf{X}\times\textsf{X}}\|u-v\|^p \check{M}_n((x,y),d(u,v)) \leq C(|x-y| + \Delta_l^{1/2})^p.\end{equation}

When $p=2$ ( $\|u-v\|$ is the $L_2-$ norm), the term $\Delta_l$ is the so-called forward strong error rate. In the proof of [Reference Jasra, Kamatani, Law and Zhou22, Theorem 4.3], it is shown that for any $n\geq 0$ there exists a $C<+\infty$ such that for $\varphi\in\mathcal{A}$ , $l\geq 0$ , one has

(7) \begin{equation}|\eta_n(\varphi)-\eta_n^l(\varphi)| \leq C\Delta_l.\end{equation}

Note that this bound can be deduced from (5), but that C may depend on n; this latter point is ignored for now.

2.3.1. Multilevel Monte Carlo

Suppose one can exactly sample from $\eta_n^l$ for any $l,n\geq 0$ . The Monte Carlo estimate of $\eta_n^l(\varphi)$ , with $\varphi\in\mathcal{A}$ , is then of course $\frac{1}{N}\sum_{i=1}^N\varphi(x_n^i)$ , where the $X_n^i$ are independent and identically distributed (i.i.d.) from $\eta_n^l$ . One has the mean square error (MSE)

\[\mathbb{E}\Bigg[\Bigg(\frac{1}{N}\sum_{i=1}^N\varphi(x_n^i) - \eta_n(\varphi)\Bigg)^2\Bigg] = \frac{\mathbb{V}\textrm{ar}_{\eta_n^l}[\varphi(X)]}{N} + |\eta_n(\varphi)-\eta_n^l(\varphi)|^2,\]

where $\mathbb{V}\textrm{ar}_{\eta_n^l}[\varphi(X)]$ is the variance of $\varphi(X)$ with respect to $\eta_n^l$ . Then, for $\varepsilon>0$ given, to target an MSE of $\mathcal{O}(\varepsilon^2)$ , by (7), one chooses $l=\mathcal{O}(|\log(\varepsilon)|)$ (as one sets $\Delta_l^2=\varepsilon^2$ ). Then one must choose $N=\mathcal{O}(\varepsilon^{-2})$ , and we suppose the cost of simulating one sample is $\mathcal{O}(\Delta_l^{-1})$ (see [Reference Jasra, Kamatani, Law and Zhou22] for a justification), again ignoring n. Then the cost of achieving an MSE of $\mathcal{O}(\varepsilon^2)$ is $\mathcal{O}(\varepsilon^{-3})$ .

Now, for any $L\geq 1$ , one has the multilevel identity ([Reference Giles13], [Reference Heinrich15])

\[\eta_n^L(\varphi) = \eta_n^0(\varphi) + \sum_{l=1}^{L}\{[\eta_n^l-\eta_n^{l-1}](\varphi)\}.\]

Suppose that, for $l\geq 1$ , it is possible to exactly sample a coupling of $(\eta_n^l,\eta_n^{l-1})$ —let us denote it $\check{\eta}_n^l$ —so that one has

\[\int_{\textsf{X}\times\textsf{X}} \|x-y\|^2 \check{\eta}_n^l(d(x,y)) \leq C\Delta_l,\]

and the cost of such a simulation is $\mathcal{O}(\Delta_l^{-1})$ . The rate $\Delta_l$ has been taken from the strong error rate in (6). Now, to estimate $[\eta_n^l-\eta_n^{l-1}](\varphi)$ , one simulates

\[((X_n^{1,l},X_n^{1,l-1}),\dots, (X_n^{N_l,l},X_n^{N_l,l-1}))\]

i.i.d. from $\check{\eta}_n^l$ for some $N_l\geq 1$ , and the approximation is

\[\frac{1}{N_l}\sum_{l=1}^{N_l}\{\varphi(x_n^{i,l})-\varphi(x_n^{i,l-1})\}.\]

For $1\leq l\leq L$ this is repeated independently for estimating $[\eta_n^l-\eta_n^{l-1}](\varphi)$ , and the case $l=0$ is performed, independently, using the Monte Carlo method above. One can then easily show that the MSE of the estimate of $\eta_n^L$ is bounded above by

\[\frac{\mathbb{V}\textrm{ar}_{\eta_n^l}[\varphi(X)]}{N_0} + C \sum_{l=1}^L\frac{\Delta_l}{N_l} + |\eta_n(\varphi)-\eta_n^l(\varphi)|^2\]

(recall $\varphi\in\mathcal{A}$ ). Then, for $\varepsilon>0$ given, to target an MSE of $\mathcal{O}(\varepsilon^2)$ , set $L=\mathcal{O}(|\log(\varepsilon)|)$ and $N_l=\mathcal{O}(\varepsilon^{-2}|\log(\varepsilon)|\Delta_l)$ (see [Reference Giles13]); then the MSE is $\mathcal{O}(\varepsilon^2)$ . The overall computational effort is $\mathcal{O}(\sum_{l=0}^L N_l \Delta_l^{-1})=\mathcal{O}(\varepsilon^{-2}|\log(\varepsilon)|^2)$ ; if $\varepsilon$ is suitably small this is a significant reduction in computational effort relative to the MC method above.

The main point here is that, at present, there are no computational methods to perform the exact simulation mentioned, and thus we focus on particle filters, which have been developed in the literature. The estimates (variance/cost) above typically depend on n, and our objective is to consider whether the variance of PF estimates can be $\mathcal{O}(\Delta_l)$ uniformly in time, so that the ML gain is retained uniformly in time. We focus on the asymptotic variance in the CLT (versus the finite sample variance) as this is more straightforward to deal with.

3.1. Independent pair resampling

The first procedure we consider is as follows. Let $n\geq 1$ , $B\in\mathcal{X}\times\mathcal{X}$ , and $\mu\in\mathscr{P}(\textsf{X}\times\textsf{X})$ , and define the probability measure

\[\check{\Phi}_n^I(\mu)(B) = \frac{\mu([G_{n-1}\otimes G_{n-1}]\check{M}_n(B)) }{\mu(G_{n-1}\otimes G_{n-1})}.\]

Set $u_n=(x_n^f,x_n^c)\in\textsf{X}\times\textsf{X}$ , then consider the joint probability measure on $(\textsf{X}\times\textsf{X})^{n+1}$ given by

(8) \begin{equation}\mathbb{P}(d(u_0,\dots,u_n)) = \check{\eta}_0(du_0) \prod_{p=1}^n \check{\Phi}_p^I(\eta_{p-1}^f\otimes\eta_{p-1}^c)(du_p).\end{equation}

Noting (3), it is easily checked that the marginal of $x_n^f$ (resp. $x_n^c$ ) is $\eta_n^f$ (resp. $\eta_n^c$ ). Denote by $\check{\eta}_n^I$ the marginal of $u_n$ induced by this joint probability measure. Thus, if one could sample a trajectory of $(u_0,\dots,u_n)$ from $\mathbb{P}(d(u_0,\dots,u_n))$ one could easily approximate quantities such as (1) or (2) using Monte Carlo methods. In most practical applications of interest, this is not possible.

The particle approximation of (8) is taken as

\[\mathbb{P}(d(u_0^{1:N},\dots,u_n^{1:N})) = \Bigg(\prod_{i=1}^N \check{\eta}_0(du_0^i)\Bigg)\Bigg(\prod_{p=1}^n\prod_{i=1}^N\check{\Phi}_p^I(\eta_{p-1}^{N,f}\otimes\eta_{p-1}^{N,c})(du_p^i)\Bigg),\]

where for $p\geq 1$ , $s\in\{\,f,c\}$ ,

\[\eta_{p-1}^{N,s}(dx) = \frac{1}{N}\sum_{i=1}^N \delta_{x_{p-1}^{i,s}}(dx).\]

The IRCPF algorithm which simulates from $\mathbb{P}(d(u_0^{1:N},\dots))$ is presented in Algorithm 1. The key point is that in this algorithm, the resampled indices for a pair of particles are generated (conditionally) independently. Set

\[\check{\eta}_{p}^{N,I}(du) = \frac{1}{N}\sum_{i=1}^N \delta_{u_{p}^{i}}(du).\]

As has been mentioned by many authors (e.g. [Reference Sen, Thiery and Jasra28]), one does not expect this procedure to be effective, in the sense that $\check{\eta}_n^I$ would not provide an appropriate dependence between $(\eta_n^f,\eta_n^c)$ .

Algorithm 1:

The IRCPF Algorithm.

3.2. Maximally coupled resampling

Let $n\geq 1$ , $B\in\mathcal{X}\times\mathcal{X}$ , and $\mu\in\mathscr{P}(\textsf{X}\times\textsf{X})$ , and define the probability measure

\begin{eqnarray*}\check{\Phi}_n^M(\mu)(B) & \,=\, & \mu\Big(\{F_{n-1,\mu,f} \wedge F_{n-1,\mu,c}\} \check{M}_n(B)\Big)+ \Big(1-\mu\Big(\{F_{n-1,\mu,f} \wedge F_{n-1,\mu,c}\}\Big)\Big) \\ & &\times (\mu\otimes\mu)\Big(\Big\{\overline{F}_{n-1,\mu,f}\otimes \overline{F}_{n-1,\mu,c}\Big\}\bar{M}_n(B)\Big),\end{eqnarray*}

where for $(x,y)\in\textsf{X}\times\textsf{X}$ ,

\begin{eqnarray*}\overline{F}_{n-1,\mu,f}(x,y) & \,=\, & \frac{F_{n-1,\mu,f}(x,y)-\{F_{n-1,\mu,f}(x,y)\wedge F_{n-1,\mu,c}(x,y)\}}{\mu(F_{n-1,\mu,f}-\{F_{n-1,\mu,f}\wedge F_{n-1,\mu,c}\})}, \\\overline{F}_{n-1,\mu,c}(x,y) & \,=\, & \frac{F_{n-1,\mu,c}(x,y)-\{F_{n-1,\mu,f}(x,y)\wedge F_{n-1,\mu,c}(x,y)\}}{\mu(F_{n-1,\mu,c}-\{F_{n-1,\mu,f}\wedge F_{n-1,\mu,c}\})}, \\F_{n-1,\mu,f}(x,y) & \,=\, & \check{G}_{n-1,\mu,f}(x)\otimes1, \\F_{n-1,\mu,c}(x,y) & \,=\, & 1\otimes \check{G}_{n-1,\mu,c}(y), \\\check{G}_{n-1,\mu,f}(x) & \,=\, & \frac{G_{n-1}(x)}{\mu(G_{n-1}\otimes 1)}, \\\check{G}_{n-1,\mu,c}(y) & \,=\, & \frac{G_{n-1}(y)}{\mu(1\otimes G_{n-1})},\end{eqnarray*}

and for $((x,y),(u,v))\in\textsf{X}^2\times \textsf{X}^2$ and $B\in\mathcal{X}\times\mathcal{X}$ ,

\[\bar{M}_n(B)((x,y),(u,v)) = \check{M}_n(B)(x,v).\]

Now set $\check{\eta}_0^M=\check{\eta}_0$ ; we define, recursively, for any $n\geq 1$ , $B\in\mathcal{X}\times\mathcal{X}$ ,

\[\check{\eta}_n^M(B) = \check{\Phi}_n^M(\check{\eta}_{n-1}^M)(B).\]

Note that by [Reference Jasra, Kamatani, Law and Zhou22, Proposition A.1] we have, for $B\in\mathcal{X}$ ,

(9) \begin{equation}\check{\eta}_n^M(B\times\textsf{X}) = \eta_n^f(B)\quad\textrm{and}\quad\check{\eta}_n^M(\textsf{X}\times B) = \eta_n^c(B).\end{equation}

To generate an MCIRCPF, we consider the following probability measure:

\[\mathbb{P}(d(u_0^{1:N},\dots,u_n^{1:N})) = \Bigg(\prod_{i=1}^N \check{\eta}_0(du_0^i)\Bigg)\Bigg(\prod_{p=1}^n\prod_{i=1}^N\check{\Phi}_p^M(\check{\eta}_{n-1}^{N,M})(du_p^i)\Bigg),\]

where for $p\geq 1$ ,

\[\check{\eta}_{p-1}^{N,M}(du) = \frac{1}{N}\sum_{i=1}^N \delta_{u_{p}^{i}}(du),\]

and for $s\in\{\,f,c\}$ we set

\[\eta_{p-1}^{N,s}(dx) = \frac{1}{N}\sum_{i=1}^N \delta_{x_{p-1}^{i,s}}(dx).\]

The way in which one can generate from $\mathbb{P}(d(u_0^{1:N},\dots))$ is given in Algorithm 2. To explain further, the simulation from $\check{\Phi}_p^M(\check{\eta}_{p-1}^{N,M})$ consists of whether the resampled indices for the particle pair are equal (second bullet in Algorithm 2) or different (third bullet in Algorithm 2). This procedure is as in [Reference Chopin and Singh6] and adopted in, for instance, [Reference Jacob, Lindsten and Schön17] and [Reference Jasra, Kamatani, Law and Zhou22]. The idea is to provide a local optimality procedure with respect to the resampling operation. This is in the sense that for any fixed N, and conditional on the information generated so far, one will maximize the probability that the resampled indices are equal.

For the MCIRCPF, we present a preliminary result which will prove to be of interest. The following assumptions are used; note that in (A3) there is a metric $\textsf{d}$ on $\textsf{X}\times\textsf{X}$ implicit in the assumption.

1. (A1) For every $n\geq 0$ , $G_n\in \mathcal{C}_b(\textsf{X})$ .

2. (A2) For every $n\geq 1$ , $\check{M}_n$ is Feller.

3. (A3) $\textsf{X}$ is a locally compact and separable metric space.

These assumptions are adopted because of the complexity of the operator $\check{\Phi}_n^M$ . As can be seen in Appendix C, where the proofs for the following result are given, it is nontrivial to work with $\check{\Phi}_n^M$ . To weaken these assumptions would lead to further calculations, which would essentially confirm the same result.

Algorithm 2:

The MCIRCPF Algorithm.

Theorem 3.1. Assume (A1)–(A3). Then for any $\varphi\in\mathcal{C}_b(\textsf{X}\times\textsf{X})$ , $n\geq 0$ , we have

\[\check{\eta}_{n}^{N,M}(\varphi) \rightarrow_{\mathbb{P}} \check{\eta}_{n}^{M}(\varphi).\]

What is interesting here is that Theorem 3.1 verifies what is expected, given the construction above: on the product space $\textsf{X}\times\textsf{X}$ the MCIRCPF approximates the target $\check{\eta}_{n}^{M}$ . As noted in (9), $\check{\eta}_{n}^{M}$ is a coupling of $(\eta_n^f,\eta_n^c)$ , but it is not the actual maximal coupling of $(\eta_n^f,\eta_n^c)$ ; this can be easily checked to be the case in most problems of practical interest. As is well known, the maximal coupling will maximize the probability that two random variables are equal, with specified marginals, and is the optimal coupling of two probability measures with respect to the $L_0-$ Wasserstein distance. If this former property is desirable from a practical perspective, then the above algorithm should not be used. The maximally coupled resampling operation is, for a finite number of samples (particles), the optimal (in the above sense) way to couple the resampling operation, but may not lead to large sample ‘good’ couplings. This is manifested in [Reference Jasra, Kamatani, Law and Zhou22], where the forward error rate (6) is lost for the diffusion problem in Section 2.3. As mentioned above, the limit is a coupling of $(\eta_n^f,\eta_n^c)$ , but in general there is no reason to suspect that it is optimal in any sense.

3.3. Maximal coupling

We now present an algorithm which can sample (in the limit) from the maximal coupling of $(\eta_n^f,\eta_n^c)$ . We will assume that for $s\in\{\,f,c\}$ the Markov kernels $M_n^s$ , as well as $\eta_0^s,$ admit a density with respect to a $\sigma-$ finite measure dx. The densities are denoted $M_n^s$ and $\eta_0^s$ , $s\in\{\,f,c\}$ , and we assume that the densities can be evaluated numerically. To remove this latter requirement is left to future work.

Let $n\geq 1$ , $B\in\mathcal{X}\times\mathcal{X}$ , and $(\mu,\nu)\in\mathscr{P}(\textsf{X})\times\mathscr{P}(\textsf{X})$ , and define the probability measure

\begin{eqnarray*}\check{\Phi}_n^{C}(\mu,\nu)(B) & \,=\, & \int_{\textsf{X}\times\textsf{X}}\mathbb{I}_B(x,y) \bigg[\int_{\textsf{X}}F_{n-1,\mu,f}(u)\wedge F_{n-1,\nu,c}(u) \delta_{\{u,u\}}(d(x,y)) du \\ & & +\,\frac{1}{1-\int_{\textsf{X}}F_{n-1,\mu,f}(u)\wedge F_{n-1,\nu,c}(u)du}\overline{F}_{n-1,\mu,\nu,f}(x)\overline{F}_{n-1,\nu,\mu,c}(y) dxdy\bigg],\end{eqnarray*}

where, for $x\in\textsf{X}$ ,

\begin{eqnarray*}\overline{F}_{n-1,\mu,\nu,f}(x) & \,=\, & F_{n-1,\mu,f}(x) - F_{n-1,\mu,f}(x)\wedge F_{n-1,\nu,c}(x),\\\overline{F}_{n-1,\nu,\mu,c}(x) & \,=\, & F_{n-1,\nu,c}(x) - F_{n-1,\mu,f}(x)\wedge F_{n-1,\nu,c}(x),\\ F_{n-1,\mu,f}(x) & \,=\, & \frac{\mu\big(G_{n-1}M_n^f(\cdot,x)\big)}{\mu(G_{n-1})},\\ F_{n-1,\nu,c}(x) & \,=\, & \frac{\nu(G_{n-1}M_n^c(\cdot,x))}{\nu(G_{n-1})}.\end{eqnarray*}

Now set $\check{\eta}_0^C$ as the maximal coupling of $(\eta_0^f,\eta_0^c)$ , and for $B\in\mathcal{X}\times\mathcal{X}$ set

\[\check{\eta}_n^C(B) = \check{\Phi}_n^{C}(\eta_{n-1}^f,\eta_{n-1}^c)(B).\]

We have, for $B\in\mathcal{X}$ ,

\[\check{\eta}_n^C(B\times\textsf{X}) = \eta_n^f(B)\quad\textrm{and}\quad\check{\eta}_n^C(\textsf{X}\times B) = \eta_n^c(B).\]

Moreover, $\check{\eta}_n^C$ is the maximal coupling of $(\eta_n^f,\eta_n^c)$ . To see why this is the case, we note that our assumptions imply that $(\eta_n^f,\eta_n^c)$ have densities with respect to a $\sigma-$ finite measure and that these densities are simply $(F_{n-1,\eta_{n-1}^f,f},F_{n-1,\eta_{n-1}^c,c})$ ; hence the expression $\check{\Phi}_n^{C}(\eta_{n-1}^f,\eta_{n-1}^c)(\cdot)$ corresponds exactly to the definition of the maximal coupling of $(\eta_n^f,\eta_n^c)$ .

To generate an MCPF, we consider the following probability measure:

\[\mathbb{P}(d(u_0^{1:N},\dots,u_n^{1:N})) = \Bigg(\prod_{i=1}^N \check{\eta}_0^C(du_0^i)\Bigg)\Bigg(\prod_{p=1}^n\prod_{i=1}^N\check{\Phi}_p^C(\eta_{p-1}^{N,f},\eta_{p-1}^{N,c})(du_p^i)\Bigg),\]

where for $p\geq 1$ , $s\in\{\,f,c\}$ ,

\[\eta_{p-1}^{N,s}(dx) = \frac{1}{N}\sum_{i=1}^N \delta_{x_{p-1}^{i,s}}(dx).\]

For $\mu\in\mathscr{P}(\textsf{X})$ , $s\in\{\,f,c\}$ , $B\in\mathcal{X}$ , $n\geq 1$ , we set

\[\Phi_n^s(\mu)(B) = \frac{\mu(G_{n-1}M_n^s(B))}{\mu(G_{n-1})}.\]

In Algorithm 3 we present how one can simulate from $\mathbb{P}(d(u_0^{1:N},\dots))$ . We remark that as $M_n^s$ , $n\geq 1$ , and $\eta_0^s$ , $s\in\{\,f,c\}$ , can be evaluated numerically, we can sample from $\check{\eta}_0^C$ and $\check{\Phi}_n^{C}(\eta_{n-1}^{N,f},\eta_{n-1}^{N,c})$ , the latter of which is the maximal coupling of $\Phi_n^f(\eta_{n-1}^{N,f})$ and $\Phi_n^c(\eta_{n-1}^{N,c})$ , using the algorithm in [Reference Thorisson29], which is presented in the bullet points in Algorithm 3. The algorithm in [Reference Thorisson29] is a rejection sampler, which would add a random running time element per time-step, which may not be desirable for some applications. As noted in [Reference Jacob, O’Leary and Atachde18], it can be shown that the expected running time is at most 2, but that the variance of the running time depends upon the closeness, in terms of the total variation distance, of the two probability measures that one is trying to couple. In general, in the limiting case, we expect that the total variation distance between $(\eta_n^f,\eta_n^c)$ will be small (in the diffusion case, under certain assumptions, it can be $\mathcal{O}(\Delta_l)$ ) and hence that the rejection sampler will perform well.

Algorithm 3:

The MCPF Algorithm.

3.4. Wasserstein coupled resampling

We describe the WCPF used in [Reference Jasra, Ballesio, von Schwerin and Tempone20]. For this case we restrict our attention to the case that $\textsf{X}=\mathbb{R}$ . The reason for this simplification is that the approach is based upon an optimal coupling result which appears to be explicit only in the one-dimensional case. In the approach to be considered, one might consider the multi-dimensional case using, for instance Hilbert space-filling curves, but we have not tested such a method empirically and it is unclear how useful it is in practice. It is explicitly assumed that the cumulative distribution function (CDF) and its inverse associated to the probability

\[\overline{\eta}_n^s(B) = \frac{\eta_n^s(G_n\mathbb{I}_B)}{\eta_n^s(G_n)},\]

where $s\in\{\,f,c\}$ , $n\geq 0$ , and $B\in\mathcal{X}$ , exist and are continuous functions. We denote the CDF (resp. inverse) of $\overline{\eta}_n^s$ as $F_{\overline{\eta}_n^s}$ (resp. $F_{\overline{\eta}_n^s}^{-1}$ ). In general we write probability measures on $\textsf{X}$ for which the CDF and inverse are well-defined as $\mathscr{P}_F(\textsf{X})$ with the associated CDF $F_{\mu}$ . Let $n\geq 1$ , $B\in\mathcal{X}\times\mathcal{X}$ , and $\mu,\nu\in\mathscr{P}_F(\textsf{X})$ , and define the probability measure

\[\check{\Phi}_n^W(\mu,\nu)(B) = \int_{\textsf{X}\times\textsf{X}} \bigg(\int_{0}^1 \delta_{\{F_{\mu}^{-1}(w),F_{\nu}^{-1}(w)\}}(du) dw \bigg) \check{M}_n(B)(u).\]

Consider the joint probability measure on $(\textsf{X}\times\textsf{X})^{n+1}$ given by

\[\mathbb{P}(d(u_0,\dots,u_n)) = \check{\eta}_0(du_0) \prod_{p=1}^n \check{\Phi}_p^W(\overline{\eta}_{p-1}^f,\overline{\eta}_{p-1}^c)(du_p).\]

It is easily checked that the marginal of $x_n^f$ (resp. $x_n^c$ ) is $\eta_n^f$ (resp. $\eta_n^c$ ). Denote by $\check{\eta}_n^W$ the marginal of $u_n$ induced by this joint probability measure.

To generate a WCPF, we consider the following:

\[\mathbb{P}(d(u_0^{1:N},\dots,u_n^{1:N})) = \Bigg(\prod_{i=1}^N \check{\eta}_0(du_0^i)\Bigg)\Bigg(\prod_{p=1}^n\prod_{i=1}^N\check{\Phi}_p^W(\overline{\eta}_{p-1}^{N,f},\overline{\eta}_{p-1}^{N,c})(du_p^i)\Bigg),\]

where for $p\geq 1$ , $s\in\{\,f,c\}$ ,

\[\overline{\eta}_{p-1}^{N,s}(dx) = \sum_{i=1}^N \frac{G_{p-1}\big(x_{p-1}^{i,s}\big)}{\sum_{j=1}^N G_{p-1}\big(x_{p-1}^{j,s}\big)}\delta_{x_{p-1}^{i,s}}(dx).\]

As before, for $p\geq 1$ , $s\in\{\,f,c\}$ ,

\[\eta_{p-1}^{N,s}(dx) = \frac{1}{N}\sum_{i=1}^N \delta_{x_{p-1}^{i,s}}(dx),\]

and for $p\geq 0$ ,

\[\check{\eta}_{p}^{N,W}(du) = \frac{1}{N}\sum_{i=1}^N \delta_{u_{p}^{i}}(du).\]

In Algorithm 4 we present how one can simulate from $\mathbb{P}(d(u_0^{1:N},\dots))$ . It should be noted that whilst the sorting operation in Algorithm 4 is $\mathcal{O}(N\log(N))$ in the worst case, it is seen in [Reference Jasra, Ballesio, von Schwerin and Tempone20] that this step is often $\mathcal{O}(N)$ in practice and is less expensive than the sampling of $\check{M}_p$ in the diffusion case.

Algorithm 4:

The WCPF Algorithm.

4.1. Central limit theorems

Define the sequence of nonnegative kernels $\{Q_n^s\}_{n\geq 1}$ , $s\in\{\,f,c\}$ , by $Q_n^s(x,dy) = G_{n-1}(x) M_n^s(x,dy)$ . For $B\in\mathcal{X}$ , $x_p\in\textsf{X}$ , define

\[Q_{p,n}^s(B)(x_p) = \int_{\textsf{X}^{n-p}} \mathbb{I}_B(x_n) \prod_{q=p}^{n-1} Q_{q+1}^s(x_q,dx_{q+1})\]

for $0\leq p<n$ ; in the case $p=n$ , let $Q_{p,n}^s$ be the identity operator. Now, for $0\leq p<n$ , $s\in\{\,f,c\}$ , $B\in\mathcal{X}$ , $x_p\in\textsf{X}$ , define

\[D_{p,n}^s(B)(x_p) = \frac{Q_{p,n}^s(\mathbb{I}_B - \eta_n^s(B))}{\eta_p^s\big(Q_{p,n}^s(1)\big)};\]

in the case $p=n$ , $D_{p,n}^s(B)(x)=\mathbb{I}_B-\eta_n^s(B)$ .

We then have our CLTs, where the proofs are given in Appendix A. Recall that for the MCPF $M_n^s$ is used as a notation for the kernel and density of $M_n^s$ .

For Wasserstein resampling, the following result is established in [Reference Jasra, Ballesio, von Schwerin and Tempone20].

Theorem 4.2. Suppose that $\check{M}_n$ , $M_n^s$ , $s\in\{\,f,c\}$ , are Feller for every $n\geq 1$ and $G_n\in\mathcal{C}_b(\textsf{X})$ for every $n\geq 0$ . Then for any $\varphi\in\mathcal{C}_b(\textsf{X})$ , $n\geq 0$ , we have

\[\sqrt{N}[\check{\eta}_{n}^{N,W}-\check{\eta}_{n}^{W}](\varphi\otimes 1 - 1 \otimes \varphi) \Rightarrow \mathcal{N}(0,\sigma_n^{2,W}(\varphi)),\]

where

\[\sigma_n^{2,W}(\varphi) = \sum_{p=0}^n \check{\eta}_p^W(\{(D_{p,n}^f(\varphi)\otimes 1- 1\otimes D_{p,n}^c(\varphi))\}^2).\]

Remark 4.1. One can also prove a multivariate CLT using the Cramér–Wold device. Consider $1\leq t <+\infty$ , $(\varphi_1,\dots,\varphi_t,\psi_1,\dots,\psi_t)\in\mathcal{B}_b(\textsf{X})^{2t}$ (or if required $\mathcal{C}_b(\textsf{X})^{2t}$ ). Consider the $t\times t$ positive definite and symmetric matrix $\Sigma_{n}^{s}(\varphi_{1:t},\psi_{1:t})$ , $s\in\{I,M,C,W\}$ , with (i,j)th entry denoted $\Sigma_{n,(ij)}^{s}(\varphi_{1:t},\psi_{1:t})$ . Using the Cramér–Wold device under the various assumptions of each the algorithms one can easily deduce that for each $s\in\{I,M,C,W\}$ ,

\[\sqrt{N}([\check{\eta}_{n}^{N,s}-\check{\eta}_{n}^{s}](\varphi_1\otimes 1 {-} 1 \otimes \psi_1),\dots,[\check{\eta}_{n}^{N,s}-\check{\eta}_{n}^{s}](\varphi_t\otimes 1 {-} 1 \otimes \psi_t) ) \Rightarrow \mathcal{N}_t(0,\Sigma_n^{s}(\varphi_{1:t},\psi_{1:t})),\]

where

\[\Sigma_{n,(ij)}^{s}(\varphi_{1:t},\psi_{1:t}) = \sum_{p=0}^n \check{\eta}_p^s([D_{p,n}^f(\varphi_i)\otimes 1- 1\otimes D_{p,n}^c(\psi_i)][D_{p,n}^f(\varphi_j)\otimes 1- 1\otimes D_{p,n}^c(\psi_j)]).\]

4.2. Asymptotic variance

We now consider bounding the asymptotic variance in the general case. The result below (Theorem 4.3), will allow one to understand when one can expect time-uniformly ‘close’ errors in approximations of $[\eta_n^f-\eta_n^c](\varphi)$ . The following assumptions are used, which are essentially those in [Reference Whiteley30] (see also [Reference Jasra19]) with some additional assumptions ((H6)–(H7) below) as we are treating a more delicate case than that of [Reference Whiteley30]. The assumptions can hold on unbounded state spaces, as will be the case for our applications.

1. (H1) There exist an unbounded $\tilde{V}\,:\,\textsf{X}\rightarrow[1,+\infty)$ and constants $\delta\in(0,1)$ and $\underline{d}\geq 1$ with the following properties. For each $d\in(\underline{d},+\infty)$ there exists a $b_d<+\infty$ such that $\forall x\in\textsf{X}$ and any $s\in\{\,f,c\}$ ,

   \[\sup_{n\geq 1} Q_n^s(e^{\tilde{V}})(x) \leq e^{(1-\delta)\tilde{V}(x) + b_d\mathbb{I}_{C_d}(x)},\]
   where $C_d=\{x\in\textsf{X}\,:\,\tilde{V}(x)\leq d\}$ .

2. (H2)

   1. 1. There exists a $C<+\infty$ such that for any $s\in\{\,f,c\}$ , $\eta_0^s(v)\leq C$ , with $v=e^{\tilde{V}}$ , with $\tilde{V}$ as in (H1).

   2. 2. For any $r>1$ there exists a $C<+\infty$ such that for any $s\in\{\,f,c\}$ , $\eta_0^s(C_r)^{-1}\leq C$ .

3. (H3) With $\underline{d}$ as in (H1), for each $d\in[\underline{d},+\infty)$ , $s\in\{\,f,c\}$ ,

   \[\int_{C_d} G_{n-1}(x)M_n^s(x,dy) > 0\ \forall x\in\textsf{X}, n\geq 1,\]
   and there exist $\tilde{\epsilon}_d^{-}>0$ , $\nu_d\in\mathcal{P}_v$ such that for each $A\in\mathcal{X}$ , $s\in\{\,f,c\}$ ,
   \[\inf_{n\geq 1} \int_{C_d\cap A} Q_n^s(x,dy) \geq \tilde{\epsilon}_d^{-}\nu_d(C_d\cap A),\ \forall x\in C_d,\]
   with $C_d$ as in (H1).

4. (H4) With $\underline{d}$ as in (H1), and $\tilde{\epsilon}_d^-$ , $\nu_d$ as in (H2), for each $d\in[\underline{d},+\infty)$ there exists $\tilde{\epsilon}_d^+\in[\tilde{\epsilon}_d^-,+\infty)$ such that for each $A\in\mathcal{X}$ , $s\in\{\,f,c\}$ ,

   \[\sup_{n\geq 1}\int_{C_d\cap A} Q_n^s(x,dy) \leq \tilde{\epsilon}_d^+\nu(C_d\cap A), \ \forall x\in C_d,\]
   with $C_d$ as in (H1).

5. (H5)

   1. 1.

      \[\sup_{n\geq 0}\sup_{x\in\textsf{X}}G_n(x) <+\infty.\]

   2. 2. For any $r>1$ there exists a $C<+\infty$ such that $\inf_{x\in C_r}G_0(x)\geq C$ , with $C_r$ as in (H1).

6. (H6) With $\underline{d}$ as in (H1), for each $d\in[\underline{d},+\infty)$ , we have for each $n\geq 0$ , $s\in\{\,f,c\}$ , that

   \[\frac{1}{G_n M_{n+1}^s(\mathbb{I}_{C_d})}\in\mathcal{L}_{v}(\textsf{X})\]
   and $\sup_{n\geq 0}\max_{s\in\{\,f,c\}}\|1/G_n M_{n+1}^s(\mathbb{I}_{C_d})\|_{v}<+\infty$ .

7. (H7) Let $\textsf{d}$ be a given metric on $\textsf{X}$ . For any $\xi\in(0,1]$ , there exists a $C<+\infty$ such that for $\varphi\in\mathcal{L}_{v^{\xi}}(\textsf{X})\cap\textrm{Lip}_{v^{\xi},\textsf{d}}(\textsf{X})$ , $(x,y)\in\textsf{X}\times\textsf{X}$ , $s\in\{\,f,c\}$ ,

   \[\sup_{n\geq 1}|Q_{n}^s(\varphi)(x)-Q_{n}^s(\varphi)(y)| \leq C\|\varphi\|_{v^{\xi}}\textsf{d}(x,y)[v(x)v(y)]^{\xi}.\]

Let $A=\{(x,y)\in\textsf{X}\times\textsf{X}\,:\,x\neq y\}$ . For $n\geq 1$ , $0\leq p \leq n$ , $s\in\{\,f,c\}$ , $x\in\textsf{X}$ , define

\[h_{p,n}^s(x) \,:\!= \frac{Q_{p,n}^s(1)(x)}{\eta_p^s\big(Q_{p,n}^s(1)\big)}\]

and, with $B\in\mathcal{X}$ ,

\[S_{p,n}^{f,c}(B)(x) \,:\!= \frac{Q_{p,n}^f(B)(x)}{Q_{p,n}^f(1)(x)} - \frac{Q_{p,n}^c(B)(x)}{Q_{p,n}^c(1)(x)}.\]

Set

\[B(n,f,c,\varphi,\xi) =\|\varphi\|\sum_{p=0}^{n-1}\rho^{n-p}\Big\{\|h_{p,n}^fS_{p,n}^{f,c}(\varphi)\|_{v^{\xi}} + |[\eta_n^f-\eta_n^c](\varphi)| + \|\varphi\|\|h_{p,n}^f-h_{p,n}^c\|_{v^{\xi}}\rho^{n-p}\Big\}.\]

Below is our main result, whose proof can be found in Appendix E.

Theorem 4.3. Assume (H1)–(H6). Then for any $\xi\in(0,1/4)$ there exist $\rho<1$ and $C<+\infty$ depending on the constants in (H1)–(H6) such that for any $\varphi\in\mathcal{B}_b(\textsf{X})$ , $n\geq 1$ , $s\in\{I,M,C,W\}$ , we have

\begin{align*}&\sigma^{2,s}_n(\varphi) \leq \check{\eta}_n^s\Big(\{\varphi\otimes 1-1\otimes\varphi - ([\eta_n^f-\eta_n^c](\varphi))\}^2\Big)\\ &\quad + C\Bigg\{B(n,f,c,\varphi,\xi) + \|\varphi\|^2\sum_{p=0}^{n-1}\rho^{n-p}\Big(\check{\eta}_p^{s}(\mathbb{I}_A(v\otimes v)^{2\xi})(\rho^{n-p}+1)\Big)\Bigg\}. \end{align*}

Additionally assume (H7). Then if $\textsf{d}^2\in\mathcal{L}_{(v\otimes v)^{\tilde{\xi}}}$ , $\tilde{\xi}\in(0,1/2)$ , for any $\xi\in(0,(1-2\tilde{\xi})/12))$ there exist $\rho<1$ and $C<+\infty$ depending on the constants in (H1)–(H7) such that for any $\varphi\in\mathcal{B}_b(\textsf{X})\cap\textrm{\textit{Lip}}_{\textsf{d}}(\textsf{X})$ , $n\geq 1$ , $s\in\{I,M,C,W\}$ , we have

\begin{align*}&\sigma^{2,s}_n(\varphi) \leq \check{\eta}_n^s\Big(\{\varphi\otimes 1-1\otimes\varphi - ([\eta_n^f-\eta_n^c](\varphi))\}^2\Big)\\ &\quad + C\Bigg\{B(n,f,c,\varphi,\xi) + \|\varphi\|^2\sum_{p=0}^{n-1}\rho^{n-p}\Big(\check{\eta}_p^{s}(\textsf{d}^2(v^{4\xi}\otimes v^{8\xi}))(\rho^{n-p}+1)\Big)\Bigg\}.\end{align*}

The main conclusion of this theorem is that one expects that the (variances of) approximations of $[\eta_n^f-\eta_n^c](\varphi)$ are uniformly ‘close’ in time if the following are satisfied:

1. 1. The quantities

   \[\check{\eta}_n^s\Big(\{\varphi\otimes 1-1\otimes\varphi - ([\eta_n^f-\eta_n^c](\varphi))\}^2\Big)\quad\textrm{and}\quad|[\eta_n^f-\eta_n^c](\varphi)|\]
   are small uniformly in time.

2. 2. The differences

   \[ \|h_{p,n}^fS_{p,n}^{f,c}(\varphi)\|_{v^{\xi}}\quad\textrm{and}\quad\|h_{p,n}^f-h_{p,n}^c\|_{v^{\xi}}\]
   are small at a rate which decays more slowly than exponentially in time.

3. 3. The quantities

   \[\check{\eta}_p^{s}(\mathbb{I}_A(v\otimes v)^{2\xi}) \quad\textrm{or}\quad\check{\eta}_p^{s}(\textsf{d}^2(v^{4\xi}\otimes v^{8\xi}))\]
   are small at a rate which decays more slowly than exponentially in time.

We note that in 1, strictly speaking, one could have $|[\eta_n^f-\eta_n^c](\varphi)|$ close at a rate which decays more slowly than exponentially in time, but as one requires the time-uniform closeness of $\check{\eta}_n^s(\{\varphi\otimes 1-1\otimes\varphi - ([\eta_n^f-\eta_n^c](\varphi))\}^2)$ , one expects this property for the former. The use of A and $\textsf{d}$ are linked to the coupling properties associated to $\check{\eta}_p^{s}$ , as we consider in the next section.

5.1. A general result

For $x\in\textsf{X}$ , define the set

\[\textsf{B}_l(x) \,:\!= \{y\in\textsf{X}\,:\,\|y-x\|>2\Delta_l\},\]

where $\|\cdot\|$ is the $L_2-$ norm. We write the transition densities with respect to Lebesgue measure dy of the diffusion as M(x, y) and the Euler approximation as $M^l(x,y)$ , $(x,y)\in\textsf{X}\times\textsf{X}$ . We note the following two equations quoted in [Reference Del Moral, Jacod and Protter11]: there exist $0<C,C'<+\infty$ (which depend on the drift and diffusion coefficients of (4)) such that for any $l\geq 0$ ,

(10) \begin{eqnarray}M(x,y) + M^l(x,y) & \,\leq\, & C\exp\{-C'\|y-x\|^2\}\quad\forall\ (x,y)\in\textsf{X}\times\textsf{X}, \end{eqnarray}

(11) \begin{eqnarray}|M(x,y) - M^l(x,y)| &\, \leq\, & C\Delta_l\exp\{-C'\|y-x\|^2\}\quad\forall\ x\in \textsf{X},y\in \textsf{B}_l(x).\end{eqnarray}

We add an additional assumption:

1. (H8) For C $^{\prime}$ as in (10)–(11), we have the following:

   1. 1. For any $\xi\in(0,1/2)$ , there exists a $C<+\infty$ such that for any $x\in\textsf{X}$ , $l\geq 0$ ,

      \[\int_{\textsf{B}_l(x)^c}v(y)^{\xi}dy \leq C v(x)^{2\xi}\int_{\textsf{B}_l(x)^c}dy.\]

   2. 2. For any $\xi\in(0,1/2)$ , there exists a $C<+\infty$ such that for any $x\in\textsf{X}$ , $l\geq 0$ ,

      \[\int_{\textsf{B}_l(x)}v(y)^{\xi}\exp\{-C'\|y-x\|^2\}dy \leq C v(x)^{2\xi}.\]

We have the following result, whose proof is in Appendix F, Section F.1.

Proposition 5.1. Assume (H1)–(H6), (H8). Then for any $(\xi,\hat{\xi})\in(0,1/8)\times(0,1/2)$ there exists a $C<+\infty$ depending on the constants in (H1)–(H6) and (H8) such that for any $\varphi\in\mathcal{B}_b(\textsf{X})$ , $n\geq 1$ , $1\leq l \leq L$ , we have

\[B(n,l,l-1,\varphi,\xi) \leq C\|\varphi\|^2\Delta_lv(x^*)^{2\xi+\hat{\xi}}.\]

The main point of this result is that now the time-uniform coupling ability of each algorithm will now rely on the properties of the limiting coupling $\check{\eta}_n^s$ .

5.2. IRCPF

We consider the term

\[\check{\eta}_n^I\big(\{\varphi\otimes 1-1\otimes\varphi - ([\eta_n^f-\eta_n^c](\varphi))\}^2\big)\]

in the upper bound of Theorem 4.3. Let us suppose that $\varphi\in\mathcal{A}$ , where $\mathcal{A}$ is defined below (5). One has that

\begin{align*}&\check{\eta}_n^I\big(\big\{\varphi\otimes 1-1\otimes\varphi - \big(\big[\eta_n^f-\eta_n^c\big](\varphi)\big)\big\}^2\big)\\[3pt] &\quad = \frac{\big(\eta_{n-1}^f\otimes\eta_{n-1}^c\big)\big([G_{n-1}\otimes G_{n-1}]\check{M}_n\big([\varphi\otimes1-1\otimes\varphi]^2\big)\big) }{\big(\eta_{n-1}^f\otimes\eta_{n-1}^c\big)(G_{n-1}\otimes G_{n-1})}- [\eta_n^f-\eta_n^c](\varphi)^2. \end{align*}

Then using $\varphi\in\mathcal{A}$ (along with the $C_2-$ inequality) and [Reference Jasra, Kamatani, Law and Zhou22, Lemma D.2], there is a finite constant $C<+\infty$ that depends on n such that

\[\check{\eta}_n^I\big(\big\{\varphi\otimes 1-1\otimes\varphi - \big(\big[\eta_n^f-\eta_n^c\big](\varphi)\big)\big\}^2\big) \leq C\Bigg(\frac{\big(\eta_{n-1}^f\otimes\eta_{n-1}^c\big)\big([G_{n-1}\otimes G_{n-1}]\check{M}_n(\tilde{\varphi})\big) }{\big(\eta_{n-1}^f\otimes\eta_{n-1}^c\big)(G_{n-1}\otimes G_{n-1})}+\Delta_l^2\Bigg),\]

where $\tilde{\varphi}(x,y) = \|x-y\|^2$ . Now, applying (6), one has the upper bound

\[C\Bigg(\frac{\big(\eta_{n-1}^f\otimes\eta_{n-1}^c\big)\big([G_{n-1}\otimes G_{n-1}]\tilde{\varphi}\big) }{\big(\eta_{n-1}^f\otimes\eta_{n-1}^c\big)\big(G_{n-1}\otimes G_{n-1}\big)}+\Delta_l+\Delta_l^2\Bigg).\]

In general, there is no reason to expect that $(\eta_{n-1}^f\otimes\eta_{n-1}^c)([G_{n-1}\otimes G_{n-1}]\tilde{\varphi})$ is small as a function of $\Delta_l$ . This is unsuprising, because one uses an independent coupling in the resampling operation. As a result, in the sense of Section 2.3.1 for ML estimation, the IRCPF is not useful.

5.3. MCIRCPF

To start our discussion, suppose first that $G_n$ is constant for each $n\geq 0$ . This represents the most favorable scenario for the MCIRCPF, because in the resampling operation, the resampled indices for each pair are equal (of course one would not use resampling in this case). For the metric in (H7), we take $\textsf{d}_1$ to be the $L_1-$ norm and set $\varphi\in\textsf{B}_b(\textsf{X})\cap\textrm{Lip}_{\textsf{d}_1}(\textsf{X})$ ( $\textsf{X}=\mathbb{R}^d$ ). Then setting $\tilde{\varphi}(x,y)=\|x-y\|^2$ , we have

\[\check{\eta}_n^M\Big(\{\varphi\otimes 1-1\otimes\varphi - ([\eta_n^f-\eta_n^c](\varphi))\}^2\Big) \leq C\|\varphi\|_{\textrm{Lip}}^2\check{\eta}_n^M(\tilde{\varphi}),\]

where $C<+\infty$ does not depend on n. Now, as $G_n$ is constant, one can deduce (e.g. by [Reference Mao25, Theorem 2.7.3]; recall (D) is assumed) that

(12) \begin{equation}\check{\eta}_n^M\Big(\{\varphi\otimes 1-1\otimes\varphi - ([\eta_n^f-\eta_n^c](\varphi))\}^2\Big) \leq Cn(n+4)\exp\{C'n^2\}\Delta_l,\end{equation}

where $C,C'<+\infty$ do not depend on n. Now (12) is not necessarily tight in n for every diffusion that satisfies (D); for instance if one has $dZ_t = a dt + b dW_t$ , for $b>0$ , $a\in({-}1,0)$ , then there is no dependence on n in the right-hand side of (12). However, it is worrying that the coupling can be exponentially bad in time, in this favorable case. The main point here is that the principal source of coupling in this algorithm is $\check{M}$ , and if this coupling deteriorates when iterating $\check{M}$ , then for the MCIRCPF one cannot hope to obtain time-uniform couplings.

More generally, if $G_n$ is nonconstant, we first remark that the upper bound (12) is likely to be $\mathcal{O}(\Delta_l^{1/2})$ , as it is the resampling operation where the forward rate is lost (see [Reference Jasra, Kamatani, Law and Zhou22]). Secondly, one expects to find examples where $\sigma^{2,s}_n(\varphi)$ is time-uniform or grows slowly as a function of n (due to the empirical results in [Reference Jasra, Kamatani, Law and Zhou22]). However, because of the highly nonlinear and complex expression for $\check{\Phi}_n^M$ , we expect a general result to be particularly arduous to obtain; this is left to future work.

5.4. MCPF

We have the following result which establishes the time-uniform coupling of the MCPF. The proof can be found in Appendix F, Section F.2.

Proposition 5.2. Assume (H1)–(H6), (H8). Then for any $(\xi,\hat{\xi})\in(0,1/32)\times(0,1/2)$ there exists a $C<+\infty$ depending on the constants in (H1)–(H6) and (H8) such that for any $\varphi\in\mathcal{B}_b(\textsf{X})$ , $n\geq 0$ , $1\leq l \leq L$ , we have

\[\sigma^{2,C}_n(\varphi) \leq C\|\varphi\|^2 \Delta_lv(x^*)^{32\xi+2\hat{\xi}}.\]

5.5. WCPF

Note that $\textsf{X}=\mathbb{R}$ and, for the metric in (H7), we take $\textsf{d}_1$ to be the $L_1-$ norm. Let $\tilde{\varphi}(x,y)=(x-y)^2$ . The proof of the following result, which establishes the time-uniform coupling of the WCPF, can be found in Appendix F, Section F.3.

Proposition 5.3. Assume (H1)–(H8). Then let $\lambda\in(0,1)$ be given; assume $\tilde{\varphi}\in\mathcal{L}_{(v\otimes v)^{\tilde{\xi}}}(\textsf{X})$ , for any $\tilde{\xi}\in(0,1/(16(1+\lambda)))$ , and set

\[(\xi,\hat{\xi})\in(0,\min\{1/32,\lambda/(16(1+\lambda)),(1-2\tilde{\xi})/12\})\times (0,1/2).\]

Then there exists a $C<+\infty$ depending on the constants in (H1)–(H8) such that for any $\varphi\in\mathcal{B}_b(\textsf{X})\cap\textrm{\textit{Lip}}_{\textsf{d}_1}(\textsf{X})$ , $n\geq 0$ , $1\leq l \leq L$ , we have

\[\sigma^{2,W}_n(\varphi) \leq C\|\varphi\|_{\textrm{\textit{Lip}}}^2\|\tilde{\varphi}\|_{v^{\tilde{\xi}}}(\Delta_l)^{1/(1+\lambda)}v(x^*)^{20\xi + (16(\lambda+1)\tilde{\xi}+2\hat{\xi})/(1+\lambda)}.\]

As $\lambda$ is greater than 0 in Proposition 5.3 (and can be made close to zero), one almost has the (time-uniform) forward error rate for the WCPF. We believe that $\lambda=0$ is the desired case and discuss strategies to establish this in Remark F.2 in Appendix F, Section F.3.

5.6. Example

We consider $\textsf{X}=\mathbb{R}$ , the diffusion $dZ_t = -\frac{3}{2}Z_tdt+dW_t$ and $\textsf{Y}=\{0,1\}$ and

\[G(x,y) = \bigg(\frac{be^x+a}{(1+e^x)}\bigg)^y\bigg(1-\bigg(\frac{be^x+a}{(1+e^x)}\bigg)\bigg)^{1-y}\]

for $y\in\textsf{Y}$ and any $0<a<b<1$ . The metric in (H7) is the $L_1-$ norm. In practice the CPFs can only be run (with currently available computational power) with $L\leq 20$ . So we will assume that there is an $L^*>L$ for which one cannot run the algorithm (say $L^*=50$ ). This will reduce the complexity of the forthcoming discussion. For the Euler discretization (although of course it is not required here) one can determine that the transition kernel $M^l(x,y)$ is the density of a $\mathcal{N}(\alpha_lx,\beta_l)$ distribution with $\alpha_l=(1-(3/2)\Delta_l)^{\Delta_l^{-1}}$ and

\[\beta_l=\Delta_l(1-(1-(3/2)\Delta_l)^{2\Delta_l^{-1}})/(1-(1-(3/2)\Delta_l)^2).\]

Thanks to our assumption concerning $L^*$ , one can easily find constants $0<\underline{\alpha}<\overline{\alpha}<1$ , $0<\underline{\beta}<\overline{\beta}\leq 1$ that do not depend on l, with $\underline{\alpha}\leq |\alpha_l|\leq \overline{\alpha}$ , $\underline{\beta}\leq\beta_l\leq\overline{\beta}$ for each $l\in\{0,\dots,L\}$ .

To verify (H1) we use, as in [Reference Whiteley, Kantas and Jasra31, Example 4.2.1], the Lyapunov function

\[\tilde{V}(x) = \frac{1}{2(1+\delta_0)}x^2 +1\]

for some $\delta_0>0$ . One can show that for any $1\leq l \leq L$ , $n\geq 1$ , $y\in\textsf{X}$ ,

\[Q_n^l(e^{\tilde{V}})(y) \leq \exp\bigg\{\frac{\alpha_l^2(1+\delta_0)}{1+\delta_0-\beta_l}(\tilde{V}(y)-1) + C\bigg\},\]

where $[\alpha_l^2(1+\delta_0)]/[1+\delta_0-\beta_l]<1$ for every $0\leq l\leq L$ , $\delta_0>0$ , and $0<C<+\infty$ is a constant that does not depend on $l,n,\delta_0$ and can be made arbitrarily large. One can easily verify (H1) for any $\delta_0>0$ with $\underline{d}>1$ large enough. (H2)–(H5) follow by elementary calculations associated to the choice of $\tilde{V}$ , $M^l$ , and $G_n$ and are omitted.

To verify (H6) as $G_n$ is bounded below, we consider $M^l(C_d)(x)$ and suppose $x>0$ , $l\geq 1$ , $d>1$ (the case $x\leq 0$ , $l=0$ can be verified in a similar manner). We have for any $x>0$ that

\[v(x) M^l(C_d)(x) \geq C \exp\bigg\{\frac{x^2}{2}\bigg(\frac{1}{\delta_0+1}-\frac{\alpha_l^2}{\beta_l}\bigg)-\frac{\alpha_l}{\beta_l}x\tilde{d}\bigg\},\]

where $\tilde{d}=\sqrt{(d-1)2(1+\delta_0)}$ . As $\delta_0$ is arbitrary one can choose $\delta_0>0$ so that

\[\frac{1}{\delta_0+1}-\frac{\alpha_l^2}{\beta_l}\geq 0\]

for any $1\leq l \leq L$ . Checking calculations on a computer for $L^*=50$ yields that $0<\delta_0\leq 5$ suffices.

For (H7) one has the following result, whose proof is in Appendix F, Section F.4.

To verify (H8), we shall suppose that C $^{\prime}$ in (11) is at least $1+1/(2(1+\delta_0))$ (i.e. $C'>1/12$ ). It is noted that it is nontrivial to check the value of C $^{\prime}$ in [Reference Bally and Talay1], but such a choice seems reasonable. In the given situation, it is then straightforward to verify (H8).