Explicit order 3/2 Runge-Kutta method for numerical solutions of stochastic differential equations by using Itô-Taylor expansion

Abstract This paper aims to present a new pathwise approximation method, which gives approximate solutions of order 32 $\begin{array}{} \displaystyle \frac{3}{2} \end{array}$ for stochastic differential equations (SDEs) driven by multidimensional Brownian motions. The new method, which assumes the diffusion matrix non-degeneracy, employs the Runge-Kutta method and uses the Itô-Taylor expansion, but the generating of the approximation of the expansion is carried out as a whole rather than individual terms. The new idea we applied in this paper is to replace the iterated stochastic integrals Iα by random variables, so implementing this scheme does not require the computation of the iterated stochastic integrals Iα. Then, using a coupling which can be found by a technique from optimal transport theory would give a good approximation in a mean square. The results of implementing this new scheme by MATLAB confirms the validity of the method.


Introduction
The purpose of this paper is to develop a new pathwise approximation to numerical solutions of SDEs driven by multidimensional Brownian motion. There is a standard approach, which described in [1], approximates solutions of SDEs to the required order using stochastic Taylor expansion at each time step. Nevertheless, applying this method would be di cult when the deriving Brownian motion dimension is greater than to get higher order than , due to hardness of generating the iterated stochastic integrals Iα. To overcome this di culty, J. G. Gaines and T. J. Lyons [2] had introduced an e ective method to deal with double stochastic integrals to get approximate solutions of order . However, unfortunately the method has not been extended to further dimensions higher than -dimensional of Brownian motion. Kloeden and Platen [1] suggest using a Fourier expansion of a Brownian bridge process to approximate stochastic integrals in any dimension. Rydén and Wiktorsson [3] have described a method that covered -dimensional of Brownian motion, uses the fact that the Itô integrals have an in nitely divisible distribution when are conditioned on the Wiener increments. There was a di erent scheme that can simulate the iterated Itô integrals for multiple independent Brownian motions. The method which derives the conditional joint characteristic function of the iterated Ito integrals given the Brownian increments, and then proposed an algorithm for the simulation of the iterated Ito integrals and the Brownian increments introduced by Wiktorsson [4]. The schemes which had introduced in [1,3,4] are useful for any dimension of stochastic integrals, but they need a computational cost.
In this paper, we aim to present a new pathwise approximation scheme that can be used to get a higherorder approximation for Brownian motion in dimension greater than . It is based on using a coupling and a version of the perturbation method. Therefore, this new scheme does not require calculating Iα, because we replace them by random variables with the same moments conditional on the linear term. Then, we get a random vector, which is a good approximation in distribution to the original Taylor expansion, and it can be generated using algorithms in MATLAB. We have used a technique from optimal transport theory [5] to give a good approximation in mean square. Other works use a coupling to approximate SDEs numerically, such that [6][7][8][9]. This paper, organized as follows: section , gives a background material, and section , shows the implementation of the scheme.

Background
Let (Ω, F, P) be a probability space. Let σ(xs : s ≤ t) denote the smallest σ-algebra such thatxs (a stochastic process) is σ(xs : is called a ltration. F t is interpreted as corresponding to the information available at time t (the amount of information increasing as time progresses). A stochastic processx t is adapted to a ltration {F t } ifx t is F t -measurable for all t ≥ .

. Stochastic Itô-Taylor expansion
Let us consider a -dimensional Ito stochastic di erential equation in the form of an integral where a(t, x) and b(t, x) are Borel measurable functions de ned on [ , ∞) × R) with values in R. Hence, for f : R → R which is a twice continuously di erentiable function, the Itô formula [1] gives where the operators are Clearly, if f (x) ≡ x we have L f = a and L f = b , in which case that there is a reduction in (2.2) to the original Itô equation forx t , that is tox Applying the Itô formula (2.2) to the function f = a and f = b in (2.5), we obtain where the remainder is For the Taylor expansion of d-dimensional Itô stochastic di erential equations, a similar expansion holds as above forx The general Itô Taylor expansion is given in (2.12).

. Perturbation method
Consider the following: suppose we wish to simulate U = X + ϵY, where X and Y are independent, X has a smooth density, and ϵ is small. Also, suppose that X is easy to generate, while Y is hard to generate. Then, generating U by generating X and Y will be hard. Alternatively, let us suppose there is another random variable that is easy to generate to overcome this dilemma, say Z, which is independent of X. The random variable Z has the same moments up to order m − as Y (i.e. E(Z k ) = E(Y k ) for k = , , , ..., m − ). Then, we can prove that V = X + ϵZ is an approximation to U with the error of order ϵ m . The justi cation for this can be seen by writing f X for the density of X etc, so we have . As Z has the same moments, we get the same expression for . It can expect from this a Wasserstein distance estimate of the same order W(U, V) = O(ϵ m ), which gives a coupling with E(U − V) = O(ϵ m ). This way is an example for the bound of densities can give a Wasserstein bound between U and V, but it is not in an immediate way, the Wasserstein bound implies the coupling between U and V [10].

