Multi-level Monte Carlo methods with the truncated Euler–Maruyama scheme for stochastic differential equations

ABSTRACT

The truncated Euler–Maruyama method is employed together with the Multi-level Monte Carlo method to approximate expectations of some functions of solutions to stochastic differential equations (SDEs). The convergence rate and the computational cost of the approximations are proved, when the coefficients of SDEs satisfy the local Lipschitz and Khasminskii-type conditions. Numerical examples are provided to demonstrate the theoretical results.

KEYWORDS:

• Truncated Euler–Maruyama method
• multi-level Monte Carlo method
• stochastic differential equations
• nonlinear coefficients
• approximation to expectation

2010 AMS SUBJECT CLASSIFICATIONS:

• 60H10
• 65C30

1. Introduction

Stochastic differential equations (SDEs) have been broadly discussed and applied as a powerful tool to capture the uncertain phenomenon in the evolution of systems in many areas [2,6,20,25,26]. However, the explicit solutions of SDEs can rarely be found. Therefore, the numerical approximation becomes an essential approach in the applications of SDEs. Monographs [18,23] provide detailed introductions and discussions to various classic methods.

Since the nonlinear coefficients have been widely adapted in SDE models [1,10,24], explicit numerical methods that have good convergence property for SDEs with non-global Lipschitz drift and diffusion coefficients are of interest to many researchers and required by practitioners. The authors in [13] developed a quite general approach to prove the strong convergence of numerical methods for nonlinear SDEs. The approach to prove the global strong convergence via the local convergence for SDEs with non-global Lipschitz coefficients was studied in [29]. More recently, the taming technique was developed to handle the non-global Lipschitz coefficients [15,16]. Simplified proof of the tamed Euler method and the tamed Milstein method can be found in [27] and [30], respectively. The truncated Euler–Maruyama (EM) method was developed in [21,22], which is also targeting on SDEs with non-global Lipschitz coefficients. Explicit methods for nonlinear SDEs that preserve positivity can be found in, for example [12,19]. A modified truncated EM method that preserves the asymptotic stability and boundedness of the nonlinear SDEs was presented in [11].

Compared to the explicit methods mentioned above, the methods with implicit term have better convergence property in approximating non-global Lipschitz SDEs with the trade-off of the relatively expensive computational cost. We just mention a few of the works [14,28,31] and the references therein.

In many situations, the expected values of some functions of the solutions to SDEs are also of interest. To estimate the expected values, the classic Monte-Carlo method is a good and natural candidate. More recently, Giles in [7,8] developed the Multi-level Monte Carlo (MLMC) method, which improves the convergence rate and reduces the computational cost of estimating expected values. A detailed survey of recent developments and applications of the MLMC method can be found in [9]. To complement [9], we only mention some new developments that are not included in [9]. Under the global Lipschitz and linear growth conditions, the MLMC method combined with the EM method applied to SDEs with small noise is often found to be the most efficient option [3]. The MLMC method with the adaptive EM method was designed for solving SDEs driven by Lévy process [4,5]. The MLMC method was applied to SDEs driven by Poisson random measures by means of coupling with the split-step implicit tau-leap at levels. However, the classic EM method with the MLMC method has been proved divergence to SDEs with non-global Lipschitz coefficients [17]. So it is interesting to investigate the combinations of the MLMC method with those numerical methods developed particularly for SDEs with non-global Lipschitz coefficients. In [17], the tamed Euler method was combined with the MLMC method to approximate expectations of some nonlinear functions of solutions to some nonlinear SDEs.

In this paper, we embed the MLMC method with the truncated EM method and study the convergence and the computational cost of this combination to approximate expectations of some nonlinear functions of solutions to SDEs with non-global Lipschitz coefficients.

In [22], the truncated EM method has been proved to converge to the true solution with the order $\frac{1}{2}$-$\epsilon$ for any arbitrarily small $\epsilon >0$. The plan of this paper is as follows. Firstly, we make some modifications of Theorem 3.1 in [8] such that the modified theorem is able to cover the truncated EM method. Then, we use the modified theorem to prove the convergence and the computational cost of the MLMC method with the truncated EM method. At last, numerical examples for SDEs with non-global Lipschitz coefficients and expectations of nonlinear functions are given to demonstrate the theoretical results.

This paper is constructed as follows. Notations, assumptions and some existing results about the truncated EM method and the MLMC method are presented in Section 2. Section 3 contains the main result on the computational complexity. A numerical example is provided in Section 4 to illustrate theoretical results. In the appendix, we give the proof of the theorem in Section 3.

2. Mathematical preliminary

