## Abstract

We prove a stochastic Taylor expansion for stochastic partial differential equations (SPDEs) and apply this result to obtain cubature methods, i.e. high-order weak approximation schemes for SPDEs, in the spirit of Lyons and Victoir (Lyons & Victoir 2004 *Proc. R. Soc. A* **460**, 169–198). We can prove a high-order weak convergence for well-defined classes of test functions if the process starts at sufficiently regular initial values. We can also derive analogous results in the presence of Lévy processes of finite type; here the results seem to be new, even in finite dimension. Several numerical examples are added.

## 1. Introduction

Let *H* be a real separable Hilbert space. We consider the stochastic (partial) differential equation for a diffusion process *X* with values in *H*(1.1)or, in the presence of jumps, the stochastic differential equation for a jump-diffusion process *X* with values in *H*(1.2)where, in general, denotes an unbounded linear operator; *α*, *β*_{1}, …, *β*_{d}, denote *C*^{∞}-bounded vector fields, i.e. the vector fields are smooth and all the derivatives (of degree ≥1) are bounded; denotes a finite-dimensional Brownian motion; and is a compound Poisson process with jump rate *μ*_{j} for *j*=1, …, *e*. The initial value *x*∈*H* appears as superscript in the notation of the solution process, i.e. . We assume that *A* is the generator of a *C*_{0} semigroup denoted by . The main reference for equations of type (1.1) is the monograph of Da Prato & Zabczyk (1992). In the (general) Lévy case, we refer to Filipović & Tappe (2006) and Peszat & Zabczyk (2007), even though we do not need results of their strength in our case since all Lévy processes under consideration are of finite type.

Owing to the unboundedness of the operator *A*, several concepts of solutions to (1.2) arise. We refer for the precise definitions to the excellent monograph by Peszat & Zabczyk (2007). The most direct analogue to the finite-dimensional setting is the concept of a *strong solution*, which is defined by(1.3)i.e. by the integrated version of equation (1.1). Besides the obvious integrability assumptions, we need in particular that for all *t*≥0 almost surely with respect to *P*. More relevant is the concept of a *mild solution*, which is related to (1.2) via variation of constants. Indeed, a mild solution is a process satisfying(1.4)given the obvious integrability assumptions. By Itô's formula, any strong solution is a mild solution, but not the other way round. Note that mild solutions need not be semi-martingales, because there is no semi-martingale decomposition if the process evolves outside of (*A*). Consequently, a Stratonovich formulation of (1.2) is, in general, not possible for mild solutions.

Of course, neither strong nor mild solutions can usually be given explicitly, which makes numerical approximation necessary. We are interested in *weak approximation* of the solution in the sense that we want to approximate the valuefor a suitable class of test functions at initial values *x*∈*H*. It is well known that the function solves the Kolmogorov equation in the weak sense, see for instance Da Prato & Zabczyk (1992) in the diffusion case.

Let us assume for a moment that there are no jump components: usually, infinite-dimensional SDEs are numerically solved by finite-element or finite-difference schemes (see, for instance, Gyöngy 1998, 1999; Hausenblas 2003; Yan 2005). This means that the original equation is projected onto some finite-dimensional subspace *H*_{h}⊂*H* and *A* is approximated by some operator *A*_{h} defined thereon. This procedure, which corresponds to a space discretization of the stochastic PDE, is followed by a conversion of the stochastic differential equation on *H*_{h} to a stochastic difference equation by discretizing in time, using the Euler method or a related scheme. Finally, the stochastic difference equation is solved by Monte Carlo simulation, which may be interpreted as a discretization on the Wiener space. For general information about approximation of finite-dimensional SDEs, see Kloeden & Platen (1992).

We want to tackle the problem in the reverse order: we want to do the discretization on the Wiener space *Ω* first, reducing the problem to a deterministic problem; that is, one replaces the *d*-dimensional Brownian motion with finitely many trajectories of bounded variation chosen with well-defined probabilities such that certain moments (of iterated Stratonovich integrals) match. The resulting deterministic problem can be solved by standard methods for the numerical treatment of deterministic PDEs, e.g. by standard finite-element or finite-difference methods. The benefit is that, once the discretization on the Wiener space has been done, we can use the well-established theory of the corresponding deterministic PDE problems, without any complications from stochasticity. Our method of choice for discretization on *Ω* is ‘*cubature on Wiener space*’, developed by Lyons & Victoir (2004) and Kusuoka (2001, 2004); see also Teichmann (2006). In the spirit of these methods, we shall obtain weak approximation schemes of any prescribed order of convergence. Note here that we discretize in the presence of the unbounded operator *A* in the drift vector field. Certainly our assumption 2.1 seems very restrictive, but this assumption is the appropriate analogues of the assumptions in finite dimension that the vector fields are bounded, *C*^{∞} bounded (see Kusuoka 2004; Lyons & Victoir 2004).

