# Higher-order spring-coupled multilevel Monte Carlo method for invariant
measures

Håkon Hoel  Department of Mathematics, University of Oslo, Norway
<EMAIL_ADDRESS>and  Sankarasubramanian Ragunathan⋆ Chair of Mathematics
for Uncertainty Quantification, RWTH Aachen University, Germany
<EMAIL_ADDRESS>
###### Abstract.

A higher-order change-of-measure multilevel Monte Carlo (MLMC) method is
developed for computing weak approximations of the invariant measures of SDE
with drift coefficients that do not satisfy the contractivity condition. This
is achieved by introducing a spring term in the pairwise coupling of the MLMC
trajectories employing the order 1.5 strong Itô–Taylor method. Through this,
we can recover the contractivity property of the drift coefficient while still
retaining the telescoping sum property needed for implementing the MLMC
method.

We show that the variance of the change-of-measure MLMC method grows linearly
in time $T$ for all $T>0$, and for all sufficiently small timestep size $h>0$.
For a given error tolerance $\epsilon>0$, we prove that the method achieves a
mean-square-error accuracy of $\mathcal{O}(\epsilon^{2})$ with a computational
cost of
$\mathcal{O}(\epsilon^{-2}\left|\log\epsilon\right|^{3/2}(\log\left|\log\epsilon\right|)^{1/2})$
for uniformly Lipschitz continuous payoff functions and
$\mathcal{O}\big{(}\epsilon^{-2}\left|\log\epsilon\right|^{5/3+\xi}\big{)}$
for discontinuous payoffs, respectively, where $\xi>0$. We also observe an
improvement in the constant associated with the computational cost of the
higher-order change-of-measure MLMC method, marking an improvement over the
Milstein change-of-measure method in the aforementioned seminal work by M.
Giles and W. Fang. Several numerical tests were performed to verify the
theoretical results and assess the robustness of the method.

###### Key words and phrases:

change of measure, ergodic measure, SDE, order 1.5 strong Itô–Taylor method,
Girsanov’s theorem, MLMC method

###### 2020 Mathematics Subject Classification:

60H10, 60H35, 65C05, 37M25

⋆Corresponding author: S. Ragunathan<EMAIL_ADDRESS>
##  1\. Introduction

In this paper, we construct a higher-order multilevel Monte Carlo (MLMC)
method for weak approximations of the stationary distribution of
$d$-dimensional stochastic differential equations of the form

$\begin{split}\operatorname{d\\!}X(t)&=a(X(t))\operatorname{d\\!}t+\operatorname{d\\!}W(t)\\\
X(0)&=x_{0}\in\mathbb{R}^{d}.\end{split}$ (1)

Here, $a\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ denotes the drift coefficient,
$W(t)$ is a $d$-dimensional standard Wiener process on a filtered, complete
probability space $\big{(}\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq
0},\mathbb{P}\big{)}$, and the deterministic initial condition $x_{0}$ is
given. More precisely, we develop a higher-order change-of-measure numerical
method for computing expected values of the form

$\mathbb{E}\big{[}\varPhi(X_{\infty})\big{]}\coloneqq\int_{\mathbb{R}^{d}}\varPhi(x)\,\pi(\operatorname{d\\!}x)=\lim\limits_{t\to\infty}\mathbb{E}\big{[}\varPhi(X_{t})\big{]},$

where $\varPhi\colon\mathbb{R}^{d}\to\mathbb{R}$ represents a quantity of
interest (QoI) and $\pi$ denotes the invariant distribution of (1). Further
details on the regularity and contractivity of the drift coefficient, ensuring
that the invariant distribution is well-defined, are provided in Section 2.

The main contribution of this work is to show that, under sufficient
regularity conditions, the higher-order change-of-measure MLMC method weakly
approximates $\mathbb{E}[\varPhi(X_{\infty})]$ with a mean-square-error
accuracy of $\mathcal{O}(\epsilon^{2})$ at the computational cost of
$\mathcal{O}(\epsilon^{-2}\left|\log\epsilon\right|^{3/2}(\log\left|\log\epsilon\right|)^{1/2})$,
and
$\mathcal{O}\big{(}\epsilon^{-2}\left|\log\epsilon\right|^{5/3+\xi}\big{)}$
for all $\xi>0$, for a payoff function $\varPhi$ that is uniformly Lipschitz
continuous, and discontinuous, respectively. The higher-order change-of-
measure scheme also has a better constant associated with the computational
cost than the Milstein scheme. This is an improvement from
$\mathcal{O}(\epsilon^{-2}\left|\log\epsilon\right|^{2})$ obtained in [9,
Theorem 5].

Computing weak approximations of the invariant measure have various
applications in physics, biology, statistical mechanics and other fields
alike. In molecular dynamics, ergodicity plays an important role in
determining the properties of the system [12]. In quantum mechanics, the
ergodic measure provides a statistical description of the spacing between the
energy levels of many-body quantum systems [5, 23]. In molecular biology, the
ergodic measure relates to the motion of bio-molecules and dictates the
essential functions of a cell [25].

The invariant measure can be computed either using Monte Carlo techniques, by
the numerical simulation of SDE paths, or by solving the associated partial
differential equation; the Fokker–Planck equation [33, equation (8.6)]. In
high-dimensional state space, computing the closed-form analytical solution to
the Fokker–Planck equation is oftentimes not tractable, if not impossible, due
to the curse of dimensionality. Since we are only interested in computing weak
approximations of the invariant measure, we use Monte Carlo simulations
because it is numerically tractable and scales well to higher dimensions.

Note that the naive Monte Carlo method can be computationally expensive as one
might need many Monte Carlo samples to achieve $\mathcal{O}(\epsilon^{2})$
mean-squared error of the estimated QoI. In [13, 14], the MLMC method is
formulated to compute weak approximations of a QoI w.r.t the measure of an
SDE. For the case of a uniformly Lipschitz continuous payoff function, and
using the Euler–Maruyama scheme, it was shown that a mean-squared error of
$\mathcal{O}(\epsilon^{-2})$ can be achieved using a computational cost of
$\mathcal{O}(\epsilon^{-2}\left|\log\epsilon\right|^{2})$, which is an
improvement from $\mathcal{O}(\epsilon^{-3})$ using the naive Monte Carlo
method. The study was later extended to the adaptive Euler-Maruyama scheme
with non-uniform and path-dependent time-step sizes in [18, 19]. The results
from these works form the basis for the development of an MLMC method to
compute weak approximations of invariant measures.

Following the idea above, the MLMC method was extended to compute weak
approximations of invariant measures in [16] for SDE with drift coefficients
that satisfy the contractivity condition, i.e. the drift coefficient is one-
sided Lipschitz continuous with a negative one-sided Lipschitz constant. The
authors show that under appropriate pairwise contracting coupling of the
trajectories of the MLMC estimator, one can obtain a uniform-in-time estimate
of the MLMC variance. The Milstein scheme is proposed in [34] to compute weak
approximations of invariant measures of SDE. Assuming that the diffusion
coefficient satisfies the commutativity condition [22, equation (3.13)] and
the drift coefficient satisfies the contractivity condition, the authors prove
a strong convergence rate of order 1. The authors also show that the
probability density function associated with the numerical solution converges
exponentially to the true invariant distribution.

Higher-order numerical methods for approximating the ergodic measure of SDE is
developed in [1]. Sufficient conditions necessary to weakly approximate
ergodic measure of SDE with a high order of accuracy, irrespective of the weak
order of the numerical scheme, are developed. The convergence properties of a
wide class of adaptive numerical approximations of SDE are studied in [11].
The authors show that as long as the adaptive numerical method does not
introduce any bias in the computation of the “Lévy area”, the numerical
approximation converges to the actual SDE dynamics. The seminal works provide
a framework for computing weak approximations of invariant measures of SDE and
to establish convergence results for methods employing adaptive numerical
integrators. This is of interest to our current work and for future
extensions.

Explicit adaptive time-stepping schemes employing the Euler–Maruyama scheme
were developed in [21, 24] to compute the expected value of QoIs that depend
on the solution of SDE with non-uniformly Lipschitz continuous drift
coefficients. The fundamental idea behind adaptive time-stepping is to use a
time-step size that depends on the drift term. Using a smaller time-step size
as the magnitude of the drift term becomes large ensures that the numerical
solution is stable and bounded. The idea was extended in [6, 7, 8] for the
MLMC method. The authors prove theoretically that using adaptive time-
stepping, under suitable pairwise contracting coupling of the MLMC
trajectories, the $p$-th moment of the difference between the true SDE
solution and the numerical solution is uniformly bounded in time $T$ for all
$T>0$, provided the drift coefficient satisfies the contractivity condition.

In all the above works, it is assumed that the drift coefficient satisfies the
contractivity condition. In practice, only a small class of SDE have drift
coefficients that satisfy the contractivity condition. This impedes the
effective implementation of the above methods for solving SDE that do not
satisfy the assumption. A change of measure technique has hence been developed
in [9] to extend the invariant measure computation to drift coefficients that
satisfy the one-sided Lipschitz continuity condition. The contractivity
property of the drift coefficient is recovered by constructing a numerical
scheme with a “spring” term. The authors also show a linear growth in the
variance of the change-of-measure MLMC estimator with respect to time $T$ for
all $T>0$. A similar approach is utilized in [3] and [10] to develop coupled
particle filters that work for Itô diffusion that have unstable dynamics, and
to compute pathwise sensitivities of SDE, respectively.

The rest of the paper is organized as follows: Sections 2 and 3 introduce the
necessary notation, regularity assumptions, and the main theoretical results
on convergence and computational cost associated with the higher-order change-
of-measure MLMC scheme, respectively. In Section 4, we will look at how the
higher-order change-of-measure algorithm performs by considering different
numerical test cases. In Section 5, we will discuss the conclusions that can
be drawn from the current work and the possible future extensions. Finally,
Section 6 contains the proofs for the theorems stated in Section 3.

##  2\. Notation and the higher-order change-of-measure method

In this section, we state the necessary notation, and assumptions on the SDE
drift coefficient and its derivatives to ensure the existence and uniqueness
of a solution to the SDE (1), and to ensure that it has an invariant measure.
We also state the higher-order change-of-measure scheme and the corresponding
theoretical results for the convergence and the computational cost.

###  2.1. Notation and assumptions

Let $\left\langle v,w\right\rangle$ denote the inner product between two
vectors $v,w\in\mathbb{R}^{d}$ defined as

$\left\langle v,w\right\rangle\coloneqq
v_{1}w_{1}+v_{2}w_{2}+\cdots+v_{d}w_{d},$

where $v_{i}$ represents the $i$-th component of the vector
$v\in\mathbb{R}^{d}$ for all $i=1,\ldots,d$. The standard Euclidean norm is
induced from the vector inner product and is defined as

$\left\|v\right\|\coloneqq\sqrt{\sum_{i=1}^{d}v_{i}^{2}},\quad\forall
v\in\mathbb{R}^{d}.$

The definition of the standard Euclidean norm can be extended similarly to a
$k$-th order tensor $A\in(\mathbb{R}^{d})^{\otimes k}$

$\left\|A\right\|\coloneqq\sqrt{\sum_{i_{1},i_{2},\cdots,i_{k}=1}^{d}\left|A_{i_{1}i_{2}\cdots
i_{k}}\right|^{2}},\qquad\forall k\in\mathbb{N}\setminus\\{1\\},$

where $\otimes$ denotes the tensor outer product [31].

For any integer $k\geq 1$, let $C^{k}_{b}(\mathbb{R}^{d})$ denote the set of
scalar-valued functions with continuous and uniformly bounded partial
derivatives up to and including order $k$. We make the following assumption on
the regularity of the drift coefficient of SDE (1),

######  Assumption 1 (Regularity of $a(x)$).

The components of the drift coefficient $a_{i}\in
C_{b}^{3}\big{(}\mathbb{R}^{d}\big{)}$, for all $i=1,\ldots,d$. We also assume
that there exists a constant $K_{b}>0$ such that

$\sup_{x\in\mathbb{R}^{d}}\max\Bigg{(}\left|\frac{\partial a_{i}}{\partial
x_{j}}(x)\right|,\left|\frac{\partial^{2}a_{i}}{\partial x_{j}\partial
x_{k}}(x)\right|,\left|\frac{\partial^{3}a_{i}}{\partial x_{j}\partial
x_{k}\partial x_{l}}(x)\right|\Bigg{)}\eqqcolon K_{b}<\infty,\quad\forall
i,j,k=1,\ldots,d.$

From Assumption 1, using the Mean value theorem, the following implications
hold for all $i=1,\ldots,d$

$\begin{split}\left\|a(x)-a(y)\right\|&=\sqrt{\sum_{i=1}^{d}\left|a_{i}(x)-a_{i}(y)\right|^{2}}\leq\sqrt{d}K_{b}\left\|x-y\right\|,\\\
~{}\left\|Da(x)-Da(y)\right\|&\leq dK_{b}\left\|x-y\right\|,\\\
~{}\left\|D^{2}a(x)-D^{2}a(y)\right\|&\leq
d^{3/2}K_{b}\left\|x-y\right\|,\end{split}$

where

$\displaystyle Da(x)$ $\displaystyle\coloneqq\frac{\partial a_{i}}{\partial
x_{j}}(x)\;e_{i}\otimes e_{j},$ $\displaystyle Da[u,v]$
$\displaystyle=\sum_{i,j=1}^{d}\frac{\partial a_{i}}{\partial
x_{j}}u_{i}v_{j},$ $\displaystyle D^{2}a(x)$
$\displaystyle\coloneqq\frac{\partial^{2}a_{i}}{\partial x_{j}\partial
x_{k}}(x)\;e_{i}\otimes e_{j}\otimes e_{k},$ $\displaystyle D^{2}a[u,v,w]$
$\displaystyle=\sum_{i,j,k=1}^{d}\frac{\partial^{2}a_{i}}{\partial
x_{j}\partial x_{k}}u_{i}v_{j}w_{k},$ $\displaystyle D^{3}a(x)$
$\displaystyle\coloneqq\frac{\partial^{3}a_{i}}{\partial x_{j}\partial
x_{k}\partial x_{l}}(x)\;e_{i}\otimes e_{j}\otimes e_{k}\otimes e_{l},$
$\displaystyle D^{3}a[u,v,w,z]$
$\displaystyle=\sum_{i,j,k,l=1}^{d}\frac{\partial^{3}a_{i}}{\partial
x_{j}\partial x_{k}\partial x_{l}}u_{i}v_{j}w_{k}z_{l}$

for all $u,v,w,z\in\mathbb{R}^{d}$, where $e_{k}$ denotes the standard
d-dimensional basis vector, i.e. a column vector with 1 in the $k$-th entry
and 0 everywhere else, and $\otimes$ denotes the tensor product. To keep
things concise, we also introduce the short-hand notation $Da(x)u\coloneqq
Da(x)[\cdot,u]$ and $(Da(x))^{\top}u\coloneqq Da(x)[u,\cdot]$ for all
$u\in\mathbb{R}^{d}$. Similarly, we also introduce the notation
$D^{2}a(x)uv\coloneqq D^{2}a(x)[.,u,v]$ and $D^{2}a(x)u^{2}\coloneqq
D^{2}a(x)[.,u,u]$ for all $u,v\in\mathbb{R}^{d}$.

Assumption 1 ensures that there exists a unique solution to the SDE (1), [30,
Theorem 5.2.1]. Assumption 1 also implies that there exists a constant
$\lambda>0$ such that the following inequality holds

$\displaystyle\left\langle x-y,a(x)-a(y)\right\rangle$
$\displaystyle\leq\lambda\left\|x-y\right\|^{2},\quad\forall
x,y\in\mathbb{R}^{d},$

where $\lambda$ denotes the one-sided Lipschitz constant.

######  Assumption 2 (uniform Lipschitz continuity of $Da(x)a(x)$).

There exists a constants $L>0$, such that

$\left|\Big{(}a_{i}\frac{\partial a_{k}}{\partial
x_{i}}\Big{)}(x)-\Big{(}a_{i}\frac{\partial a_{k}}{\partial
x_{i}}\Big{)}(y)\right|\leq L\left\|x-y\right\|,\qquad\forall
x,y\in\mathbb{R}^{d}$

for all $i,k=1,\ldots,d$.

Let $\mathcal{A}$ denote the infinitesimal generator of the Itô process (1)
defined as

$\displaystyle\mathcal{A}\coloneqq\sum_{i=1}^{d}a_{i}(x)\frac{\partial}{\partial
x_{i}}+\frac{1}{2}\sum_{i=1}^{d}\frac{\partial^{2}}{\partial
x_{i}^{2}}=a(x)\cdot D+\frac{1}{2}\Delta,$ (2)

where $\Delta$ represents the Laplacian operator. The Laplacian applied to the
vector-valued drift coefficient $a$ is defined as the following component-wise
operation

$\displaystyle\Delta
a(x)\coloneqq\Bigg{(}\sum_{i=1}^{d}\frac{\partial^{2}a_{1}}{\partial
x_{i}^{2}}(x),\cdots,\sum_{i=1}^{d}\frac{\partial^{2}a_{d}}{\partial
x_{i}^{2}}(x)\Bigg{)}^{\top}.$

The infinitesimal generator is needed to construct the order 1.5 strong
Itô–Taylor scheme and occurs naturally due to the Itô–Taylor expansion.

######  Proposition 1 (uniformly Lipschitz continuity of $\mathcal{A}$).

Let Assumptions 1 and 2 hold. Then, there exists a constant $\widehat{L}>0$
such that

$\left\|\mathcal{A}a(x)-\mathcal{A}a(y)\right\|\leq\widehat{L}\left\|x-y\right\|,\quad\forall
x,y\in\mathbb{R}^{d}.$

###### Proof.

Consider the operator $\mathcal{A}$ applied to the $k$-th component of the
drift coefficient $a$. We then obtain

$\begin{split}\left|\mathcal{A}a_{k}(x)-\mathcal{A}a_{k}(y)\right|&=\left|\sum_{i=1}^{d}\bigg{(}a_{i}(x)\frac{\partial
a_{k}(x)}{\partial x_{i}}-a_{i}(y)\frac{\partial a_{k}(y)}{\partial
x_{i}}\bigg{)}+\frac{1}{2}\sum_{i=1}^{d}\bigg{(}\frac{\partial^{2}a_{k}(x)}{\partial
x_{i}^{2}}-\frac{\partial^{2}a_{k}(y)}{\partial x_{i}^{2}}\bigg{)}\right|\\\
&\leq\sum_{i=1}^{d}\left|a_{i}(x)\frac{\partial a_{k}(x)}{\partial
x_{i}}-a_{i}(y)\frac{\partial a_{k}(y)}{\partial
x_{i}}\right|+\frac{1}{2}\sum_{i=1}^{d}\left|\frac{\partial^{2}a_{k}(x)}{\partial
x_{i}^{2}}-\frac{\partial^{2}a_{k}(y)}{\partial x_{i}^{2}}\right|\\\
&\leq\sum_{i=1}^{d}L\left\|x-y\right\|+\frac{1}{2}\sum_{i=1}^{d}K_{b}\left\|x-y\right\|\\\
&=\widehat{L}_{k}\left\|x-y\right\|,\end{split}$

where $\widehat{L}_{k}\coloneqq d\Big{(}L+\frac{1}{2}K_{b}\Big{)}$. We then
have