Throughout this paper, unless otherwise specified, we let $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ be a complete probability space with a filtration $\{{\mathcal{F}}_t{\}}_{t\ge 0}$ satisfying the usual condition (that is, it is right continuous and increasing while ${\mathcal{F}}_0$ contains all $\mathbb{P}-$null sets). Let $\mathbb{E}$ denote the expectation corresponding to $\mathbb{P}$. Let $B\left(t\right)$ be an m-dimensional Brownian motion defined on the space. If A is a vector or matrix, its transpose is denoted by $A^T$. If $x\in {\mathbb{R}}^d$, then $|x|$ is the Euclidean norm. If A is a matrix, we let $|A|=\sqrt{\text{trace}\left(A^{T}A\right)}$ be its trace norm. If A is a symmetric matrix, denote by ${\lambda }_{max}\left(A\right)$ and ${\lambda }_{min}\left(A\right)$ its largest and smallest eigenvalue, respectively. Moreover, for two real numbers a and b, set $a\vee b=max(a,b)$ and $a\wedge b=min(a,b)$. If G is a set, its indicator function is denoted by $I_G\left(x\right)=1$ if $x\in G$ and 0 otherwise.

Here we consider an SDE (1) $\mathrm{d}X\left(t\right)=\mu \left(X\right(t\left)\right)\mathrm{d}t+\sigma \left(X\right(t\left)\right)\mathrm{d}B\left(t\right)$(1) on $t\ge 0$ with the initial value $X\left(0\right)=X_0\in {\mathbb{R}}^d$, where $\mu :{\mathbb{R}}^d\to {\mathbb{R}}^d\mathrm{and}\sigma :{\mathbb{R}}^d\to {\mathbb{R}}^{d\times m}.$

When the coefficients obey the global Lipschitz condition, the strong convergence of numerical methods for SDEs has been well studied [18]. When the coefficients μ and σ are locally Lipschitz continuous without the linear growth condition, Mao [21,22] recently developed the truncated EM method. To make this paper self-contained, we give a brief review of this method firstly.