Before going into detail, let us motivate the use of cubature formulae in the present context. Let denote the solution of (1.1), formally rewritten in Stratonovich form, if each ‘’ is replaced by ‘’, i.e.(1.5)for a curve function of bounded variation. Roughly speaking the idea of cubature on the Wiener space is to construct short-time asymptotics (for some given degree of accuracy *m*≥2)with some positive, time-independent weights *λ*_{1}, …, *λ*_{N} satisfying *λ*_{1}+⋯+*λ*_{N}=1 and some *d*-dimensional paths *ω*_{1}, …, *ω*_{N} of bounded variation. Of course, the weights and paths are chosen in a specific way, which will be described later in more detail, and the asymptotics will hold only for some class of test functions *f*. Note in particular that the cubature paths *ω*_{1}, …, *ω*_{N} depend on the interval [0,*t*]—they become rougher as *t* approaches 0. The aforementioned procedure replaces the SDE by *N* deterministic well-defined PDEs that have unique mild solutions. The iteration of the short-time asymptotics due to the Markov property then yields a weak high-order approximation scheme.

Here also the main advantage of cubature methods in contrast to Taylor methods becomes visible. The time discretization in the realm of cubature methods always leads to reasonable expressions, namely to reasonable partial differential equations of type (1.5). If we wanted to apply the usual discretization methods in time like the Euler–Maruyama method, we might run into problems. Indeed, the naive Euler scheme is well suited for the differential formulation (1.1) of the problemfor *n*≥1, but it might immediately lead to some . Even in the case of a well-defined strong solution, there is no reason why the discrete approximation should always stay in (*A*). Hence, the naive implementation of the Euler–Maruyama method does not work.

Only formulation (1.4) seems to be suitable for using a Euler-like method. If one understands the semigroup *S* well, one can approximate by expressing the integrals in (1.4) as Riemannian sums, involving evaluations of , which yields a sort of strong Euler method (see, for instance, the very interesting book Prévôt & Röckner (2007) for strong convergence theorems in this direction).

We do not discretize the integral in formulation (1.4), but (weakly) approximate the Brownian motion's paths by paths of bounded variation such that we obtain a weak approximation of .

The presence of jumps does not lead to more complicated expressions since the short-time asymptotics of a jump diffusion can be easily derived from a diffusion's short-time asymptotics by conditioning on the jumps. The arising picture is the following. Discretizing equation (1.2) means to allow a certain number of jumps between consecutive points in the time grid. Between two jumps, we apply a diffusion cubature formula to express the short-time asymptotics. This yields as a corollary of the diffusion theory also the jump-diffusion theory.

A direct application of the cubature on the Wiener space technique to jump diffusions driven by Lévy processes of infinite activity is not possible. However, note that any Lévy process can be approximated by processes with finite activity in a weak sense (see, for instance, Cont & Tankov 2004; Peszat & Zabczyk 2007), and therefore the solutions of the corresponding Hilbert space-valued SDEs converge weakly. In this sense, cubature methods can be used for jump diffusions driven by Lévy processes of infinite activity, too.

The article is organized as follows. In §2, we describe the analytic setting for a stochastic Taylor expansion to work. This is a delicate question since we deal with one unbounded vector field. In §§3 and 4, we work out the cubature method from scratch and prove the relevant convergence results in the diffusion case. In §5, we allow for jumps and prove the associated short-time asymptotics that is relevant to set up a weak approximation scheme. In §6, we apply our method to several examples to demonstrate the results.

## 2. Setting and assumptions

Let be a filtered probability space with a filtration satisfying the usual conditions. Let be a *d*-dimensional Brownian motion and , *j*=1, …, *e* be *e* independent compound Poisson processes given bywhere denotes a Poisson process with jump rate *μ*_{j}>0 and is an i.i.d. sequence of random variables with distribution *ν*_{j} for *j*=1, …, *e*, such that each *ν*_{j} admits all moments. All notions are with respect to the given filtration. We assume furthermore that all sources of randomness are mutually independent.

Let *H* be a separable Hilbert space. We furthermore fix a strongly continuous semigroup *S* on *H* with generator *A*. Let *α*, *β*_{1}, …, *β*_{d}, the diffusion-vector fields, and *δ*_{1}, …, *δ*_{e}, the jump-vector fields, be *C*^{∞} bounded on *H*; that is, the vector fields are infinitely often differentiable with bounded partial derivatives of all proper orders *n*≥1. We consider the mild càdlàg solution of a stochastic differential equation(2.1)(2.2)See Filipović & Tappe (2006) for all necessary details on existence and uniqueness for the previous equation. Note furthermore the decomposition theorem in Forster *et al*. (submitted), which states that we do not need any further existence and uniqueness results in this case: in particular, we do not need to impose further (contractivity) conditions on *A* as in Filipović & Tappe (2006) in our finite activity case.