$\begin{split}\left\|\mathcal{A}a(x)-\mathcal{A}a(y)\right\|=\sqrt{\sum_{k=1}^{d}\left|\mathcal{A}a_{k}(x)-\mathcal{A}a_{k}(y)\right|^{2}}\leq\underbrace{d^{3/2}\Big{(}L+\frac{1}{2}K_{b}\Big{)}}_{\eqqcolon\widehat{L}}\left\|x-y\right\|\end{split}$

∎

The one-sided Lipschitz property of the operator $\mathcal{A}$ acting on the
drift coefficient $a$ follows from the uniform Lipschitz continuity of
$\mathcal{A}$, i.e. there exists a constant $\widehat{\lambda}>0$ such that

$\left\langle
x-y,\mathcal{A}a(x)-\mathcal{A}a(y)\right\rangle\leq\widehat{\lambda}\left\|x-y\right\|^{2}.$

The regularity assumptions on the drift coefficient and its derivatives ensure
that the order 1.5 strong Itô–Taylor scheme is well defined and is required
for the implementation of the higher-order change-of-measure scheme [22,
Theorem 10.6.3].

Note that the regularity assumption on the drift coefficient alone does not
guarantee the existence of an invariant measure. To this end, we assume that
the drift coefficient also satisfies the dissipativity condition given below

######  Assumption 3 (Dissipativity condition).

The drift coefficient satisfies the dissipativity condition, i.e. there exists
constants $\alpha,\beta>0$ such that

$\left\langle
x,a(x)\right\rangle\leq-\alpha\left\|x\right\|^{2}+\beta,\qquad\forall
x\in\mathbb{R}^{d}$

Assumptions 1 and 3 are sufficient to guarantee that the SDE (1) has an
invariant measure $\pi(x)$ and that the solution to the SDE converges
exponentially to $\pi(x)$ [26, 27, 28, 29], i.e.

$\left|\mathbb{E}[\varPhi(X_{T})]-\mathbb{E}[\varPhi(X_{\infty})]\right|\leq\mu^{\star}\operatorname{e}^{-\lambda^{\star}T},$
(3)

for some constants $\mu^{\star},\lambda^{\star}>0$. More details regarding the
sufficient conditions for an Itô process to be ergodic can be found in [17,
22].

###  2.2. The higher-order change-of-measure method

Having made the necessary assumptions needed to ensure the well-posedness of a
solution to the SDE (1), and the existence of an invariant measure $\pi(x)$,
the expected value of a QoI under the $\mathbb{Q}$ measure, given by
$\mathbb{E}^{\mathbb{Q}}[\varPhi(X_{T})]$, can be approximated using the
standard MLMC method as

$\mathbb{E}^{\mathbb{Q}}[\varPhi(\overline{X}_{T}^{L})]=\mathbb{E}^{\mathbb{Q}}[\varPhi(\overline{X}_{T}^{0})]+\sum_{\ell=1}^{L}\mathbb{E}^{\mathbb{Q}}[\varPhi(\overline{X}_{T}^{f,\ell})-\varPhi(\overline{X}_{T}^{c,\ell-1})]$

where $L$ denotes the total number of levels in the MLMC estimator, $T>0$
represents the terminal time, and $\overline{X}^{f}$ and $\overline{X}^{c}$
denote the fine and the coarse trajectories, respectively. Both the fine and
the coarse trajectories are numerically simulated using the order 1.5 strong
Itô–Taylor scheme given by
Fine trajectory, $\mathbb{Q}$-measure:

$\displaystyle\begin{split}\overline{X}^{f}(t_{2n+1})&=\overline{X}^{f}(t_{2n})+ha(\overline{X}^{f}(t_{2n}))+\Delta
W_{2n}^{\mathbb{Q}}+Da(\overline{X}^{f}(t_{2n}))\Delta
Z_{2n}^{\mathbb{Q}}+\frac{h^{2}}{2}\mathcal{A}a(\overline{X}^{f}(t_{2n})),\end{split}$
(4)
$\displaystyle\begin{split}\overline{X}^{f}(t_{2n+2})&=\overline{X}^{f}(t_{2n+1})+ha(\overline{X}^{f}(t_{2n+1}))+\Delta
W_{2n+1}^{\mathbb{Q}}+Da(\overline{X}^{f}(t_{2n+1}))\Delta
Z_{2n+1}^{\mathbb{Q}}\\\
&\quad+\frac{h^{2}}{2}\mathcal{A}a(\overline{X}^{f}(t_{2n+1})),\end{split}$
(5)

Coarse trajectory, $\mathbb{Q}$-measure:

$\displaystyle\begin{split}\overline{X}^{c}(t_{2n+2})&=\overline{X}^{c}(t_{2n})+2ha(\overline{X}^{c}(t_{2n}))+\Big{(}\Delta
W_{2n+1}^{\mathbb{Q}}+\Delta W_{2n}^{\mathbb{Q}}\Big{)}\\\
&\quad+Da(\overline{X}^{c}(t_{2n}))\Big{(}\Delta Z_{2n+1}^{\mathbb{Q}}+\Delta
Z_{2n}^{\mathbb{Q}}+h\Delta
W_{2n}^{\mathbb{Q}}\Big{)}+2h^{2}\mathcal{A}a(\overline{X}^{c}(t_{2n})),\end{split}$
(6)

for $n=0,1,\ldots,N/2-1$ where $t_{n}=nh$, $N=T/h\in\mathbb{N}$,
$h=2^{-\ell}h_{0}$ denotes the uniform time-step size at level $\ell$ of the
multilevel estimator, $h_{0}$ denotes the time-step size at level $\ell=0$,
$\overline{X}^{f}(0)=\overline{X}^{c}(0)=x_{0}$, and $\Delta Z_{n}$ is a
random variable defined as

$\Delta
Z_{n}^{k}\coloneqq\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{s}\operatorname{d\\!}W_{s}^{k}\operatorname{d\\!}s,\qquad\forall
k=1,\ldots,d.$

The tuple of correlated random variables $(\Delta W_{n}^{\mathbb{Q}},\Delta
Z_{n}^{\mathbb{Q}})$, under the $\mathbb{Q}$ measure, can be generated by

$\Delta W_{n}^{\mathbb{Q}}\coloneqq\Delta U_{1,n},\qquad\Delta
Z_{n}^{\mathbb{Q}}\coloneqq\frac{1}{2}h\Big{(}\Delta
U_{1,n}+\frac{1}{\sqrt{3}}\Delta U_{2,n}\Big{)},$

where $\Delta U_{1,n},\Delta U_{2,n}$ are independent identically distributed
(i.i.d.) $\mathcal{N}^{\mathbb{Q}}(0,h\operatorname{I_{d}})$ samples, and
$\operatorname{I_{d}}$ denotes the $d\times d$ identity matrix. Note that

$\mathbb{E}\big{[}\left\langle\Delta W_{i}^{\mathbb{Q}},\Delta
W_{j}^{\mathbb{Q}}\right\rangle\big{]}=dh\delta_{ij},\qquad\mathbb{E}\big{[}\left\langle\Delta
Z_{i}^{\mathbb{Q}},\Delta
Z_{j}^{\mathbb{Q}}\right\rangle\big{]}=\frac{dh^{3}}{3}\delta_{ij},\qquad\mathbb{E}\big{[}\left\langle\Delta
W_{i}^{\mathbb{Q}},\Delta
Z_{j}^{\mathbb{Q}}\right\rangle\big{]}=\frac{dh^{2}}{2}\delta_{ij},$

where

$\delta_{ij}=\begin{cases}1,&i=j\\\ 0,&i\neq j\end{cases}$

represents the Kronecker delta function.

We cannot use the standard MLMC algorithm for computing the weak
approximations of the ergodic measure as the fine and the coarse trajectories
of the MLMC estimator may drift away from each other. This can lead to an
exponential growth in the multilevel variance w.r.t. the terminal time $T$,
which can lead to an increase in the computational cost needed to achieve good
variance estimates [8, 9]. To circumvent this issue, we introduce a spring
term in both the fine and the coarse trajectories.

Let $S>0$ denote the spring constant. The order 1.5 strong Itô–Taylor scheme,
after the introduction of the spring term, for both the fine and the coarse
trajectories can be written as
Fine trajectory, $\mathbb{P}$-measure:

$\displaystyle\begin{split}\widehat{Y}^{f}(t_{2n+1})&=\widehat{Y}^{f}(t_{2n})+Sh\big{(}\widehat{Y}^{c}(t_{2n})-\widehat{Y}^{f}(t_{2n})\big{)}+ha(\widehat{Y}^{f}(t_{2n}))+\Delta
W_{2n}^{\mathbb{P}},\\\ &\quad+Da(\widehat{Y}^{f}(t_{2n}))\Delta
Z_{2n}^{\mathbb{P}}+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}^{f}(t_{2n})),\end{split}$
(7)
$\displaystyle\begin{split}\widehat{Y}^{f}(t_{2n+2})&=\widehat{Y}^{f}(t_{2n+1})+Sh\big{(}\widehat{Y}^{c}(t_{2n+1})-\widehat{Y}^{f}(t_{2n+1})\big{)}+ha(\widehat{Y}^{f}(t_{2n+1}))+\Delta
W_{2n+1}^{\mathbb{P}},\\\ &\qquad+Da(\widehat{Y}^{f}(t_{2n+1}))\Delta
Z_{2n+1}^{\mathbb{P}}+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}^{f}(t_{2n+1})),\end{split}$
(8)

Coarse trajectory, $\mathbb{P}$-measure:

$\displaystyle\begin{split}\widehat{Y}^{c}(t_{2n+1})&=\widehat{Y}^{c}(t_{2n})+Sh\big{(}\widehat{Y}^{f}(t_{2n})-\widehat{Y}^{c}(t_{2n})\big{)}+ha(\widehat{Y}^{c}(t_{2n}))+\Delta
W_{2n}^{\mathbb{P}}\\\ &\quad+Da(\widehat{Y}^{c}(t_{2n}))\Delta
Z_{2n}^{\mathbb{P}}+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}^{c}(t_{2n})),\end{split}$
(9)
$\displaystyle\begin{split}\widehat{Y}^{c}(t_{2n+2})&=\widehat{Y}^{c}(t_{2n})+2Sh\big{(}\widehat{Y}^{f}(t_{2n})-\widehat{Y}^{f}(t_{2n})\big{)}+2ha(\widehat{Y}^{c}(t_{2n}))+(\Delta
W_{2n}^{\mathbb{P}}+\Delta W_{2n+1}^{\mathbb{P}})\\\
&\qquad+Da(\widehat{Y}^{c}(t_{2n}))(\Delta Z_{2n+1}^{\mathbb{P}}+\Delta
Z_{2n}^{\mathbb{P}}+h\Delta
W_{2n}^{\mathbb{P}})+2h^{2}\mathcal{A}a(\widehat{Y}^{c}(t_{2n})),\end{split}$
(10)

for $n=0,1,\ldots,N/2-1$. The tuple of random variables $(\Delta
W_{n}^{\mathbb{P}},\Delta Z_{n}^{\mathbb{P}})$ under the $\mathbb{P}$ measure
can be generated as

$\Delta W_{n}^{\mathbb{P}}\coloneqq\Delta V_{1,n},\qquad\Delta
Z_{n}^{\mathbb{P}}\coloneqq\frac{1}{2}h\Big{(}\Delta
V_{1,n}+\frac{1}{\sqrt{3}}\Delta V_{2,n}\Big{)},$ (11)

where $\Delta V_{1,n},\Delta
V_{2,n}\sim\mathcal{N}^{\mathbb{P}}(0,h\operatorname{I_{d}})$ are i.i.d.
samples. Note that both the fine and the coarse trajectories have the same
initial value, i.e. $\widehat{Y}^{f}(0)=\widehat{Y}^{c}(0)=x_{0}$.

Note that the order 1.5 strong Itô–Taylor scheme with the spring term
included, defined by the equations (7)-(10), does not solve the original SDE
(1). Using Girsanov’s theorem, we can compute weak approximations of the
invariant measure of the original SDE with numerical schemes that include the
“spring” term. Girsanov’s theorem states that there indeed exists a
transformation from the $\mathbb{P}$ measure, associated with the numerical
process with the “spring” term, to the $\mathbb{Q}$ measure, associated with
the original SDE (1) [33, Theorem 9.2].

Consider the fine trajectory of the numerical scheme with and without the
“spring” term, sampled from the measures $\mathbb{P}$ and $\mathbb{Q}^{f}$,
respectively, where $\mathbb{Q}^{f}$ denotes the measure associated with the
fine trajectory of the numerical scheme without the “spring” term. Under the
transformation

$\begin{split}\Delta U_{1,n}^{\mathbb{Q}^{f}}&=\widehat{S}^{f}_{n}h+\Delta
V_{1,n}^{\mathbb{P}}\\\ \Delta
U_{2,n}^{\mathbb{Q}^{f}}&=-\sqrt{3}\widehat{S}^{f}_{n}h+\Delta
V_{2,n}^{\mathbb{P}},\end{split}$ (12)

we have that

$\widehat{Y}_{n}^{f}(\omega)=\overline{X}_{n}^{f}(\omega),\quad\forall\omega\in\Omega$

for all $n=0,1,\ldots,N$, where

$\begin{split}\widehat{S}^{f}_{n}\coloneqq
S\Big{(}\widehat{Y}_{n}^{c}-\widehat{Y}_{n}^{f}\Big{)}.\end{split}$ (13)

The Radon–Nikodym derivative for the fine trajectory, according to Girsanov’s
theorem, is given by

$R^{f}_{T}=\frac{\operatorname{d\\!}\mathbb{Q}^{f}}{\operatorname{d\\!}\mathbb{P}}=\prod_{n=0}^{N-1}\exp\Bigg{(}-\left\langle\widehat{S}^{f}_{n},\Delta
V_{1,n}\right\rangle+\sqrt{3}\left\langle\widehat{S}^{f}_{n},\Delta
V_{2,n}\right\rangle-2h\left\|\widehat{S}^{f}_{n}\right\|^{2}\Bigg{)}.$ (14)

Similarly, consider the coarse trajectory of the numerical scheme with and
without the “spring” term, sampled from the measures $\mathbb{P}$ and
$\mathbb{Q}^{c}$, respectively, where $\mathbb{Q}^{c}$ denotes the measure
associated with the coarse trajectory of the numerical scheme without the
“spring” term. Under the transformation

$\begin{split}\Delta U_{1,2n}^{\mathbb{Q}^{c}}&=\widehat{S}^{c}_{2n}h+\Delta
V_{1,2n}^{\mathbb{P}}\\\ \Delta
U_{2,2n}^{\mathbb{Q}^{c}}&=-2\sqrt{3}\widehat{S}^{c}_{2n}h+\Delta
V_{2,2n}^{\mathbb{P}},\end{split}$ (15)

we obtain

$\widehat{Y}_{2n}^{c}(\omega)=\overline{X}_{2n}^{c}(\omega),\quad\forall\omega\in\Omega$

for all $n=0,\ldots,N/2$, where

$\begin{split}\widehat{S}^{c}_{2n}\coloneqq
S\Big{(}\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\Big{)}.\end{split}$ (16)

The Radon–Nikodym derivative for the coarse trajectory is then given by

$\begin{split}R^{c}_{T}=\frac{\operatorname{d\\!}\mathbb{Q}^{c}}{\operatorname{d\\!}\mathbb{P}}&=\prod_{n=0}^{N/2-1}\exp\Bigg{(}-\left\langle\widehat{S}^{c}_{2n},\Delta
V_{1,2n}+\Delta
V_{1,2n+1}\right\rangle+2\sqrt{3}\left\langle\widehat{S}^{c}_{2n},\Delta
V_{2,2n}+\Delta V_{2,2n+1}\right\rangle\\\
&\qquad-13h\left\|\widehat{S}^{c}_{2n}\right\|^{2}\Bigg{)}.\end{split}$ (17)

Using the Radon–Nikodym derivatives given by the equations (14) and (17), the
MLMC estimator under the change-of-measure can be written as

$\mathbb{E}^{\mathbb{Q}}\Big{[}\varPhi(\overline{X}_{T}^{L})\Big{]}=\mathbb{E}^{\mathbb{Q}}\Big{[}\varPhi(\overline{X}_{T}^{0})\Big{]}+\sum_{\ell=1}^{L}\mathbb{E}^{\mathbb{P}}\Big{[}\varPhi(\widehat{Y}_{T}^{f,\ell})R_{T}^{f,\ell}-\varPhi(\widehat{Y}_{T}^{c,\ell})R_{T}^{c,\ell}\Big{]}.$
(18)

Theorems 5 and 6 show that $\varPhi$ is $L^{p}(\Omega)$-stable under the
$\mathbb{P}$ measure. We also make the assumption that $\varPhi$ is atleast
$L_{1}(\Omega)$-stable under the $\mathbb{Q}$ measure, i.e.

$\mathbb{E}^{\mathbb{Q}}\Big{[}\left|\varPhi(\overline{X}_{T})\right|\Big{]}<\infty,\quad\forall
T>0.$

From the above assumption, and the results from Theorems 5 and 6, we have that
equation 18 is well-defined for all $T>0$.

Having defined the higher-order change-of-measure method, we will discuss the
corresponding main theoretical results in Section 3.

##  3\. Theoretical results

In this section, we state the key theoretical results that establish the
stability of the fine and the coarse trajectories of the MLMC estimator. We
also prove $L^{p}(\Omega)$-norm convergence of the difference between the
pairwise coupled numerical paths, the boundedness of the Radon–Nikodym
derivative, and the resulting computational complexity of the higher-order
change-of-measure MLMC method.

######  Theorem 1 (Stability of numerical solution).

Consider the SDE (1) and let Assumptions 1, 2 and 3 hold. Then for any spring
constant $S>0$, there exist constants $\kappa_{1},\kappa_{2}>0$, independent
of $T$ and $p$, such that for any $T>0$, $p\geq 1$, and for all
$0<h<\kappa_{1}$, it holds that

$\displaystyle\sup_{0\leq n\leq
N}\mathbb{E}\Big{[}\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{p}\Big{]}^{1/p}\leq\kappa_{2}p^{1/2},\qquad\sup_{0\leq
n\leq
N}\mathbb{E}\Big{[}\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{p}\Big{]}^{1/p}\leq\kappa_{2}p^{1/2}.$

We defer the proof to Section 6. Theorem 1 states that the trajectories of the
higher-order change-of-measure MLMC method are $L^{p}(\Omega)$ stable for all
sufficiently small $h>0$, and that the moments are uniformly bounded in $T$.
Note that the stability of the trajectories is not affected by the spring
term. A similar theoretical result is obtained for the numerical scheme
without the spring term

######  Theorem 2 (Stability of numerical solution without the spring term).

Consider the SDE (1) and let Assumptions 1, 2, and 3 hold. Then there exist
constants $C_{1},C_{2}>0$, independent of $T$ and $p$, such that for any
$T>0$, $p\geq 1$, and for all $0<h<C_{1}$, it holds that

$\sup_{0\leq n\leq
N}\mathbb{E}\Big{[}\left\|\overline{X}_{t_{n}}^{f}\right\|^{p}\Big{]}^{1/p}\leq
C_{2}p^{1/2},\qquad\sup_{0\leq n\leq
N/2}\mathbb{E}\Big{[}\left\|\overline{X}_{t_{2n}}^{c}\right\|^{p}\Big{]}^{1/p}\leq
C_{2}p^{1/2}.$

Using the results from Theorem 1, we obtain the $L_{p}(\Omega)$-norm bounds
for the pairwise coupled trajectories difference.

