Multilevel MC method for weak approximation of stochastic differential equation with the exact coupling scheme

1 Introduction

A stochastic differential equation is a mixture of deterministic and probabilistic elements. It is fundamental in many fields of application, and therefore, we find it with great interest in financial and actuarial mathematics, engineering applications, and insurance companies. In recent years, specialists in stochastic differential equations have been trying to find approximate solutions to it, either by using strong approximations or weak approximations. Strong approximation methods have had the most significant use in many scientific papers with a clear difference in applications and conditions. For example, when a stochastic differential equation is one-dimensional, there are often no difficulties in its application. When the Wiener process is multidimensional, the task gets more difficult. The previous articles [1,2,3] explore methods for estimating double integrals in any dimension using Fourier expansion. For the interested reader to learn more about the accomplishment of the simulation of the stochastic differential equation, we refer to [4,5]. Due to its superior properties for many problems, solutions to fractional stochastic differential equations (FSDEs) driven by the Brownian motion have recently received much attention from scientific researchers, see [6,7]. In [8], the authors look at a stochastic viscoelastic wave equation with nonlinear damping and logarithmic nonlinear source terms that established a blow-up result. On the other hand, a strong approximation of an Ito process is not always needed. It is possible that you are only interested in a function of the value of the Ito process at a certain point in time T , one of the first two moments, for example. This method and the use of expectation for some functions give an excellent approximation to the probability distribution of random variables. As a result, the kind of approximation that is needed in this case is significantly weaker than that provided by the pathwise convergence method. In this paper, we will apply the process of weak approximations using the technique of coupling developed by Davie [9].

This paper is structured as follows: Section 2 summarizes the recent stochastic differential equation (SDE) studies and addresses Davie’s methodology [9]. Section 3 discusses the integrated correct relation and the system of Euler. Section 4, offers a numerical example of convergence behavior. In this work, for the exact coupling, we use the multilevel Monte Carlo (MC) technique [10,11] on the two-dimensional SDEs.

2 Stochastic differential equations

Assume the SDE is as follows:

where i = 1 , … , n on [ 0 , T ] , Y ( t ) is a n -dimensional vector, and V ( t ) is a n -dimensional path. Further, the coefficients B i k ( t , Y ( t ) ) satisfy a global Lipschitz condition ∣ O ( t , Y ) − O ( t , y ) ∣ ≤ A ∣ Y − y ∣ , and ∣ B ( t , Y ) − B ( t , y ) ∣ ≤ A ∣ Y − y ∣ , for all t ∈ [ t 0 , T ] and Y , y ∈ R , with A > 0 is a constant. Now assuming A and B are continuously on t for each Y , then there is a unique solution Y ( t ) to equation (2.1). To find an approximate solution on this interval [ 0 , T ] , we divided this interval into positive integer N equal length of intervals, i.e., S = T / N . Adding the following quadratic terms to the Euler scheme yields the Milstein scheme: ∑ k , l = 1 n χ i k l ( j S , y ( j ) ) P k l ( j ) . This will give the following scheme:

where Δ V k ( j ) = V k ( ( j + 1 ) S ) − V k ( j S ) , P k l ( j ) = ∫ j S ( j + 1 ) S { V k ( t ) − V k ( j S ) } d V l ( t ) , and χ i k l ( t , y ) = ∑ m = 1 q B m k ( t , y ) ∂ B i l ∂ y m ( t , y ) . If the commutativity condition

is satisfied for all y ∈ R d , t ∈ [ 0 , T ] , and all i , k , l , After that, the Milstein method is reduced to

It should be noted that the previous strategy is solely dependent on the Brownian motion Δ V k ( j ) being implemented . We may use Δ V k ( j ) only and apply the special equations for the Milstein method. This can be implemented from: the observation that P k l ( j ) + P l k ( j ) = 2 F k l ( j ) , where F k l ( j ) = 1 2 Δ V k ( j ) Δ V l ( j ) for k ≠ l and F k k ( j ) = 1 2 { ( Δ V k ( j ) ) 2 − S } . Also note that scheme (2.4) has the order 1 for n = 1 , but if n > 1 , we obtain the order 1 2 . As delineated in the previous paper [9], we tend to modify to scheme (2.4), which can provide the order 1 with the invertible diffusion. As a result, we will have the exact coupling strategy as follows.