We first choose a strictly increasing continuous function $\omega :{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ such that $\omega \left(r\right)\to \mathrm{\infty }$ as $r\to \mathrm{\infty }$ and (2) $\underset{|x|\le u}{sup}\left(|\mu \right(x)|\vee |\sigma (x\left)|\right)\le \omega \left(u\right),\mathrm{\forall }u\ge 1.$(2) Denote by ${\omega }^{-1}$ the inverse function of ω and we see that ${\omega }^{-1}$ is a strictly increasing continuous function from $\left[\omega \right(0),\mathrm{\infty })$ to $R_{+}$. We also choose a number $s_l^{\ast }\in (0,1]$ and a strictly decreasing function $h:(0,s_l^{\ast }]\to (0,\mathrm{\infty })$ such that (3) $h\left(s_l^{\ast }\right)\ge \omega \left(2\right),\underset{s_l\to 0}{lim}h\left(s_l\right)=\mathrm{\infty }\mathrm{and}{s}_l^{1/4}h\left(s_l\right)\le 1,\mathrm{\forall }{s}_l\in (0,s_l^{\ast }].$(3) For a given stepsize $s_l\in (0,1)$, let us define the truncated functions ${\mu }_{s_l}\left(x\right)=\mu \left(\right(|x|\wedge {\omega }^{-1}\left(h\right(s_l\left)\right)\left)\frac{x}{|x|}\right)\text{and}{\sigma }_{s_l}\left(x\right)=\sigma \left(\right(|x|\wedge {\omega }^{-1}\left(h\right(s_l\left)\right)\left)\frac{x}{|x|}\right)$ for $x\in {\mathbb{R}}^d$, where we set $x/|x|=0$ when x=0. Moreover, let ${\overline{X}}_{s_l}\left(t\right)$ denote the approximation to $X\left(t\right)$ using the truncated EM method with time step size $s_l=M^{-l}T$ for $l=0,1,\dots ,L$. The numerical solutions ${\overline{X}}_{s_l}\left(t_k\right)$ for $t_k=k{s}_l$ are formed by setting ${\overline{X}}_{s_l}\left(0\right)=X_0$ and computing (4) ${\overline{X}}_{s_l}\left(t_{k+1}\right)={\overline{X}}_{s_l}\left(t_k\right)+{\mu }_{s_l}\left({\overline{X}}_{s_l}\right(t_k\left)\right)s_l+{\sigma }_{s_l}\left({\overline{X}}_{s_l}\right(t_k\left)\right)\mathrm{\Delta }{B}_k$(4) for $k=0,1,\dots ,$ where $\mathrm{\Delta }{B}_k=B\left(t_{k+1}\right)-B\left(t_k\right)$ is the Brownian motion increment.

Now we give some assumptions to guarantee that the truncated EM solution (4(4) ${\overline{X}}_{s_l}\left(t_{k+1}\right)={\overline{X}}_{s_l}\left(t_k\right)+{\mu }_{s_l}\left({\overline{X}}_{s_l}\right(t_k\left)\right)s_l+{\sigma }_{s_l}\left({\overline{X}}_{s_l}\right(t_k\left)\right)\mathrm{\Delta }{B}_k$(4) ) will converge to the true solution to the SDE (1(1) $\mathrm{d}X\left(t\right)=\mu \left(X\right(t\left)\right)\mathrm{d}t+\sigma \left(X\right(t\left)\right)\mathrm{d}B\left(t\right)$(1) ) in the strong sense.

Assumption 2.1

The coefficients μ and σ satisfy the local Lipschitz condition that for any real number R>0, there exists a $K_R>0$ such that (5) $|\mu \left(x\right)-\mu \left(y\right)|\vee |\sigma \left(x\right)-\sigma \left(y\right)|\le K_R|x-y|$(5) for all $x,y\in R^d$ with $|x|\vee |y|\le R$.

Assumption 2.2

The coefficients μ and σ satisfy the Khasminskii-type condition that there exists a pair of constants p>2 and K>0 such that (6) $x^T\mu \left(x\right)+\frac{p-1}{2}|\sigma \left(x\right){|}^2\le K(1+|x{|}^2)$(6) for all $x\in {\mathbb{R}}^d$.

Assumption 2.3

There exists a pair of constants $q\ge 2$ and $H_1>0$ such that (7) $(x-y{)}^T(\mu \left(x\right)-\mu \left(y\right))+\frac{q-1}{2}|\sigma (x)-\sigma (y){|}^2\le H_1|x-y{|}^2$(7) for all $x,y\in {\mathbb{R}}^d$.

Assumption 2.4

There exists a pair of positive constants ρ and $H_2$ such that (8) $|\mu \left(x\right)-\mu \left(y\right){|}^2\vee |\sigma \left(x\right)-\sigma \left(y\right){|}^2\le H_2(1+|x{|}^{\rho }+|y{|}^{\rho })|x-y{|}^2$(8) for all $x,y\in {\mathbb{R}}^d$.

Let $f\left(X\right(t\left)\right)$ denote a payoff function of the solution to some SDE driven by a given Brownian path $B\left(t\right)$. In this paper, we need f satisfies the following assumption.

Assumption 2.5

There exists a constant c>0 such that (9) $|f\left(x\right)-f\left(y\right)|\le c(1+|x{|}^c+|y{|}^c)|x-y|$(9) for all $x,y\in {\mathbb{R}}^d$.

Using the idea in [7,8], the expected value of $f\left({\overline{X}}_{s_l}\right(t\left)\right)$ can be decomposed in the following way (10) $\mathbb{E}\left[f\right({\overline{X}}_{s_L}\left(T\right)\left)\right]=\mathbb{E}\left[f\right({\overline{X}}_{s_0}\left(T\right)\left)\right]+\sum _{l=1}^L\mathbb{E}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f({\overline{X}}_{s_{l-1}}\left(T\right)\left)\right].$(10)

Let $Y_0$ be an estimator for $\mathbb{E}\left[f\right({\overline{X}}_{s_0}\left(T\right)\left)\right]$ using $N_0$ samples. Let $Y_l$ be an estimator for $\mathbb{E}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f({\overline{X}}_{s_{l-1}}\left(T\right)\left)\right]$ using $N_l$ paths such that $Y_l=\frac{1}{N_l}\sum _{i=1}^{N_l}\left[f\right({\overline{X}}_{s_l}^{\left(i\right)}\left(T\right))-f({\overline{X}}_{s_{l-1}}^{\left(i\right)}\left(T\right)\left)\right].$ The multi-level method independently estimates each of the expectations on the right-hand side of Equation (10(10) $\mathbb{E}\left[f\right({\overline{X}}_{s_L}\left(T\right)\left)\right]=\mathbb{E}\left[f\right({\overline{X}}_{s_0}\left(T\right)\left)\right]+\sum _{l=1}^L\mathbb{E}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f({\overline{X}}_{s_{l-1}}\left(T\right)\left)\right].$(10) ) such that the computational complexity can be minimized, see [8] for more details.

3. Main results

In this section, Theorem 3.1 in [8] is slightly generalized. Then the convergence rate and computational complexity of the truncated EM method combined with the MLMC method are studied.

3.1. Generalized theorem for the MLMC method

Theorem 3.1

If there exist independent estimators $Y_l$ based on $N_l$ Monte Carlo samples, and positive constants $\alpha ,\beta ,c_1,c_2,c_3$ such that