######  Theorem 3 ($L^{p}(\Omega)$-norm pairwise coupled trajectories
difference).

Consider the SDE (1) and let Assumptions 1, 2 and 3 hold. Then for any
$S>\lambda/2$, there exist constants $\kappa_{3},\kappa_{4},\kappa_{5}>0$,
independent of $T$ and $p$, such that for any $T>0$, $p\geq 1$, and for all
$0<h<\min(\kappa_{3},\kappa_{4}/\sqrt[3]{p})$, it holds that

$\sup_{0\leq n\leq
N}\Bigg{(}\mathbb{E}\Bigg{[}\left\|\widehat{Y}_{t_{n}}^{f}-\widehat{Y}_{t_{n}}^{c}\right\|^{p}\Bigg{]}\Bigg{)}^{1/p}\leq\kappa_{5}\min\Big{(}p^{1/2}h^{1/2},\;ph,\;p^{3/2}h^{3/2}\Big{)}.$

We defer the proof to Section 6. Theorem 3 states that the
$L^{p}(\Omega)$-norm of the pairwise coupled difference of the MLMC
trajectories converges uniformly in time $T$ for all sufficiently small $h>0$.
By choosing $S>\lambda/2$, we can achieve uniform-in-time bounds obtained in
Theorem 3.

######  Remark.

Note that for all $p\geq 1$ and $h\in(0,1)$, the term
$ph\geq\min(p^{1/2}h^{1/2},p^{3/2}h^{3/2})$ and hence can be removed from
Theorem 3. We have included the intermediate term as we use it recursively to
obtain the higher-order convergence rate of $\mathcal{O}(p^{3/2}h^{3/2})$.

Making use of Theorems 1 and 3, we can bound the $L^{p}(\Omega)$-norm of the
Radon–Nikodym derivatives.

######  Theorem 4 ($L^{p}(\Omega)$-norm Radon–Nikodym derivative).

Consider the SDE (1) and let Assumptions 1, 2 and 3 hold. Then for all
$S>\lambda/2$, $T\geq 1$, $p\geq 1$, and for all $0<h\leq h_{\max}(T)$, it
holds that

$\mathbb{E}\Big{[}\left|R^{f}_{T}\right|^{p}\Big{]}\leq\eta_{1},\qquad\mathbb{E}\Big{[}\left|R^{c}_{T}\right|^{p}\Big{]}\leq\eta_{2},$

where $h_{\max}(T)$ denotes the maximum permissible time-step size as a
function of the terminal time $T$ given by

$h_{\max}(T)\coloneqq\min\big{(}1,\;C_{0}/\sqrt{T\log T}\big{)},$

$C_{0},\eta_{1},\eta_{2}>0$ are constants that are independent of $T$ and $p$.

We defer the proof to Section 6. Theorem 4 shows that the Radon–Nikodym
derivative is uniformly bounded in time $T$ for all sufficiently small $h>0$.
Using Theorems 1, 3 and 4, we can obtain the bounds for the
$L^{p}(\Omega)$-norm of the higher-order change-of-measure MLMC estimator.
Before proceeding with the analysis, we consider two cases for the regularity
of the payoff function $\varPhi\colon\mathbb{R}^{d}\to\mathbb{R}$. The first
case we consider is that of a uniformly Lipschitz continuous payoff function

######  Assumption 4 (uniform Lipschitz continuity of $\varPhi(x)$).

There exists a constant $K>0$, such that

$\left|\varPhi(x)-\varPhi(y)\right|\leq K\left\|x-y\right\|,\qquad\forall
x,y\in\mathbb{R}^{d},$

where $K$ is the uniform Lipschitz constant.

Using Assumption 4, we can bound the moments of the MLMC estimator.

######  Theorem 5 (Multilevel estimator moments - uniformly Lipschitz
continuous payoff).

Consider the SDE (1) and let Assumptions 1, 2, and 3 hold. Let
$\varPhi\colon\mathbb{R}^{d}\to\mathbb{R}$ denote a function that satisfies
Assumption 4. Then for any $S>\lambda/2$, for all $T>1$, $p\geq 1$, and for
all $h\in\big{(}0,h_{\max}(T)\big{]}$, it holds that

$\mathbb{E}\Bigg{[}\left|\varPhi(\widehat{Y}_{T}^{f})R_{T}^{f}-\varPhi(\widehat{Y}_{T}^{c})R_{T}^{c}\right|^{p}\Bigg{]}^{1/p}\leq\kappa_{7}\sqrt{T}p^{5/2}h^{3/2},$

where $\kappa_{7}>0$ is a constant that is independent of $T$ and $p$.

We defer the proof to Section 6. The second case we consider is a class of
discontinuous payoff functions of the following form
$\varPhi(x)=\mathbbm{1}_{x\in G}$, for a Borel set
$G\in\mathcal{B}\big{(}\mathbb{R}^{d}\big{)}$ and

$\begin{split}\mathbbm{1}_{x\in G}\coloneqq\begin{cases}1,&\text{if }x\in G\\\
0,&\text{otherwise.}\end{cases}\end{split}$ (19)

Let $\overline{G}$ denote the closure of the set $G$, and $\partial G$ denote
the boundary of the set $G$. Note that $\partial
G=\overline{G}\setminus\text{int}\big{(}G\big{)}$, where $\text{int}(A)$
denotes the interior of the set $A\subset\mathbb{R}^{d}$. We define the
distance measure from a point $x\in\mathbb{R}^{d}$ to the boundary $\partial
G$ as

$d_{\partial G}(x)\coloneqq\inf_{y\in\partial
G}\left\|y-x\right\|,\qquad\forall x\in\mathbb{R}^{d}.$ (20)

Using the distance measure, we define the set

$\begin{split}G_{\delta}\coloneqq\Big{\\{}x\in\mathbb{R}^{d}\;\big{|}\;d_{\partial
G}(x)\leq\delta\Big{\\}}\end{split}$ (21)

that denotes the set of points $x\in\mathbb{R}^{d}$ that are distance $\delta$
away from the boundary $\partial G$. We make the following assumption on the
probability density function of the numerical process

######  Assumption 5 (Existence and boundedness of density).

For any $T>1$, $\xi>0$, and for all sufficiently small $h\in(0,h_{\max}(T)]$,
there exists a probability density function $\rho_{\widehat{Y}_{T}^{f}}(x)$
associated to the numerical process $\widehat{Y}_{T}^{f}$, such that

$\sup_{T>1}\sup_{h\in(0,h_{\max}(T)]}\sup_{x\in\mathbb{R}^{d}}\rho_{\widehat{Y}_{T}^{f}}(x)\leq
C_{\rho}<\infty,$

where $C_{\rho}>0$ is a constant that is independent of $h$ and $T$.

Having made the necessary assumptions, we now state the following lemma that
ensures that the distance measure is bounded

######  Lemma 1 (Probability of distance measure).

For all $\delta>0$, for any $T>1$, and for all sufficiently small
$h\in\big{(}0,h_{\max}(T)\big{]}$, assuming that there exists a constant
$C_{v}>0$ such that

$\operatorname{vol}(G_{\delta})\coloneqq\int_{G_{\delta}}\operatorname{d\\!}x\leq
C_{v}\delta$

holds, then there exists another constant $C_{b}>0$ that is independent of
$\delta$, $T$, and $h$ such that

$\mathbb{P}\Big{(}d_{\partial G}(\widehat{Y}_{T}^{f})\leq\delta\Big{)}\leq
C_{b}\delta.$

###### Proof.

Note that the event

$\Big{\\{}\omega\in\Omega\;\big{|}\;d_{\partial
G}(\widehat{Y}_{T}^{f}(\omega))\leq\delta\Big{\\}}=\Big{\\{}\omega\in\Omega\;\big{|}\;\widehat{Y}_{T}^{f}(\omega)\in
G_{\delta}\Big{\\}}$

Using the assumption that the volume of the set $G_{\delta}$ is bounded, we
obtain

$\begin{split}\mathbb{P}\Big{(}\widehat{Y}_{T}^{f}\in
G_{\delta}\Big{)}=\mathbb{E}\Big{[}\mathbbm{1}_{\widehat{Y}_{T}^{f}\in
G_{\delta}}\Big{]}&=\int_{\mathbb{R}^{d}}\mathbbm{1}_{x\in
G_{\delta}}\rho_{\widehat{Y}_{T}^{f}}(x)\operatorname{d\\!}x\\\
&=\int_{G_{\delta}}\rho_{\widehat{Y}_{T}^{f}}(x)\operatorname{d\\!}x\\\
&\leq\int_{G_{\delta}}\sup_{T>1}\sup_{h\in(0,h_{\max}(T)]}\sup_{x\in\mathbb{R}^{d}}\rho_{\widehat{Y}_{T}^{f}}(x)\operatorname{d\\!}x\\\
&\leq
C_{\rho}\Big{(}\underbrace{\int_{G_{\delta}}\operatorname{d\\!}x}_{\eqqcolon\operatorname{vol}(G_{\delta})}\Big{)}\leq\underbrace{C_{\rho}C_{v}}_{\eqqcolon
C_{b}}\delta.\end{split}$

∎

We now extend Theorem 5 to the case of a discontinuous payoff function defined
by equation (19).

######  Theorem 6 (Multilevel estimator moments - discontinuous payoff).

Consider the SDE (1) and let Assumptions 1, 2, and 3 hold. Let
$\varPhi\colon\mathbb{R}^{d}\to\mathbb{R}$ denote the payoff function given by
equation (19). Then for all $S>\lambda/2$, $T>1$, $p,q\geq 1$,
$h\in\big{(}0,h_{\max}(T)\big{]}$, and $\xi>0$, it holds that

$\mathbb{E}\Bigg{[}\left|\varPhi(\widehat{Y}_{T}^{f})R_{T}^{f}-\varPhi(\widehat{Y}_{T}^{c})R_{T}^{c}\right|^{p}\Bigg{]}\leq\kappa_{8}T^{p/2}p^{2p}(1+q)q^{\frac{1}{2}\frac{q}{q+1}}h^{\frac{3}{2}-\xi},$

where $\kappa_{8}>0$ is a constant that is independent of $T$, $p$, and $q$.

We defer the proof to Section 6. From Theorems 5 and 6 we have that the
variance of the MLMC estimator is $\mathcal{O}(Th^{3})$ and
$\mathcal{O}(Th^{3/2-\xi})$ for uniformly Lipschitz continuous and
discontinuous payoff functions, respectively. Having established all the
necessary theoretical results, we now determine the computational cost of
implementing the higher-order change-of-measure MLMC method.

######  Theorem 7 (MLMC computational complexity).

Let $\widehat{\varPhi}$ denote the MLMC estimator corresponding to the higher-
order change-of-measure scheme

$\begin{split}\widehat{\varPhi}&\coloneqq\frac{1}{N_{0}}\sum_{k=1}^{N_{0}}\big{(}\varPhi(\overline{X}_{T}^{0})\big{)}_{k}+\sum_{\ell=1}^{L}\frac{1}{N_{\ell}}\sum_{k=1}^{N_{\ell}}\big{(}\varPhi(\widehat{Y}_{T}^{f,\ell})R_{T}^{f,\ell}-\varPhi(\widehat{Y}_{T}^{c,\ell})R_{T}^{c,\ell}\big{)}_{k},\end{split}$

where $N_{\ell}$ denotes the number of Monte Carlo samples on each level of
the MLMC estimator. Let Assumptions 1, 2, and 3 hold. Then for any
sufficiently small error tolerance $\epsilon>0$, and for a suitable choice of
$L,T,N_{\ell}$, and $h_{0}$, the mean square error (MSE) of the estimator
$\widehat{\varPhi}$ is bounded as

$\mathbb{E}\Big{[}\big{(}\widehat{\varPhi}-\mathbb{E}[\varPhi(X_{\infty})]\big{)}^{2}\Big{]}\leq\epsilon^{2}$

and the expected total computational cost of the MLMC estimator is

$\operatorname{Cost}=\mathcal{O}\big{(}\epsilon^{-2}\left|\log\epsilon\right|^{3/2}(\log\left|\log\epsilon\right|)^{1/2}\big{)},$

for all QoI that satisfy Assumption 4. The expected total computational cost
is

$\operatorname{Cost}=\mathcal{O}\big{(}\epsilon^{-2}\left|\log\epsilon\right|^{5/3+\xi}\big{)},$

for all $\xi>0$, QoI of the form given by equation (19), and satisfying
Assumption 5.

We defer the proof to Section 6. The computational cost argument similarly
applies to the Milstein scheme. Since the Radon–Nikodym term for the Milstein
scheme [9, equation (13)] shares a similar structure with the higher-order
scheme given by equations (14) and (17), and considering that we have not
utilized the $\mathcal{O}(h^{3/2})$ convergence rate in deriving Theorem 4, we
believe the results can be extended to the Milstein scheme, yielding the
condition $h_{0}=\mathcal{O}(1/\sqrt{T\log T})$. Assuming that this is the
case, for a uniformly Lipschitz continuous payoff function, recall that

$\begin{split}\left|\mathbb{E}\big{[}\varPhi(\widehat{Y}_{T}^{f,\ell})R_{T}^{f,\ell}-\varPhi(\widehat{Y}_{T}^{c,\ell})R_{T}^{c,\ell}\big{]}\right|&\leq
C_{1}\sqrt{T}(2^{-\ell}h_{0})\\\
\mathbb{V}\Big{(}\varPhi(\widehat{Y}_{T}^{f,\ell})R_{T}^{f,\ell}-\varPhi(\widehat{Y}_{T}^{c,\ell})R_{T}^{c,\ell}\Big{)}&\leq
c_{1}T(2^{-\ell}h_{0})^{2}\end{split},\qquad\forall\ell\geq 1,$

where $c_{1},C_{1}>0$ are constants. Then for the choice of
$h_{0}=\mathcal{O}(1/\sqrt{T\log T})$ and

$L=\left\lceil\log_{2}\Big{(}\epsilon^{-1}\left|\log
T\right|^{-1/2}\Big{)}+\log_{2}(\sqrt{6}\widetilde{C}_{1})\right\rceil,$

where $\widetilde{C}_{1}>0$ is a constant, the computational cost is

$\operatorname{Cost}=\mathcal{O}\Big{(}\epsilon^{-2}\left|\log\epsilon\right|^{3/2}(\log\left|\log\epsilon\right|)^{1/2}\Big{)}.$

For a discontinuous payoff function of the form given by equation (19) and for
all $\xi>0$, we can show that the strong error rate and the variance decay
rate for the Milstein scheme are bounded as

$\begin{split}\left|\mathbb{E}\big{[}\varPhi(\widehat{Y}_{T}^{f,\ell})R_{T}^{f,\ell}-\varPhi(\widehat{Y}_{T}^{c,\ell})R_{T}^{c,\ell}\big{]}\right|&\leq
C_{1}\sqrt{T}(2^{-\ell}h_{0})^{1-\xi}\\\
\mathbb{V}\Big{(}\varPhi(\widehat{Y}_{T}^{f,\ell})R_{T}^{f,\ell}-\varPhi(\widehat{Y}_{T}^{c,\ell})R_{T}^{c,\ell}\Big{)}&\leq
c_{1}T(2^{-\ell}h_{0})^{1-\xi}\end{split},\qquad\forall\ell\geq 1,$

where $c_{1},C_{1}>0$ are constants. Then for the choice of
$h_{0}=\mathcal{O}\Big{(}1/\sqrt{T\log T}\Big{)}$ and

$L=\left\lceil\frac{1}{1-\xi}\log_{2}\Big{(}\epsilon^{-1}\sqrt{T}h_{0}^{1-\xi}\Big{)}+\frac{1}{1-\xi}\log_{2}(\sqrt{6}\widetilde{C}_{1})\right\rceil,$

where $\widetilde{C}_{1}>0$ is a constant, all $\xi>0$, and
$\hat{\xi}=\frac{\xi}{1-\xi}$, we obtain

$\operatorname{Cost}=\mathcal{O}\Big{(}\epsilon^{-(2+\hat{\xi})}\left|\log\epsilon\right|^{2+\hat{\xi}/2}\Big{)}.$

##  4\. Numerical results

In this section, we run several numerical simulations to verify the results
obtained in Section 2. Algorithm 1 describes the implementation of the higher-
order change-of-measure scheme to compute the weak approximations of the
invariant measure.

Initial condition $x_{0}$, Terminal time $T$, # of multilevel estimator levels
$L$, # of Monte Carlo samples $M$, Time-step size on level–0 $h_{0}$.
$i=1,\ldots,M$

1:Initialize the fine and the coarse trajectories at time $t=0$, i.e.
$\displaystyle\widehat{Y}_{0}^{f}=\widehat{Y}_{0}^{c}=x_{0}$.

2:Set the initial values for the Radon–Nikodym derivatives, i.e.
$R_{0}^{f}=R_{0}^{c}=1$. $n=0,\ldots,N/2-1$ Updating the state of the process
from $t_{2n}$ to $t_{2n+1}$

3:Using equation (11), generate the tuple of correlated random increments
$(\Delta W_{2n},\Delta Z_{2n})$.

4:Update the fine and the coarse trajectories using equations (7) and (9).

5:Compute the Radon–Nikodym derivative $R^{f}_{2n}$ for the fine trajectory
using equation (14), i.e.
$\begin{multlined}R_{2n+1}^{f}=R_{2n}^{f}\cdot\exp\Bigg{(}-\left\langle\widehat{S}^{f}_{2n},\Delta
V_{1,2n}\right\rangle+\sqrt{3}\left\langle\widehat{S}^{f}_{2n},\Delta
V_{2,2n}\right\rangle-2h\left\|\widehat{S}^{f}_{2n}\right\|^{2}\Bigg{)}\end{multlined}R_{2n+1}^{f}=R_{2n}^{f}\cdot\exp\Bigg{(}-\left\langle\widehat{S}^{f}_{2n},\Delta
V_{1,2n}\right\rangle+\sqrt{3}\left\langle\widehat{S}^{f}_{2n},\Delta
V_{2,2n}\right\rangle-2h\left\|\widehat{S}^{f}_{2n}\right\|^{2}\Bigg{)}$
Updating the state of the process from $t_{2n+1}$ to $t_{2n+2}$

6:Using equation (11), generate the tuple of correlated random increments
$(\Delta W_{2n+1},\Delta Z_{2n+1})$.

7:Update the fine and the coarse trajectories using equations (8) and (10).

8:Compute the Radon–Nikodym derivatives, $R_{2n+1}^{f}$ and $R_{2n+1}^{c}$,
for the fine and the coarse trajectories, respectively, using equations (14)
and (17), i.e.
$\begin{split}R_{2n+2}^{f}&=R_{2n+1}^{f}\cdot\exp\bigg{(}-\left\langle\widehat{S}^{f}_{2n+1},\Delta
V_{1,2n+1}\right\rangle\\\
&\qquad\qquad\qquad\qquad+\sqrt{3}\left\langle\widehat{S}^{f}_{2n+1},\Delta
V_{2,2n+1}\right\rangle-2h\left\|\widehat{S}^{f}_{2n+1}\right\|^{2}\bigg{)}\\\
~{}R_{2n+2}^{c}&=R_{2n}^{c}\cdot\exp\bigg{(}-\left\langle\widehat{S}^{c}_{2n},\Delta
V_{1,2n}+\Delta
V_{1,2n+1}\right\rangle+2\sqrt{3}\left\langle\widehat{S}^{c}_{2n},\Delta
V_{2,2n}+\Delta V_{2,2n+1}\right\rangle\\\
&\qquad\qquad\qquad\qquad-13h\left\|\widehat{S}^{c}_{2n}\right\|^{2}\bigg{)}\end{split}$
For a given payoff function, $\varPhi\colon\mathbb{R}^{d}\to\mathbb{R}$, the
expected value of the payoff w.r.t. the invariant measure, i.e.
$\mathbb{E}\Bigg{[}\left|\varPhi\big{(}\widehat{Y}_{T}^{f}\big{)}R_{T}^{f}-\varPhi\big{(}\widehat{Y}_{T}^{c}\big{)}R_{T}^{c}\right|^{p}\Bigg{]}^{1/p}$