The previous conditions are slightly more than standard for existence and uniqueness of mild solutions, i.e. in Filipović & Tappe (2006) the authors need Lipschitz conditions on the vector fields, whereas we assume them to be *C*^{∞} bounded. In order to formulate a stochastic Taylor expansion, we shall need one main assumption that we formulate in the sequel. This assumption has already been successfully applied in several circumstances (e.g. Filipović & Teichmann (2003) and Baudoin & Teichmann (2005)).

We apply the following notations for the Hilbert space (*A*^{k}):which we need in order to specify the main analytic condition for our considerations:

*We assume that* *α, β*_{1}, …, *β*_{d}*, the diffusion-vector fields, and* *δ*_{1}, …, *δ*_{e}*, the jump-vector fields, map* *and are* *C*^{∞} *bounded thereon for each* *k*≥0*; that is, the vector fields are infinitely often differentiable with bounded partial derivatives of all orders* *n*≥1 *on the Hilbert space* (*A*^{k}) *for each* *k*≥0.

These assumptions are the appropriate analogues of the assumptions in finite dimension that the vector fields are bounded, *C*^{∞} bounded (see Kusuoka 2004; Lyons & Victoir 2004). In order to establish a true convergence rate, one needs an additional cut-off argument, which is outlined in remark 4.9. This can, like in the finite-dimensional setting, certainly be improved.

In order to show examples of vector fields, which are *C*^{∞} bounded on (*A*^{k}) consider the following structure. Let *H* be a separable Hilbert space and *A* the generator of a strongly continuous semigroup. We know that (*A*^{∞}) is a Fréchet space and an injective limit of the Hilbert space (*A*^{k}) for *k*≥0. Following the analysis (as developed in Filipović & Teichmann (2003); see also Hamilton 1982; Kriegl & Michor 1997) in which the analytic concepts were originally developed, we can consider a vector field . If *V* is smooth in the sense explained in Filipović & Teichmann (2003) and has the property that its derivatives of order *n*≥1 are bounded on *U*⊂*H*, then *V* is obviously a *C*^{∞} bounded vector field and additionally is a Banach map-vector field in the sense of Filipović & Teichmann (2003). Such vector fields constitute a class, where the above assumptions can be readily checked.

## 3. The case when *A* is bounded linear

We shall assume in this section that there are no jumps, i.e. we consider(3.1)In order to give an introduction to cubature on the Wiener space, we consider the problem for a bounded operator *A*. In this case, there are virtually no differences to the finite-dimensional setting, except the fact that the drift vector field does not need to be bounded by some constant on the whole Hilbert space (due to the presence of one linear operator in it). Remember that, in Lyons & Victoir (2004) and Kusuoka (2001, 2004), one deals with globally bounded vector fields. We shall circumvent this problem by a small refinement of the arguments.

Since mild and strong solutions coincide, we can always work with strong solutions, which are semi-martingales. Consequently, we can rewrite (3.1) into its Stratonovich form(3.2)where denotes the Stratonovich-corrected drift, i.e.(3.3)wheredenotes the Fréchet derivative of a function or vector field *F*. This notation enables us to writewhere we use the convention that ‘’.

The following notions form the core of cubature on the Wiener space. We give only a short description and refer the reader to Lyons & Victoir (2004) and Teichmann (2006) for more details. Let denote the set of all multi-indices in {0, …, *d*}. We define a degree on by setting*k*∈, . We have to count all the zeros twice owing to the different scalings for and the Brownian motion.

Recall that any vector field *β* can be interpreted as a first-order differential operator on test functions *f* byFor a multi-index , *k*∈, letdenote the corresponding iterated Stratonovich integral. The iterated integrals form the building blocks of the stochastic Taylor formula (see Baudoin 2004).

*Let* *and fix* 0<*t*<1, *m*∈, *m*≥1 *and x*∈*H*. *Then we have**with*

Note that under the assumptions of proposition 3.1, the bound for the remainder term can be infinite. However, if *f* is bounded, *C*^{∞} bounded and has bounded support, then we can guarantee thatIn the unbounded case, this question is more subtle (see example 4.8 and remark 4.9).

Proposition 3.1 shows that iterated Stratonovich integrals play the role of polynomials in the stochastic Taylor expansion. Consequently, it is natural to use them in order to define cubature formulae. Let denote the space of paths of bounded variation taking values in ^{d}. As for the Brownian motion, we append a component for any . Furthermore, we establish the following convention: whenever is the solution to some stochastic differential equation of type (1.1) driven by *B*, whether on a finite- or infinite-dimensional space, and , we denote by the solution of the deterministic differential equation given by formally replacing all the occurrences of ‘’ with ‘’ (with the same initial values; see equation (1.5) as compared with (1.1)). Note that it is necessary that the SDE for *X* is formally formulated in the Stratonovich sense (recall that the Stratonovich formulation does not necessarily make sense for (1.1)).