. Construction of pathwise approximation of SDEs using stochastic Taylor expansions
Consider an Itô SDE on an interval [ , T], for a q-dimensional vectorx t , with a d-dimensional driving Brownian path W t . Assume (2.9) satis es an existence and uniqueness theorem. The basic idea to get pathwise approximation is to divide [ , T] into a nite number N of subintervals, where the length of the step equals h = T N . We have used a stochastic Taylor expansion in order to approximate the equation on each subinterval.
The basic scheme which we have is Euler-Maruyamâ In order to get the Milstein scheme, we add the quadratic termŝ [10]. Our constructed method, which approximates SDEs at order , assumes nondegeneracy of b in (2.9) and uses Itô-Taylor expansion and Runge-Kutta method. It is needed to use the following notations, which had mentioned in this reference [1] to construct the scheme for approximate solutions of SDEs at a higher order. The following discussion is derived from [10]. Let M be the set of all multi-indices, α = (j , ..., j l ) of length where n(α) is the number of zero indices in α.
The case l(α) = n(α) = m+ is needed (for example when α = ( )) because Iα has non-zero mean, so we need a higher order for the local error. The order γ = m Taylor approximation to (2.9) is given bŷ where v is the multi-index of zero length, and f α,i denotes the ith component of fα.
The main point which will be discussed in this paper; how do we generate RHS of (2.12) easily? Actually, this research has developed a scheme where to generate RHS of (2.12) by using perturbation and coupling method; how?
Let us apply this method in (2.12) in order to approximate SDEs at higher order. In order to simplify application, we will be deal with the rst iteration from to h of (2.12) As Iα is not independent of ∆W, so it would be di cult to generate these stochastic integrals when we want to nd approximate solutions. By following the construction in [10], we can overcome this issue using the relation For α = (j , ..., j l ) we can replace the iterated stochastic integrals Iα using random variables with the same moments conditional on the linear term as follows where we calculate the sum over all β = (i , ..., i l ) such that for each k ∈ { , ..., l} we have either i k = j k or i k = < j k .
(For examples I =h (K + K V + K V + K V V )). The random variables K β and V k ∼ N( , ) are generated independently.
By setting ϵ = h , and de ning Mm = α ∈ M : ≤ l(α) ≤ m , we can write this where V k is independent N( , ) of a random variable Kα and Q i is a polynomial, each monomial of which has overall order at least 2 in ϵ.
We have found it hard to generate Y directly by using (2.14). To tackle this problem, we use a version of the perturbation method by assuming random variables L β and using them rather than K β which easy to generate for α ∈ M such that ≤ l(α) ≤ m.
It will be discussed the way of construction L β later in this section with details. Then, (2.14) can be modi ed as follows where we calculate the sum over all β = (i , ..., i l ) such that for each k ∈ { , ..., l} we have either i k = j k or i k = < j k . (For exampleĪ =h (L + L V + L V + L V V )). The random variables L β and V k ∼ N( , ) are generated independently.
Then, we can be rede ne the formula (2.15) as follows whereV k are independent N( , ) and independent of the random variable Lα.
The following theorem is essential to give the error bound in the approximation by the modi ed Y by considering equation (2.9) and using the notation of section (2.12). [10]. Proof. We refer the reader to [11].

Theorem 1. This theorem is L p version of
How do we apply theorem 1 to generate approximate solution for (2.12)? Basically, we assume coe cients a k and b ik are su ciently regular (uniform bounds for the coe cients and their derivatives up to order two will certainly su ce) and the matrix (b ik ) has rank q everywhere with uniformly bounded right inverse. Then, the iterated integrals I α,jh,(j+ )h in (2.12) are hard to generate as it has been mentioned so we should use modi cation of iterated integrals which are easy to generateĪ α,jh,(j+ )h [10]. The modi ed formula which results is the followinĝ We now return to think about the application of coupling to the simulations of SDEs in the rst step between Y andȲ. We have assumed the L β satisfy the hypothesis of the theorem 1 such that E(Kα ...Kα r ) = E(Lα ...Lα r ) whenever α , ..., αr ∈ Mm satisfy r k= (l(α k )− ) ≤ m − . Therefore, Wasserstein bound between Y andȲ implies the existence of a suitable coupling between Y andȲ for which the following relation (2.19) holds. The coupling which has been deduced between Y andȲ can be extended to be between the random variables (V k , Kα) and (V k ,Lα), which leads to deduce the following relation, . (2.19) For the bound at j step, we can use the same argument as used in (2.19) to deduce the bound for (2.21), except that the coe cient functions evaluated at random variablesx (j) i such as ((V (j) , Y (j) andȲ (j) ), so the extended coupling between (V (j) k , K (j) α ) and (V (j) k ,L (j) α )) are conditional on the σ-algebra F j , where F j = σ{V i , K i , L i .i < j}. Hence, if we take an expectation w.r.t σ-algebra and conditional on (F j ), then we get rst a conditional bound as follows . (2.20) Therefore, by nested expectation, we have It is required to know the calculations of the relevant expectation of products of K β for generating a suitable random variables L β . However, before doing that, we would clarify how to construct proper random variables. First, the new random variables L β must satisfy the moment conditions in the statement of theorem 1. Thus, once the relevant moments of the K β are calculated, then it would be taken any choice of L β which satis es these moment conditions of the theorem 1. To calculate the moments, we have used the following three lemmas. We refer the reader to [12] for more details of how we compute the non-deterministic moments.