Algorithm 1 Higher-order change-of-measure algorithm

###  4.1. 1D Triple-well potential – Discontinuous payoff

We start with the numerical tests by considering the 1D SDE

$\begin{split}\operatorname{d\\!}X&=\frac{X^{3}(2-X^{2})(X^{8}+2X^{6}+4X^{2}-4)}{2(X^{6}+1)^{2}}\operatorname{d\\!}t+\operatorname{d\\!}W\\\
X(0)&=1,\end{split}$

corresponding to the potential function
$f(X)\coloneqq(X^{4}-2X^{2})^{2}/(4(X^{6}+1))$. The invariant probability
density of the SDE is

$\rho(X)\propto\exp(-2f(x))=\exp\Big{(}-\frac{(X^{4}-2X^{2})^{2}}{2(X^{6}+1)}\Big{)},$

and is as depicted in Figure 1. Note that the drift coefficient satisfies
Assumptions 1, 2, and 3. We are interested in computing the weak approximation
of the payoff function $\varPhi(x)\coloneqq\mathbbm{1}_{x\in[0,2]}$ w.r.t. the
invariant measure of the SDE for the Double-well potential. The reference
solution to the weak approximation of the payoff function w.r.t the invariant
measure was computed by solving the elliptic Fokker–Planck PDE [33, equation
(8.24)], using the continuous, piecewise quadratic finite elements for the
spatial discretization. For this purpose, we use the Fenics library in Python
programming language. The pseudo-reference solution is
$\mathbb{E}[\varPhi(X_{\infty})]=0.42863$, rounded to 5 significant digits. We
run the simulation for $M=10^{8}$ Monte Carlo samples, until the terminal time
$T=10$, to compute the convergence rate for the higher-order numerical scheme.
We use an initial time-step size $h_{0}=2^{-5}$.

Figure 1. Plots associated to the SDE considered in Section 4.1: Probability
density function of the triple-well potential (top left), Curve fit to compute
the parameters $\mu^{\star}$ and $\lambda^{\star}$ (top right), Variance
growth rate of the MLMC estimator w.r.t terminal time on level $\ell=1$
(bottom).

Figure 1 shows that the variance of the MLMC estimator increases linearly with
time $T$ for a fixed time-step size $h$, as shown in Theorem 6. Figure 2 shows
that the convergence rates obtained from the numerical simulation are
consistent with theory. We observe that both the strong error and the variance
rate for the higher-order change-of-measure scheme are $\mathcal{O}(h^{3/2})$
as proved in Theorem 6. Recall that the set
$\Theta_{T}\coloneqq\big{\\{}\omega\in\Omega\;\big{|}\;\left\|\widehat{Y}_{T}^{f}-\widehat{Y}_{T}^{c}\right\|\geq\nu_{1}\left|\log
h\right|\big{\\}}$ is as defined in Theorem 4. From Figure 2, we also observe
that the probability of the fine and the coarse trajectories diverging
$\mathbb{P}(\Theta_{T})=0$ for $\nu_{1}=1$. We observe that the kurtosis
increases with the level of the multilevel estimator. This can be attributed
to our computing the weak approximation of a discontinuous payoff function.

Figure 2. Numerical results for the problem considered in Section 4.1 at the
terminal time, $T=10$: strong error (top left), variance rate (top right),
kurtosis (bottom left), divergence probability (bottom right)

The numerical experiment matches the theoretical results as shown in Figure 3.
For the choice of $h_{0}=\mathcal{O}(T^{-1/(3/2-\xi)})$, for a discontinuous
QoI, and using the higher-order change-of-measure scheme, we can achieve a
computational cost of
$\mathcal{O}\big{(}\left|\log\epsilon\right|^{5/3+\xi_{1}}\big{)}$, for all
$\xi_{1}>0$, while still achieving an MSE of $\mathcal{O}(\epsilon^{2})$. As
for the Milstein scheme developed in [9], we observe that to achieve an MSE of
$\mathcal{O}(\epsilon^{2})$, the computational cost is
$\mathcal{O}(\epsilon^{-\xi_{2}}\left|\log\epsilon\right|^{2+\xi_{2}/2})$ for
the choice of $h_{0}=\mathcal{O}(1/\sqrt{T\log T})$ and all $\xi_{2}>0$.
Although we use a larger initial time-step size for the Milstein scheme, the
associated computational cost is lower for the higher-order scheme compared to
the Milstein scheme.

Figure 3. Change-of-measure MLMC scheme results for the problem considered in
Section 4.1: computational cost (top left), mean squared error (top right),
asymptotic computational cost (bottom)

###  4.2. 2D potential well – Discontinuous payoff

We will now consider the SDE

$\begin{split}\operatorname{d\\!}X_{1}&=\Bigg{(}X_{1}\bigg{(}\frac{4}{(1+X_{1}^{2})^{2}}-2\bigg{)}+\frac{1}{2}\operatorname{sech}^{2}(X_{1})\tanh(X_{2})\Bigg{)}\operatorname{d\\!}t+\operatorname{d\\!}W^{1}\\\
\operatorname{d\\!}X_{2}&=\Bigg{(}X_{2}\bigg{(}\frac{4}{(1+X_{2}^{2})^{2}}-2\bigg{)}+\frac{1}{2}\operatorname{sech}^{2}(X_{2})\tanh(X_{1})\Bigg{)}\operatorname{d\\!}t+\operatorname{d\\!}W^{2}\\\
X(0)&=(0,0)^{\top}\end{split}$

corresponding to the potential function

$f(X_{1},X_{2})\coloneqq\frac{X_{1}^{4}+1}{X_{1}^{2}+1}+\frac{X_{2}^{4}+1}{X_{2}^{2}+1}-\frac{1}{2}\tanh(X_{1})\tanh(X_{2})-2$

that models the dynamics of a particle under the influence of a potential
well. The drift coefficient satisfies Assumptions 1, 2, and 3. The drift
coefficient exhibits non-linearity close to the origin. We are interested in
computing the weak approximation of the payoff function

$\varPhi(X_{1},X_{2})\coloneqq\begin{cases}1,&\text{if }0\leq X_{1}+X_{2}\leq
1.4\text{ and }-0.75\leq X_{2}-X_{1}\leq 0.75\\\
0,&\text{otherwise}\end{cases},$

w.r.t. the invariant measure of the SDE. The pseudo-reference solution is
$\mathbb{E}[\varPhi(X_{\infty})]=0.1674$, rounded to 4 significant digits. We
run the simulation until the terminal time $T=10$, using an initial time-step
size $h_{0}=2^{-4}$, and $M=10^{7}$ Monte Carlo samples.

Figure 4. The logarithm of the probability density function of the SDE
considered in Section 4.2.

Figure 4 depicts the logarithmic probability density of the 2D potential well
and shows that the variance of the MLMC estimator increases linearly with time
$T$ for a fixed time-step size $h$, as shown in Theorem 6. Figure 5 shows that
the convergence rates obtained from the numerical simulation are consistent
with theory. We observe that both the strong error and the variance rate for
the higher-order change-of-measure scheme are $\mathcal{O}(h^{3/2})$ as proved
in Theorem 6. From Figure 2, we also observe that the probability of the fine
and the coarse trajectories diverging $\mathbb{P}(\Theta_{T})=0$ for
$\nu_{1}=1$. We observe that the kurtosis increases with the level of the
multilevel estimator. This can be attributed to computing the weak
approximation of a discontinuous payoff function.

Figure 5. Numerical results for the problem considered in Section 4.2 at the
terminal time, $T=10$: strong error (top left), variance rate (top right),
kurtosis (bottom left), divergence probability (bottom right)

From Figure 6, it is observed that the numerical results match the theoretical
results. By choosing $h_{0}=\mathcal{O}(1/\sqrt{T\log T})$ and
$h_{0}=\mathcal{O}\big{(}T^{-1/(3/2-\xi_{1})}\big{)}$ for the Milstein and the
higher-order schemes, respectively, for all $\xi_{1},\xi_{2}>0$, it is noted
that the numerical results match the theoretical results. Furthermore, the
computational cost of the higher-order scheme is lower than the Milstein
scheme in terms of both the rate and the associated constant. We also note
that the results extend to 2D SDE as shown in theory.

Figure 6. Change-of-measure MLMC scheme results for the problem considered in
Section 4.2: computational cost (top left), mean squared error (top right),
asymptotic computational cost (bottom)

###  4.3. Thomas’ cyclically symmetric attractor – Lipschitz payoff

In [32], a cyclically symmetric attractor system is formulated, which like the
Lorenz system, is a chaotic attractor system. To this end, we consider the SDE

$\begin{split}\operatorname{d\\!}X_{1}&=\big{(}\sin(X_{2})-0.18X_{1}\big{)}\operatorname{d\\!}t+\operatorname{d\\!}W^{1}\\\
\operatorname{d\\!}X_{2}&=\big{(}\sin(X_{3})-0.18X_{2}\big{)}\operatorname{d\\!}t+\operatorname{d\\!}W^{2}\\\
\operatorname{d\\!}X_{3}&=\big{(}\sin(X_{1})-0.18X_{3}\big{)}\operatorname{d\\!}t+\operatorname{d\\!}W^{3}\\\
X(0)&=(1,2,2)^{\top}.\end{split}$

Note that the drift term given above satisfies Assumptions 1, 3, and 2. We are
interested in computing the weak approximation of the uniformly Lipschitz
continuous payoff function $\varPhi(x)\coloneqq\left\|x\right\|$ w.r.t. the
invariant measure of the SDE. The pseudo-reference solution is
$\mathbb{E}[\varPhi(X_{\infty})]=3.9664$, rounded to 5 significant digits. We
run $M=10^{6}$ Monte Carlo simulations until the terminal time $T=40$ using an
initial time-step size $h_{0}=2^{-4}$.

Figure 7. Numerical results for the problem considered in Section 4.3 at the
terminal time, $T=40$: Strong error (top left), Variance rate (top right),
Kurtosis (bottom left), and the expected total computational cost (bottom
right)

From Figure 7, we observe that the numerical rates are consistent with the
theoretical rates for a uniformly Lipschitz continuous payoff function. Unlike
for the case of a discontinuous payoff, we observe that the kurtosis of both
the numerical schemes is bounded when considering a uniformly Lipschitz
continuous payoff. This leads to a better estimation of the sample variance on
all levels of the multilevel estimator. We also observe the probability that
the fine and the coarse trajectories diverge $\mathbb{P}(\Theta_{T})=0$ for
$\nu_{1}=1$.

Figure 8. Change-of-measure MLMC scheme results for the problem considered in
Section 4.3: computational cost (top left), mean squared error (top right),
asymptotic computational cost (bottom left), variance growth rate of the MLMC
estimator w.r.t terminal time on level $\ell=1$ (bottom right).

From Figure 8 we observe that for $h_{0}=2^{-4}$, the variance grows linearly
w.r.t. the terminal time $T$ which is consistent with the theoretical results
in Theorem 5. The numerical experiments show that the computational cost of
both the higher-order and the Milstein schemes are bounded from above by
$\mathcal{O}\big{(}\epsilon^{-2}\left|\log\epsilon\right|^{3/2}(\log\left|\log\epsilon\right|)^{1/2}\big{)}$.
This is because the majority of the computational cost is attributed to
generating samples on level zero of the MLMC estimator for both numerical
schemes. Although both numerical schemes have the same computational cost
order, the computational cost constant associated with the higher-order
numerical scheme is lower than the Milstein scheme to achieve an MSE of
$\mathcal{O}(\epsilon^{2})$.

##  5\. Conclusion and future scope

In this paper, we developed a tractable higher-order numerical scheme to
compute weak approximations of the invariant measure of an SDE using the MLMC
method, that does not satisfy the contractivity condition. The fundamental
idea is similar to one in [9], where we introduce a spring term in the 1.5
strong Itô–Taylor scheme that pulls the fine and the coarse trajectories of
the multilevel estimator close to each other. This ensures that the fine and
the coarse trajectories do not drift far away from each other while recovering
the contractivity property of the SDE needed for the efficient implementation
of MLMC method. We then use Girsanov’s transformation to compute the weak
approximations w.r.t. the new measure.

We theoretically proved that the fine and the coarse trajectories of the
multilevel estimator are $L^{p}(\Omega)$ stable. We proved that the pairwise
coupled difference between the fine and the coarse trajectories of the
multilevel estimator has $L^{p}(\Omega)$ convergence rate
$\mathcal{O}\big{(}\min(p^{1/2}h^{1/2},ph,p^{3/2}h^{3/2})\big{)}$ for all
$p\geq 1$, and for all sufficiently small $h>0$. We proved the
$L^{p}(\Omega)$-norm convergence of the multilevel estimator to be
$\mathcal{O}\big{(}T^{1/2}p^{5/2}h^{3/2}\big{)}$ for a uniformly Lipschitz
continuous payoff function. For a discontinuous payoff function defined by
equation 19, the $p$-th moment of the multilevel estimator was proved to be
$\mathcal{O}\Big{(}T^{p/2}p^{2p}h^{\frac{3}{2}-\xi}\Big{)}$ for all $p\geq 1$
and for all $\xi>0$. Provided that the SDE under consideration satisfies
Assumptions 1, 2, and 3, we proved that the computational cost involved in the
implementation of the higher-order MLMC scheme is
$\mathcal{O}\big{(}\epsilon^{-2}\left|\log\epsilon\right|^{3/2}(\log\left|\log\epsilon\right|)^{1/2}\big{)}$,
and
$\mathcal{O}\big{(}\epsilon^{-2}\left|\log\epsilon\right|^{5/3+\xi}\big{)}$,
where $\xi>0$, for all uniformly Lipschitz continuous payoffs, and all
discontinuous payoffs of the form given by equation (19), respectively.

There are several ways in which the current work can be extended. Firstly, we
can extend the analysis of the higher-order change-of-measure scheme to SDE
with state-dependent diffusion coefficients. Another possible extension would
be to use adaptive time-stepping based on the idea devised in [6, 7, 8]. By
choosing time-steps adapted to the state of the process, and using the new
higher-order change-of-measure scheme, we can extend the computation of weak
approximations of invariant measure of SDE with non-uniformly Lipschitz
continuous drift coefficients. It is important to note that the coupling
between the coarse and the fine trajectories and the computation of the
corresponding Radon–Nikodym derivatives for the higher-order change-of-measure
scheme are not straightforward unlike the Milstein scheme due to the higher-
order noise terms. Further analysis is needed to develop an optimal coupling
and to compute the corresponding Radon–Nikodym derivatives. This poses an
interesting problem for future research.

The authors of [20] construct a higher-order adaptive method for strong
approximations of exit times of Itô stochastic differential equations (SDE).
The fundamental idea is to use a smaller time-step size as the numerical
solution gets closer to the boundary of the domain, and to use higher-order
integration schemes for better approximation of the state of the Itô-
diffusion. The authors prove that by using a higher-order order numerical
integrator along with adaptive time-steps leads to a better strong error
convergence rate of exit-time computations. Combining the ideas of adaptive
timestepping for exit-time problems with the higher-order change-of-measure
scheme may help in computing weak approximations of quantities of interest
(QoI) that depend on the exit-time and the state of the process. This makes an
interesting problem for future research.

##  6\. Theoretical proofs

#### Proof of Theorem 1

We adopt a methodology similar to one in [9] to prove the stability of the
numerical trajectories with the spring constant. To keep the notation concise,
we denote the tuple of random variables $(\Delta W_{n}^{\mathbb{P}},\Delta
Z_{n}^{\mathbb{P}})=(\Delta W_{n},\Delta Z_{n})$ for all $n=1,\ldots,N-1$. We
introduce the constants $C,\widetilde{\alpha},\widetilde{\beta}>0$ that depend
only on $S,d$, and the constants from Assumptions 1, 2, and 3. The constants
$C,\widetilde{\alpha},\widetilde{\beta}$ may change from line to line in the
proof argument.

Computing the square of the Euclidean norm of the fine and the coarse
trajectories at the odd time-step, we obtain

$\begin{split}\left\|\widehat{Y}_{2n+1}^{f}\right\|^{2}&=\left\|Sh\widehat{Y}_{2n}^{c}+(1-Sh)\bigg{\\{}\widehat{Y}_{2n}^{f}+\frac{A_{2n}^{f}}{1-Sh}\bigg{\\}}\right\|^{2},\\\
\left\|\widehat{Y}_{2n+1}^{c}\right\|^{2}&=\left\|Sh\widehat{Y}_{2n}^{f}+(1-Sh)\bigg{\\{}\widehat{Y}_{2n}^{c}+\frac{A_{2n}^{c}}{1-Sh}\bigg{\\}}\right\|^{2},\end{split}$

where

$\begin{split}A_{2n}^{f}\coloneqq ha(\widehat{Y}_{2n}^{f})+\Delta
W_{2n}+Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}_{2n}^{f}),\\\
A_{2n}^{c}\coloneqq ha(\widehat{Y}_{2n}^{c})+\Delta
W_{2n}+Da(\widehat{Y}_{2n}^{c})\Delta
Z_{2n}+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}_{2n}^{c}).\end{split}$

Note that for the given spring constant $S>0$, and for all sufficiently small
$h>0$,

$\frac{1}{1-Sh}\leq 2.$

Since the Euclidean norm squared is a convex function, from Jensen’s
inquality, we have

$\displaystyle\left\|\theta
x_{1}+(1-\theta)x_{2}\right\|^{2}\leq\theta\left\|x_{1}\right\|^{2}+(1-\theta)\left\|x_{2}\right\|^{2},\qquad\forall
x_{1},x_{2}\in\mathbb{R}^{d},\;\;\theta\in(0,1).$

Using this inequality, we obtain

$\begin{split}\left\|\widehat{Y}_{2n+1}^{f}\right\|^{2}&\leq
Sh\left\|\widehat{Y}_{2n}^{c}\right\|^{2}+(1-Sh)\left\|\widehat{Y}_{2n}^{f}+\frac{A_{2n}^{f}}{1-Sh}\right\|^{2}\\\
&\leq
Sh\left\|\widehat{Y}_{2n}^{c}\right\|^{2}+(1-Sh)\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+2\left\|A_{2n}^{f}\right\|^{2}+2\left\langle\widehat{Y}_{2n}^{f},A_{2n}^{f}\right\rangle.\end{split}$
(22)

Let $b_{i}\in\mathbb{R}^{d}$ for $i=1,2,\ldots,m$. From Hölder’s inequality,
we then have

$\left\|\sum_{i=1}^{m}b_{i}\right\|^{p}\leq
m^{p-1}\Biggl{(}\sum_{i=1}^{m}\left\|b_{i}\right\|^{p}\Biggr{)},\qquad\forall
p\geq 1.$

Using the above inequality, we can bound $\left\|A_{2n}^{f}\right\|^{2}$ as