Fix *t*>0 and *m*≥1. Positive weights *λ*_{1}, …, *λ*_{N} summing up to 1 and paths form a *cubature formula on Wiener space* of degree *m* if for all multi-indices with , *k*∈, we have thatwhere we used the convention in line with the previous one, namely

Lyons & Victoir (2004) show the existence of cubature formulae on the Wiener space for any *d* and some size by applying Chakalov's theorem on cubature formulae and Chow's theorem for nilpotent Lie groups. Moreover, because of the scaling properties of the Brownian motion (and its iterated Stratonovich integrals), i.e.it is sufficient to construct cubature paths for *t*=1.

*Once and for all, we fix one cubature formula* *with weights* *λ*_{1}, …, *λ*_{N} *of degree* *m*≥1 *on the interval* [0,1]*. By abuse of notation, for any* *t*>0*, we will denote* *,* *s*∈[0,*t*]*,* *l*=1, …, *N**, the corresponding cubature formula for* [0,*t*].

For *d*=1 Brownian motions, a cubature formula on the Wiener space of degree *m*=3 is given by *N*=2 pathsfor fixed time horizon *t*. The corresponding weights are given by *λ*_{1}=*λ*_{2}=1/2.

Combining the stochastic Taylor expansion, the deterministic Taylor expansion for solutions of ODEs on *H* and a cubature formula on the Wiener space yields a one-step scheme for weak approximation of the SDE (1.1) for bounded operators *A* in precisely the same way as in Lyons & Victoir (2004). Indeed, we get(3.4)for 0<*t*<1, .

Divide the interval [0,*T*] into *p* subintervals according to the partition . For a multi-index , consider the path defined by concatenating the paths , i.e. for andfor *r* such that , where is scaled to be a cubature path on the interval .

*Fix T*>0 *and m*∈*, a cubature formula of degree m as in* *definition* 3.3 *and a partition of* [0,*T*] *as above. For every* *with**for all* *with* , *k*∈*, there is a constant D independent of the partition such that*

Note that the assumption on *f* in proposition 3.6 is always satisfied if *f* is bounded, *C*^{∞} bounded and has bounded support and if all vector fields *α*, *β*_{1}, …, *β*_{d} have bounded support (compare with remark 3.2). In this case, *P*_{t}*f* has bounded support, too, and all derivatives are bounded on the whole of *H* (see the discussion after theorem 4.4 in the next section).

Proposition 3.6 also yields *deterministic a priori* estimates for the weak rate of convergence, which hold true if we are able to evaluate the respective tree deterministically (see also §6). In principle, there are methods to do so; however, they might be more cumbersome to implement than a Monte Carlo evaluation of the tree.

There is one general case where proposition 3.6 can be applied: namely for stochastic differential equations of type (1.1). If we replace the unbounded operator *A* with a bounded operator Ã, which is close to *A* for a large enough set of values *x*∈(*A*), then we can apply the previous result for bounded linear operators. One candidate for this procedure is the Yosida approximation of *A*.

## 4. The case when *A* is unbounded linear

We shall also assume in this section that there are no jumps and refer to equation (1.1). We recall the analytic background: let us denote by (*A*^{k}) the Hilbert space given by (*A*^{k}) equipped with the graph norm and *k*≥1. Furthermore, we introduce the space(*A*^{∞}) is topologized as a projective limit of the Hilbert space (*A*^{k}), *k*≥0, i.e. the topology on (*A*^{∞}) is the initial topology of the maps , *k*∈. Note that (*A*^{∞}) is no longer a Hilbert space, but only a Fréchet space, i.e. a locally convex vector space that is completely metrizable by a translation invariant metric. We impose our main assumption 2.1.

The following proposition collects a few easy, but interesting corollaries from the existence and uniqueness theorem for equation (1.1) applied to the situation specified in assumption 2.1.

*Fix k*∈*. For any x*∈(*A*^{k})*, there is a unique continuous mild solution of* (3.1) *interpreted as an SDE in the Hilbert space* (*A*^{k})*. If* *, then the mild solution in* (*A*^{k}) *coincides with the mild solution in* (*A*^{k+1})*. Consequently, for x*∈(*A*^{k+1})*, the mild solution of* (1.1) *in* (*A*^{k}) *is a strong solution and, in particular, a semi-martingale*.

If we start in *x*∈(*A*^{∞}), then we get a continuous process , such that is the (mild and strong) solution of (3.1) in any Hilbert space (*A*^{k}), *k*∈. Furthermore, proposition 4.1 allows us to avoid any problems due to the topological structure of (*A*^{∞}) by reverting to an appropriate Hilbert space (*A*^{k}) and interpreting the results again in (*A*^{∞}). The meaning of (*A*^{∞}) is that it is the largest subspace of the Hilbert space *H*, where we can innocently do the necessary analysis on differential operators. Note that there are subtle phenomena of explosion which can occur in this setting: for instance, the support of the law of a strong solution process *X* solving equation (1.1) might be bounded in *H* but unbounded in (*A*), where it is a mild solution. Owing to such phenomena, example 4.8 and remark 4.9 after theorem 4.4 are in fact quite subtle.