Implementation of the method
In the preceding section, it has introduced the method that can be used to get approximate solutions of higherorder for any dimension of Brownian motions. Therefore, to get numerical solutions for SDEs of order , the formula (2.12) is used by putting m = , then using theorem 1 to determine the moments that we need forĪα and then using lemmas 1, 2 and 3 to compute the moments. For β in the de nition ofĪα in (2.16), we need all indices of length 2 and 3. The random variables L β in (2.16) must satisfy the moment condition in the statement of theorem 1, that E(Kα ...Kα r ) = E(Lα ...Lα r ) whenever α , ..., αr ∈ Mm satisfy r k= (l(α k ) − ) ≤ m − . Hence, in order to get numerical solutions of order , L β must satisfy the moments E(Lα) for all α of length , α of length with non zero-indices and the expectation of products E(Lα L β ) with α, β have length . As a result of that, we have the following moments K β which the random variables L β must satisfy Deterministic K α

. Construction of L β
In this section, we discuss the choices that we make for Lα. Firstly, for all the deterministic α, the random variables Lα have the same values as E(Kα). For α of length , we only need the expectations, so we can choose them to be deterministic by making them equal to E(Kα). For the cases of the products, we can use the fact of the antisymmetry K lk = −K kl for k < l, so we choose the random variables L lk with k < l which satisfy it. Therefore, the expectations for the products (Lα L β ) (α and β are distinct and not deterministic) are zeros such as E(L lk L l k ) = for k < l, k < l and k, l ≠ k , l .
Hence, the random variables L kl with k < l are required to have all moments nite, mean zero, variance and uncorrelated.Hence, we can take any random variable which satis es the assumption of nite moments and has the mean zero, variance , and generate independent copies. We generate L kl with k < l using normal distribution with mean zero, variance and uncorrelated, but this is not the only way of generating L β .
As a result , we set L lk = −L kl for ≤ k < , and by using lemma 1,2 and 3, they give L = and L kk = for k > . By setting L = , L kk = L k k = L kk = − for k > , and all other L β of length to . As a result of these and by using (2.16) we have the followinḡ for k, l, n > .
gives approximate solutions of order [1]. where . Simulation result The following system of SDEs had simulated by using (3.4) for the number of steps N = . In Table 1, we calculate the error of the method at a time interval [ , T], where T = by using the iteration numbers N = , and , step sizes h = / , / and / and the simulation number R = .   Figure 2 shows the con dence intervals along the line of the least square for the expectations against step sizes. The blue line is the least square. The level of the con dence interval is %. The green line is the upper con dence bounds, while the red line is the lower con dence bounds. The red line has just two points because the lower bound of one of the con dence intervals is negative. Therefore, we swap the negative value in the con dence interval with the zero, so when we take the log for this con dence interval, then we get −∞ as the lower bound for this con dence interval, but it does not show in the graph. Therefore, we have used this method to handle the negative value in the con dence interval because if we do not do this, then the least square regression line does not pass through con dence intervals as it should be. This scheme aims to give approximation at order , where the slope of the line is . .

Conclusion
The paper presents a numerical method that gives approximate solutions of stochastic di erential equations with a strong error rate of O(h ). It is not the only method that can give this result. As mentioned in the introduction, a quiet few methods are existing, which can give approximate solutions for SDEs of order , but they have a computational cost such as the one which was developed by Kloeden and Platen [1] in page using Fourier series expansion. Therefore, our method does not involve a computational cost due to not requiring simulation of iterated stochastic integrals Iα. We can extend the method by adding additional terms from the formula (2.12) to achieve a higher-order of numerical solutions.