$\begin{split}\left\|A_{2n}^{f}\right\|^{2}&\leq
4\bigg{\\{}h^{2}\left\|a(\widehat{Y}_{2n}^{f})\right\|^{2}+\left\|\Delta
W_{2n}\right\|^{2}+\left\|Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}\right\|^{2}+\frac{h^{4}}{4}\left\|\mathcal{A}a(\widehat{Y}_{2n}^{f})\right\|^{2}\bigg{\\}}\\\
~{}&\leq
C\bigg{\\{}h^{2}\Big{(}1+\left\|\widehat{Y}_{2n}^{f}\right\|^{2}\Big{)}+\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}\bigg{\\}}.\end{split}$
(23)

Using Young’s inequality, Cauchy-Schwarz inequality, and Assumption 3, the
inner product term $\left\langle\widehat{Y}_{2n}^{f},A_{2n}^{f}\right\rangle$
can be bounded as

$\begin{split}\left\langle\widehat{Y}_{2n}^{f},A_{2n}^{f}\right\rangle&=h\left\langle\widehat{Y}_{2n}^{f},a(\widehat{Y}_{2n}^{f})\right\rangle+\left\langle\widehat{Y}_{2n}^{f},\Delta
W_{2n}\right\rangle+\left\langle\widehat{Y}_{2n}^{f},Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}\right\rangle\\\
&\quad+\frac{h^{2}}{2}\left\langle\widehat{Y}_{2n}^{f},\mathcal{A}a(\widehat{Y}_{2n}^{f})\right\rangle\\\
~{}&\leq-\widetilde{\alpha}h\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\langle\widehat{Y}_{2n}^{f},\Delta
W_{2n}\right\rangle+\left\langle\widehat{Y}_{2n}^{f},Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}\right\rangle+\widetilde{\beta}h.\end{split}$ (24)

Combining equations (23) and (24) with (22), we obtain

$\begin{split}\left\|\widehat{Y}_{2n+1}^{f}\right\|^{2}&\leq
Sh\left\|\widehat{Y}_{2n}^{c}\right\|^{2}+(1-Sh-2\widetilde{\alpha}h)\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+C\Big{\\{}\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}\Big{\\}}\\\
&\quad+2\left\langle\widehat{Y}_{2n}^{f},\Delta
W_{2n}\right\rangle+2\left\langle\widehat{Y}_{2n}^{f},Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}\right\rangle+2\widetilde{\beta}h.\end{split}$ (25)

Similarly, we have the following upper bound for the coarse path at the odd
time-step

$\begin{split}\left\|\widehat{Y}_{2n+1}^{c}\right\|^{2}&\leq
Sh\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+(1-Sh-2\widetilde{\alpha}h)\left\|\widehat{Y}_{2n}^{c}\right\|^{2}+C\Big{\\{}\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}\Big{\\}}\\\
&\quad+2\left\langle\widehat{Y}_{2n}^{c},\Delta
W_{2n}\right\rangle+2\left\langle\widehat{Y}_{2n}^{c},Da(\widehat{Y}_{2n}^{c})\Delta
Z_{2n}\right\rangle+2\widetilde{\beta}h.\end{split}$ (26)

Using the similar approach as above, we can write the square of the Euclidean
norm for the fine and the coarse paths at the even time-step as

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}\right\|^{2}&=\left\|Sh\widehat{Y}_{2n+1}^{c}+(1-Sh)\bigg{\\{}\widehat{Y}_{2n+1}^{f}+\frac{A_{2n+1}^{f}}{1-Sh}\bigg{\\}}\right\|^{2},\\\
~{}\left\|\widehat{Y}_{2n+2}^{c}\right\|^{2}&=\left\|2Sh\widehat{Y}_{2n}^{f}+(1-2Sh)\bigg{\\{}\widehat{Y}_{2n}^{c}+\frac{A_{2n+1}^{c}}{1-2Sh}\bigg{\\}}\right\|^{2},\end{split}$

where

$\begin{split}A_{2n+1}^{f}&\coloneqq ha(\widehat{Y}_{2n+1}^{f})+\Delta
W_{2n+1}+Da(\widehat{Y}_{2n+1}^{f})\Delta
Z_{2n+1}+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}_{2n+1}^{f}),\\\
~{}A_{2n+1}^{c}&\coloneqq 2ha(\widehat{Y}_{2n}^{c})+(\Delta W_{2n}+\Delta
W_{2n+1})\\\ &\quad+Da(\widehat{Y}_{2n}^{c})(\Delta Z_{2n+1}+\Delta
Z_{2n}+h\Delta W_{2n})+2h^{2}\mathcal{A}a(\widehat{Y}_{2n}^{c}).\end{split}$

Using Young’s inequality, Cauchy-Schwarz inequality, and Assumption 3, we
obtain

$\displaystyle\begin{split}\left\|\widehat{Y}_{2n+2}^{f}\right\|^{2}&\leq
Sh\left\|\widehat{Y}_{2n+1}^{c}\right\|^{2}+(1-Sh-2\widetilde{\alpha}h)\left\|\widehat{Y}_{2n+1}^{f}\right\|^{2}\\\
&\quad+C\Big{\\{}\left\|\Delta W_{2n+1}\right\|^{2}+\left\|\Delta
Z_{2n+1}\right\|^{2}\Big{\\}}+2\left\langle\widehat{Y}_{2n+1}^{f},\Delta
W_{2n+1}\right\rangle\\\
&\quad+2\left\langle\widehat{Y}_{2n+1}^{f},Da(\widehat{Y}_{2n+1}^{f})\Delta
Z_{2n+1}\right\rangle+2\widetilde{\beta}h,\\\ \end{split}$ (27)
$\displaystyle\begin{split}\left\|\widehat{Y}_{2n+2}^{c}\right\|^{2}&\leq
2Sh\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+(1-2Sh-4\widetilde{\alpha}h)\left\|\widehat{Y}_{2n}^{c}\right\|^{2}+C\Big{\\{}\left\|\Delta
Z_{2n+1}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}\\\ &\quad+\left\|\Delta
W_{2n+1}\right\|^{2}+\left\|\Delta
W_{2n}\right\|^{2}\Big{\\}}+2\left\langle\widehat{Y}_{2n}^{c},\Delta
W_{2n}\right\rangle+2\left\langle\widehat{Y}_{2n}^{c},\Delta
W_{2n+1}\right\rangle\\\
&\quad+2\left\langle\widehat{Y}_{2n}^{c},Da(\widehat{Y}_{2n}^{c})\Delta
Z_{2n+1}\right\rangle+2\left\langle\widehat{Y}_{2n}^{c},Da(\widehat{Y}_{2n}^{c})\Delta
Z_{2n}\right\rangle+4\widetilde{\beta}h.\end{split}$ (28)

Substituting equations (25) and (26) into (27) and simplifying the inequality,
we obtain

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}\right\|^{2}&\leq
2Sh(1-Sh-2\widetilde{\alpha}h)\left\|\widehat{Y}_{2n}^{c}\right\|^{2}+(1-2Sh-4\widetilde{\alpha}h+Ch^{2})\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+C\Big{\\{}\left\|\Delta
W_{2n}\right\|^{2}\\\ &\quad+\left\|\Delta Z_{2n}\right\|^{2}+\left\|\Delta
W_{2n+1}\right\|^{2}+\left\|\Delta Z_{2n+1}\right\|^{2}\Big{\\}}\\\
&\quad+2\left\langle
Sh\widehat{Y}_{2n}^{c}+(1-Sh-2\widetilde{\alpha}h)\widehat{Y}_{2n}^{f},\Delta
W_{2n}\right\rangle\\\ &\quad+2\left\langle
Sh(Da(\widehat{Y}_{2n}^{c}))^{\top}\widehat{Y}_{2n}^{c}+(1-Sh-2\widetilde{\alpha}h)(Da(\widehat{Y}_{2n}^{f}))^{\top}\widehat{Y}_{2n}^{f},\Delta
Z_{2n}\right\rangle\\\ &\quad+2\left\langle\widehat{Y}_{2n+1}^{f},\Delta
W_{2n+1}\right\rangle+2\left\langle(Da(\widehat{Y}_{2n+1}^{f}))^{\top}\widehat{Y}_{2n+1}^{f},\Delta
Z_{2n+1}\right\rangle+4\widetilde{\beta}h-4\widetilde{\alpha}\widetilde{\beta}h^{2}.\end{split}$
(29)

Combining equations (29) and (28), note that for any
$\gamma\in(0,\widetilde{\alpha})$, and for all sufficiently small $h>0$, the
following inequality holds

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n+2}^{c}\right\|^{2}&\leq(1-4\gamma
h)\Bigg{(}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{f}\right\|^{2}\Bigg{)}+C\Big{\\{}\left\|\Delta
Z_{2n+1}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}\\\ &\quad+\left\|\Delta
W_{2n+1}\right\|^{2}+\left\|\Delta W_{2n}\right\|^{2}\Big{\\}}\\\
&\quad+2\left\langle(1+Sh)\widehat{Y}_{2n}^{c}+(1-Sh-2\widetilde{\alpha}h)\widehat{Y}_{2n}^{f},\Delta
W_{2n}\right\rangle\\\
&\quad+2\left\langle\widehat{Y}_{2n}^{c}+\widehat{Y}_{2n+1}^{f},\Delta
W_{2n+1}\right\rangle\\\
&\quad+2\left\langle(Da(\widehat{Y}_{2n}^{c}))^{\top}\widehat{Y}_{2n}^{c}+(Da(\widehat{Y}_{2n+1}^{f}))^{\top}\widehat{Y}_{2n+1}^{f},\Delta
Z_{2n+1}\right\rangle\\\
&\quad+2\left\langle(1+Sh)(Da(\widehat{Y}_{2n}^{c}))^{\top}\widehat{Y}_{2n}^{c},\Delta
Z_{2n}\right\rangle\\\
&\quad+2\left\langle(1-Sh-2\widetilde{\alpha}h)(Da(\widehat{Y}_{2n}^{f}))^{\top}\widehat{Y}_{2n}^{f},\Delta
Z_{2n}\right\rangle+8\widetilde{\beta}h.\end{split}$ (30)

Note that for all sufficiently small $h>0,\;(1-4\gamma
h)\leq\operatorname{e}^{-4\gamma h}$, and $\operatorname{e}^{4\gamma h}\leq
2$. Multiplying both sides of equation (30) by $\operatorname{e}^{2\gamma
t_{2n+2}}$, we obtain

$\begin{split}\operatorname{e}^{2\gamma
t_{2n+2}}\Bigg{\\{}\left\|\widehat{Y}_{2n+2}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n+2}^{c}\right\|^{2}\Bigg{\\}}&\leq\operatorname{e}^{2\gamma
t_{2n}}\Bigg{\\{}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\Bigg{\\}}+C\operatorname{e}^{2\gamma
t_{2n}}\Bigg{\\{}\left\|\Delta Z_{2n+1}\right\|^{2}\\\ &\quad+\left\|\Delta
Z_{2n}\right\|^{2}+\left\|\Delta W_{2n+1}\right\|^{2}+\left\|\Delta
W_{2n}\right\|^{2}\Bigg{\\}}\\\ &\quad+2\operatorname{e}^{2\gamma
t_{2n}}\left\langle\varphi_{1,2n},\Delta
W_{2n}\right\rangle+2\operatorname{e}^{2\gamma
t_{2n}}\left\langle\varphi_{1,2n+1},\Delta W_{2n+1}\right\rangle\\\
&\quad+2\operatorname{e}^{2\gamma t_{2n}}\left\langle\varphi_{2,2n},\Delta
Z_{2n}\right\rangle+2\operatorname{e}^{2\gamma
t_{2n}}\left\langle\varPhi_{2,2n+1},\Delta Z_{2n+1}\right\rangle\\\
&\quad+16\widetilde{\beta}h\operatorname{e}^{2\gamma t_{2n}},\end{split}$

where

$\begin{split}\varphi_{1,2n}&\coloneqq\operatorname{e}^{4\gamma
h}\Big{\\{}(1+Sh)\widehat{Y}_{2n}^{c}+(1-Sh-2\widetilde{\alpha}h)\widehat{Y}_{2n}^{f}\Big{\\}},\\\
\varphi_{1,2n+1}&\coloneqq\operatorname{e}^{4\gamma
h}\Big{\\{}\widehat{Y}_{2n}^{c}+\widehat{Y}_{2n+1}^{f}\Big{\\}},\\\
\varphi_{2,2n}&\coloneqq\operatorname{e}^{4\gamma
h}\Big{\\{}(1+Sh)(Da(\widehat{Y}_{2n}^{c}))^{\top}\widehat{Y}_{2n}^{c}+(1-Sh-2\widetilde{\alpha}h)(Da(\widehat{Y}_{2n}^{f}))^{\top}\widehat{Y}_{2n}^{f}\Big{\\}},\\\
\varphi_{2,2n+1}&\coloneqq\operatorname{e}^{4\gamma
h}\Big{\\{}(Da(\widehat{Y}_{2n}^{c}))^{\top}\widehat{Y}_{2n}^{c}+(Da(\widehat{Y}_{2n+1}^{f}))^{\top}\widehat{Y}_{2n+1}^{f}\Big{\\}}.\end{split}$

Summing the above inequality over multiple time-steps, we obtain

$\begin{split}\operatorname{e}^{2\gamma
t_{2n}}\Bigg{\\{}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\Bigg{\\}}&\leq
2\left\|x_{0}\right\|^{2}+C\Bigg{\\{}\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta
Z_{k}\right\|^{2}+\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta W_{k}\right\|^{2}\Bigg{\\}}\\\
&\quad+2\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\langle\varphi_{1,k},\Delta
W_{k}\right\rangle+2\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\langle\varphi_{2,k},\Delta Z_{k}\right\rangle\\\
&\quad+\sum_{k=0}^{n-1}16\widetilde{\beta}h\operatorname{e}^{2\gamma
t_{2k}}.\end{split}$ (31)

Note that for all sufficiently small
$h>0,\;\frac{1}{1-2\widetilde{\alpha}h}\leq 2$. Combining equations (25), (26)
and using Young’s inequality, we have

$\displaystyle\left\|\widehat{Y}_{2n+1}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n+1}^{c}\right\|^{2}$
$\displaystyle\leq(1-2\widetilde{\alpha}h)C\Bigg{\\{}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}+\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta
Z_{2n}\right\|^{2}+8\widetilde{\beta}h\Bigg{\\}}$

Multiplying both sides of the above inequality by the
$\operatorname{e}^{2\gamma t_{2n+1}}$ and substituting equation (31) into the
above inequality, we obtain

$\begin{split}\operatorname{e}^{2\gamma
t_{2n+1}}\Bigg{\\{}\left\|\widehat{Y}_{2n+1}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n+1}^{c}\right\|^{2}\Bigg{\\}}&\leq
C\Bigg{\\{}2\left\|x_{0}\right\|^{2}+\sum_{k=0}^{2n}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta Z_{k}\right\|^{2}+\sum_{k=0}^{2n}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta W_{k}\right\|^{2}\\\
&\quad+\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\langle\varphi_{1,k},\Delta
W_{k}\right\rangle+\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\langle\varPhi_{2,k},\Delta Z_{k}\right\rangle\\\
&\quad+\sum_{k=0}^{n-1}16\widetilde{\beta}h\operatorname{e}^{2\gamma
t_{2k}}\Bigg{\\}}.\end{split}$ (32)

Combining equations (31) and (32), raising both sides of the inequality to the
power $p/2$, taking the supremum over $n\in[0,N]$ and taking the expectation
on both sides of the inequality, we obtain

$\displaystyle\mathbb{E}\Bigg{[}\sup_{0\leq n\leq N}\operatorname{e}^{\gamma
pt_{n}}\Bigg{\\{}\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{2}+\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{2}\Bigg{\\}}^{p/2}\Bigg{]}$
$\displaystyle\leq\widetilde{C}\Big{(}J_{1}+J_{2}+J_{3}+J_{4}+J_{5}+J_{6}\Big{)},$
(33)

where $\widetilde{C}\coloneqq\max(6^{p/2-1}C^{p/2},1)$, and

$\displaystyle J_{1}$ $\displaystyle\coloneqq
2^{p/2}\left\|x_{0}\right\|^{p},$ $\displaystyle J_{4}$
$\displaystyle\coloneqq\mathbb{E}\Bigg{[}\sup_{1\leq n\leq
N/2}\Bigg{|}\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\langle\varphi_{1,k},\Delta
W_{k}\right\rangle\Bigg{|}^{p/2}\Bigg{]},$ $\displaystyle~{}J_{2}$
$\displaystyle\coloneqq\mathbb{E}\Bigg{[}\Bigg{|}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta Z_{k}\right\|^{2}\Bigg{|}^{p/2}\Bigg{]},$ $\displaystyle
J_{5}$ $\displaystyle\coloneqq\mathbb{E}\Bigg{[}\sup_{1\leq n\leq
N/2}\Bigg{|}\sum_{k=0}^{2n-1}\operatorname{e}^{2\gamma
t_{k}}\left\langle\varphi_{2,k},\Delta
Z_{k}\right\rangle\Bigg{|}^{p/2}\Bigg{]},$ $\displaystyle~{}J_{3}$
$\displaystyle\coloneqq\mathbb{E}\Bigg{[}\Bigg{|}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta W_{k}\right\|^{2}\Bigg{|}^{p/2}\Bigg{]},$ $\displaystyle
J_{6}$
$\displaystyle\coloneqq\left|\sum_{k=0}^{N-1}16\widetilde{\beta}\operatorname{e}^{2\gamma
t_{k}}h\right|^{p/2}.$

Recall that

$\Delta W_{k}=\Delta V_{1},\qquad\Delta Z_{k}=\frac{h}{2}\bigg{(}\Delta
V_{1}+\frac{\Delta V_{2}}{\sqrt{3}}\bigg{)}.$

Note that the term $J_{1}$ is a constant and is bounded. To bound the terms
$J_{2},J_{3},J_{4}$ and $J_{6}$, we make use of the fact that

$\displaystyle\mathbb{E}[\left\|\Delta W_{k}\right\|^{p}/h^{p/2}]$
$\displaystyle=\mathbb{E}[\left\|\Delta V_{1,k}\right\|^{p}/h^{p/2}]\leq
d^{p/2}p!!\leq d^{p/2}p^{p/2},$ $\displaystyle\mathbb{E}[\left\|\Delta
Z_{k}\right\|^{p}/h^{3p/2}]$
$\displaystyle=\frac{1}{2^{p}h^{p/2}}\mathbb{E}\Bigg{[}\left\|\Delta
V_{1,k}+\frac{1}{\sqrt{3}}\Delta V_{2,k}\right\|^{p}\Bigg{]}$
$\displaystyle\leq\frac{1}{2h^{p/2}}\Bigg{(}\mathbb{E}\Big{[}\left\|\Delta
V_{1,k}\right\|^{p}\Big{]}+\frac{1}{3^{p/2}}\mathbb{E}\Big{[}\left\|\Delta
V_{2,k}\right\|^{p}\Big{]}\Bigg{)}\leq 2d^{p/2}p^{p/2}.$

The terms $J_{3}$ and $J_{6}$ can be bounded in a similar way as in [9]
yielding the bounds