As in §3, we introduce the vector field *β*_{0} defined by(4.1)*β*_{0} is defined for *x*∈(*A*). As a vector field taking values in (*A*^{k}), it is only well defined on (*A*^{k+1}). Consequently, for *x*∈(*A*^{k+1}), we may reformulate the SDE (1.1)—understood as equation in (*A*^{k})—in Stratonovich form (3.2).

Now we formulate the stochastic Taylor expansion in some (*A*^{r(m)}) with a degree of regularity *r*(*m*)≥0 depending on *m*≥1. For the estimation of the error term, we will use the *extended support* defined by(4.2)where *t*>0, *x*∈*H* and *ω*_{1}, …, *ω*_{N} are paths of bounded variation. Here, means the support of the law of in (*A*^{r(m)}). Despite assumption 3.4, let us, for one moment, enter the dependence of the cubature formula on the interval [0,*t*] explicitly into the notation, in the sense that are the paths of bounded variation scaled in such a way that they, together with the weights, form a cubature formula on [0,*t*]. Then we denote

If a general support theorem holds in infinite dimensions, we can replace the extended support by the ordinary support of , since the solution of the corresponding ODE driven by paths of bounded variation lie in the support of the solution of the SDE according to the support theorem. However, to our knowledge, no general support theorem has been established for our setting so far.

*Let m*≥1 *be fixed, then there is r*(*m*)≥0 *such that for any* *,* *and* 0<*t*<1 *we have**with**We can choose* *, where* *is the largest integer smaller than m*/2.

The proof is similar to the finite-dimensional situation, but one has to switch between different spaces on the way.

Fix *m* and *f* as above and . We interpret the equation in . By the above remarks, we can express the SDE in its Stratonovich form (3.2). By Itô's formula,(4.3)The idea is to express again by Itô's formula and insert it in equation (4.3). This is completely unproblematic for . For *i*=0, recall thatBy re-expressing (4.3) in Itô formulation, applying Itô's formula and re-expressing it back to Stratonovich formulation, we see thatwhere(4.4)provided that all the new vector fields are well defined and the processes are still semi-martingales. Both conditions are satisfied if —note that the maps , are *C*^{∞}, *k*∈. By induction, we finally getwithNote that *R*_{m} is well defined for because integration of non-semi-martingales with respect to *dt* is possible, which corresponds to the index *i*_{0}=0.

As in the finite-dimensional case, we re-express *R*_{m} in terms of Itô integrals and use the (one-dimensional) Itô isometry several times, until we arrive at the desired estimate. ▪

We recall the notation for bounded measurable functions . Analogous to proposition 3.6, we immediately get the following theorem.

*Fix T*>0*,* *m*≥1*,* *r*(*m*) *and* *as in* *theorem* 4.3*, a cubature formula on the Wiener space of degree m as in* *definition* 3.3 *and a partition* . *Under* *assumption* 2.1*, for any* *with**for all* *with* , *k*∈, *there is a constant D independent of the partition such that**where* *is, again, understood as the mild solution to an ODE in* *for any path ω of bounded variation*.

The proof follows Kusuoka (2001, 2004); see also Kloeden & Platen (1992). For and *x*∈*H* letwhere are scaled to form a cubature formula on [0,*t*]. Denote , *r*=1,…,*p*, the increments of the time partition given in the statement of the theorem. By iterating the operators (and the semigroup property of ODEs), we immediately obtain

(4.5)By ordinary Taylor expansion, keeping in mind the degree function deg, we note thatwhere andIn the sequel *C* denotes a constant independent of the partition and *x*, which may change from line to line. We can estimate the approximation error byCombining this result with theorem 4.3, we may conclude that(4.6)By telescopic sums,For the estimation of the rear termwe may use (4.6) with *f*(*x*) being replaced by , giving usfrom which we may easily conclude the theorem. ▪

Gyöngy & Shmatkov (2006) show a strong Wong–Zakai-type approximation result, where they also need to impose smoothness assumptions on the initial value *x*. Otherwise, the assumptions in Gyöngy & Shmatkov (2006) are different from ours. They allow linear, densely defined vector fields and general adapted coefficients; on the other hand, the generator *A* needs to be elliptic.

Under the previous assumptions, we can also prove a Donsker-type result on the weak convergence of the ‘cubature tree’ to the diffusion. This result will be presented elsewhere.

If *f* is smooth then we can show by (first and higher) variation processes, as introduced for instance in Peszat & Zabczyk (2007), that is smooth on (*A*^{k}).