1. $\mathbb{E}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f(X\left(T\right)\left)\right]\le c_1{s}_l^{\alpha },$

2. $\mathbb{E}\left[Y_l\right]=\{\begin{array}{ll}\mathbb{E}\left[f\right({\overline{X}}_{s_0}\left(T\right)\left)\right],& l=0,\\ \mathbb{E}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f({\overline{X}}_{s_{l-1}}\left(T\right)\left)\right],& l>0,\end{array}$

3. $\mathrm{Var}\left[Y_l\right]\le c_2{N}_l^{-1}{s}_l^{\beta },$

4. the computational complexity of $Y_l$, denoted by $C_l,$ is bounded by $C_l\le c_3{N}_l{s}_l^{-1},$

then there exists a positive constant $c_4$ such that for any $\epsilon <e^{-1}$ the multi-level estimator $Y=\sum _{l=0}^L{Y}_l$ has a mean square error $($MSE$)$ $\mathrm{MSE}\equiv \mathbb{E}\left[\right(Y-\mathbb{E}\left[f\right(X\left(T\right)\left)\right]{)}^2]<{\epsilon }^2.$

Furthermore, the upper bound of computational complexity of Y, denoted by C, is given by $C\le \{\begin{array}{ll}{c}_3(2c_5^2{c}_2+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }){\epsilon }^{-1/\alpha },& \alpha \le (-\log \epsilon )/\log \left[\right(\log \epsilon /\epsilon {)}^2],\\ c_3(2c_5^2{c}_2+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }){\epsilon }^{-2}(\log \epsilon {)}^2,& \alpha >(-log\epsilon )/\log \left[\right(\log \epsilon /\epsilon {)}^2]\end{array}$ for $\beta =1$, $C\le \{\begin{array}{ll}{c}_3\left[2c_2{T}^{\beta -1}\right(1-M^{-(\beta -1)/2}{)}^{-2}+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }\right]{\epsilon }^{-2},& \alpha \ge \frac{1}{2},\\ c_3\left[2c_2{T}^{\beta -1}\right(1-M^{-(\beta -1)/2}{)}^{-2}+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }\right]{\epsilon }^{-1/\alpha },& \alpha <\frac{1}{2}\end{array}$ for $\beta >1$, and $C\le \{\begin{array}{ll}{c}_3\left[2c_2\right(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }(1-M^{-(1-\beta )/2}{)}^{-2}+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }]{\epsilon }^{-2-(1-\beta )/\alpha },& \beta \le 2\alpha ,\\ c_3\left[2c_2\right(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }(1-M^{-(1-\beta )/2}{)}^{-2}+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }]{\epsilon }^{-1/\alpha },& \beta >2\alpha \end{array}$ for $0<\beta <1$.

The proof is in the appendix.

Remark 3.1

The main difference of Theorem 3.1 and Theorem 3.1 in [8] lies in the first condition. In [8], one needs $\alpha \ge \frac{1}{2}$. In this paper, this requirement is weaken by any $\alpha >0$.

3.2. Specific theorem for truncated Euler with the MLMC

Next we consider the MLMC path simulation with truncated EM method and discuss their computational complexity using Theorem 3.1.

From Theorem 3.8 in [22], under Assumptions 2.1–2.4, for every small $s_l\in (0,s_l^{\ast })$, where $s_l^{\ast }\in (0,1)$ and for any real number T>0, we have (11) $\mathbb{E}|X\left(T\right)-{\overline{X}}_{s_l}\left(T\right){|}^{\overline{q}}\le c{s}_l^{\overline{q}/2}\left(h\right(s_l){)}^{\overline{q}},$(11) for $\overline{q}\ge 2$. If $\overline{q}=1$, by using the Holder inequality, we also know that $\mathbb{E}|X\left(T\right)-{\overline{X}}_{s_l}\left(T\right)|\le {\left(\mathbb{E}|X\right(T)-{\overline{X}}_{s_l}(T\left){|}^2\right)}^{1/2}\le \left(c{s}_l\right(h\left(s_l\right){)}^2{)}^{1/2}=c{s}_l^{1/2}h\left(s_l\right),$ so we can obtain (12) $\begin{array}{rl}& \mathbb{E}\left[|f\right({\overline{X}}_{s_l}\left(T\right))-f(X\left(T\right)\left)|\right]\\ & \le \mathbb{E}\left[c\right(1+|{\overline{X}}_{s_l}\left(T\right){|}^c+|X\left(T\right){|}^c\left)|{\overline{X}}_{s_l}\right(T)-X(T\left)|\right]\le c\left(\mathbb{E}|{\overline{X}}_{s_l}\right(T)-X(T){|}^2{)}^{1/2}\\ & \le c{s}_l^{1/2}h\left(s_l\right)\end{array}$(12) with the polynomial growth condition (9(9) $|f\left(x\right)-f\left(y\right)|\le c(1+|x{|}^c+|y{|}^c)|x-y|$(9) ). This implies that $\alpha =\frac{1}{4}$ for the truncated EM scheme.

