We consider a one-dimensional continuous stochastic process (Xt)t∈[0,T] that satisfies the stochastic differential equation (1.1) dXt=b(t,Xt)dt+σ(t,Xt)dWtwithX0=x0, where (Wt)t∈[0,T] is a Brownian motion defined on a filtered probability space (Ω,F,(Ft)t∈[0,T],P) . Various numerical approximation schemes for the SDE (1.1) have been proposed and studied in the literature in the past decades. The probably most prominent numerical approximation for the solution of (1.1) is the Euler scheme, which can be described as follows. Let φn:R+→R+ be the function defined by φn(t)=i/n when t∈[i/n,(i+1)/n) . The continuous Euler approximation scheme is described by (1.2) dXtn=b(φn(t),Xφn(t)n)dt+σ(φn(t),Xφn(t)n)dWtwithX0n=x0. The probabilistic properties of the Euler approximation scheme have been investigated in numerous papers. We refer to the classical works [2, 3, 14, 16, 17, 23] for the studies on weak and strong approximation errors among many others. Asymptotic results in the framework of non-regular coefficients can be found in e.g. [1, 8, 9, 27].

In this paper we take a different route and propose to use the Wiener chaos expansion (also called polynomial chaos in the literature) to approximate the solution of the SDE (1.1). To explain ideas let us fix an orthonormal basis (ei)i≥1 of the separable Hilbert space L2([0,T]) . It is a well-known statement (cf. [6]) that if Xt∈L2(Ω,F,P) for all t∈[0,T] , where F is generated by the Brownian motion (Wt)t∈[0,T] , it admits the chaotic expansion (1.3) Xt=∑α∈Ixα(t)Ψα, where xα are deterministic functions, Ψα are mutually orthogonal random projections that are associated to the basis (ei)i≥1 , and the index set I is defined via (1.4) I:={α=(αi)i≥1:αi∈N0and almost allαi’s are0}. Such orthogonal expansions have been successfully applied in numerous fields of stochastic and numerical analysis. We refer e.g. to [11–13, 25] for applications of the Wiener chaos expansion in the context of SDEs and to [20–22, 26] for applications of polynomial expansion in modelling, simulation and filtering of stochastic partial differential equations. The aim of this work is to use the chaos expansion (1.3) to numerically approximate the solution of the SDE (1.1). For this purpose we need to truncate the infinite sum in (1.3). An obvious approach is to consider the approximation (1.5) Xtp,k=∑α∈Ip,kxα(t)Ψα, where the subset Ip,k⊂I refers to using orthogonal projections Ψα only with respect to the first k basis elements (ei)1≤i≤k and only up to the pth order Wiener chaos. This method is mostly related to the articles [11, 12, 22, 21]. More specifically, the works [21, 22] study the L2 -error associated with the approximation (1.5), but only for a particular choice of the basis (ei)i≥1 . In this paper we will study the decay rate of E[(Xt−Xtp,k)2] when k,p→∞ for a general basis (ei)i≥1 of L2([0,T]) applying methods from Malliavin calculus.

The paper is structured as follows. In Section 2 we present the elements of Malliavin calculus. The main results of the paper are demonstrated in Section 3. Section 4 is devoted to proofs. In Section 5 we illustrate our approach with exemplary numerical results for the Haar and a trigonometric basis, and propose a heuristic based on sparse indices for computational speedup.