$\begin{split}J_{3}&=\mathbb{E}\Bigg{[}\Bigg{|}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta W_{k}\right\|^{2}\Bigg{|}^{p/2}\Bigg{]}\leq
d^{p/2}p^{p/2}\frac{\operatorname{e}^{\gamma pT}}{(2\gamma)^{p/2}},\\\
J_{6}&=\left|\sum_{k=0}^{N-1}16\widetilde{\beta}\operatorname{e}^{2\gamma
t_{k}}h\right|^{p/2}\leq
4^{p}\widetilde{\beta}^{p/2}\frac{\operatorname{e}^{\gamma
pT}}{(2\gamma)^{p/2}}.\end{split}$

Similarly, the term $J_{2}$ yields the bounds

$\displaystyle
J_{2}=\mathbb{E}\Bigg{[}\Bigg{|}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta Z_{k}\right\|^{2}\Bigg{|}^{p/2}\Bigg{]}$
$\displaystyle=\mathbb{E}\Bigg{[}\Bigg{|}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h\frac{\left\|\Delta Z_{k}\right\|^{2}}{h}\Bigg{|}^{p/2}\Bigg{]}$
$\displaystyle\leq\left|\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h\right|^{p/2-1}\mathbb{E}\Bigg{[}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h^{1+p}\frac{\left\|\Delta Z_{k}\right\|^{p}}{h^{3p/2}}\Bigg{]}$
$\displaystyle\leq 2d^{p/2}p^{p/2}h^{p}\frac{\operatorname{e}^{\gamma
pT}}{(2\gamma)^{p/2}}.f$

To bound the terms $J_{4}$ and $J_{5}$, we rewrite the summations in the terms
$J_{4}$ and $J_{5}$ as Itô integrals, and remind ourselves that

$dW_{s}=dV_{1},\qquad
dZ_{s}=\frac{h}{2}\bigg{(}dV_{1}+\frac{dV_{2}}{\sqrt{3}}\bigg{)}.$

where $dV_{1}$ and $dV_{2}$ are independent standard Wiener processes.
Applying the Burkholder–Davis–Gundy (BDG) inequality as stated in [4], we then
have

$\begin{split}J_{4}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq t\leq
T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor s/h\right\rfloor
h}\left\langle\varphi_{1,\left\lfloor
s/h\right\rfloor},\operatorname{d\\!}W_{s}\right\rangle\right|^{p/2}\Bigg{]}\leq(C_{\text{BDG}}p)^{p/4}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{4\gamma
t_{k}}h\left\|\varphi_{1,k}\right\|^{2}\right|^{p/4}\Bigg{]},\\\ \end{split}$

where $C_{\text{BDG}}>0$ is a constant that is independent of $p$. Using
Young’s inequality, for all sufficiently small $h>0$, we have

$\left\|\varphi_{1,2n}\right\|^{2}\leq
8\Big{(}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\Big{)},\qquad\left\|\varphi_{1,2n+1}\right\|^{2}\leq
4\Big{(}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n+1}^{c}\right\|^{2}\Big{)}.$

Using the above inequalities, we can write

$\begin{split}J_{4}&\leq(12C_{\text{BDG}}p)^{p/4}\left|\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h\right|^{p/4-1}\mathbb{E}\Bigg{[}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h\operatorname{e}^{\gamma
pt_{k}/2}\Big{(}\left\|\widehat{Y}_{t_{k}}^{c}\right\|^{2}+\left\|\widehat{Y}_{t_{k}}^{f}\right\|^{2}\Big{)}^{p/4}\Bigg{]}\\\
~{}&\leq(12C_{\text{BDG}}p)^{p/4}\Bigl{(}\frac{\operatorname{e}^{2\gamma
T}}{2\gamma}\Bigr{)}^{p/4-1}\mathbb{E}\Bigg{[}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h\cdot\sup_{0\leq n\leq N}\operatorname{e}^{\gamma
pt_{n}/2}\Big{(}\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{2}+\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{2}\Big{)}^{p/4}\Bigg{]}\\\
~{}&\leq(12C_{\text{BDG}}p)^{p/4}\mathbb{E}\Bigg{[}\frac{\operatorname{e}^{\gamma
pT/2}}{(2\gamma)^{p/4}}\cdot\sup_{0\leq n\leq N}\operatorname{e}^{\gamma
pt_{n}/2}\Big{(}\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{2}+\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{2}\Big{)}^{p/4}\Bigg{]}.\end{split}$

Using Young’s inequality, and Jensen’s inequality, for any $\zeta>0$, we can
write

$\begin{split}J_{4}\leq\frac{1}{4\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N}\operatorname{e}^{\gamma
pt_{n}}\Big{(}\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{2}+\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{2}\Big{)}^{p/2}\Bigg{]}+\zeta\Big{(}\frac{6C_{\text{BDG}}}{\gamma}\Big{)}^{p/2}p^{p/2}\operatorname{e}^{\gamma
pT}.\end{split}$

From Young’s inequality we have

$\left\|\varphi_{2,2n}\right\|^{2}\leq
8dK_{b}\Big{(}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\Big{)},\qquad\left\|\varphi_{2,2n+1}\right\|^{2}\leq
4dK_{b}\Big{(}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\Big{)}.$

Using the above inequality, and Burkholder–Davis–Gundy’s inequality, we can
bound $J_{5}$ as

$\begin{split}J_{5}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq t\leq
T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor s/h\right\rfloor
h}\left\langle\varphi_{2,\left\lfloor
s/h\right\rfloor},\operatorname{d\\!}Z_{s}\right\rangle\right|^{p/2}\Bigg{]}\\\
~{}&\leq(4dK_{b}C_{\text{BDG}}p)^{p/4}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{4\gamma
t_{k}}h^{3}\left\|\varphi_{2,k}\right\|^{2}\right|^{p/4}\Bigg{]}\\\
~{}&\leq(4dK_{b}C_{\text{BDG}}p)^{p/4}\left|\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h\right|^{p/4-1}\mathbb{E}\Bigg{[}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h^{1+\frac{p}{2}}\operatorname{e}^{\frac{\gamma
pt_{k}}{2}}\Big{(}\left\|\widehat{Y}_{t_{k}}^{c}\right\|^{2}+\left\|\widehat{Y}_{t_{k}}^{f}\right\|^{2}\Big{)}^{\frac{p}{4}}\Bigg{]}\\\
~{}&\leq(4dK_{b}C_{\text{BDG}}p)^{p/4}\Bigl{(}\frac{\operatorname{e}^{2\gamma
T}}{2\gamma}\Bigr{)}^{p/4-1}\mathbb{E}\Bigg{[}\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}h^{1+p/2}\cdot\sup_{0\leq n\leq N}\operatorname{e}^{\gamma
pt_{n}/2}\Big{(}\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{2}+\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{2}\Big{)}^{p/4}\Bigg{]}.\end{split}$

Using Young’s inequality, and Jensen’s inequality, for any $\zeta>0$, we have

$\begin{split}J_{5}\leq\frac{1}{4\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N}\operatorname{e}^{\gamma
pt_{n}}\Big{(}\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{2}+\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{2}\Big{)}^{p/2}\Bigg{]}+\zeta\bigg{(}\frac{2dK_{b}C_{\text{BDG}}}{\gamma}\bigg{)}^{p/2}p^{p/2}h^{p}\operatorname{e}^{\gamma
pT}.\end{split}$

Combining all the terms together, choosing $\zeta=\widetilde{C}$, using
Jensen’s inequality, there exists a constant $\kappa_{2}>0$ independent of $T$
and $p$, for all $T>0$ and $p\geq 1$, such that

$\begin{split}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq N}\operatorname{e}^{\gamma
pt_{n}}\Bigg{(}\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{2}+\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{2}\Bigg{)}^{p/2}\Bigg{]}&\leq\kappa_{2}^{p}p^{p/2}\operatorname{e}^{\gamma
pT}.\end{split}$

Since the above inequality holds for all $T>0$, there exists a constant
$\kappa_{1}>0$ such that for any $0<h<\kappa_{1}$,

$\sup_{0\leq n\leq
N}\mathbb{E}\Big{[}\left\|\widehat{Y}_{t_{n}}^{f}\right\|^{p}\Big{]}\leq\kappa_{2}^{p}p^{p/2},\qquad\sup_{0\leq
n\leq
N}\mathbb{E}\Big{[}\left\|\widehat{Y}_{t_{n}}^{c}\right\|^{p}\Big{]}\leq\kappa_{2}^{p}p^{p/2}.$
(34)

#### Proof of Theorem 2

We obtain the bound only for the fine trajectory to keep the argument brief.
The proof argument can be extended similarly for the coarse trajectory. We
also introduce the constants $C,\widetilde{\alpha},\widetilde{\beta}>0$ that
depend only on $d$, and the constants from Assumptions 1, 2, and 3. The
constants $C,\widetilde{\alpha},\widetilde{\beta}$ may change from line to
line in the proof argument.

From equations (4), and (5), computing the square of the Euclidean norm of the
fine trajectory, we obtain

$\begin{split}\left\|\overline{X}_{n+1}^{f}\right\|^{2}&=\left\|X_{n}\right\|^{2}+\left\|\widehat{I}_{n}\right\|^{2}+2\left\langle\overline{X}_{n}^{f},\widehat{I}_{n}\right\rangle,\end{split}$
(35)

where

$\widehat{I}_{n}\coloneqq ha(\overline{X}_{n}^{f})+\Delta
W_{n}+Da(\overline{X}_{n}^{f})\Delta
Z_{n}+\frac{h^{2}}{2}\mathcal{A}a(\overline{X}_{n}^{f}).$

Using Hölder’s inequality and Assumptions 1, and 2, the term
$\left\|\widehat{I}_{n}\right\|^{2}$ yields the bound

$\begin{split}\left\|\widehat{I}_{n}\right\|^{2}&\leq
4\bigg{(}h^{2}\left\|a(\overline{X}_{n}^{f})\right\|^{2}+\left\|\Delta
W_{n}\right\|^{2}+\left\|Da(\overline{X}_{n}^{f})\Delta
Z_{n}\right\|^{2}+\frac{h^{4}}{4}\left\|\mathcal{A}a(\overline{X}_{n}^{f})\right\|^{2}\bigg{)}\\\
~{}&\leq
C\bigg{(}h^{2}\big{(}1+\left\|\overline{X}_{n}^{f}\right\|^{2}\big{)}+\left\|\Delta
W_{n}\right\|^{2}+\left\|\Delta Z_{n}\right\|^{2}\bigg{)}.\end{split}$ (36)

Using Assumption 3, and Young’s inequality, the term
$\left\langle\overline{X}_{n}^{f},\widehat{I}_{n}\right\rangle$ yields

$\begin{split}\left\langle\overline{X}_{n}^{f},\widehat{I}_{n}\right\rangle&=h\left\langle\overline{X}_{n}^{f},a(\overline{X}_{n}^{f})\right\rangle+\left\langle\overline{X}_{n}^{f},\Delta
W_{n}\right\rangle+\left\langle\overline{X}_{n}^{f},Da(\overline{X}_{n}^{f})\Delta
Z_{n}\right\rangle+\frac{h^{2}}{2}\left\langle\overline{X}_{n}^{f},\mathcal{A}a(\overline{X}_{n}^{f})\right\rangle\\\
~{}&\leq-\alpha h\left\|\overline{X}_{n}^{f}\right\|^{2}+\beta
h+\left\langle\overline{X}_{n}^{f},\Delta
W_{n}\right\rangle+\left\langle\overline{X}_{n}^{f},Da(\overline{X}_{n}^{f})\Delta
Z_{n}\right\rangle+Ch^{2}\Big{(}1+\left\|\overline{X}_{n}^{f}\right\|^{2}\Big{)}\\\
~{}&\leq-\widetilde{\alpha}h\left\|\overline{X}_{n}^{f}\right\|^{2}+\left\langle\overline{X}_{n}^{f},\Delta
W_{n}\right\rangle+\left\langle\overline{X}_{n}^{f},Da(\overline{X}_{n}^{f})\Delta
Z_{n}\right\rangle+\widetilde{\beta}h.\end{split}$ (37)

Substituting equations (36) and (37) in (35), we obtain

$\begin{split}\left\|\overline{X}_{n+1}^{f}\right\|^{2}&\leq(1-2\widetilde{\alpha}h+Ch^{2})\left\|\overline{X}_{n}^{f}\right\|^{2}+C(\left\|\Delta
W_{n}\right\|^{2}+\left\|\Delta
Z_{n}\right\|^{2})+2\left\langle\overline{X}_{n}^{f},\Delta
W_{n}\right\rangle\\\
&\quad+2\left\langle\overline{X}_{n}^{f},Da(\overline{X}_{n}^{f})\Delta
Z_{n}\right\rangle+2\widetilde{\beta}h.\end{split}$

Note that for any $\gamma\in(0,\widetilde{\alpha}-Ch/2)$, and for all
sufficiently small $h>0$, the following inequality holds

$\begin{split}\left\|\overline{X}_{n+1}^{f}\right\|^{2}&\leq(1-2\gamma
h)\left\|\overline{X}_{n}^{f}\right\|^{2}+C(\left\|\Delta
W_{n}\right\|^{2}+\left\|\Delta
Z_{n}\right\|^{2})+2\left\langle\overline{X}_{n}^{f},\Delta
W_{n}\right\rangle\\\
&\quad+2\left\langle\overline{X}_{n}^{f},Da(\overline{X}_{n}^{f})\Delta
Z_{n}\right\rangle+2\widetilde{\beta}h.\end{split}$

Multiplying both sides of the above inequality by $\operatorname{e}^{2\gamma
t_{n+1}}$, noting that $\operatorname{e}^{2\gamma h}\leq 2$ and $(1-2\gamma
h)\leq\operatorname{e}^{-2\gamma h}$, for all sufficiently small $h>0$, we
obtain

$\begin{split}\operatorname{e}^{2\gamma
t_{n+1}}\left\|\overline{X}_{n+1}^{f}\right\|^{2}&\leq\operatorname{e}^{2\gamma
t_{n}}\left\|\overline{X}_{n}^{f}\right\|^{2}+C\operatorname{e}^{2\gamma
t_{n}}(\left\|\Delta W_{n}\right\|^{2}+\left\|\Delta
Z_{n}\right\|^{2})+4\operatorname{e}^{2\gamma
t_{n}}\left\langle\overline{X}_{n}^{f},\Delta W_{n}\right\rangle\\\
&\quad+4\operatorname{e}^{2\gamma
t_{n}}\left\langle\overline{X}_{n}^{f},Da(\overline{X}_{n}^{f})\Delta
Z_{n}\right\rangle+4\widetilde{\beta}\operatorname{e}^{2\gamma
t_{n}}h.\end{split}$

Recursively summing the above inequality over $n=0,1,\ldots,N-1$, taking the
power $p/2$ on both sides of the inequality, and taking the supremum over all
$n$, and taking the expected value on both sides of the inequality, we obtain

$\begin{split}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq N-1}\operatorname{e}^{\gamma
pt_{n+1}}\left\|\overline{X}_{n+1}^{f}\right\|^{p}\Bigg{]}&\leq\widetilde{C}\Big{(}J_{1}+J_{2}+J_{3}+J_{4}+J_{5}+J_{6}\Big{)},\end{split}$

where $\widetilde{C}\coloneqq\max(6^{p/2-1}C,\;1)$, and the terms $J_{1}$
through $J_{6}$ are defined and bounded below.

The terms $J_{1}$ and $J_{6}$ yield the bounds

$\begin{split}J_{1}\coloneqq\left\|x_{0}\right\|^{p}<\infty,\qquad
J_{6}\coloneqq\left|\sum_{k=0}^{N-1}4\widetilde{\beta}\operatorname{e}^{2\gamma
t_{k}}h\right|^{p/2}\leq
2^{p}\widetilde{\beta}^{p/2}\frac{\operatorname{e}^{\gamma
pT}}{(2\gamma)^{p/2}}.\end{split}$

The terms $J_{2},J_{3},J_{4}$, and $J_{5}$ can be bounded in a similar way as
in Theorem 1 to yield the bounds

$\begin{split}J_{2}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta W_{k}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
d^{p/2}p^{p/2}\frac{\operatorname{e}^{\gamma pT}}{(2\gamma)^{p/2}},\\\
~{}J_{3}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{2\gamma
t_{k}}\left\|\Delta Z_{k}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
2d^{p/2}p^{p/2}h^{p}\frac{\operatorname{e}^{\gamma pT}}{(2\gamma)^{p/2}},\\\
~{}J_{4}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq t\leq
T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor s/h\right\rfloor
h}\left\langle\overline{X}_{\left\lfloor
s/h\right\rfloor}^{f},\operatorname{d\\!}W_{s}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{4\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N}\operatorname{e}^{\gamma
pt_{n}}\left\|\overline{X}_{n}^{f}\right\|^{p}\Bigg{]}+\zeta\bigg{(}\frac{C_{\text{BDG}}}{\gamma}\bigg{)}^{p/2}p^{p/2}\operatorname{e}^{\gamma
pT},\\\ ~{}J_{5}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq t\leq
T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor s/h\right\rfloor
h}\left\langle\overline{X}_{\left\lfloor
s/h\right\rfloor}^{f},Da(\overline{X}_{\left\lfloor
s/h\right\rfloor}^{f})\operatorname{d\\!}Z_{s}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{4\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N}\operatorname{e}^{\gamma
pt_{n}}\left\|\overline{X}_{n}^{f}\right\|^{p}\Bigg{]}+\zeta\bigg{(}\frac{dK_{b}C_{\text{BDG}}}{\gamma}\bigg{)}^{p/2}p^{p/2}h^{p}\operatorname{e}^{\gamma
pT},\end{split}$

where $C_{\text{BDG}}>0$ is a constant independent of $p$.

Combining all the terms together, choosing $\zeta=\widetilde{C}$, using
Jensen’s inequality, there exists a constant $C_{2}>0$ independent of $T$ and
$p$, for all $T>0$ and $p\geq 1$, such that

$\begin{split}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq N}\operatorname{e}^{\gamma
pt_{n}}\left\|\overline{X}_{t_{n}}^{f}\right\|^{p}\Bigg{]}&\leq
C_{2}^{p}p^{p/2}\operatorname{e}^{\gamma pT}.\end{split}$

Since the above inequality holds for all $T>0$, there exists a constant
$C_{1}>0$ such that for any $0<h<C_{1}$,

$\sup_{0\leq n\leq
N}\mathbb{E}\Big{[}\left\|\overline{X}_{t_{n}}^{f}\right\|^{p}\Big{]}\leq
C_{2}^{p}p^{p/2}.$

#### Proof of Theorem 3

Using equations (7) to (10), we can write the difference between the fine and
the coarse trajectories at the even time-step as

$\begin{split}\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}&=(1-4Sh+2S^{2}h^{2})(\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c})+(2-Sh)h\big{(}a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\big{)}\\\
&\quad+h\big{(}a(\widehat{Y}_{2n+1}^{f})-a(\widehat{Y}_{2n}^{f})\big{)}+(1-Sh)\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n}\\\
&\quad+\big{(}Da(\widehat{Y}_{2n+1}^{f})-Da(\widehat{Y}_{2n}^{f})\big{)}\Delta
Z_{2n+1}+\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n+1}\\\ &\quad-hDa(\widehat{Y}_{2n}^{c})\Delta
W_{2n}+\frac{h^{2}}{2}(2-Sh)\big{(}\mathcal{A}a(\widehat{Y}_{2n}^{f})-\mathcal{A}a(\widehat{Y}_{2n}^{c})\big{)}\\\
&\quad+\frac{h^{2}}{2}\big{(}\mathcal{A}a(\widehat{Y}_{2n+1}^{f})-\mathcal{A}a(\widehat{Y}_{2n}^{f})\big{)}-h^{2}\mathcal{A}a(\widehat{Y}_{2n}^{c}).\end{split}$
(38)