Next we consider the variance of $Y_l$. It follows that (13) $\mathrm{Var}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f(X\left(T\right)\left)\right]\le \mathbb{E}\left[\right(f\left({\overline{X}}_{s_l}\right(T\left)\right)-f\left(X\right(T\left){)}^2\right]\le c{s}_l\left(h\right(s_l){)}^2$(13) using Equations (9(9) $|f\left(x\right)-f\left(y\right)|\le c(1+|x{|}^c+|y{|}^c)|x-y|$(9) ) and (11(11) $\mathbb{E}|X\left(T\right)-{\overline{X}}_{s_l}\left(T\right){|}^{\overline{q}}\le c{s}_l^{\overline{q}/2}\left(h\right(s_l){)}^{\overline{q}},$(11) ).

In addition, it can be noted that $f\left({\overline{X}}_{s_l}\right(T\left)\right)-f\left({\overline{X}}_{s_{l-1}}\right(T\left)\right)=\left[f\right({\overline{X}}_{s_l}\left(T\right))-f(X\left(T\right)\left)\right]-\left[f\right({\overline{X}}_{s_{l-1}}\left(T\right))-f(X\left(T\right)\left)\right],$

thus we have $\begin{array}{rl}& \mathrm{Var}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f({\overline{X}}_{s_{l-1}}\left(T\right)\left)\right]\\ & \le (\sqrt{\mathrm{Var}\left[f\right({\overline{X}}_{s_l}\left(T\right))-f(X\left(T\right)\left)\right]}+\sqrt{\mathrm{Var}\left[f\right({\overline{X}}_{s_{l-1}}\left(T\right))-f(X\left(T\right)\left)\right]}{)}^2\\ & \le c{s}_l\left(h\right(s_l){)}^2+c{s}_{l-1}(h\left(s_{l-1}\right){)}^2\\ & \le c{s}_l^{1/2},\end{array}$ where the fact $s_l^{1/4}h\left(s_l\right)\le 1$ from Equation (3(3) $h\left(s_l^{\ast }\right)\ge \omega \left(2\right),\underset{s_l\to 0}{lim}h\left(s_l\right)=\mathrm{\infty }\mathrm{and}{s}_l^{1/4}h\left(s_l\right)\le 1,\mathrm{\forall }{s}_l\in (0,s_l^{\ast }].$(3) ) is used.

Now we have $\mathrm{Var}\left[Y_l\right]=N_l^{-1}\mathrm{Var}\left[f\right({\overline{X}}_{s_l}^{\left(i\right)}\left(T\right))-f({\overline{X}}_{s_{l-1}}^{\left(i\right)}\left(T\right)\left)\right]\le c{N}_l^{-1}{s}_l^{1/2}.$ So we have $\beta =\frac{1}{2}$ for the truncated EM method.

According to the Theorem 3.1, it is easy to find that the upper bound of the computational complexity of Y is $\left[4c_1^2{c}_2{c}_3\sqrt{M}\right(1-M^{-1/4}{)}^{-2}+\frac{4M^2}{M-1}{c}_1^4{c}_3]{\epsilon }^{-4}.$

4. Numerical simulations

To illustrate the theoretical results, we consider a nonlinear scalar SDE (14) $\mathrm{d}x\left(t\right)=\left(x\right(t)-x^3(t\left)\right)\mathrm{d}t+|x\left(t\right){|}^{3/2}\mathrm{d}B\left(t\right),t\ge 0,x\left(0\right)=x_0\in \mathbb{R},$(14) where $B\left(t\right)$ is a scalar Brownian motion. This is a specified Lewis stochastic volatility model. According to Examples 3.5 and 3.9 in [22], we sample over 1000 discretized Brownian paths and use stepsizes $s_l=T/2^l$ for $l=1,2,\dots ,5$ in the truncated EM method. Let $\widehat{Y_l}$ denote the sample value of $Y_l$. Here we set T=1 and $h\left(s_l\right)=s_l^{-1/4}$.

Firstly, we show some computational results of the classic EM method with the MLMC method.

It can be seen from Table 1 that the simulation result of (14(14) $\mathrm{d}x\left(t\right)=\left(x\right(t)-x^3(t\left)\right)\mathrm{d}t+|x\left(t\right){|}^{3/2}\mathrm{d}B\left(t\right),t\ge 0,x\left(0\right)=x_0\in \mathbb{R},$(14) ) computed by the MLMC approach together with the classic EM method is divergent.