In this section we introduce some basic concepts of Malliavin calculus. The interested readers are referred to [24] for more thorough exposition on this subject. Set H=L2([0,T]) and let ⟨·,·⟩H denote the scalar product on H . We note that H is a separable Hilbert space and denote by (ei)i≥1 an orthonormal basis of H . We consider the isonormal Gaussian family W={W(h):h∈H} indexed by H defined on a probability space (Ω,F,P) , i.e. the random variables W(h) are centered Gaussian with a covariance structure determined via E[W(g)W(h)]=⟨g,h⟩H. Here and throughout the paper we assume that F=σ(W) . In our setting we consider W(h)=∫0ThsdWs , where W is a standard Brownian motion. We define the normalised Hermite polynomials through the identities (2.1) H0(x):=1,Hn(x):=(−1)nn!exp(x22)dndxn(exp(−x22)),n≥1. The nth Wiener chaos Hn is the closed linear subspace of L2(Ω,F,P) generated by the family of random variables {Hn(W(h)):‖h‖H=1} . The vector spaces Hn , n≥0 , are orthogonal and we have the Wiener chaos expansion (2.2) L2(Ω,F,P)= ⨁n=0∞Hn. (See [24, Theorem 1.1.1] for more details.) Next, for α∈I , where I has been introduced in (1.4), we define the random variable (2.3) Ψα:= ∏i=1∞Hαi(W(ei)), where (ei)i≥1 is a fixed orthonormal basis of H . We define |α|=∑i=1∞αi and α!=∏i=1∞αi! for α∈I . We deduce that the set {Ψα:α∈Iwith|α|=n} forms a complete orthonormal basis of the nth Wiener chaos Hn and consequently {Ψα:α∈I} is a complete orthonormal basis of L2(Ω,F,P) (cf. [24, Proposition 1.1.1]).

Now, we introduce multiple stochastic integrals of order n∈N , which are denoted by In . For an element h⊗n:=h⊗⋯⊗h in H⊗n with ‖h‖H=1 , we define (2.4) In(h⊗n):=n!Hn(W(h)),n≥1. Assuming that the mapping In is linear, the definition (2.4) can be extended to all symmetric elements h∈H⊗n by polarisation identity. Finally, for an arbitrary function h∈H⊗n we set In(h):=In(h˜), where h˜ denotes the symmetrised version of h∈H⊗n . By definition In maps H into Hn , so the multiple integrals of different orders are orthogonal. In particular, they satisfy the isometry property (2.5) E[Im(g)In(h)]=n!⟨g˜,h˜⟩H⊗n1{n=m},h∈H⊗n,g∈H⊗m. Furthermore, for any symmetric h∈H⊗n,g∈H⊗m , the following multiplication formula holds: (2.6) Im(g)In(h)= ∑r=0min(m,n)r!mrnrIm+n−2r(g⊗rh), where the rth contraction g⊗rh is defined by g⊗rh(t1,…,tm+n−2r):=∫[0,T]rg(u1,…,ur,t1,…tm−r)×h(u1,…,ur,tm−r+1,…tm+n−2r)du1…dur. (See [24, Proposition 1.1.3].) The Wiener chaos expansion (2.2) transfers to the context of multiple integrals as follows. For each random variable F∈L2(Ω,F,P) we obtain the orthogonal decomposition (2.7) F= ∑n=0∞In(gn),gn∈H⊗n, where I0=E[F] and the decomposition is unique when gn , n≥1 , are assumed to be symmetric (cf. [24, Theorem 1.1.2]).

Next, we introduce the notion of Malliavin derivative and its adjoint operator. We define the set of smooth random variables via S={F=f(W(h1),…,W(hn)):n≥1,hi∈H}, where f∈Cp∞(Rn) (i.e. the space of infinitely differentiable functions such that all derivatives exhibit polynomial growth). The kth order Malliavin derivative of F∈S , denoted by DkF , is defined by (2.8) DkF= ∑i1,…,ik=1n∂k∂xi1⋯∂xikf(W(h1),…,W(hn))hi1⊗⋯⊗hik. Notice that DkF is a H⊗k -valued random variable, and we write DxkF for the realisation of the function DkF at the point x∈[0,T]k . The space Dk,q denotes the completion of the set S with respect to the norm ‖F‖k,q:=(E[|F|q]+ ∑m=1kE[‖DmF‖H⊗mq])1/q. We define Dk,∞=∩q>1Dk,q }1}}{\mathbb{D}^{k,q}}$]]>. The Malliavin derivative of the random variable Ψα∈S , α∈I , can be easily computed using the definition (2.8) and the formula Hn′(x)=nHn−1(x) : (2.9) DsΨα= ∑i=1∞αiei(s)Ψα−(i),where (2.10) α−(i):=(α1,…,αi−1,αi−1,αi+1,…)ifαi≥1. Higher order Malliavin derivatives of Ψα are computed similarly.