Fix *k*≥0. Let denote the first variation process of in direction , i.e. is the mild solution to an SDE of the type (3.1). Consequently, it is bounded in and we may conclude thatexists and is bounded by boundedness of D*f* and integrability of the first variation. Similarly, we get existence and continuity of higher order derivatives on (*A*^{k}).

We shall provide examples in which the assumptions of theorem 4.4 are satisfied, i.e. where we obtain high-order convergence of the respective cubature methods. The conditions seem at first sight restrictive (see the following remark 4.9 for a concrete example under assumption 2.1); however, the conditions are parallel to those obtained in Kusuoka (2004) and Lyons & Victoir (2004), where the functions and vector fields have to be bounded and *C*^{∞} bounded. Here, we have an additional complication of one certainly unbounded, but not even continuous drift vector field, which leads to the following set of assumptions.

The vector fields *α*, *β*_{1}, …, *β*_{d} have the following property (compare also with tame maps in Hamilton 1982): there exist smooth mapssuch that for *k*≥0 the restrictions to the respective subspaces take values in , i.e.and such that(4.7)for *i*=1, …, *d*. Here *R*(*λ*, *A*) denotes the resolvent map for some *λ*∈*ρ*(*A*). We assume that have bounded support on . The function *f* is of the same type for a bounded, *C*^{∞} bounded function .

Under these assumptions, we can readily check that the law of the mild solution starting at the initial value has bounded support in *H*: outside a ball of radius *R*>0 in *H* the solution process becomes deterministic, on some interval, hence by the uniform boundedness theorem there is a large *R*′ such that the image of the ball with radius *R*>0 under the maps *S*_{t} lies in a ball with radius *R*′>0 on [0,*T*].

For smooth functions *f* of the stated type, we then haveSince we take only the supremum over bounded sets, namely the extended supports of , this implies the assumption of theorem 4.4.

The previous assumptions (4.7) on the vector fields are not too restrictive since we can always obtain them by a linear isomorphism and (smoothly) cutting off outside a ball in . Both operations are numerically innocent. Under assumption 2.1, we can apply the following isomorphism to the solution of our SDE (1.1):This isomorphism transforms the solution on to an *H*-valued processwith . The transformed process satisfies an SDE, where the transformed vector fields (if well defined) factor over such as previously assumed in assumption (4.7), namely(4.8)The assumptions (4.7) mean that we must (smoothly) cut off the vector fields *α*, *β*_{1}, …, *β*_{d} outside sets of large norm , which is an event—under assumption 2.1—of small probability (recall that the vector fields are Lipschitz on and therefore second moments with respect to the norm exist). Note that the cut-off vector fields do not have an extension to *H* since continuous functions with bounded support on do generically not have a continuous extension onto *H*. For , we can take initial values *y*∈*H*; however, those initial values correspond to quite regular initial values for the original process .

From the point of view of the process *Y*, we have hence proved thatis weakly approximated by evaluating *f* on the cubature tree for *Y*. This is equivalent to evaluating *g* on the cubature tree for *X* in order to approximate the expected value .

## 5. The cubature method in the presence of jumps

The extension of cubature formulae to jump diffusions seems to be new even in the finite-dimensional case. We shall heavily use the fact that only finitely many jumps occur in compact time intervals almost surely.

We shall first prove an asymptotic result on jump diffusions.

*Consider equation* (*1.2*). *Let* *be a bounded measurable function, then we obtain*(5.1)*for t*≥0.

We condition on the jump times and read off the results by inserting the probabilities for a Poisson process with intensity *μ*_{j} to reach level *n*_{j} at time *t*. ▪

This result gives us the time asymptotics with respect to the jump structure. It can now be combined with the original cubature result for the diffusion between the jumps, in order to obtain a result for jump diffusions. We denote by the jump time of the Poisson process *N*^{j} for the *n*th jump. We know that for each Poisson process, the vector is uniformly distributed if conditioned on the event that . The uniform distribution is on the *n*_{j} simplex . This allows us to apply an original cubature formula between two jumps of order , since we gain, for each jump, one order of time asymptotics from the jump structure.

Assume now that the jump distributions *ν*_{j} are concentrated at one point *Z*_{j}≠0, i.e. for *j*=1, …, *e* and *k*≥1. If we want to consider a general jump structure, this amounts to an additional integration with respect to the jump distribution *ν*_{j}.

Now we define a short-time approximation for the conditional expectationsof order with . Expressed in words, we are going to do the following: starting from the initial value we solve the stochastic differential equation (1.2) along the cubature paths *ω*_{l} with probability *λ*_{l}>0, *l*=1, …, *N* until the first jump appears. We collect the endpoints of the trajectories, add the jump size at these points and start a new cubature method from the resulting points on. Note that we can take a cubature method of considerably lower degree since every jump increases the local order of time asymptotics by 1. The jump times are chosen independent and uniformly distributed on simplices of certain dimension *n*_{j} such that . We denote the cubature trajectory between the jump times and for with . If *n*_{j}=0 no trajectories are associated. Hence, we obtain the following theorem.

