We consider a one-dimensional continuous stochastic process ${({X_{t}})}_{t\in [0,T]}$ that satisfies the stochastic differential equation where ${({W_{t}})}_{t\in [0,T]}$ is a Brownian motion defined on a filtered probability space $(\varOmega ,\mathcal{F},{({\mathcal{F}_{t}})}_{t\in [0,T]},\mathbb{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 ${\varphi _{n}}:{\mathbb{R}_{+}}\to {\mathbb{R}_{+}}$ be the function defined by ${\varphi _{n}}(t)=i/n$ when $t\in [i/n,(i+1)/n)$. The continuous Euler approximation scheme is described by 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 ${({e_{i}})}_{i\ge 1}$ of the separable Hilbert space ${L^{2}}([0,T])$. It is a well-known statement (cf. [6]) that if ${X_{t}}\in {L^{2}}(\varOmega ,\mathcal{F},\mathbb{P})$ for all $t\in [0,T]$, where $\mathcal{F}$ is generated by the Brownian motion ${({W_{t}})}_{t\in [0,T]}$, it admits the chaotic expansion where ${x_{\alpha }}$ are deterministic functions, ${\varPsi ^{\alpha }}$ are mutually orthogonal random projections that are associated to the basis ${({e_{i}})}_{i\ge 1}$, and the index set $\mathcal{I}$ is defined via 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 where the subset ${\mathcal{I}_{p,k}}\subset \mathcal{I}$ refers to using orthogonal projections ${\varPsi ^{\alpha }}$ only with respect to the first k basis elements ${({e_{i}})}_{1\le i\le 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 ${L^{2}}$-error associated with the approximation (1.5), but only for a particular choice of the basis ${({e_{i}})}_{i\ge 1}$. In this paper we will study the decay rate of $\mathbb{E}[{({X_{t}}-{X_{t}^{p,k}})^{2}}]$ when $k,p\to \infty $ for a general basis ${({e_{i}})}_{i\ge 1}$ of ${L^{2}}([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 $\mathbb{H}={L^{2}}([0,T])$ and let ${\langle \cdot ,\cdot \rangle _{\mathbb{H}}}$ denote the scalar product on $\mathbb{H}$. We note that $\mathbb{H}$ is a separable Hilbert space and denote by ${({e_{i}})}_{i\ge 1}$ an orthonormal basis of $\mathbb{H}$. We consider the isonormal Gaussian family $W=\{W(h):h\in \mathbb{H}\}$ indexed by $\mathbb{H}$ defined on a probability space $(\varOmega ,\mathcal{F},\mathbb{P})$, i.e. the random variables $W(h)$ are centered Gaussian with a covariance structure determined via Here and throughout the paper we assume that $\mathcal{F}=\sigma (W)$. In our setting we consider $W(h)={\int _{0}^{T}}{h_{s}}d{W_{s}}$, where W is a standard Brownian motion. We define the normalised Hermite polynomials through the identities The nth Wiener chaos ${\mathcal{H}_{n}}$ is the closed linear subspace of ${L^{2}}(\varOmega ,\mathcal{F},\mathbb{P})$ generated by the family of random variables $\{{H_{n}}(W(h)):\hspace{2.5pt}\| h{\| _{\mathbb{H}}}=1\}$. The vector spaces ${\mathcal{H}_{n}}$, $n\ge 0$, are orthogonal and we have the Wiener chaos expansion (See [24, Theorem 1.1.1] for more details.) Next, for $\alpha \in \mathcal{I}$, where $\mathcal{I}$ has been introduced in (1.4), we define the random variable where ${({e_{i}})}_{i\ge 1}$ is a fixed orthonormal basis of $\mathbb{H}$. We define $|\alpha |={\sum _{i=1}^{\infty }}{\alpha _{i}}$ and $\alpha !={\prod _{i=1}^{\infty }}{\alpha _{i}}!$ for $\alpha \in \mathcal{I}$. We deduce that the set $\{{\varPsi ^{\alpha }}:\hspace{2.5pt}\alpha \in \mathcal{I}\hspace{2.5pt}\text{with}\hspace{2.5pt}|\alpha |=n\}$ forms a complete orthonormal basis of the nth Wiener chaos ${\mathcal{H}_{n}}$ and consequently $\{{\varPsi ^{\alpha }}:\hspace{2.5pt}\alpha \in \mathcal{I}\}$ is a complete orthonormal basis of ${L^{2}}(\varOmega ,\mathcal{F},\mathbb{P})$ (cf. [24, Proposition 1.1.1]).

Now, we introduce multiple stochastic integrals of order $n\in \mathbb{N}$, which are denoted by ${I_{n}}$. For an element ${h^{\otimes n}}:=h\otimes \cdots \otimes h$ in ${\mathbb{H}^{\otimes n}}$ with $\| h{\| _{\mathbb{H}}}=1$, we define Assuming that the mapping ${I_{n}}$ is linear, the definition (2.4) can be extended to all symmetric elements $h\in {\mathbb{H}^{\otimes n}}$ by polarisation identity. Finally, for an arbitrary function $h\in {\mathbb{H}^{\otimes n}}$ we set where $\widetilde{h}$ denotes the symmetrised version of $h\in {\mathbb{H}^{\otimes n}}$. By definition ${I_{n}}$ maps $\mathbb{H}$ into ${\mathcal{H}_{n}}$, so the multiple integrals of different orders are orthogonal. In particular, they satisfy the isometry property Furthermore, for any symmetric $h\in {\mathbb{H}^{\otimes n}},g\in {\mathbb{H}^{\otimes m}}$, the following multiplication formula holds: where the rth contraction $g{\otimes _{r}}h$ is defined by (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\in {L^{2}}(\varOmega ,\mathcal{F},\mathbb{P})$ we obtain the orthogonal decomposition where ${I_{0}}=\mathbb{E}[F]$ and the decomposition is unique when ${g_{n}}$, $n\ge 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 where $f\in {C_{p}^{\infty }}({\mathbb{R}^{n}})$ (i.e. the space of infinitely differentiable functions such that all derivatives exhibit polynomial growth). The kth order Malliavin derivative of $F\in \mathcal{S}$, denoted by ${D^{k}}F$, is defined by Notice that ${D^{k}}F$ is a ${\mathbb{H}^{\otimes k}}$-valued random variable, and we write ${D_{x}^{k}}F$ for the realisation of the function ${D^{k}}F$ at the point $x\in {[0,T]^{k}}$. The space ${\mathbb{D}^{k,q}}$ denotes the completion of the set $\mathcal{S}$ with respect to the norm We define ${\mathbb{D}^{k,\infty }}={\cap _{q\mathrm{>}1}}{\mathbb{D}^{k,q}}$. The Malliavin derivative of the random variable ${\varPsi ^{\alpha }}\in \mathcal{S}$, $\alpha \in \mathcal{I}$, can be easily computed using the definition (2.8) and the formula ${H^{\prime }_{n}}(x)=\sqrt{n}{H_{n-1}}(x)$: Higher order Malliavin derivatives of ${\varPsi ^{\alpha }}$ are computed similarly.

The operator ${D^{k}}$ possesses an unbounded adjoint denoted by ${\delta ^{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\in \text{Dom}({\delta ^{k}})$ and $F\in {\mathbb{D}^{k,2}}$, where $\text{Dom}({\delta ^{k}})$ consists of all elements $u\in {L^{2}}(\varOmega ;{\mathbb{H}^{\otimes k}})$ such that the inequality $|\mathbb{E}[{\langle {D^{k}}F,u\rangle _{{\mathbb{H}^{\otimes k}}}}]|\le c{(E[{F^{2}}])^{1/2}}$ holds for some $c\mathrm{>}0$ and all $F\in {\mathbb{D}^{k,2}}$, then we have the identity In case the random variable $F\hspace{0.1667em}\in \hspace{0.1667em}{L^{2}}(\varOmega ,\mathcal{F},\mathbb{P})$ has a chaos decomposition as displayed in (2.7), the statement $F\in {\mathbb{D}^{k,2}}$ is equivalent to the condition ${\sum _{n=1}^{\infty }}{n^{k}}n!\| {g_{n}}{\| _{{\mathbb{H}^{\otimes n}}}^{2}}\mathrm{<}\infty $. In particular, when $F\in {\mathbb{D}^{1,2}}$ we deduce an explicit chaos representation of the derivative $DF={({D_{t}}F)}_{t\in [0,T]}$ where ${g_{n}}(\cdot ,t):{[0,T]^{n-1}}\to \mathbb{R}$ is obtained from the function ${g_{n}}$ 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 ${({X_{t}})}_{t\in [0,T]}$ is a solution of a stochastic differential equation (1.1) and $b,\sigma \in {C^{1}}(\mathbb{R})$. In this setting $D{X_{t}}={({D_{s}}{X_{t}})}_{s\in [0,T]}$ is given as the solution of the SDE for $s\le t$, and ${D_{s}}{X_{t}}=0$ if $s\mathrm{>}t$ (see [24, Theorem 2.2.1]). Throughout the paper $\partial /\partial x$ denotes the derivative with respect to the space variable and $\partial /\partial 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 In order to study the strong approximation error we require a good control of the coefficients ${x_{\alpha }}(t)$. When ${X_{t}}$ is sufficiently smooth in the Malliavin sense, we deduce the identity applying the duality formula (2.11). In fact, the coefficients ${x_{\alpha }}(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 $\mathcal{I}$ in the expansion (1.3). We introduce the index set The approximation of ${X_{t}}$ is now defined via (1.5): We remark that the quantity ${X_{t}^{p,k}}$ is more useful than the Euler approximation ${X_{t}^{n}}$ introduced in (1.2) if we are interested in the approximation of the first two moments of ${X_{t}}$. Indeed, the first two moments of ${X_{t}^{p,k}}$ are given explicitly by 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 ${({e_{i}})}_{i\ge 1}$ of ${L^{2}}([0,T])$ whereas we are interested in the asymptotic analysis for a general basis ${({e_{i}})}_{i\ge 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 ${L^{2}}$-error $\mathbb{E}[{({X_{t}^{p,k}}-{X_{t}})^{2}}]$.

Let us give some remarks about the statement (3.6). First of all, recalling the Karhunen–Loéve expansion ${W_{t}}={\sum _{l=1}^{\infty }}{E_{l}}(t){\int _{0}^{T}}{e_{l}}(s)d{W_{s}}$, we readily deduce that Hence, the upper bound on the right-hand side of (3.6) indeed converges to 0 when $p,k\to \infty $. We also note that the error associated with the truncation of the chaos order does not depend on the particular choice of the basis ${({e_{i}})}_{i\ge 1}$, while the error associated with truncation of basis strongly depends on ${({e_{i}})}_{i\ge 1}$ (which is not really surprising). Furthermore, we remark that it is extremely computationally costly to compute all coefficients ${x_{\alpha }}(t)$, $\alpha \in {\mathcal{I}_{p,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\le 4$).

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 ${L^{2}}$-norm of the Malliavin derivatives of ${X_{t}}$.

Now, we proceed with the proof of the main result of Theorem 3.3, which follows the ideas of [21, Proof of Theorem 2.2]. Observe the decomposition where $d(\alpha ):=\max \{i\ge 1:\hspace{2.5pt}{\alpha _{i}}\mathrm{>}0\}$. To determine an upper bound for ${A_{1}}(p)$, we use the chaos expansions (1.3) and (2.7) of ${X_{t}}$, i.e., with ${\mathbf{t}^{n}}=({t_{1}},\dots ,{t_{n}})$ and symmetric kernel functions ${\xi _{n}}({\mathbf{t}^{n}};t)$ being defined by for all $n\in {\mathbb{N}_{0}}$. We conclude that where we abbreviate Therefore, we deduce via (4.7) that Finally, using (4.8) and applying Proposition 4.1 we obtain The treatment of the term ${A_{2}}(p,k)$ is more involved. Recalling the definition of $d(\alpha )=\max \{i\ge 1:\hspace{2.5pt}{\alpha _{i}}\mathrm{>}0\}$, we introduce the characteristic set $({i_{1}},\dots ,{i_{n}})$ of $\alpha \in \mathcal{I}$ for ${i_{1}}\le {i_{2}}\le \cdots \le {i_{n}}$ and $|\alpha |=n$. It is defined as follows: ${i_{1}}$ is the index of the first non-zero component of α. If ${\alpha _{{i_{1}}}}=1$ then ${i_{2}}$ is the index of the second non-zero component of α; otherwise ${i_{1}}={i_{2}}=\cdots ={i_{{\alpha _{{i_{1}}}}}}$ and ${i_{{\alpha _{{i_{1}}}}+1}}$ is the second non-zero component of α. The same operation is repeated for the index ${i_{{\alpha _{{i_{1}}}}+1}}$ and further non-zero components of α. In this fashion, the characteristic set is constructed, resulting in the observation that $d(\alpha )={i_{n}}$. To give a simple example, consider the multiindex $\alpha =(2,0,1,4,0,0,\dots )\in \mathcal{I}$. Here the characteristic set of α is given by For any $\alpha \in \mathcal{I}$ with $|\alpha |=n$ and a basis ${({e_{i}})}_{i\ge 1}$ of $\mathbb{H}={L^{2}}([0,T])$, we denote by ${\widetilde{e}_{\alpha }}({\mathbf{t}^{n}})$ the (scaled) symmetrised form of $\bigotimes {e_{i}^{\otimes {\alpha _{i}}}}$ defined via where $({i_{1}},\dots ,{i_{n}})$ is the characteristic set of α and the sum runs over all permutations π within the permutation group ${\mathfrak{P}^{n}}$ of $\{1,\dots ,n\}$. From (3.1) we know that ${x_{\alpha }}(t)=\mathbb{E}[{X_{t}}{\varPsi ^{\alpha }}]$. On the other hand, we have the representation Now, for any $\alpha \in \mathcal{I}$ with $|\alpha |=n$, we obtain the identity by (2.11). Since where ${\mathbf{t}_{j}^{n}}$ is obtained from ${\mathbf{t}^{n}}$ by omitting ${t_{j}}$ and ${\alpha ^{-}}(\cdot )$ denotes the diminished multi-index as defined in (2.10), we deduce that with ${t_{0}}:=0$ and ${t_{n+1}}:=t$, by changing the order of integration. Then for any ${i_{n}}=l\ge 1$ we integrate by parts to deduce where Now, we use substitution in (4.12) by renaming ${\mathbf{t}_{j}^{n}}$ in the following way for each j: With ${s_{i}}={t_{i}}$ for all $i\le j-1$ and ${s_{i}}={t_{i+1}}$ for all $i\mathrm{>}j-1$, we have ${\mathbf{s}^{n-1}}:={\mathbf{t}_{j}^{n}}$ by setting ${s_{0}}=0$ and ${s_{n}}=t$. Moreover, we denote with ${\mathbf{s}^{n-1,r}}$, $r=1,\dots ,n-1$, the set that is generated from ${\mathbf{s}^{n-1}}$ by taking ${s_{r}}$ twice. To finalize this notation, we set ${\mathbf{s}^{n-1,0}}=({s_{0}},{s_{1}},\dots ,{s_{n-1}})$ and ${\mathbf{s}^{n-1,n}}=({s_{1}},\dots ,{s_{n-1}},{s_{n}})$. Then Because ${E_{l}}({s_{0}})={E_{l}}(0)=0$ and ${E_{l}}({s_{n}})={E_{l}}(t)$, from (4.13) we see that by summing over j all terms except one cancel out. Hence, setting we obtain from (4.12) since $|{\alpha ^{-}}({i_{|\alpha |}})|=|\alpha |-1$ and $\alpha !\ge {\alpha ^{-}}({i_{|\alpha |}})!$. In order to interpret the last sum, consider the random variable ${I_{n-1}}({\psi _{l}}({\mathbf{s}^{n-1}};t))$. Then, following (4.7) and the identity (4.11), we conclude that Thus, we get the estimate Now, using the Cauchy–Schwarz inequality we deduce that Recall the identity ${\xi _{n}}(\cdot ;t)={(n!)^{-1}}\mathbb{E}[{D^{n}}{X_{t}}]$. Using similar arguments as in the proof of Proposition 4.1, we obtain the inequality Since $|{s_{j}}-{s_{j-1}}|\le t$, we finally conclude that Consequently, we deduce Putting things together we conclude Combining (4.9) and (4.15) we obtain the result of Theorem 3.3.

To get an idea of the performance of the polynomial chaos expansion for different bases and truncations, we apply the trigonometric basis from Example 3.4 (i) and the Haar basis from Example 3.4 (ii) to the geometric Brownian motion from Example 3.2, on the horizon $t\in [0,1]$ with $\mu =\sigma ={x_{0}}=1$. In this setting we have the identity which can be used to compare the absolute errors of approximated variances using different bases and values for k and p. We used the ode45 solver with standard settings within Matlab 2017b (64bit, single core) on an Intel Core i5-6300U (2.4GHz, 8GB RAM).