3 Multilevel method for weak order of SDEs

Giles is developing the multilevel MC technique and based on using the MC technique many times with various time increments. So, let us assume that we have the SDE.

where i = 1 , … , q on [ 0 , T ] , x ( t ) is a n -dimensional vector, and W ( t ) is a n -dimensional Brownian path. We will estimate E f ( x ( T ) ) , by the mean of N successive estimates, i.e., 1 N ∑ i = 1 N f ( x S ( i ) ) , where x S ( i ) is a solution of (3.1) and f is globally Lipschitz. The multilevel approach has the benefit of reducing the estimation’s computational load. In the multilevel MC simulations, we examine MC simulations with various time steps S ( r ) = T 2 r , while the standard method depends on a fixed time step. Now, we assume that μ r = E f ( x S ( r ) ) . Then, we will estimate μ n , where n is large enough.

From Giles [10], we have

The multilevel MC approach depends on approximating all terms in the right-hand side of (3.2) independently.

First, for μ 0 , we have

and

Next, for μ k − μ k − 1 , we have

Thus, the right-hand side of (3.2) becomes

It is obvious that the quantity ( f ( x S ( k ) ) − f ( x S ( k − 1 ) ) ) is the outcome of two approximations with successive time steps S ( k ) and S ( k − 1 ) , but the same Brownian path. The variance may be reduced by using the same Brownian path, so that a smaller value for N k can be obtained.

Now, V k represents the variance of a single sample ( f ( x S ( k ) ) − f ( x S ( k − 1 ) ) ) . Y ˆ k simple estimator’s variance is expressed as follows:

As a result, if the variance of the starting level is

then the variance of Y ˆ = ∑ k = 1 n Y ˆ k is expressed as follows:

Furthermore, the computational load being proportional to

Therefore, by choosing N k to be proportional to V k S k , the variance, for a fixed computational load, is minimized.

The following theorem presents the general application of the multilevel MC approach, see [10]. By applying this theorem to the exact coupling, we obtain the computational load results.

The previous theorem may be used to obtain β , which is twice the scheme’s strong order.

We’ll examine the order of the variance in the following applications, i.e., V k = O ( S k β ) , and use Theorem 3.1 to calculate the computational load.

4 The multilevel MC method for the exact coupling method using scheme (2.7)

In [13], Giles and Szpruch have worked on the multilevel approach and obtained O ( ε − 2 ) computational load although they employ a more stringent demand for the function’s regularity f than Lipschitz. So, we showed in [14] that under nondegeneracy, the order of method (2.7) with the exact coupling is O ( S ) , and a single sample’s variance is

As a result of (4.1) and because f is a Lipschitz function, the vector and matrix functions fulfill the Lipschitz conditions with uniformly bounded derivatives, and we have the weak order O ( S ) for (2.7) with the invertibility condition for matrix b i k ( t , X ( t ) ) . So that

Thus, the single sample variance V k is O ( S k 2 ) for the standard estimator (3.5). In addition, the best option for N k is proportional to ( S k 3 / 2 ) asymptotically. As a result, since the variance is β = 2 > 1 , the computational load will be O ( ε − 2 ) , according to Theorem 3.1.

As a consequence, the number of simulations N k will drop rapidly as the variance of the multilevel MC for scheme (2.7) with the exact coupling decreases. So the main order of the convergence of this method will be obtained from the first level of the simulation, which will be the dominant term with the computational load O ( ε − 2 ) .

5 Conclusion

We established a weak convergence method to numerically solve stochastic differential equations with a possibly degenerate diffusion coefficient. We used an exact coupling in a two-dimensional case, which has a condition that the SDE should be invertible.