Table 1. Numerical results using the MLMC with the classic EM method.

l	1	2	3	4	5
${\hat{Y}}_l$	1.00	$2.59\mathrm{e}+102$	$-2.94\mathrm{e}+159$	–	–

The simulation results using the MLMC method combined with the truncated EM method is presented in Table 2. It is clear that some convergent trend is displayed.

Table 2. Numerical results using the MLMC with the truncated EM method.

l	1	2	3	4	5
${\hat{Y}}_l$	0.39	−0.18	−0.024	−0.003	−0.0006

Next, it is noted that compared with the standard Monte Carlo method the computational cost can be saved by using MLMC method. From Figure 1, we can see that the MLMC method is approximately 10 times more efficient than the standard Monte Carlo method when $\epsilon$ is sufficient small.

Figure 1. Computational cost.

Acknowledgements

The authors would like to thank the referee and editor for their very useful comments and suggestions, which have helped to improve this paper a lot.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

X. Mao  http://orcid.org/0000-0002-6768-9864

Additional information

Funding

The authors would also like to thank the Natural Science Foundation of Shanghai [14ZR1431300], the Royal Society [Wolfson Research Merit Award WM160014], the Leverhulme Trust [RF-2015-385], the Royal Society of London [IE131408], the EPSRC [EP/E009409/1], the Royal Society of Edinburgh [RKES115071], the London Mathematical Society [11219], the Edinburgh Mathematical Society [RKES130172], the Ministry of Education (MOE) of China [MS2014DHDX020], Shanghai Pujiang Program [16PJ1408000], the Natural Science Fund of Shanghai Normal University [SK201603], Young Scholar Training Program of Shanghai's Universities for their financial support.

Appendix

Proof

Proof of Theorem 3.1

Using the notation $\lceil x\rceil$ to denote the unique integer n satisfying the inequalities $x\le n<x+1$, we start by choosing L to be $L=\lceil \frac{\log \left(\sqrt{2}{c}_1{T}^{\alpha }{\epsilon }^{-1}\right)}{\alpha \log M}\rceil ,$ so that $\frac{1}{\sqrt{2}}{M}^{-\alpha }\epsilon <c_1{s}_L^{\alpha }\le \frac{1}{\sqrt{2}}\epsilon .$ Hence, by the condition 1 and 2 we have (A1) $\begin{array}{rl}& \left(\mathbb{E}\right[Y]-\mathbb{E}[f\left(X\right(T\left)\right)]{)}^2\\ & ={(\mathbb{E}\left[\sum _{l=0}^L{Y}_l\right]-\mathbb{E}[f\left(X\right(T\left)\right)\left]\right)}^2\\ & =\left(\mathbb{E}\right[f\left({\overline{X}}_{s_L}\right(T\left)\right)-f\left(X\right(T\left)\right)]{)}^2\\ & \le (c_1{s}_L^{\alpha }{)}^2\triangleq \frac{1}{2}{\epsilon }^2.\end{array}$(A1) Therefore, we have $\left(\mathbb{E}\right[Y]-\mathbb{E}[f\left(X\right)]{)}^2\le \frac{1}{2}{\epsilon }^2.$ This upper bound on the square of bias error together with the upper bound of $\frac{1}{2}{\epsilon }^2$ on the variance of the estimator, which will be proved later, gives a upper bound of ${\epsilon }^2$ to the MSE.

Noting $\sum _{l=0}^L{s}_l^{-1}=s_L^{-1}\sum _{i=0}^L{M}^{-i}<\frac{M}{M-1}{s}_L^{-1},$ using the standard result for a geometric series and the inequality $\left(1/\sqrt{2}\right)M^{-\alpha }\epsilon <c_1{s}_L^{\alpha }$, we can obtain $s_L^{-1}<M{\left(\frac{\epsilon }{\sqrt{2}{c}_1}\right)}^{-1/\alpha }.$ Then, we have (A2) $\sum _{l=0}^L{s}_l^{-1}<\frac{M}{M-1}{s}_L^{-1}<\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }.$(A2)

We now consider the different possible values of β and to compare them to the α.

(a) If $\beta =1$, we set $N_l=\lceil 2{\epsilon }^{-2}(L+1)c_2{s}_l\rceil$ so that $V\left[Y\right]=\sum _{l=0}^{L}V\left[Y_l\right]\le \sum _{l=0}^L{c}_2{N}_l^{-1}{s}_l\le \frac{1}{2}{\epsilon }^2,$ which is the required.