As a first step, we aim to obtain a loose upper-bound on the difference
between the fine and the coarse trajectories

$\begin{split}\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}&=(1-4Sh+2S^{2}h^{2})(\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c})+(2-Sh)h\big{(}a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\big{)}+\widehat{I}_{2n+1},\end{split}$
(39)

where

$\begin{split}\widehat{I}_{2n+1}&\coloneqq
h\big{(}a(\widehat{Y}_{2n+1}^{f})-a(\widehat{Y}_{2n}^{f})\big{)}+(1-Sh)\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n}\\\
&\quad+\big{(}Da(\widehat{Y}_{2n+1}^{f})-Da(\widehat{Y}_{2n}^{f})\big{)}\Delta
Z_{2n+1}+\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n+1}\\\ &\quad-hDa(\widehat{Y}_{2n}^{c})\Delta
W_{2n}+\frac{h^{2}}{2}(2-Sh)\big{(}\mathcal{A}a(\widehat{Y}_{2n}^{f})-\mathcal{A}a(\widehat{Y}_{2n}^{c})\big{)}\\\
&\quad+\frac{h^{2}}{2}\big{(}\mathcal{A}a(\widehat{Y}_{2n+1}^{f})-\mathcal{A}a(\widehat{Y}_{2n}^{f})\big{)}-h^{2}\mathcal{A}a(\widehat{Y}_{2n}^{c}).\end{split}$

Taking the square of the Euclidean norm on both sides of equation (39), using
the Cauchy–Schwarz inequality and the one-sided Lipschitz property of the
drift coefficient, we can for sufficiently small $h>0$ write

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}&=\left\|(1-4Sh+2S^{2}h^{2})(\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c})+(2-Sh)h\big{(}a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\big{)}+\widehat{I}_{2n+1}\right\|^{2}\\\
~{}&\leq(1-4Sh+2S^{2}h^{2})^{2}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+4h^{2}\left\|a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\right\|^{2}\\\
&\quad+\left\|\widehat{I}_{2n+1}\right\|^{2}+2h(1-4Sh+2S^{2}h^{2})(2-Sh)\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\right\rangle\\\
&\quad+2(1+2\sqrt{d}K_{b}h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\left\|\widehat{I}_{2n+1}\right\|\\\
~{}&\leq(1-4(2S-\lambda)h+Ch^{2})\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+\left\|\widehat{I}_{2n+1}\right\|^{2}\\\
&\quad+2(1+2\sqrt{d}K_{b}h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\left\|\widehat{I}_{2n+1}\right\|,\end{split}$

where $C>0$ is a constant that depends only on $S,d$, and the constants from
Assumptions 1, 2 and 3. The constant $C$ may change from line to line in the
proof argument. From the above inequality, we observe that for all
sufficiently small $h>0$, there exists a $\gamma\in(0,2S-\lambda)$ such that
the following holds

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}&\leq(1-4\gamma
h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+3\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\left\|\widehat{I}_{2n+1}\right\|+\left\|\widehat{I}_{2n+1}\right\|^{2}.\end{split}$

Further applying Young’s inequality to

$\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\left\|\widehat{I}_{2n+1}\right\|\leq\frac{2\gamma
h}{3}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+\frac{3}{8\gamma
h}\left\|\widehat{I}_{2n+1}\right\|^{2},$

we obtain

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}&\leq(1-2\gamma
h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+\frac{C}{h}\left\|\widehat{I}_{2n+1}\right\|^{2}.\end{split}$
(40)

Using Hölder’s inequality, and Assumptions 1, 3, and 2, the term
$\widehat{I}_{2n+1}$ can be bounded as follows

$\begin{split}\left\|\widehat{I}_{2n+1}\right\|^{2}&\leq
8C\Bigg{\\{}h^{2}\left\|\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\right\|^{2}+h^{2}\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}+\left\|\Delta
Z_{2n+1}\right\|^{2}\\\
&\quad+h^{4}\Big{(}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\Big{)}+h^{4}\Bigg{\\}}.\end{split}$
(41)

From equation (7), we have

$\begin{split}\left\|\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\right\|^{p}&\leq
6^{p-1}C\Bigg{\\{}h^{p}\bigg{(}1+\left\|\widehat{Y}_{2n}^{f}\right\|^{p}+\left\|\widehat{Y}_{2n}^{c}\right\|^{p}\bigg{)}+\left\|\Delta
W_{2n}\right\|^{p}+\left\|\Delta Z_{2n}\right\|^{p}\Bigg{\\}}.\end{split}$
(42)

Substituting equation (42) in (41), for $p=2$, we obtain

$\begin{split}\left\|\widehat{I}_{2n+1}\right\|^{2}&\leq
C\Bigg{\\{}h^{4}\bigg{(}1+\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\bigg{)}+h^{2}\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}+\left\|\Delta
Z_{2n+1}\right\|^{2}\Bigg{\\}}.\end{split}$ (43)

Substituting equation (43) in (40), we obtain

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}&\leq(1-2\gamma
h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+C\Bigg{\\{}h^{3}\bigg{(}\left\|\widehat{Y}_{2n}^{f}\right\|^{2}+\left\|\widehat{Y}_{2n}^{c}\right\|^{2}\bigg{)}+h\left\|\Delta
W_{2n}\right\|^{2}\\\ &\quad+\frac{1}{h}\Big{(}\left\|\Delta
Z_{2n}\right\|^{2}+\left\|\Delta
Z_{2n+1}\right\|^{2}\Big{)}+h^{3}\Bigg{\\}}.\end{split}$

Multiplying both sides of the above inequality by $\operatorname{e}^{\gamma
t_{2n}}$, noting that $1-2\gamma h\leq e^{2\gamma h}\leq 2$ holds for all
sufficiently small $h>0$, and from recursive summation of the inequality over
multiple time-steps, we obtain

$\begin{split}\operatorname{e}^{\gamma
t_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}&\leq
C\Bigg{\\{}h^{3}\sum_{k=0}^{2n-1}\operatorname{e}^{\gamma
t_{k}}\bigg{(}\left\|\widehat{Y}_{k}^{f}\right\|^{2}+\left\|\widehat{Y}_{k}^{c}\right\|^{2}\bigg{)}+h\sum_{k=0}^{2n-1}\operatorname{e}^{\gamma
t_{k}}\left\|\Delta W_{k}\right\|^{2}\\\
&\quad+\frac{1}{h}\sum_{k=0}^{2n-1}\operatorname{e}^{\gamma
t_{k}}\left\|\Delta
Z_{k}\right\|^{2}+\sum_{k=0}^{2n-1}\operatorname{e}^{\gamma
t_{k}}h^{3}\Bigg{\\}}.\end{split}$

Raising both sides of the above inequality to the power $p/2$, taking the
supremum over $n\in[0,N/2]$, and taking the expectation on both sides of the
inequality, we obtain

$\begin{split}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq N/2}\operatorname{e}^{\gamma
pt_{2n}/2}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}&\leq\widetilde{C}\Bigg{\\{}J_{1}+J_{2}+J_{3}+J_{4}\Bigg{\\}},\end{split}$

where $\widetilde{C}\coloneqq\max(4^{p/2-1}C^{p/2},1)$, and

$\begin{aligned}
J_{1}&=h^{3p/2}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}\Big{(}\left\|\widehat{Y}_{k}^{f}\right\|^{2}+\left\|\widehat{Y}_{k}^{c}\right\|^{2}\Big{)}\right|^{p/2}\Bigg{]},\\\
J_{3}&=h^{3p/2}\left|\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}\right|^{p/2},\end{aligned}\qquad\begin{aligned}
J_{2}&=h^{p/2}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}\left\|\Delta W_{k}\right\|^{2}\right|^{p/2}\Bigg{]},\\\
J_{4}&=h^{-p/2}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}\left\|\Delta Z_{k}\right\|^{2}\right|^{p/2}\Bigg{]}.\end{aligned}$

To keep the argument concise, we introduce the constants $C_{J_{i}}>0$ for
$i=1,2,3,4$ that depend only on $S,d$, $\gamma$, and the constants from
Assumptions 1, 2 and 3. Using Theorem 1, term $J_{1}$ yields

$\begin{split}J_{1}&=h^{p}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}h\Big{(}\left\|\widehat{Y}_{k}^{f}\right\|^{2}+\left\|\widehat{Y}_{k}^{c}\right\|^{2}\Big{)}\right|^{p/2}\Bigg{]}\\\
&\leq h^{p}\left|\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}h\right|^{p/2-1}\mathbb{E}\Bigg{[}\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}h\Big{(}\left\|\widehat{Y}_{k}^{f}\right\|^{2}+\left\|\widehat{Y}_{k}^{c}\right\|^{2}\Big{)}^{p/2}\Bigg{]}\\\
&\leq 2^{p/2-1}h^{p}\left|\int_{0}^{T}\operatorname{e}^{\gamma
t}\operatorname{d\\!}t\right|^{p/2-1}\sum_{k=0}^{N-1}\operatorname{e}^{\gamma
t_{k}}h\mathbb{E}\Big{[}\left\|\widehat{Y}_{k}^{f}\right\|^{p}+\left\|\widehat{Y}_{k}^{c}\right\|^{p}\Big{]}\\\
&\leq C_{J_{1}}^{p}p^{p/2}h^{p}\operatorname{e}^{\gamma pT/2}.\end{split}$

The terms $J_{2}$, $J_{3}$, and $J_{4}$ can be bounded in a similar way as in
Theorem 1 yielding the bounds

$J_{2}\leq C_{J_{2}}^{p}p^{p/2}h^{p/2}\operatorname{e}^{\gamma pT/2},\qquad
J_{3}\leq C_{J_{3}}^{p}h^{p}\operatorname{e}^{\gamma pT/2},\qquad J_{4}\leq
C_{J_{4}}^{p}p^{p/2}h^{p/2}\operatorname{e}^{\gamma pT/2}.$

Combining all the estimates together, there exists a constant
$\widetilde{C}_{1}>0$ that is independent of $T$ and $p$, such that

$\begin{split}\mathbb{E}\Big{[}\sup_{0\leq n\leq N/2}\operatorname{e}^{\gamma
pt_{2n}/2}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Big{]}&\leq\widetilde{C}_{(1)}^{p}p^{p/2}h^{p/2}\operatorname{e}^{\gamma
pT/2}.\end{split}$

Since the above inequality holds for any arbitrary $T>0$, there exists a
constant $\kappa_{3}>0$, such that for any $0<h<\kappa_{3}$,

$\sup_{0\leq n\leq
N/2}\mathbb{E}\Big{[}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Big{]}\leq\widetilde{C}_{1}^{p}p^{p/2}h^{p/2}.$
(44)

Having obtained a loose upper bound on the difference between the fine and the
coarse trajectories, we aim to sharpen the bound obtained in equation (44).
Taking the square of the Euclidean norm on both sides of equation (39) and
using the one-sided Lipschitz property of the drift coefficient, we can write

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}&=\left\|(1-4Sh+2S^{2}h^{2})(\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c})+(2-Sh)h\big{(}a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\big{)}+\widehat{I}_{2n+1}\right\|^{2}\\\
~{}&\leq(1-4Sh+2S^{2}h^{2})^{2}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+6h^{2}\left\|a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\right\|^{2}\\\
&\quad+3\left\|\widehat{I}_{2n+1}\right\|^{2}+2h(1-4Sh+2S^{2}h^{2})(2-Sh)\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},a(\widehat{Y}_{2n}^{f})-a(\widehat{Y}_{2n}^{c})\right\rangle\\\
&\quad+2(1-4Sh+2S^{2}h^{2})\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\widehat{I}_{2n+1}\right\rangle\\\
~{}&\leq(1-4(2S-\lambda)h+Ch^{2})\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+3\left\|\widehat{I}_{2n+1}\right\|^{2}\\\
&\quad+2(1-4Sh+2S^{2}h^{2})\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\widehat{I}_{2n+1}\right\rangle.\end{split}$

As in the first approach, there exists a $\gamma\in(0,2S-\lambda)$, for all
sufficiently small $h>0$, such that

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}&\leq(1-4\gamma
h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+3\left\|\widehat{I}_{2n+1}\right\|^{2}\\\
&\qquad+2(1-4Sh+2S^{2}h^{2})\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\widehat{I}_{2n+1}\right\rangle.\end{split}$
(45)

Using the Mean value theorem, we can write

$\begin{split}a(\widehat{Y}_{2n+1}^{f})-a(\widehat{Y}_{2n}^{f})&=Da(\widehat{Y}_{2n}^{f})\Big{(}\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\Big{)}+R_{2},\end{split}$
(46)

where $R_{2}$ denotes the second-order remainder from the Mean value theorem.
From Assumption 1, it follows that $\left\|R_{2}\right\|\leq
d^{3/2}K_{b}\left\|\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\right\|^{2}$.

Substituting equation (46) in (45), using Hölder’s inequality, Young’s
inequality, and equation (42), we have

$\begin{split}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}&\leq(1-4\gamma
h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+C\Bigg{\\{}\Big{(}h^{2}\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}\\\ &\quad+\left\|\Delta
Z_{2n+1}\right\|^{2}\Big{)}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+h^{2}\left\|\Delta
Z_{2n}\right\|^{2}+h^{4}\Big{(}1+\left\|\widehat{Y}_{2n}^{f}\right\|^{2}\Big{)}\\\
&\quad+\Big{(}h^{4}+\left\|\Delta
Z_{2n+1}\right\|^{2}\Big{)}\left\|\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\right\|^{2}+h^{2}\left\|R_{2}\right\|^{2}\Bigg{\\}}\\\
&\quad+2(1-4Sh+2S^{2}h^{2})\Bigg{\\{}h\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},R_{2}\right\rangle\\\
&\quad+h\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
W_{2n}\right\rangle\\\
&\quad+h\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},Da(\widehat{Y}_{2n}^{f})Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}\right\rangle\\\
&\quad+(1-Sh)\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n}\right\rangle\\\
&\quad+\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n+1}^{f})-Da(\widehat{Y}_{2n}^{f})\big{)}\Delta
Z_{2n+1}\right\rangle\\\
&\quad+\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n+1}\right\rangle\\\
&\quad+\frac{h^{2}}{2}\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\mathcal{A}a(\widehat{Y}_{2n+1}^{f})-\mathcal{A}a(\widehat{Y}_{2n}^{f})\right\rangle-\frac{h^{2}}{2}\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\Delta
a(\widehat{Y}_{2n}^{f})\right\rangle\Bigg{\\}}.\end{split}$

Multiplying both sides of the above inequality by $\operatorname{e}^{2\gamma
t_{2n}}$, noting that $\operatorname{e}^{4\gamma h}\leq 2$ for all
sufficiently small $h>0$, summing over multiple time-steps, raising both sides
of the inequality to the power $p/2$, taking supremum over $n\in[0,N/2]$,
taking the expectation on both sides of the inequality, we obtain

$\displaystyle\mathbb{E}\Big{[}\sup_{0\leq n\leq N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{p}\Big{]}$
$\displaystyle\leq\widetilde{C}\sum_{m=1}^{16}J_{m},$

where $\widetilde{C}\coloneqq\max(16^{p/2-1}C^{p/2},1)$ and the terms $J_{m}$
for $m=1,\ldots,16$ are defined and bounded below. For the term $J_{1}$, we
have

$\begin{split}J_{1}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{2}\left\|\Delta
W_{2k}\right\|^{2}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{2}\right|^{p/2}\Bigg{]}\\\
&\leq\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h\right|^{p/2-1}\mathbb{E}\Bigg{[}\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{1+p}\frac{\left\|\Delta
W_{2k}\right\|^{p}}{h^{p/2}}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{p}\Bigg{]}\\\
&\leq d^{p/2}p^{p/2}h^{p}\left|\int_{0}^{T}\operatorname{e}^{2\gamma
t}\operatorname{d\\!}t\right|^{p/2-1}\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h\mathbb{E}\Big{[}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{p}\Big{]}.\end{split}$

Using equation (44), we can write

$\begin{split}J_{1}&\leq C_{J_{1}}^{p}p^{p}h^{3p/2}\operatorname{e}^{\gamma
pT}.\end{split}$

Similarly, the terms $J_{2},J_{3}$, and $J_{4}$, yield the bounds

$\begin{split}J_{2}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}\left\|\Delta
Z_{2k}\right\|^{2}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
C_{J_{2}}^{p}p^{p}h^{3p/2}\operatorname{e}^{\gamma pT},\\\
J_{3}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}\left\|\Delta
Z_{2k+1}\right\|^{2}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
C_{J_{3}}^{p}p^{p}h^{3p/2}\operatorname{e}^{\gamma pT},\\\
J_{4}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{2}\left\|\Delta Z_{2k}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
C_{J_{4}}^{p}p^{p/2}h^{2p}\operatorname{e}^{\gamma pT}.\end{split}$

The term $J_{5}$, using the result from Theorem 1, yields the bound

$\begin{split}J_{5}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{4}\Big{(}1+\left\|\widehat{Y}_{2k}^{f}\right\|^{2}\Big{)}\right|^{p/2}\Bigg{]}\leq
C_{J_{5}}^{p}p^{p/2}h^{3p/2}\operatorname{e}^{\gamma pT}.\end{split}$

For the term $J_{6}$, using equation (47), we obtain

$\begin{split}J_{6}&\coloneqq
h^{3p/2}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
C_{J_{6}}^{p}p^{p/2}h^{2p}\operatorname{e}^{\gamma pT},\end{split}$

where, using the result from Theorem 1, we have

$\mathbb{E}\bigg{[}\left\|\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\right\|^{p}\bigg{]}\leq
C_{f}^{p}p^{p/2}h^{p/2},$ (47)

for all $p\geq 1$, where $C_{f}>0$ is a constant that independent of $T$ and
$p$. Similarly, using equation (46) we have

$\begin{split}J_{7}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}\left\|\Delta
Z_{2k+1}\right\|^{2}\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
C_{J_{7}}^{p}p^{p}h^{3p/2}\operatorname{e}^{\gamma pT},\\\
~{}J_{8}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{2}\left\|R_{2}\right\|^{2}\right|^{p/2}\Bigg{]}\leq
C_{J_{8}}^{p}p^{p}h^{3p/2}\operatorname{e}^{\gamma pT}.\end{split}$

For the term $J_{9}$, using Cauchy–Schwarz’s inequality, Jensen’s inequality,
and Young’s inequality for any $\zeta>0$, we have

$\begin{split}J_{9}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},R_{2}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq\mathbb{E}\Bigg{[}\left|\sup_{0\leq n\leq N/2}\operatorname{e}^{\gamma
t_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\right|^{p/2}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|^{2}\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+2\zeta\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|^{2}\right|^{p}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+\zeta
C_{J_{9}}^{p}p^{p}h^{p}\operatorname{e}^{\gamma pT}.\end{split}$

For the term $J_{10}$, using the Burkholder–Davis–Gundy inequality as stated
in [4], and equation (44), we have

$\begin{split}J_{10}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},\big{(}Da(\widehat{Y}_{2k}^{f})-Da(\widehat{Y}_{2k}^{c})\big{)}\Delta
W_{2k}\right\rangle\right|^{p/2}\Bigg{]}\\\ &=\mathbb{E}\Bigg{[}\sup_{0\leq
t\leq T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor
s/2h\right\rfloor h}h\left\langle\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{f}-\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{c},\big{(}Da(\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{f})-Da(\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{c})\big{)}\operatorname{d\\!}W_{s}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq(C_{\text{BDG}}p)^{p/4}d^{3p/4}K_{b}^{p/2}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{4\gamma
t_{2k}}h^{3}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{4}\right|^{p/4}\Bigg{]}\\\
&\leq(C_{\text{BDG}}p)^{p/4}d^{3p/4}K_{b}^{p/2}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{4\gamma
t_{2k}}h\right|^{p/4-1}\sum_{k=0}^{N/2-1}\operatorname{e}^{4\gamma
t_{2k}}h^{1+p/2}\mathbb{E}\Bigg{[}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{p}\Bigg{]}\\\
&\leq C_{J_{10}}^{p}p^{3p/4}h^{p}\operatorname{e}^{\gamma pT}.\end{split}$