The operator Dk possesses an unbounded adjoint denoted by δk , which is often referred to as multiple Skorokhod integral. The following integration by parts formula holds (see [24, Exercise 1.3.7]): if u∈Dom(δk) and F∈Dk,2 , where Dom(δk) consists of all elements u∈L2(Ω;H⊗k) such that the inequality |E[⟨DkF,u⟩H⊗k]|≤c(E[F2])1/2 holds for some c>0 }0$]]> and all F∈Dk,2 , then we have the identity (2.11) E[Fδk(u)]=E[⟨DkF,u⟩H⊗k]. In case the random variable F∈L2(Ω,F,P) has a chaos decomposition as displayed in (2.7), the statement F∈Dk,2 is equivalent to the condition ∑n=1∞nkn!‖gn‖H⊗n2<∞ . In particular, when F∈D1,2 we deduce an explicit chaos representation of the derivative DF=(DtF)t∈[0,T] (2.12) DtF= ∑n=1∞nIn−1(gn(·,t)),t∈[0,T], where gn(·,t):[0,T]n−1→R is obtained from the function gn by setting the last argument equal to t (see [24, Proposition 1.2.7 ]). Finally, we present an explicit formula for the Malliavin derivative of a solution of a stochastic differential equation. Assume that (Xt)t∈[0,T] is a solution of a stochastic differential equation (1.1) and b,σ∈C1(R) . In this setting DXt=(DsXt)s∈[0,T] is given as the solution of the SDE (2.13) DsXt=σ(s,Xs)+∫st∂∂xb(u,Xu)Ds(Xu)du+∫st∂∂xσ(u,Xu)Ds(Xu)dWu, for s≤t , and DsXt=0 if s>t }t$]]> (see [24, Theorem 2.2.1]). Throughout the paper ∂/∂x denotes the derivative with respect to the space variable and ∂/∂t denotes the derivative with respect to the time variable.

We start with the analysis of the Wiener chaos expansion introduced in (1.3). Under square integrability assumption on the solution X of the SDE (1.1) we obtain the Wiener chaos expansion Xt=∑α∈Ixα(t)ΨαwithΨα= ∏i=1∞Hαi(W(ei)). In order to study the strong approximation error we require a good control of the coefficients xα(t) . When Xt is sufficiently smooth in the Malliavin sense, we deduce the identity (3.1) xα(t)=E[XtΨα]=1α!E[⟨D|α|Xt, ⨂i=1∞ei⊗αi⟩H⊗|α|] applying the duality formula (2.11). In fact, the coefficients xα(t) satisfy a system of ordinary differential equations. The following propagator system has been derived in [12, 21]. We state the proof for completeness.

We remark that the propagator system (3.2) is recursive. Let us give some simple examples to illustrate how (3.2) can be solved explicitly.

For a general specification of the drift coefficient b and diffusion coefficient σ in model (1.1) the propagator system (3.2) cannot be solved explicitly. Thus, precise infinite dimensional Wiener chaos expansion (1.3) is out of reach. For simulation purposes it is an obvious idea to consider a finite subset of I in the expansion (1.3). We introduce the index set (3.3) Ip,k:={α∈I:|α|≤pandαi=0for alli>k}. }k\big\}.\]]]> The approximation of Xt is now defined via (1.5): (3.4) Xtp,k=∑α∈Ip,kxα(t)Ψα. We remark that the quantity Xtp,k is more useful than the Euler approximation Xtn introduced in (1.2) if we are interested in the approximation of the first two moments of Xt . Indeed, the first two moments of Xtp,k are given explicitly by E[Xtp,k]=x0(t)andE[(Xtp,k)2]=∑α∈Ip,kxα2(t), while higher order moments can be computed via an application of the multiplication formula (2.6).