For the bound of the computational complexity C, we have $\begin{array}{rl}C& =\sum _{l=0}^L{C}_l\le c_3\sum _{l=0}^L{N}_l{s}_l^{-1}\\ & \le c_3\sum _{l=0}^L\left(2{\epsilon }^{-2}\right(L+1)c_2{s}_l+1)s_l^{-1}\\ & \le c_3\left(2{\epsilon }^{-2}\right(L+1{)}^2{c}_2+\sum _{l=0}^L{s}_l^{-1})\\ & \le c_3\left(2{\epsilon }^{-2}\right(L+1{)}^2{c}_2+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }\right).\end{array}$ According to the definition of L, we have $L\le \frac{\log {\epsilon }^{-1}}{\alpha \log M}+\frac{\log \left(\sqrt{2}{c}_1{T}^{\alpha }\right)}{\alpha \log M}+1.$ Given that $1<\log {\epsilon }^{-1}$ for $\epsilon <e^{-1}$, we have $L+1\le c_5\log {\epsilon }^{-1},$ where $c_5=\frac{1}{\alpha \log M}+max(0,\frac{log\left(\sqrt{2}{c}_1{T}^{\alpha }\right)}{\alpha \log M})+2.$ Hence, the computation complexity is bounded by $\begin{array}{rl}C& \le c_3\left(2{\epsilon }^{-2}{c}_5^2\right(\log {\epsilon }^{-1}{)}^2{c}_2+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }\right)\\ & =c_3\left(2{\epsilon }^{-2}{c}_5^2\right(\log \epsilon {)}^2{c}_2+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }\right).\end{array}$ So if $\alpha \le (-log\epsilon )/log\left[\right(log\epsilon /\epsilon {)}^2]$, we have $C\le c_3(2c_5^2{c}_2+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }){\epsilon }^{-1/\alpha }.$ If $\alpha >(-\log \epsilon )/\log \left[\right(\log \epsilon /\epsilon {)}^2]$, we have $C\le c_3(2c_5^2{c}_2+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }){\epsilon }^{-2}(\log \epsilon {)}^2.$

(b) For $\beta >1$, setting $N_l=\lceil 2{\epsilon }^{-2}{c}_2{T}^{(\beta -1)/2}(1-M^{-(\beta -1)/2}{)}^{-1}{s}_l^{(\beta +1)/2}\rceil ,$ then we have $\begin{array}{rl}V\left[Y\right]& =\sum _{l=0}^{L}V\left[Y_l\right]\le \sum _{l=0}^L{c}_2{N}_l^{-1}{s}_l^{\beta }\\ & \le \frac{1}{2}{\epsilon }^2{T}^{-(\beta -1)/2}(1-M^{-(\beta -1)/2})\sum _{l=0}^L{s}_l^{(\beta -1)/2}.\end{array}$

Using the stand result for a geometric series (A3) $\begin{array}{rl}\sum _{l=0}^L{s}_l^{(\beta -1)/2}& =T^{(\beta -1)/2}\sum _{l=0}^L(M^{-(\beta -1)/2}{)}^l\\ & <T^{(\beta -1)/2}(1-M^{-(\beta -1)/2}{)}^{-1},\end{array}$(A3) we obtain that the upper bound of variance is $\frac{1}{2}{\epsilon }^2$. So the computation complexity is bounded by $\begin{array}{rl}C& \le c_3\sum _{l=0}^L{N}_l{s}_l^{-1}\\ & \le c_3\sum _{l=0}^L\left(2{\epsilon }^{-2}{c}_2{T}^{(\beta -1)/2}\right(1-M^{-(\beta -1)/2}{)}^{-1}{s}_l^{(\beta +1)/2}+1)s_l^{-1}\\ & =c_3\left[2{\epsilon }^{-2}{c}_2{T}^{(\beta -1)/2}\right(1-M^{-(\beta -1)/2}{)}^{-1}\sum _{l=0}^L{s}_l^{(\beta -1)/2}+\sum _{l=0}^L{s}_l^{-1}]\\ & \le c_3\left[2{\epsilon }^{-2}{c}_2{T}^{(\beta -1)/2}\right(1-M^{-(\beta -1)/2}{)}^{-1}{T}^{(\beta -1)/2}(1-M^{-(\beta -1)/2}{)}^{-1}\\ & +\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }\right]\\ & =c_3\left[2{\epsilon }^{-2}{c}_2{T}^{\beta -1}\right(1-M^{-(\beta -1)/2}{)}^{-2}+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }\right].\end{array}$ So when $\alpha \ge \frac{1}{2}$, we have $C\le c_3\left[2c_2{T}^{\beta -1}\right(1-M^{-(\beta -1)/2}{)}^{-2}+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }\right]{\epsilon }^{-2},$ When $\alpha <\frac{1}{2}$, we have $C\le c_3\left[2c_2{T}^{\beta -1}\right(1-M^{-(\beta -1)/2}{)}^{-2}+\frac{M^2}{M-1}\left(\sqrt{2}{c}_1{)}^{1/\alpha }\right]{\epsilon }^{-1/\alpha }.$

