Rooted Tree Analysis for Order Conditions of Stochastic Runge-Kutta Methods for the Weak Approximation of Stochastic Differential Equations

# Rooted Tree Analysis for Order Conditions of Stochastic Runge-Kutta Methods for the Weak Approximation of Stochastic Differential Equations

###### Abstract

A general class of stochastic Runge-Kutta methods for the weak approximation of Itô and Stratonovich stochastic differential equations with a multi-dimensional Wiener process is introduced. Colored rooted trees are used to derive an expansion of the solution process and of the approximation process calculated with the stochastic Runge-Kutta method. A theorem on general order conditions for the coefficients and the random variables of the stochastic Runge-Kutta method is proved by rooted tree analysis. This theorem can be applied for the derivation of stochastic Runge-Kutta methods converging with an arbitrarily high order.

s

tochastic Runge-Kutta method \sepstochastic differential equation \sepweak approximation \seprooted tree analysis \seporder condition
MSC 2000: 65C30 \sep60H35 \sep65L05 \sep60H10 \sep34F05

## 1 Introduction

In recent years many numerical methods have been proposed for the approximation of stochastic differential equations (SDEs), see e.g. [7], [9], [10], [11], [14], [19] and [20]. Mainly, numerical methods for strong and for weak approximations can be distinguished. While strong approximations focus on a good approximation of the path of a solution, weak approximations are applied if a good distributional approximation is needed. In Section 2 of the present paper, a class of stochastic Runge-Kutta (SRK) methods for the weak approximation of Itô and Stratonovich SDEs is introduced. As in the deterministic setting, order conditions for SRK methods are calculated by comparing the numerical solution with the exact solution over one step assuming exact initial values. Therefore, the actual solution of the SDE and the numerical approximation process have to be expanded by a stochastic Taylor series. However, even for low orders such expansions become much more complex than in the deterministic setting where it is already a lengthy task. In order to handle this task in an easy way, a rooted tree theory based on three different kinds of colored nodes is established in Section 3, which is a generalization of the rooted tree theory due to Butcher [3]. Thus, colored trees are applied in Section 4 and 5 to give a representation of the solution and the approximation process calculated with the SRK method in order to allow a rooted tree analysis of order conditions. A similar approach with two different kinds of nodes has been introduced by Burrage & Burrage [1], [2] for a SRK method converging in the strong sense as well as in Komori et al. [8] for ROW-type schemes for Stratonovich SDEs. Finally, the main Theorem 6.4 presented in Section 6 immediately yields all order conditions for the coefficients and the random variables of the introduced SRK method such that it converges with an arbitrarily given order in the weak sense. As a result of this theorem, the lengthy calculation and comparison of Taylor expansions can be avoided.

Let be a probability space with a filtration and let for some . We consider the solution of either a -dimensional Itô stochastic differential equation system

 dXt=a(t,Xt)dt+b(t,Xt)dWt (1)

or a -dimensional Stratonovich stochastic differential equation system

 dXt=a(t,Xt)dt+b(t,Xt)∘dWt. (2)

Let be the -measurable initial condition such that for some holds where denotes the Euclidean norm if not stated otherwise. Here, is an -dimensional Wiener process w.r.t. . SDE (1) and (2) can be written in integral form

 Xt=x0+∫tt0a(s,Xs)ds+m∑j=1∫tt0bj(s,Xs)∗dWjs (3)

for , where the th column of the -matrix function is denoted by for . Here, the second integral w.r.t. the Wiener process has to be interpreted either as an Itô integral in case of SDE (1) or as a Stratonovich integral in case of SDE (2), which is indicated by the asterisk.

The solution of a Stratonovich SDE with drift and diffusion is also a solution of an Itô SDE as in (1) and therefore also a diffusion process, however with the modified drift

 ~ai(t,x)=ai(t,x)+12d∑j=1m∑k=1bj,k(t,x)∂bi,k∂xj(t,x) (4)

for and provided that is sufficiently differentiable, i.e.

 Xt=Xt0+∫tt0a(s,Xs)ds+∫tt0b(s,Xs)∘dWs=Xt0+∫tt0~a(s,Xs)ds+∫tt0b(s,Xs)dWs. (5)