*Fix m*≥1. *Consider the stochastic differential equation* (*1.2*) *under the condition* *for j*=1, …, *m with* *along concatenated trajectories of type ω*_{l,j,q}. *Choose a cubature method of degree**Concatenation is performed only with increasing q index and a typical concatenated trajectory is denoted by* . *Here we have in mind that the intervals, where the chosen path is ω*_{l,j,q}*, come from a jump of N*^{j} *and have length* . *Then there is* *such that*(5.2)*where* *means the solution of stochastic differential equation* (*1.2*) *in Stratonovich form*(5.3)*along the trajectory ω*.

By our main assumption 2.1, we know that the linkage operators are *C*^{∞} bounded on each (*A*^{k}); hence, through concatenation the errors, which appear on each subinterval , are of the type for some 1≤*q*≤*n*. Taking the supremum yields the result. ▪

Combining the previous result with proposition 5.1 yields under certain conditions on the vector fields (see the discussion after theorem 4.3 in the previous section) by the triangle inequality that there is a constant *D*>0 such that(5.4)

By iteration of the previous result, we obtain in precisely the same manner as in §4 a cubature method of order *m* by applying several cubature methods of order *m*′≤*m* between the jumps.

The only random element in the expectation is given by the jump times , which vary on certain simplices. For the implementation, one has to simulate the uniform distributions on the simplices *tΔ*^{k}. Since the integrals on the simplices *tΔ*^{k} only have continuous integrands, we cannot hope for other methods than Monte Carlo. The evaluation by a Monte Carlo algorithm can also be seen as a random choice of the concatenation grid for the constructed . Another view could be to see a deterministic grid for the diffusion, which is saturated by points where jumps occur. Implementation will be done elsewhere.

## 6. Numerical examples

Owing to the previous results, a numerical scheme for equation (1.1) can be set up by the following steps. Note that we have a weak order of approximation (*m*−1)/2 only under the assumptions of the previous section. In order to obtain those assumptions, one has to modify a general equation of type (1.1) by smoothing procedures. These modifications can be done in a controlled way; more precisely, for each modification we have a rate of convergence to the unmodified object.

Approximate the vector fields

*α*,*β*_{1}, …,*β*_{d}by those satisfying assumption 2.1. If the original vector fields are globally Lipschitz, one can do this approximation with a rate of convergence for the*L*^{2}distance of the original solution process and its approximation.Choose a degree of accuracy

*m*≥2, which determines the weak order of convergence (*m*−1)/2 in the sequel. Associated to*m*the number*r*(*m*) can be identified, which tells us about the degree of regularity of the initial value*x*∈*H*, which one needs for the assertions of theorem 4.4.Identify as a result of the previous specifications a radius

*R*>0 such that the norm of is rarely beyond*R*. Cut off the vector fields smoothly on and verify assumptions (4.7), maybe after smoothing, for the transformed process*Y*such as exercised in remark 4.9.The resulting tree of trajectories yields a finite number of non-autonomous PDEs (1.5), which have to be evaluated. In the implementation, one calculates with the smoothened vector fields satisfying assumption 2.1, but not with the cut-off vector fields, since a large norm value

*R*is reached with small probability because of its very choice.

We test the cubature method for two concrete examples: one toy example, where explicit solutions of the stochastic partial differential equations (SPDE) are readily available, and the other, more interesting but still very easy example. Since cubature on the Wiener space is a weak method, we calculate the expected value of a functional of the solution to the SPDE in both cases; that is, the outputs of our computations are real numbers.

The results presented here are calculated in Matlab using the built-in PDE solver `pdepe` for solving the deterministic PDEs given by inserting the cubature paths into the SPDE under consideration. This PDE solver depends on a space grid given by the user as well as on a time grid, which is not very critical because it is adaptively refined by the program.

We use the simplest possible cubature formula for *d*=1 Brownian motion:with weights *λ*_{1}=*λ*_{2}=1/2 define a cubature formula of degree *m*=3 on [0,*T*]. Consequently, solving an SDE on a Hilbert space with cubature on the Wiener space for the above cubature formula and *p* iterations means solving 2^{p} PDEs. This starts to get restrictive even for a very simple problem for, say, *p*=10, where one already has to solve more than 1000 PDEs. One possibility to overcome these tight limitations is to use ‘a Monte Carlo simulation on the tree’. Recall that a *p*-step cubature method approximatesSince , we can interpret the r.h.s. as the expectation of a random variable on the tree . Therefore, we can approximate the r.h.s. by picking tuples at random—according to their probabilities —and calculating the average of the corresponding outcomes . Of course, by following this strategy, we have to replace the deterministic error estimates by a stochastic one, which heavily depends on the standard deviation of understood as a random variable on the tree. Note, however, that this is the usual situation for weak approximation methods for SDEs.