The strong approximation error associated with the truncation (3.4) has been previously studied in [21, 22] in a slightly different context. In particular, in both papers the authors consider one specific basis (ei)i≥1 of L2([0,T]) whereas we are interested in the asymptotic analysis for a general basis (ei)i≥1 . While [22] mostly uses methods from analysis, our approach is based upon Malliavin calculus and is close in spirit to [21]. The main result of our paper gives an upper bound on the L2 -error E[(Xtp,k−Xt)2] .

Let us give some remarks about the statement (3.6). First of all, recalling the Karhunen–Loéve expansion Wt=∑l=1∞El(t)∫0Tel(s)dWs , we readily deduce that E[Wt2]= ∑l=1∞El2(t)<∞and∫0tE[Ws2]ds= ∑l=1∞∫0tEl2(τ)dτ<∞. Hence, the upper bound on the right-hand side of (3.6) indeed converges to 0 when p,k→∞ . We also note that the error associated with the truncation of the chaos order does not depend on the particular choice of the basis (ei)i≥1 , while the error associated with truncation of basis strongly depends on (ei)i≥1 (which is not really surprising). Furthermore, we remark that it is extremely computationally costly to compute all coefficients xα(t) , α∈Ip,k , for a large chaos order p. Thanks to the first bound ((p+1)!)−1 in (3.6) it is sufficient to use small values of p in practical situations (usually p≤4 ). Example 3.4.

Let us explicitly compute the last terms of the upper bound (3.6) for two prominent bases of L2([0,1]) .

(i) (Trigonometric basis) Consider the orthonormal basis (ei)i≥1 given by e1(t)=1,e2j(t)=2sin(2πjt),e2j+1(t)=2cos(2πjt),j≥1. In this setting we obtain that E2j2(t)=12π2j2(1−cos(2πjt))2,E2j+12(t)=12π2j2(1−sin(2πjt))2, for any j≥1 . Consequently, we deduce that (3.8) ∑l=k+1∞(El2(t)+∫0tEl2(τ)dτ)≤Ck for some C>0 }0$]]>.

(ii) (Haar basis) The Haar basis is a collection of functions {e0,ej,n:j=1,…,2n−1,n≥1} defined as follows: e0(t)=1,ej,n(t)=2(n−1)/2:t∈[2−n+1(j−1),2−n(2j−1)],−2(n−1)/2:t∈[2−n(2j−1),2−n+1j],0:else. In this case we deduce the following representation for Ej,n(t) : Ej,n(t)=2(n−1)/2(t−2−n+1(j−1)):t∈[2−n+1(j−1),2−n(2j−1)],2(n−1)/2(2−n+1j−t):t∈[2−n(2j−1),2−n+1j],0:else. Since maxt∈[0,1]Ej,n2(t)=2−(n+1) and Ej,n(t)≠0 only for t∈[2−n+1(j−1),2−n+1j] , we finally obtain that (3.9) ∑l=n+1∞ ∑j=12n(Ej,l2(t)+∫0tEj,l2(τ)dτ)≤C2−n. Note that the exponential asymptotic rate for a fixed pth order Wiener chaos is equivalent to O(k−1) , if we choose k=2n basis elements (ei) . We will come back to these two bases in Section 5 where we provide exemplary numerical results that highlight this identical asymptotic rate, but also the different error distributions over time.

Throughout the proofs C=C(t,K) denotes a generic positive constant, which might change from line to line. We start with a proposition that gives an upper bound for the L2 -norm of the Malliavin derivatives of Xt .