The solution of the stochastic differential equation (3) is sometimes denoted by in order to emphasize the initial condition. We suppose that the drift and the diffusion are measurable functions satisfying a linear growth and a Lipschitz condition

 ∥a(t,x)∥+∥b(t,x)∥≤C(1+∥x∥) (6) ∥a(t,x)−a(t,y)∥+∥b(t,x)−b(t,y)∥≤C∥x−y∥ (7)

for all and all with some constant . Then the conditions of the Existence and Uniqueness Theorem are fulfilled for the Itô SDE (1) (see, e.g., [6]). If the conditions also hold with replaced by the modified drift in the Itô SDE, then the Existence and Uniqueness Theorem also applies to the Stratonovich SDE (2).

In the following, let denote the space of times continuously differentiable functions for which all partial derivatives up to order have polynomial growth. That is, for which there exist constants and depending on , such that holds for all and any partial derivative of order .

Let be a discretization of the time interval such that

 0≤t0

and define for with the maximum step size

 h=max0≤n≤N−1hn.

In the following, we consider a class of approximation processes of the type where is a random variable or in general a vector of random variables, with moments of sufficiently high order, and is a vector valued function of dimension . We write and we construct the sequence

 Y0=Xt0Yn+1=A(tn,Yn,hn;ξn),n=0,1,…,N−1, (9)

where is independent of , while for is independent of and . Then we can define weak convergence with some order of an approximation process.

###### Definition 1.1

A time discrete approximation process converges weakly with order  to the solution process of SDE (1) or SDE (2) as if for each there exists a constant , which does not depend on , and a finite such that

 |E(f(Xt))−E(f(Y(t)))|≤Cfhp (10)

holds for each and .

Since we are interested in calculating a global approximation converging in the weak sense with some desired order , we make use of the following theorem due to Milstein (1986) [13] which is stated with an appropriate notation.

###### Theorem 1.2

Let be the solution of SDE (1) or of SDE (2). Suppose the following conditions hold:

1. the coefficients in the case of SDE (1), in the case of SDE (2) and are continuous, satisfy a Lipschitz condition (7) and belong to with respect to for , ,

2. for sufficiently large (specified below) the moments do exist and are uniformly bounded with respect to and ,

3. assume that for all the following local error estimation

 |E(f(Xt,x(t+h)))−E(f(Yt,x(t+h)))|≤K(x)hp+1 (11)

is valid for , any with and .

Then for all and all the following global error estimation

 |E(f(Xt0,Xt0(tn)))−E(f(Yt0,Xt0(tn)))|≤Chp (12)

holds for all , where is a constant and where is the maximum step size of the corresponding discretization , i.e. the method (9) has order of accuracy in the sense of weak approximation.

A proof of Theorem 1.2 can be found in [11], [12], [13] and [17]. Lemma 1.3 gives sufficient conditions such that condition (ii) of Theorem 1.2 (see also [11], [12], [13]) holds.

###### Lemma 1.3

Suppose that for given by (9) and the conditions

 ∥E(A(tn,x,h;ξn)−x)∥≤C1(1+∥x∥)h, (13) ∥A(tn,x,h;ξn)−x∥≤M(ξn)(1+∥x∥)h1/2 (14)

hold where has moments of all orders, i.e. , , with constants and independent of . Then for every even number the expectations exist and are uniformly bounded with respect to and , if only exists.

## 2 A Class of Stochastic Runge-Kutta Methods

In the following a class of stochastic Runge-Kutta methods is introduced for the approximation of both Itô and Stratonovich stochastic differential equation systems w.r.t. an -dimensional Wiener process. In order to preserve the most possible generality, the considered class of stochastic Runge-Kutta methods is of type (9) and has the following structure

 Y0=x0Yn+1=A(tn,Yn,hn;θν(hn):ν∈M) (15)

where is an arbitrary finite set of multi-indices with elements and , , are some suitable random variables. For the weak approximation of the solution of the -dimensional SDE system (3), considered either with respect to Itô or Stratonovich calculus, the general class of -stage stochastic Runge-Kutta methods is given by

 Y0 = x0 Yn+1 = Yn +s∑i=1z(0,0)ia(tn+c(0,0)ihn,H(0,0)i) (16) +s∑i=1m∑k=1∑ν∈Mz(k,ν)ibk(tn+c(k,ν)ihn,H(k,ν)i)

