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.


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,21,26,27].However, the explicit solutions of SDEs can rarely be found.Therefore, the numerical approximation becomes an essential approach in the applications of SDEs.Monographs [19,24] provide detailed introductions and discussions to various classic methods.
Since the non-linear coefficients have been widely adapted in SDE models [1,10,25], 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 [14] 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 [30].More recently, the taming technique was developed to handle the nonglobal Lipschitz coefficients [16,17].Simplified proof of the tamed Euler method and the tamed Milstein method can be found in [28] and [31], respectively.The truncated Euler-Maruyama (EM) method was developed in [22,23], 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,20].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 [15,29,32] 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 [18].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 [18], 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 [23], the truncated EM method has been proved to converge to the true solution with the order 1/2 − ε for any arbitrarily small ε > 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.

Mathematical Preliminary
Throughout this paper, unless otherwise specified, we let (Ω , F , P) be a complete probability space with a filtration {F t } t≥0 satisfying the usual condition (that is, it is right continuous and increasing while F 0 contains all P−null sets).Let E denote the expectation corresponding to P. Let B(t) 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 ∈ R d , then |x| is the Euclidean norm.If A is a matrix, we let |A| = trace(A T A) be its trace norm.If A is a symmetric matrix, denote by λ max (A) and λ min (A) its largest and smallest eigenvalue, respectively.Moreover, for two real numbers a and b, set Here we consider an SDE on t ≥ 0 with the initial value X(0) = X 0 ∈ R d , where When the coefficients obey the global Lipschitz condition, the strong convergence of numerical methods for SDEs has been well studied [19].When the coefficients µ and σ are locally Lipschitz continuous without the linear growth condition, Mao [22,23] 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 ω : ( Denote by ω −1 the inverse function of ω and we see that ω −1 is a strictly increasing continuous function from [ω(0), ∞) to R + .We also choose a number s * l ∈ (0, 1] and a strictly decreasing function h : (0, For a given stepsize s l ∈ (0, 1), let us define the truncated functions x |x| for x ∈ R d , where we set x/|x| = 0 when x = 0.Moreover, let X s l (t) denote the approximation to X(t) using the truncated EM method with time step size s l = M −l T for l = 0, 1, . . ., L. The numerical solutions X s l (t k ) for t k = ks l are formed by setting X s l (0) = X 0 and computing for k = 0, 1, . . ., where ∆B k = B(t k+1 ) − B(t k ) is the Brownian motion increment.Now we give some assumptions to guarantee that the truncated EM solution (4) will converge to the true solution to the SDE (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 The coefficients µ and σ satisfy the Khasminskii-type condition that there exists a pair of constants p > 2 and K > 0 such that for all x ∈ R d .
Assumption 2.3 There exists a pair of constants q ≥ 2 and H 1 > 0 such that for all x, y ∈ R d .
Assumption 2.4 There exists a pair of positive constants ρ and H 2 such that for all x, y ∈ R d .
Let f (X(t)) denote a payoff function of the solution to some SDE driven by a given Brownian path B(t).In this paper, we need f satisfies the following assumption.
Assumption 2.5 There exists a constant c > 0 such that for all x, y ∈ R d .
Using the idea in [7,8], the expected value of f (X s l (t)) can be decomposed in the following way Let Y 0 be an estimator for E[f (X s0 (T ))] using N 0 samples.Let Y l be an estimator for The multi-level method independently estimates each of the expectations on the right-hand side of (10) such that the computational complexity can be minimized, see [8] for more details.

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

Generalised theorem for the MLMC method
then there exists a positive constant c 4 such that for any ε < e −1 the multi-level estimator has a mean square error (MSE) Furthermore, the upper bound of computational complexity of Y , denoted by C, is given by for β > 1, and The proof is in Appendix.

Remark 3.2
The main difference of Theorem 3.1 and Theorem 3.1 in [8] lies in the first condition.In [8], one needs α ≥ 1 2 .In this paper, this requirement is weaken by any α > 0.

Specific theorem for truncated Euler with the MLMC
Next we consider the multi-level Monte Carlo path simulation with truncated EM method and discuss their computational complexity using Theorem 3.1.
From Theorem 3.8 in [23], under Assumptions 2.1-2.4,for every small s l ∈ (0, s * l ), where s * l ∈ (0, 1) and for any real number T > 0, we have for q ≥ 2. If q = 1, by using the Holder inequality, we also know that with the polynomial growth condition (9).This implies that α = 1/4 for the truncated EM scheme.
Next we consider the variance of Y l .It follows that using ( 9) and (11).In addition, it can be noted that thus we have where the fact s So we have β = 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

Numerical Simulations
To illustrate the theoretical results, we consider a non-linear scalar SDE where B(t) is a scalar Brownian motion.This is a specified Lewis stochastic volatility model.According to Examples 3.5 and 3.9 in [23], we sample over 1000 discretized Brownian paths and use stepsizes s l = T /2 l for l = 1, 2, . . ., 5 in the truncated EM method.Let Ŷl denote the sample value of Y l .Here we set T = 1 and h(s l ) = s −1/4 l .Firstly, we show some computational results of the classic EM method with the MLMC method.  1 that the simulation result of ( 14) computed by the MLMC approach together with the classic EM method is divergent.
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 Ŷ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 ε is sufficient small.

Appendix
Proof of Theorem 3.1.
Using the notation ⌈x⌉ to denote the unique integer n satisfying the inequalities x ≤ n < x + 1, we start by choosing L to be Hence, by the condition 1 and 2 we have Therefore, we have This upper bound on the square of bias error together with the upper bound of 1 2 ε 2 on the variance of the estimator, which will be proved later, gives a upper bound of ε 2 to the MSE.Noting using the standard result for a geometric series and the inequality 1 Then, we have We now consider the different possible values of β and to compare them to the α.
which is the required.
For the bound of the computational complexity C, we have According to the definition of L, we have Given that 1 < logε −1 for ε < e −1 , we have where Hence, the computation complexity is bounded by then we have Using the stand result for a geometric series we obtain that the upper bound of variance is 1 2 ε 2 .So the computation complexity is bounded by So when α ≥ 1 2 , we have then we have we obtain the upper bound on the variance of the estimator to be 1 2 ε 2 .Finally, using the upper bound of N l , the computational complexity is where (18) is used in the last inequality.Moreover, because of the inequality 1 If β ≤ 2α, then ε −2−(1−β)/α > ε − 1 α , so we have If β > 2α, then ε −2−(1−β)/α < ε − 1 α , so we have

Theorem 3 . 1
If there exist independent estimators Y l based on N l Monte Carlo samples, and positive constants α, β, c 1 , c 2 , c 3 such that

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