Consider the Ornstein–Uhlenbeck process, defined as the solution to the equation(6.1)on the Hilbert space . *Δ* denotes the Dirichlet Laplacian on ]0,1[, i.e. *Δ* is a negative definite self-adjoint operator on *H* with extending the classical Laplace operator defined on . It is easy to see that *Δ* is dissipative and therefore, by the Lumer–Phillips theorem, it is the generator of a *C*_{0} contraction semigroup on *H*. The coefficient *ϕ*∈*H* is some fixed vector.

In this case, the definition of a mild solution(6.2)already gives a representation of the solution provided that the heat semigroup *S*_{t} applied to the starting vector *x* and to *ϕ* is available. We choose *x*(*u*)=sin(*πu*), u∈]0,1[, and may conclude thatbecause *x* is an eigenvector of *Δ* with eigenvalue −*π*^{2}. Consider the linear functional given by(6.3)We want to compute

In table 1, the error, i.e. the output of the method minus the true value given above, is presented. *p* is the number of cubature steps, i.e. the number of iterations of the one-step cubature method. The discretization in space, i.e. of ]0,1[, used by the PDE solver contains 50 uniform points and the discretization of the time interval—additional to the one induced by the cubature method—contains 500 points. The stochastic perturbation factor *ϕ* is chosen to be *ϕ*(*u*)=sin(*πu*), i.e. even. We see a very fast decrease of the error in this simple situation. On the other hand, the variance of the random variable on the tree considered before is clearly too high for the Monte Carlo simulation on the tree to work. Indeed, has true standard deviation of(6.4)Assuming that the central limit theorem applies, confidence intervals around the solution given by a Monte Carlo method are proportional to the standard deviation divided by the square root of the number of trajectories. Consequently, we would roughly need to calculate 10^{12} paths on the tree in order to achieve a similar level of exactness as in table 1! Indeed, note that the standard deviation of the solution is of order 10^{−1}, while the error in the last row of table 1 is of order 10^{−7}. The equationthen gives *m*≈10^{12}. Note that this heuristics is also confirmed by our experiments, where Monte Carlo simulation on the tree clearly fails. The data are also shown in figure 1.

The failure of Monte Carlo simulation on the tree also applies to any other (naive) Monte Carlo approach to problem (6.1), including the usual finite-element or finite-difference approaches.

As a more realistic example, we consider the heat equation with a stochastic perturbation involving a Nemicky operator. More precisely, consider(6.5)with *x*(*u*)=sin(*πu*). Even though we do not know the law of the solution of (6.5), we are still able to calculate explicitly because *Φ* is a linear functional. Indeed, is given by(6.6)and, consequently,

The expectation of the (one-dimensional) Itô integral is 0 and we get the same result as before, i.e.for *x*(*u*)=sin(*πu*). Nevertheless, we believe that this example is already quite difficult, especially since the cubature method actually has to work with the Stratonovich formulation(6.7)In particular, the equation (in Stratonovich form) has a nonlinear drift and a nonlinear volatility.

Note that we expect the standard deviation of the solution of the above equation to be smaller than before, because decreases as decreases in *t*. The space discretization used by the PDE solver has size 50, which already seems to be sufficient, because using a finer discretization (100 grid points) does not change the results significantly. Table 3 shows the results using Monte Carlo simulation on the tree. *m* denotes the number of trajectories followed, while the ‘statistical error’ in the table is an indicator for the error of the Monte Carlo simulation. More precisely, the values in the last column are the empirical standard deviations of the result divided by the square root of the number of trajectories. Comparable to the Ornstein–Uhlenbeck process, the convergence of the pure cubature method is very fast (table 2). The (empirical) variance is, however, quite large such that the Monte Carlo-aided method does not work at all. Note that the statistical error in table 3 is of the order of the total computational error, which can be almost completely attributed to the Monte Carlo simulation.

To test the method further, we also try more irregular data. Let(6.8)The exact value of the quantity of interest is calculated by solving the corresponding heat equation numerically. This gives the value . The initial vector *x* given in (6.8) is in but its derivative is no longer square integrable. Consequently, and the theory does not provide an order of approximation. Nevertheless, probably because of the smoothing properties of the Laplace operator, numerical results show the same behaviour as before (figure 1).

If we replace heat equation (6.1) by an evolution equation of the form(6.9)then we still see the same kind of behaviour if we fix the space discretization for the PDE solver. This time, the PDEs require a much finer space resolution in order to give reliable numbers.

## Acknowledgments

The authors gratefully acknowledge the support from the FWF grant Y328.

## Footnotes

- Received January 14, 2008.
- Accepted April 23, 2008.

- © 2008 The Royal Society