for with

 H(0,0)i = Yn +s∑j=1Z(0,0),(0,0)ija(tn+c(0,0)jhn,H(0,0)j) +s∑j=1m∑r=1∑μ∈MZ(0,0),(r,μ)ijbr(tn+c(r,μ)jhn,H(r,μ)j) H(k,ν)i = Yn +s∑j=1Z(k,ν),(0,0)ija(tn+c(0,0)jhn,H(0,0)j) +s∑j=1m∑r=1∑μ∈MZ(k,ν),(r,μ)ijbr(tn+c(r,μ)jhn,H(r,μ)j)

for , and , where

 z(0,0)i =αihn z(k,ν)i =∑ι∈Mγ(ι)i(k,ν)θι(hn) Z(0,0),(0,0)ij =A(0,0),(0,0)ijhn Z(0,0),(r,μ)ij =∑ι∈MB(ι)ij(0,0),(r,μ)θι(hn) Z(k,ν),(0,0)ij =A(k,ν),(0,0)ijhn Z(k,ν),(r,μ)ij =∑ι∈MB(ι)ij(k,ν),(r,μ)θι(hn)

for . Here are the coefficients of the SRK method and as usual the weights can be defined by

 c(0,0)=A(0,0),(0,0)e,c(k,ν)=A(k,ν),(0,0)e, (17)

with . If for and , , then (2) is called an explicit SRK method, otherwise it is called implicit. The class of SRK methods introduced above can be characterized by an extended Butcher array

(18)

for and for . We assume that the random variables satisfy the moment condition

 E(θp1ν1(hn)⋅…⋅θpκνκ(hn))=O(h(p1+…+pκ)/2n) (19)

for all and , . The moment condition ensures a contribution of each random variable having an order of magnitude . This condition is in accordance with the order of magnitude of the increments of the Wiener process. Further, the moment condition is necessary for the estimates of the reminder terms of the Taylor expansion of the SRK approximation presented in Section 6.

Some SRK schemes which belong to the introduced general class of SRK methods can be found in [15], [16] and [17]. Further, many Runge-Kutta type schemes proposed in recent literature like in [7], [8], [10] or [21] are covered. Usually, the set may consist of some multi-indices with for and the random variables may be chosen as multiple Itô or Stratonovich integrals of type or , depending on the calculus that is used.

For example, the SRK scheme RI1WM due to Rößler [17] for the Itô SDE (1) in the case of and with is defined by (2) with

 z(0,0)i =αi⋅hn z(k,l)i =γ(k)i(k,l)^I(k)+γ(k,l)i(k,l)^I(k,l)√hn Z(0,0),(0,0)ij =A(0,0),(0,0)ij⋅hn Z(0,0),(r,s)ij =B(r)ij(0,0),(r,s)^I(r) Z(k,l),(0,0)ij =A(k,l),(0,0)ij⋅hn Z(k,l),(r,s)ij =B(0)ij(k,l),(r,s)√hn

for . Further, we define in the case of , in the case of or and in the case of for . Here, , , are independent random variables defined by and . The , , are defined by with independent random variables such that for , and for and . Thus, we can characterize the SRK method (2) by the following Butcher array for with :

The coefficients of the order SRK scheme RI1WM are given in Table 1. For detailed calculations of the order conditions and the corresponding coefficients we refer to [17].

As an example for a SRK scheme due to Rößler applicable to the Stratonovich SDE (2) with and fulfilling a commutativity condition (see [16], [17] for details) we choose now and

 z(0,0)i =αi⋅hn z(k,k)i =γ(k)i(k,k)^I(k)(hn) Z(0,0),(0,0)ij =A(0,0),(0,0)ij⋅hn Z(0,0),(k,k)ij =B(k)ij(0,0),(k,k)^I(k)(hn) Z(k,k),(0,0)ij =A(k,k),(0,0)ij⋅hn Z(k,k),(l,l)ij =B(l)ij(k,k),(l,l)^I(l)(hn)

for and . The coefficients of such a method can be represented by the Butcher array taking for the form

For detailed calculations of the order conditions we refer to [16] and [17]. The coefficients of the order 2.0 SRK scheme RS1WM are presented in Table 2.

## 3 Stochastic Rooted Tree Theory