(c) For $0<\beta <1$, setting $N_l=\lceil 2{\epsilon }^{-2}{c}_2{s}_L^{-(1-\beta )/2}(1-M^{-(1-\beta )/2}{)}^{-1}{s}_l^{(\beta +1)/2}\rceil ,$ then we have $\begin{array}{rl}V\left[Y\right]& =\sum _{l=0}^{L}V\left[Y_l\right]\le \sum _{l=0}^L{c}_2{N}_l^{-1}{s}_l^{\beta }\\ & \le \frac{1}{2}{\epsilon }^2{s}_L^{(1-\beta )/2}(1-M^{-(1-\beta )/2})\sum _{l=0}^L{s}_l^{-(1-\beta )/2}.\end{array}$ Because (A4) $\begin{array}{rl}\sum _{l=0}^L{s}_l^{-(1-\beta )/2}& =s_L^{-(1-\beta )/2}\sum _{l=0}^L(M^{-(1-\beta )/2}{)}^l\\ & <s_L^{-(1-\beta )/2}(1-M^{-(1-\beta )/2}{)}^{-1},\end{array}$(A4) we obtain the upper bound on the variance of the estimator to be $\frac{1}{2}{\epsilon }^2$.

Finally, using the upper bound of $N_l$, the computational complexity is $\begin{array}{rl}C& \le c_3\sum _{l=0}^L{N}_l{s}_l^{-1}\\ & \le c_3\sum _{l=0}^L\left(2{\epsilon }^{-2}{c}_2{s}_L^{-(1-\beta )/2}\right(1-M^{-(1-\beta )/2}{)}^{-1}{s}_l^{(\beta +1)/2}+1)s_l^{-1}\\ & =c_3\left[2{\epsilon }^{-2}{c}_2{s}_L^{-(1-\beta )/2}\right(1-M^{-(1-\beta )/2}{)}^{-1}\sum _{l=0}^L{s}_l^{-(1-\beta )/2}+\sum _{l=0}^L{s}_l^{-1}]\\ & \le c_3\left[2{\epsilon }^{-2}{c}_2{s}_L^{-(1-\beta )}\right(1-M^{-(1-\beta )/2}{)}^{-2}+\sum _{l=0}^L{s}_l^{-1}],\end{array}$ where (A4) is used in the last inequality.

Moreover, because of the inequality $\left(1/\sqrt{2}\right)M^{-\alpha }\epsilon <c_1{s}_L^{\alpha }$, we have $s_L^{-(1-\beta )}<(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }{\epsilon }^{-(1-\beta )/\alpha },$ then $\begin{array}{rl}C& \le c_3\left[2{\epsilon }^{-2}{c}_2{s}_L^{-(1-\beta )}\right(1-M^{-(1-\beta )/2}{)}^{-2}+\sum _{l=0}^L{s}_l^{-1}]\\ & \le c_3\left[2{\epsilon }^{-2}{c}_2\right(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }{\epsilon }^{-(1-\beta )/\alpha }(1-M^{-(1-\beta )/2}{)}^{-2}+\sum _{l=0}^L{s}_l^{-1}]\\ & \le c_3\left[2{\epsilon }^{-2}{c}_2\right(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }{\epsilon }^{-(1-\beta )/\alpha }(1-M^{-(1-\beta )/2}{)}^{-2}+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }]\\ & =c_3\left[2c_2\right(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }(1-M^{-(1-\beta )/2}{)}^{-2}{\epsilon }^{-2-(1-\beta )/\alpha }+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }{\epsilon }^{-1/\alpha }].\end{array}$ If $\beta \le 2\alpha$, then ${\epsilon }^{-2-(1-\beta )/\alpha }>{\epsilon }^{-1/\alpha }$, so we have $C\le c_3\left[2c_2\right(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }(1-M^{-(1-\beta )/2}{)}^{-2}+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }]{\epsilon }^{-2-(1-\beta )/\alpha }.$ If $\beta >2\alpha$, then ${\epsilon }^{-2-(1-\beta )/\alpha }<{\epsilon }^{-1/\alpha }$, so we have $C\le c_3\left[2c_2\right(\sqrt{2}{c}_1{)}^{(1-\beta )/\alpha }{M}^{1-\beta }(1-M^{-(1-\beta )/2}{)}^{-2}+\frac{M^2}{M-1}(\sqrt{2}{c}_1{)}^{1/\alpha }]{\epsilon }^{-1/\alpha }.◼$