The terms $J_{12}$ and $J_{14}$ can be bounded in a similar way as the term
$J_{10}$ to yield the bounds

$\begin{split}J_{12}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},\big{(}Da(\widehat{Y}_{2k}^{f})-Da(\widehat{Y}_{2k}^{c})\big{)}\Delta
Z_{2k}\right\rangle\right|^{p/2}\Bigg{]}\\\ &=\mathbb{E}\Bigg{[}\sup_{0\leq
t\leq T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor
s/2h\right\rfloor h}\left\langle\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{f}-\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{c},\big{(}Da(\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{f})-Da(\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{c})\big{)}\operatorname{d\\!}Z_{s}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq C_{J_{12}}^{p}p^{3p/4}h^{p}\operatorname{e}^{\gamma pT},\\\
~{}J_{14}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},\big{(}Da(\widehat{Y}_{2k}^{f})-Da(\widehat{Y}_{2k}^{c})\big{)}\Delta
Z_{2k+1}\right\rangle\right|^{p/2}\Bigg{]}\\\ &\leq
C_{J_{14}}^{p}p^{3p/4}h^{p}\operatorname{e}^{\gamma pT}.\end{split}$

The term $J_{11}$, using the Burkholder–Davis–Gundy inequality, Jensen’s
inquality, and Young’s inquality for any $\zeta>0$, can be bounded as

$\begin{split}J_{11}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},Da(\widehat{Y}_{2k}^{f})Da(\widehat{Y}_{2k}^{f})\Delta
Z_{2k}\right\rangle\right|^{p/2}\Bigg{]}\\\ &=\mathbb{E}\Bigg{[}\sup_{0\leq
t\leq T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor
s/2h\right\rfloor h}h\left\langle\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{f}-\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{c},Da(\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{f})Da(\widehat{Y}_{(2\left\lfloor
s/2h\right\rfloor)}^{f})\operatorname{d\\!}Z_{s}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+2\zeta(C_{\text{BDG}}p)^{p/2}d^{2p}K_{b}^{2p}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{5}\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+\zeta
C_{J_{11}}^{p}p^{p/2}h^{2p}\operatorname{e}^{\gamma pT}.\end{split}$

The term $J_{13}$ yields the bound

$\begin{split}J_{13}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},(Da(\widehat{Y}_{2k+1}^{f})-Da(\widehat{Y}_{2k}^{f}))\Delta
Z_{2k+1}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq(C_{\text{BDG}}p)^{p/4}d^{3p/4}K_{b}^{p/2}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{4\gamma
t_{2k}}h^{5}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{2}\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|^{2}\right|^{p/4}\Bigg{]}\\\
&\leq C_{J_{13}}^{p}p^{3p/4}h^{3p/2}\operatorname{e}^{\gamma pT}.\end{split}$

The term $J_{15}$, using the Cauchy–Schwarz inequality, Jensen’s inequality,
and Young’s inequality for any $\zeta>0$, can be bounded as

$\begin{split}J_{15}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{2}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|\left\|\mathcal{A}a(\widehat{Y}_{2k+1}^{f})-\mathcal{A}a(\widehat{Y}_{2k}^{f})\right\|\right|^{p/2}\Bigg{]}\\\
&\leq\widehat{L}^{p/2}\mathbb{E}\Bigg{[}\left|\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
t_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\right|^{p/2}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h^{2}\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+2\zeta\widehat{L}^{p}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h^{2}\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|\right|^{p}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+\zeta
C_{J_{15}}^{p}p^{p/2}h^{3p/2}\operatorname{e}^{\gamma pT}.\end{split}$

The term $J_{16}$, using the Cauchy–Schwarz inequality, Jensen’s inequality,
and Young’s inequality for any $\zeta>0$, yields the bound

$\begin{split}J_{16}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h^{2}\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},\Delta
a(\widehat{Y}_{2k}^{f})\right\rangle\right|^{p/2}\Bigg{]}\\\ &\leq
d^{3p/4}K_{b}^{p/2}\mathbb{E}\Bigg{[}\left|\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
t_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\right|^{p/2}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h^{2}\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+2\zeta
d^{3p/2}K_{b}^{p}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h^{2}\right|^{p}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+\zeta
C_{J_{16}}^{p}h^{p}\operatorname{e}^{\gamma pT}.\end{split}$

Combining all the above estimates, choosing $\zeta=\widetilde{C}$ there exists
a constant $\widetilde{C}_{2}>0$ that is independent of $T$ and $p$, such that

$\begin{split}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq N/2}\operatorname{e}^{\gamma
pt_{2n}}\Bigg{(}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{)}\Bigg{]}\leq\widetilde{C}_{2}^{p}p^{p}h^{p}\operatorname{e}^{\gamma
pT}.\end{split}$

Since the above inequality holds for all $T>0$, we have for all
$0<h<\kappa_{3}$,

$\sup_{0\leq n\leq
N/2}\mathbb{E}\Big{[}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Big{]}\leq\widetilde{C}_{2}^{p}p^{p}h^{p}.$
(48)

Note that since we have used equation (44) to obtain equation (48). The bound
obtained in equation (48) can be further sharpened through argumentum in
recursum.

Using the Mean value theorem, we can write

$a(\widehat{Y}_{2n+1}^{f})-a(\widehat{Y}_{2n}^{f})=Da(\widehat{Y}_{2n}^{f})\Big{(}\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\Big{)}+\frac{1}{2}D^{2}a(\widehat{Y}_{2n}^{f})\Big{(}\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\Big{)}^{2}+R_{3}$
(49)

where $\left\|R_{3}\right\|\leq
d^{2}K_{b}\left\|\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\right\|^{3}$.
Substituting equation (49) in (45), using Hölder’s inequality, Young’s
inequality, and equation (42), we obtain

$\displaystyle\left\|\widehat{Y}_{2n+2}^{f}-\widehat{Y}_{2n+2}^{c}\right\|^{2}$
$\displaystyle\leq(1-4\gamma
h)\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+C\Bigg{\\{}\Big{(}h^{2}\left\|\Delta
W_{2n}\right\|^{2}+\left\|\Delta Z_{2n}\right\|^{2}$
$\displaystyle\quad+\left\|\Delta
Z_{2n+1}\right\|^{2}\Big{)}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}+h^{2}\left\|\Delta
Z_{2n}\right\|^{2}+h^{4}\Big{(}1+\left\|\widehat{Y}_{2n}^{f}\right\|^{2}\Big{)}$
$\displaystyle\quad+\Big{(}h^{4}+\left\|\Delta
Z_{2n+1}\right\|^{2}\Big{)}\left\|\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\right\|^{2}+h^{2}\left\|R_{3}\right\|^{2}\Bigg{\\}}$
$\displaystyle\quad+2(1-4Sh+2S^{2}h^{2})\Bigg{\\{}h\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},R_{3}\right\rangle$
$\displaystyle\quad+h\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
W_{2n}\right\rangle$
$\displaystyle\quad+h\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},Da(\widehat{Y}_{2n}^{f})Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}\right\rangle$
$\displaystyle\quad+(1-Sh)\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n}\right\rangle$
$\displaystyle\quad+\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n+1}^{f})-Da(\widehat{Y}_{2n}^{f})\big{)}\Delta
Z_{2n+1}\right\rangle$
$\displaystyle\quad+\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\big{(}Da(\widehat{Y}_{2n}^{f})-Da(\widehat{Y}_{2n}^{c})\big{)}\Delta
Z_{2n+1}\right\rangle$
$\displaystyle\quad+\frac{h^{2}}{2}\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},\mathcal{A}a(\widehat{Y}_{2n+1}^{f})-\mathcal{A}a(\widehat{Y}_{2n}^{f})\right\rangle+$
$\displaystyle\quad+\frac{h}{2}\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},D^{2}a(\widehat{Y}_{2n}^{f})\big{(}\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\big{)}^{2}-h\Delta
a(\widehat{Y}_{2n}^{f})\right\rangle\Bigg{\\}}.$

Observe that the order of convergence of the terms
$J_{9},J_{10},J_{12},J_{14}$, and $J_{16}$ can be improved using equation (48)
and (49). The term $J_{9}$ using the Cauchy–Schwarz inequality, and Young’s
inequality for any $\zeta>0$, then yields the bound

$\begin{split}J_{9}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},R_{3}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq\mathbb{E}\Bigg{[}\left|\sup_{0\leq n\leq N/2}\operatorname{e}^{\gamma
t_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\right|^{p/2}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h\left\|\widehat{Y}_{2k+1}^{f}-\widehat{Y}_{2k}^{f}\right\|^{3}\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+\zeta
C_{J_{9}}^{p}p^{3p/2}h^{3p/2}\operatorname{e}^{\gamma pT}.\end{split}$

Using equation (48), the bounds for the terms $J_{10},J_{12}$ and $J_{14}$ can
be improved to

$J_{10}\leq C_{J_{10}}^{p}p^{5p/4}h^{3p/2}\operatorname{e}^{\gamma pT},\qquad
J_{12}\leq C_{J_{12}}^{p}p^{5p/4}h^{3p/2}\operatorname{e}^{\gamma pT},\quad
J_{14}\leq C_{J_{14}}^{p}p^{5p/4}h^{3p/2}\operatorname{e}^{\gamma pT}.$

The term $J_{16}$, using equation (49), can be redefined as

$J_{16}\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c},D^{2}a(\widehat{Y}_{2n}^{f})\Big{(}\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\Big{)}^{2}-h\Delta
a(\widehat{Y}_{2n}^{f})\right\rangle\right|^{p/2}\Bigg{]}.$

We can express the term
$D^{2}a(\widehat{Y}_{2n}^{f})\Big{(}\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\Big{)}^{2}$
as

$\begin{split}D^{2}a(\widehat{Y}_{2n}^{f})\Big{(}\widehat{Y}_{2n+1}^{f}-\widehat{Y}_{2n}^{f}\Big{)}^{2}&=D^{2}a(\widehat{Y}_{2n}^{f})\Delta
W_{2n}^{2}+D^{2}a(\widehat{Y}_{2n}^{f})Q_{2n}^{2}+2D^{2}a(\widehat{Y}_{2n}^{f})Q_{2n}\Delta
W_{2n},\end{split}$ (50)

where

$Q_{2n}\coloneqq
Sh(\widehat{Y}_{2n}^{c}-\widehat{Y}_{2n}^{f})+ha(\widehat{Y}_{2n}^{f})+Da(\widehat{Y}_{2n}^{f})\Delta
Z_{2n}+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}_{2n}^{f}),$

and

$\mathbb{E}\Big{[}\left\|Q_{2n}\right\|^{p}\Big{]}\leq C_{Q}^{p}p^{p/2}h^{p}$
(51)

where $C_{Q}>0$ is a constant that depends only on $S,d$, and the constants
from Assumptions 1, 2, and 3. Using equation (50), we can express $J_{16}$ as

$J_{16}\leq
3^{p/2-1}2^{p/2}(\widehat{J}_{1}+\widehat{J}_{2}+\widehat{J}_{3}),$

where $\widehat{J}_{1},\widehat{J}_{2}$ and $\widehat{J}_{3}$ are defined as

$\begin{split}\widehat{J}_{1}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},D^{2}a(\widehat{Y}_{2k}^{f})\Big{(}\Delta
W_{2k}^{2}-\mathbb{E}\Big{[}\Delta
W_{2k}^{2}\Big{]}\Big{)}\right\rangle\right|^{p/2}\Bigg{]},\\\
\widehat{J}_{2}&\coloneqq\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|\left\|D^{2}a(\widehat{Y}_{2k}^{f})Q_{2k}^{2}\right\|\right|^{p/2}\Bigg{]},\\\
\widehat{J}_{3}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},D^{2}a(\widehat{Y}_{2k}^{f})Q_{2k}\Delta
W_{2k}\right\rangle\right|^{p/2}\Bigg{]}.\end{split}$

To bound the term $\widehat{J}_{1}$, we define the Itô process
$\operatorname{d\\!}P_{1}(s)$ in $\mathbb{R}^{d\times d}$, for all
$s\in[0,T]$, such that

$\begin{split}\operatorname{d\\!}P_{1}^{ij}(s)&=\begin{cases}2\Big{(}V_{1}^{i}(s)-V_{1}^{i}\big{(}2\left\lfloor
s/2h\right\rfloor h\big{)}\Big{)}\operatorname{d\\!}V_{1}^{i}(s),&\text{if
}i=j\\\ \Big{(}V_{1}^{i}(s)-V_{1}^{i}\big{(}2\left\lfloor s/2h\right\rfloor
h\big{)}\Big{)}\operatorname{d\\!}V_{1}^{j}(s)+\Big{(}V_{1}^{j}(s)-V_{1}^{j}\big{(}2\left\lfloor
s/2h\right\rfloor h\big{)}\Big{)}\operatorname{d\\!}V_{1}^{i}(s),&\text{if
}i\neq j\end{cases},\\\ P_{1}^{ij}(0)&=0,\end{split}$

for all $i,j=1,2,\ldots,d$. The process $P_{1}$ is a Martingale process [2,
Definition 4.75]. Using the Burkholder–Davis–Gundy inequality, Jensen’s
inequality, and Young’s inequality for any $\zeta>0$, we have

$\begin{split}\widehat{J}_{1}&=\mathbb{E}\Bigg{[}\sup_{0\leq t\leq
T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor s/2h\right\rfloor
h}h\left\langle\widehat{Y}_{2\left\lfloor
s/2h\right\rfloor}^{f}-\widehat{Y}_{2\left\lfloor
s/2h\right\rfloor}^{c},D^{2}a(\widehat{Y}_{2\left\lfloor
s/2h\right\rfloor}^{f})\operatorname{d\\!}P_{1}(s)\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq(2C_{\text{BDG}}p)^{p/4}d^{3p/4}K_{b}^{p/2}\mathbb{E}\Bigg{[}\left|\sup_{0\leq
n\leq N/2}\operatorname{e}^{2\gamma
t_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}\right|^{p/4}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{4}\right|^{p/4}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+C_{\widehat{J}_{1}}^{p}p^{p/2}h^{3p/2}\operatorname{e}^{\gamma
pT}.\end{split}$

For the term $\widehat{J}_{2}$, we have

$\begin{split}\widehat{J}_{2}&\leq
d^{3p/4}K_{b}^{p/2}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|\left\|Q_{2k}\right\|^{2}\right|^{p/2}\Bigg{]}\\\
&\leq d^{3p/4}K_{b}^{p/2}\mathbb{E}\Bigg{[}\left|\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
t_{2k}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|\right|^{p/2}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{\gamma
t_{2k}}h\left\|Q_{2k}\right\|^{2}\right|^{p/2}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2k}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+\zeta
C_{\widehat{J}_{2}}^{p}p^{p}h^{2p}\operatorname{e}^{\gamma pT}.\end{split}$

To bound the term $\widehat{J}_{3}$, we rewrite it as

$\begin{split}\widehat{J}_{3}\leq
2^{p/2-1}\Big{(}\widetilde{J}_{1}+\widetilde{J}_{2}\Big{)},\end{split}$

where

$\begin{split}\widetilde{J}_{1}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},D^{2}a(\widehat{Y}_{2k}^{f})\widetilde{Q}_{2k}\Delta
W_{2k}\right\rangle\right|^{p/2}\Bigg{]},\\\
\widetilde{J}_{2}&\coloneqq\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2-1}\left|\sum_{k=0}^{n}\operatorname{e}^{2\gamma
t_{2k}}h\left\langle\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c},D^{2}a(\widehat{Y}_{2k}^{f})\big{(}Da(\widehat{Y}_{2k}^{f})\Delta
Z_{2k}\big{)}\Delta W_{2k}\right\rangle\right|^{p/2}\Bigg{]},\end{split}$

and

$\widetilde{Q}_{2k}\coloneqq
Sh(\widehat{Y}_{2k}^{c}-\widehat{Y}_{2k}^{f})+ha(\widehat{Y}_{2k}^{f})+\frac{h^{2}}{2}\mathcal{A}a(\widehat{Y}_{2k}^{f}),$

where

$\mathbb{E}\bigg{[}\left\|\widetilde{Q}_{2k}\right\|^{p}\bigg{]}\leq
C_{Q}^{p}p^{p/2}h^{p}.$

For the term $\widetilde{J}_{1}$, using the Burkholder–Davis–Gundy inequality,
Jensen’s inequality, and Young’s inequality for any $\zeta>0$, we obtain

$\begin{split}\widetilde{J}_{1}&=\mathbb{E}\Bigg{[}\sup_{0\leq t\leq
T}\left|\int_{0}^{t}\operatorname{e}^{2\gamma\left\lfloor s/2h\right\rfloor
h}h\left\langle\widehat{Y}_{2\left\lfloor
s/2h\right\rfloor}^{f}-\widehat{Y}_{2\left\lfloor
s/2h\right\rfloor}^{c},D^{2}a(\widehat{Y}_{2\left\lfloor
s/2h\right\rfloor}^{f})\widetilde{Q}_{2\left\lfloor
s/2h\right\rfloor}\operatorname{d\\!}W_{s}\right\rangle\right|^{p/2}\Bigg{]}\\\
&\leq(C_{\text{BDG}}p)^{p/4}\mathbb{E}\Bigg{[}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{4\gamma
t_{2k}}h^{3}\left\|\widehat{Y}_{2k}^{f}-\widehat{Y}_{2k}^{c}\right\|^{2}\left\|D^{2}a(\widehat{Y}_{2k}^{f})\right\|^{2}\left\|\widetilde{Q}_{2k}\right\|^{2}\right|^{p/4}\Bigg{]}\\\
&\leq(C_{\text{BDG}}p)^{p/4}d^{3p/4}K_{b}^{p/2}\mathbb{E}\Bigg{[}\left|\sup_{0\leq
n\leq N/2}\operatorname{e}^{2\gamma
t_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{2}\right|^{p/4}\left|\sum_{k=0}^{N/2-1}\operatorname{e}^{2\gamma
t_{2k}}h^{3}\left\|\widetilde{Q}_{2k}\right\|^{2}\right|^{p/4}\Bigg{]}\\\
&\leq\frac{1}{8\zeta}\mathbb{E}\Bigg{[}\sup_{0\leq n\leq
N/2}\operatorname{e}^{\gamma
pt_{2n}}\left\|\widehat{Y}_{2n}^{f}-\widehat{Y}_{2n}^{c}\right\|^{p}\Bigg{]}+\zeta
C_{\widetilde{J}_{1}}^{p}p^{p}h^{2p}\operatorname{e}^{\gamma pT}.\end{split}$

Using, the Hölder’s inequality, Jensen’s inequality, Young’s inequality for
any $\zeta>0$, we obtain