The SDE system (3) can be represented by an autonomous SDE system

 Xt=x0+∫tt0a(Xs)ds+m∑j=1∫tt0bj(Xs)∗dWjs (20)

with one additional equation representing time. Hence, it is sufficient to treat autonomous SDE systems in the following. First of all, we give a definition of colored graphs which will be suitable in the rooted tree theory for SDEs w.r.t. a multi-dimensional Wiener process (see [18]).

###### Definition 3.1

Let be a positive integer.

1. A monotonically labelled S-tree (stochastic tree) t with nodes is a pair of maps with

so that for . Unless otherwise noted, we choose the set where is a variable index with .

2. denotes the set of all monotonically labelled S-trees w.r.t. . Here two trees and just differing by their colors and are considered to be identical if there exists a bijective map with and so that holds for .

So defines a father son relation between the nodes, i.e.  is the father of the son . Furthermore the color , which consists of one element of the set , is added to the node for . Here,  \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0.1cm, radius=1.6mm, treefit=loose] \Tn \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0cm, radius=1.6mm, treefit=loose] \TC* [tnpos=r] is a deterministic node and  \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0.1cm, radius=1.6mm, treefit=loose] \Tn \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0cm, radius=1.6mm, treefit=loose] \TC [tnpos=r] is a stochastic node with a variable index . In the case of the node of type  \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0.1cm, radius=1.6mm, treefit=loose] \Tn \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0cm, radius=1.6mm, treefit=loose] \Tdot [tnpos=r] is denoted as the root and always sketched as the lowest node of the graph. However, in the case of , the nodes and may also serve as the root of the tree. The variable index is associated with the th component of the corresponding -dimensional Wiener process of the considered SDE. In case of a one-dimensional Wiener process one can omit the variable indices since we have for all (see also [17]). As an example Figure 1 presents two elements of .

For the labelled S-tree in Figure 1 we have and . The color of the nodes is given by , , and .

###### Definition 3.2

Let . We denote by the number of deterministic nodes, by the number of stochastic nodes and by the number of pairs of stochastic nodes with the same variable index. The order of the tree t is defined as with .

The order of the trees and presented in Figure 1 can be calculated as . Every labelled S-tree can be written as a combination of three different brackets defined as follows.

###### Definition 3.3

If are colored trees then we denote by

 (t1,…,tk),[t1,…,tk] and {t1,…,%tk}j

the tree in which are each joined by a single branch to \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0.1cm, radius=1.6mm, treefit=loose] \Tn \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0cm, radius=1.6mm, treefit=loose] \Tdot , \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0.1cm, radius=1.6mm, treefit=loose] \Tn \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0cm, radius=1.6mm, treefit=loose] \TC* and \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0.1cm, radius=1.6mm, treefit=loose] \Tn \pstree[treemode=U, dotstyle=otimes, dotsize=3.2mm, levelsep=0cm, radius=1.6mm, treefit=loose] \TC [tnpos=r] , respectively (see Figure 2).

Therefore proceeding recursively, for the two examples and in Figure 1 we obtain and .

Due to the fact that we are interested in calculating weak approximations, it will turn out that we can concentrate our considerations to one representative tree of each equivalence class.

###### Definition 3.4

Let and be elements of . Then the trees t and u are equivalent, i.e. , if the following hold:

1. There exist two bijective maps

 ψ:{1,…,l(t)}→{1,…,l(t)} with ψ(1)=1,π:A→A % with π(γ)=γ and π(τ)=τ,

so that the following diagram commutes

The set of all equivalence classes under the relation is denoted by . We denote by the cardinality of t, i.e. the number of possibilities of monotonically labelling the nodes of t with numbers .

Thus, a monotonically labelled S-tree u is equivalent to t, if each label is replaced by and if each stochastic node with variable index is replaced by an other stochastic node . Thus, all trees in Figure 3 belong to the same equivalence class as in the example above, since the indices and are just renamed either by and or and , respectively. Finally the graphs differ only in the labelling of their number indices.

For every rooted tree , there exists a corresponding elementary differential which is a direct generalization of the differential in the deterministic case (see, e.g., [3]). For , the elementary differential is defined recursively by

 F(γ)(x)=f(x),F(τ)(x)=a(x),F(σj)(x)=bj(x),

for single nodes and by

 F(t)(x)=⎧⎪ ⎪⎨⎪ ⎪⎩f(