The asymptotic error of chaos expansion approximations for stochastic differential equations

TL;DR

This paper introduces a Wiener chaos expansion-based numerical scheme for stochastic differential equations, providing explicit error bounds and leveraging Malliavin calculus for analysis.

Contribution

It presents a novel chaos expansion approximation method with explicit error bounds for SDEs, advancing numerical analysis in stochastic calculus.

Findings

Derived explicit upper bounds for $L^2$ approximation error.

Applied Malliavin calculus to analyze the chaos expansion approximation.

Demonstrated the effectiveness of the method for square integrable SDEs.

Abstract

In this paper we present a numerical scheme for stochastic differential equations based upon the Wiener chaos expansion. The approximation of a square integrable stochastic differential equation is obtained by cutting off the infinite chaos expansion in chaos order and in number of basis elements. We derive an explicit upper bound for the $L^2$ approximation error associated with our method. The proofs are based upon an application of Malliavin calculus.

Tables

Table 1. Table 1: Computational times in seconds and error estimates for trigonometric and Haar bases from Example  3.4 and different values of p 𝑝 p and k 𝑘 k . The number of coefficients n coeff subscript 𝑛 coeff n_{\text{coeff}} depends also on the type of truncation, compare Table  2 for details. error t=1 is the absolute error of the expansion’s variance | V ​ a ​ r ​ [ X t ] − V ​ a ​ r ​ [ X t p , k ] | 𝑉 𝑎 𝑟 delimited-[] subscript 𝑋 𝑡 𝑉 𝑎 𝑟 delimited-[] subscript superscript 𝑋 𝑝 𝑘 𝑡 |Var[X_{t}]-Var[X^{p,k}_{t}]| at t = 1 𝑡 1 t=1 while error max max {}_{\text{max}} is the maximum value over t ∈ [ 0 , 1 ] 𝑡 0 1 t\in[0,1]

				Trigonometric basis			Haar basis
$k$	$p$	$n_{\text{coeff}}$	type	time	${\text{error}}_{t=1}$	${\text{error}}_{\text{max}}$	time	${\text{error}}_{t=1}$	${\text{error}}_{\text{max}}$
2	1	3	full	0.04	6.04	6.04	0.03	5.31	5.31
4	1	5	full	0.02	5.68	5.68	0.02	5.31	5.31
8	1	9	full	0.22	5.49	5.49	0.04	5.30	5.30
16	1	17	full	0.05	5.40	5.40	0.19	5.31	5.31
32	1	33	full	0.17	5.35	5.35	1.31	5.30	5.30
64	1	65	full	1.28	5.33	5.33	10.12	5.31	5.31
2	2	6	full	0.02	3.04	3.04	0.02	1.61	1.83
4	2	15	full	0.03	2.35	2.35	0.05	1.61	1.76
8	2	45	full	0.20	1.98	1.98	0.82	1.61	1.69
8	2	41	sp1	0.31	1.99	1.99	0.24	1.61	1.69
8	2	19	sp2	0.04	2.16	2.16	0.07	1.61	1.72
16	2	153	full	3.85	1.80	1.80	18.61	1.61	1.65
16	2	141	sp3	1.00	1.80	1.80	4.60	1.61	1.65
16	2	27	sp4	0.05	2.07	2.07	0.19	1.61	1.67
32	2	561	full	86.83	1.71	1.71	554.61	1.61	1.63
32	2	537	sp5	22.80	1.71	1.71	143.90	1.61	1.63
32	2	69	sp6	0.31	1.84	1.84	3.22	1.61	1.64
64	2	2145	full	2189.70	1.66	1.66	17234.51	1.61	1.62
2	3	10	full	0.02	2.15	2.15	0.02	0.38	1.35
4	3	35	full	0.10	1.29	1.29	0.19	0.38	1.01
8	3	165	full	2.76	0.84	0.84	10.58	0.38	0.76
8	3	127	sp7	0.53	0.85	0.85	1.77	0.37	0.76
8	3	37	sp8	0.06	1.11	1.11	0.13	0.38	0.86
16	3	969	full	200.34	0.61	0.61	1070.17	0.38	0.58
16	3	763	sp9	45.61	0.62	0.62	169.16	0.38	0.59
16	3	45	sp10	0.14	1.02	1.02	0.39	0.38	0.76
2	4	15	full	0.02	1.94	1.94	0.02	0.07	1.32
4	4	70	full	0.30	1.04	1.04	0.79	0.07	0.87
8	4	495	full	26.99	0.57	0.57	117.55	0.07	0.56
8	4	303	sp11	4.34	0.57	0.57	12.30	0.07	0.52
8	4	32	sp12	0.05	0.96	0.96	0.12	0.07	0.75
16	4	4845	full	5986.44	0.32	0.32	30591.94	0.07	0.33
16	4	40	sp13	0.08	0.87	0.87	0.36	0.07	0.67
32	4	92	sp14	0.50	0.59	0.59	3.56	0.07	0.45
2	5	21	full	0.03	1.91	1.91	0.03	0.01	1.27
4	5	126	full	0.92	1.00	1.00	2.04	0.01	0.87
8	5	1287	full	192.98	0.51	0.51	855.07	0.01	0.53
8	5	599	sp15	10.04	0.51	0.51	50.03	0.01	0.55
8	5	36	sp16	0.03	0.92	0.92	0.11	0.01	0.74
16	5	20349	full	120469.96	0.26	0.26	600591.17	0.01	0.28
16	5	44	sp17	0.06	0.83	0.83	0.28	0.01	0.65
32	5	98	sp18	0.51	0.55	0.55	3.60	0.01	0.42

Table 2. Table 2: List of first and second order sparse indices used in Section  5.2 together with the number of resulting coefficient functions. The reference numbers coincide with those in Table  1

symbol	$p$	$k$	order	index $\text{mathbh}r$/$\left(\text{mathbh}r\right)$	$n_{\text{coeff}}$
sp1	2	8	1	$\text{mathbh}r=(2,2,2,2,1,1,1,1)$	41
sp2	2	8	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	19
				$\text{mathbh}{r}^2=(2,2,2,2,0,0,0,0)$
sp3	2	16	1	$\text{mathbh}r=(2,2,2,2,1,\dots ,1)$	141
sp4	2	16	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	27
				$\text{mathbh}{r}^2=(2,2,2,2,0,\dots ,0)$
sp5	2	32	1	$\text{mathbh}r=(2,2,2,2,2,2,2,2,1,\dots ,1)$	537
sp6	2	32	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	69
				$\text{mathbh}{r}^2=(2,2,2,2,2,2,2,2,0,\dots ,0)$
sp7	3	8	1	$\text{mathbh}r=(3,3,2,2,1,1,1,1)$	127
sp8	3	8	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	37
				$\text{mathbh}{r}^2=(2,2,2,2,0,0,0,0)$
				$\text{mathbh}{r}^3=(3,3,2,2,0,0,0,0)$
sp9	3	16	1	$\text{mathbh}r=(3,3,2,2,1,\dots ,1)$	763
sp10	3	16	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	45
				$\text{mathbh}{r}^2=(2,2,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^3=(3,3,2,2,0,\dots ,0)$
sp11	4	8	1	$\text{mathbh}r=(4,4,2,2,1,1,1,1)$	303
sp12	4	8	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	32
				$\text{mathbh}{r}^2=(2,2,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^3=(3,3,2,,0,\dots ,0)$
				$\text{mathbh}{r}^4=(4,3,0,\dots ,0)$
sp13	4	16	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	40
				$\text{mathbh}{r}^2=(2,2,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^3=(3,3,2,0,\dots ,0)$
				$\text{mathbh}{r}^4=(4,3,0,\dots ,0)$
sp14	4	32	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	92
				$\text{mathbh}{r}^2=(2,2,2,2,2,2,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^3=(3,3,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^4=(4,4,0,\dots ,0)$
sp15	5	8	1	$\text{mathbh}r=(5,5,2,2,1,1,1,1)$	599
sp16	5	8	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	36
				$\text{mathbh}{r}^2=(2,2,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^3=(3,3,2,,0,\dots ,0)$
				$\text{mathbh}{r}^4=(4,3,0,\dots ,0)$
				$\text{mathbh}{r}^5=(5,3,0,\dots ,0)$
sp17	5	16	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	44
				$\text{mathbh}{r}^2=(2,2,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^3=(3,3,2,0,\dots ,0)$
				$\text{mathbh}{r}^4=(4,3,0,\dots ,0)$
				$\text{mathbh}{r}^5=(5,3,0,\dots ,0)$
sp18	5	32	2	$\text{mathbh}{r}^1=(1,\dots ,1)$	98
				$\text{mathbh}{r}^2=(2,2,2,2,2,2,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^3=(3,3,2,2,0,\dots ,0)$
				$\text{mathbh}{r}^4=(4,4,0,\dots ,0)$
				$\text{mathbh}{r}^5=(5,5,0,\dots ,0)$

Equations

$$d{X}_t=b(t,X_t)dt+\sigma (t,X_t)d{W}_t\text{with}{X}_0=x_0,$$

$$\displaystyle dX_{t}^{n}=b\bigl{(}\varphi_{n}(t),X_{\varphi_{n}(t)}^{n}\bigr{)}dt+\sigma\bigl{(}\varphi_{n}(t),X_{\varphi_{n}(t)}^{n}\bigr{)}dW_{t}\quad\text{with}\ X_{0}^{n}=x_{0}.$$

$$X_t=\sum _{\alpha \in \mathcal{I}}{x}_{\alpha }(t){\Psi }^{\alpha },$$

$$\displaystyle\mathcal{I}:=\bigl{\{}\alpha=(\alpha_{i})_{i\geq 1}:~{}\alpha_{i}\in\mathbh{N}_{0}\text{ and almost all $\alpha_{i}$'s are }0\bigr{\}}.$$

$$X_t^{p,k}=\sum _{\alpha \in {\mathcal{I}}_{p,k}}{x}_{\alpha }(t){\Psi }^{\alpha },$$

$$\mathbh{E}\bigl{[}W(g)W(h)\bigr{]}=\langle g,h\rangle_{\mathbh{H}}.$$

$$\displaystyle H_{0}(x):=1,\qquad H_{n}(x):=\frac{(-1)^{n}}{\sqrt{n!}}\exp\biggl{(}\frac{x^{2}}{2}\biggr{)}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}\biggl{(}\exp\biggl{(}-\frac{x^{2}}{2}\biggr{)}\biggr{)},\quad n\geq 1.$$

$$L^2(\Omega ,\mathcal{F},\text{mathbh}P)=\underset{n=0}{\overset{\infty }{⨁}}{\mathcal{H}}_n.$$

$$\displaystyle\varPsi^{\alpha}:=\prod_{i=1}^{\infty}H_{\alpha_{i}}\bigl{(}W(e_{i})\bigr{)},$$

$$\displaystyle I_{n}\bigl{(}h^{\otimes n}\bigr{)}:=\sqrt{n!}~{}H_{n}\bigl{(}W(h)\bigr{)},\quad n\geq 1.$$

$$I_n(h):=I_n(\tilde{h}),$$

$$\displaystyle\mathbh{E}\bigl{[}I_{m}(g)I_{n}(h)\bigr{]}=n!\langle\widetilde{g},\widetilde{h}\rangle_{\mathbh{H}^{\otimes n}}1_{\{n=m\}},\quad h\in\mathbh{H}^{\otimes n},g\in\mathbh{H}^{\otimes m}.$$

$$I_m(g)I_n(h)=\sum _{r=0}^{\min (m,n)}r!\left(\genfrac{}{}{0ex}{}{m}{r}\right)\left(\genfrac{}{}{0ex}{}{n}{r}\right)I_{m+n-2r}(g{\otimes }_{r}h),$$

$$g{\otimes }_{r}h(t_1,\dots ,t_{m+n-2r})$$

$$\times h(u_1,\dots ,u_r,t_{m-r+1},\dots t_{m+n-2r})d{u}_1\dots d{u}_r.$$

$$F=\sum _{n=0}^{\infty }{I}_n(g_n),g_n\in \text{mathbh}{H}^{\otimes n},$$

$$\mathcal{S}=\bigl{\{}F=f\bigl{(}W(h_{1}),\ldots,W(h_{n})\bigr{)}:~{}n\geq 1,h_{i}\in\mathbh{H}\bigr{\}},$$

$$\displaystyle D^{k}F=\sum_{i_{1},\ldots,i_{k}=1}^{n}\frac{\partial^{k}}{\partial x_{i_{1}}\cdots\partial x_{i_{k}}}f\bigl{(}W(h_{1}),\ldots,W(h_{n})\bigr{)}h_{i_{1}}\otimes\cdots\otimes h_{i_{k}}.$$

$$\|F\|_{k,q}:=\Biggl{(}\mathbh{E}\bigl{[}|F|^{q}\bigr{]}+\sum_{m=1}^{k}\mathbh{E}\bigl{[}\bigl{\|}D^{m}F\bigr{\|}_{\mathbh{H}^{\otimes m}}^{q}\bigr{]}\Biggr{)}^{1/q}.$$

$$D_s{\Psi }^{\alpha }$$

$${\alpha }^{-}(i)$$

$$\displaystyle\mathbh{E}\bigl{[}F\delta^{k}(u)\bigr{]}=\mathbh{E}\bigl{[}\bigl{\langle}D^{k}F,u\bigr{\rangle}_{\mathbh{H}^{\otimes k}}\bigr{]}.$$

$$\displaystyle D_{t}F=\sum_{n=1}^{\infty}nI_{n-1}\bigl{(}g_{n}(\cdot,t)\bigr{)},\quad t\in[0,T],$$

$$D_s{X}_t=\sigma (s,X_s)+{\int }_s^t\frac{\partial }{\partial x}b(u,X_u)D_s(X_u)du+{\int }_s^t\frac{\partial }{\partial x}\sigma (u,X_u)D_s(X_u)d{W}_u,$$

$$\displaystyle X_{t}=\sum_{\alpha\in\mathcal{I}}x_{\alpha}(t)\varPsi^{\alpha}\quad\text{with}\ \varPsi^{\alpha}=\prod_{i=1}^{\infty}H_{\alpha_{i}}\bigl{(}W(e_{i})\bigr{)}.$$

$$\displaystyle x_{\alpha}(t)=\mathbh{E}\bigl{[}X_{t}\varPsi^{\alpha}\bigr{]}=\frac{1}{\sqrt{\alpha!}}\mathbh{E}\Biggl{[}\Biggl{\langle}D^{|\alpha|}X_{t},\bigotimes_{i=1}^{\infty}e_{i}^{\otimes\alpha_{i}}\Biggr{\rangle}_{\mathbh{H}^{\otimes|\alpha|}}\Biggr{]}$$

$$x_{\alpha }^{\prime }(t)$$

$$x_{\alpha }(0)$$

$$\displaystyle x_{\alpha}(t)=x_{0}1_{\{\alpha=0\}}+\int_{0}^{t}\mathbh{E}\bigl{[}\varPsi^{\alpha}b(s,X_{s})\bigr{]}ds+\mathbh{E}\Biggl{[}\varPsi^{\alpha}\int_{0}^{t}\sigma(s,X_{s})dW_{s}\Biggr{]}.$$

$$b_{\alpha}(s,X_{s})=\mathbh{E}\bigl{[}\varPsi^{\alpha}b(s,X_{s})\bigr{]}.$$

$$\displaystyle\mathbh{E}\Biggl{[}\varPsi^{\alpha}\int_{0}^{t}\sigma(s,X_{s})dW_{s}\Biggr{]}$$

$$=\sum _{j=1}^{\infty }{\int }_0^t\sqrt{{\alpha }_j}{e}_j(s){\sigma }_{{\alpha }^{-}(j)}(s,X_s)ds.$$

Full text

Chaos expansion,

Malliavin calculus,

numerical approximation,

stochastic differential equations,

65C30,

60H10,

60H07,

††volume: 6††issue: 2††articletype: research-article

\aid

VMSTA133 \startlocaldefs

\urlstylerm \endlocaldefs

\pretitle

Research Article

[type=corresp,id=cor1]Corresponding author.

\publishedonline\sday

23 4 2019

1 Introduction

We consider a one-dimensional continuous stochastic process $(X_{t})_{t\in[0,T]}$ that satisfies the stochastic differential equation

[TABLE]

where $(W_{t})_{t\in[0,T]}$ is a Brownian motion defined on a filtered probability space $(\varOmega,\mathcal{F},\allowbreak(\mathcal{F}_{t})_{t\in[0,T]},\mathbh{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}:\mathbh{R}_{+}\to\mathbh{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

[TABLE]

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\geq 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},\mathbh{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

[TABLE]

where $x_{\alpha}$ are deterministic functions, $\varPsi^{\alpha}$ are mutually orthogonal random projections that are associated to the basis $(e_{i})_{i\geq 1}$, and the index set $\mathcal{I}$ is defined via

[TABLE]

Such orthogonal expansions have been successfully applied in numerous fields of stochastic and numerical analysis. We refer e.g. to [11, 12, 13, 25] for applications of the Wiener chaos expansion in the context of SDEs and to [20, 21, 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

[TABLE]

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\leq i\leq k}$ and only up to the $p$th 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\geq 1}$. In this paper we will study the decay rate of $\mathbh{E}[(X_{t}-X_{t}^{p,k})^{2}]$ when $k,p\to\infty$ for a general basis $(e_{i})_{i\geq 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.

2 Background on Malliavin calculus

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 $\mathbh{H}=L^{2}([0,T])$ and let $\langle\cdot,\cdot\rangle_{\mathbh{H}}$ denote the scalar product on $\mathbh{H}$. We note that $\mathbh{H}$ is a separable Hilbert space and denote by $(e_{i})_{i\geq 1}$ an orthonormal basis of $\mathbh{H}$. We consider the isonormal Gaussian family $W=\{W(h):h\in\mathbh{H}\}$ indexed by $\mathbh{H}$ defined on a probability space $(\varOmega,\mathcal{F},\mathbh{P})$, i.e. the random variables $W(h)$ are centered Gaussian with a covariance structure determined via

[TABLE]

Here and throughout the paper we assume that $\mathcal{F}=\sigma(W)$. In our setting we consider $W(h)=\int_{0}^{T}h_{s}dW_{s}$, where $W$ is a standard Brownian motion. We define the normalised Hermite polynomials through the identities

[TABLE]

The $n$th Wiener chaos $\mathcal{H}_{n}$ is the closed linear subspace of $L^{2}(\varOmega,\mathcal{F},\mathbh{P})$ generated by the family of random variables $\{H_{n}(W(h)):~{}\|h\|_{\mathbh{H}}=1\}$. The vector spaces $\mathcal{H}_{n}$, $n\geq 0$, are orthogonal and we have the Wiener chaos expansion

[TABLE]

(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

[TABLE]

where $(e_{i})_{i\geq 1}$ is a fixed orthonormal basis of $\mathbh{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}:~{}\alpha\in\mathcal{I}\text{ with }|\alpha|=n\}$ forms a complete orthonormal basis of the $n$th Wiener chaos $\mathcal{H}_{n}$ and consequently $\{\varPsi^{\alpha}:~{}\alpha\in\mathcal{I}\}$ is a complete orthonormal basis of $L^{2}(\varOmega,\mathcal{F},\mathbh{P})$ (cf. [24, Proposition 1.1.1]).

Now, we introduce multiple stochastic integrals of order $n\in\mathbh{N}$, which are denoted by $I_{n}$. For an element $h^{\otimes n}:=h\otimes\cdots\otimes h$ in $\mathbh{H}^{\otimes n}$ with $\|h\|_{\mathbh{H}}=1$, we define

[TABLE]

Assuming that the mapping $I_{n}$ is linear, the definition (2.4) can be extended to all symmetric elements $h\in\mathbh{H}^{\otimes n}$ by polarisation identity. Finally, for an arbitrary function $h\in\mathbh{H}^{\otimes n}$ we set

[TABLE]

where $\widetilde{h}$ denotes the symmetrised version of $h\in\mathbh{H}^{\otimes n}$. By definition $I_{n}$ maps $\mathbh{H}$ into $\mathcal{H}_{n}$, so the multiple integrals of different orders are orthogonal. In particular, they satisfy the isometry property

[TABLE]

Furthermore, for any symmetric $h\in\mathbh{H}^{\otimes n},g\in\mathbh{H}^{\otimes m}$, the following multiplication formula holds:

[TABLE]

where the $r$th contraction $g\otimes_{r}h$ is defined by

[TABLE]

(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},\mathbh{P})$ we obtain the orthogonal decomposition

[TABLE]

where $I_{0}=\mathbh{E}[F]$ and the decomposition is unique when $g_{n}$, $n\geq 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

[TABLE]

where $f\in C_{p}^{\infty}(\mathbh{R}^{n})$ (i.e. the space of infinitely differentiable functions such that all derivatives exhibit polynomial growth). The $k$th order Malliavin derivative of $F\in\mathcal{S}$, denoted by $D^{k}F$, is defined by

[TABLE]

Notice that $D^{k}F$ is a $\mathbh{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 $\mathbh{D}^{k,q}$ denotes the completion of the set $\mathcal{S}$ with respect to the norm

[TABLE]

We define $\mathbh{D}^{k,\infty}=\cap_{q>1}\mathbh{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)$:

[TABLE]

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\mathbh{D}^{k,2}$, where $\text{Dom}(\delta^{k})$ consists of all elements $u\in L^{2}(\varOmega;\mathbh{H}^{\otimes k})$ such that the inequality $|\mathbh{E}[\langle D^{k}F,u\rangle_{\mathbh{H}^{\otimes k}}]|\leq c(E[F^{2}])^{1/2}$ holds for some $c>0$ and all $F\in\mathbh{D}^{k,2}$, then we have the identity

[TABLE]

In case the random variable $F\,{\in}\,L^{2}(\varOmega,\mathcal{F},\mathbh{P})$ has a chaos decomposition as displayed in (2.7), the statement $F\in\mathbh{D}^{k,2}$ is equivalent to the condition $\sum_{n=1}^{\infty}n^{k}n!\|g_{n}\|^{2}_{\mathbh{H}^{\otimes n}}<\infty$. In particular, when $F\in\mathbh{D}^{1,2}$ we deduce an explicit chaos representation of the derivative $DF=(D_{t}F)_{t\in[0,T]}$

[TABLE]

where $g_{n}(\cdot,t):[0,T]^{n-1}\to\mathbh{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}(\mathbh{R})$. In this setting $DX_{t}=(D_{s}X_{t})_{s\in[0,T]}$ is given as the solution of the SDE

[TABLE]

for $s\leq t$, and $D_{s}X_{t}=0$ if $s>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.

3 Main results

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

[TABLE]

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

[TABLE]

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.

Theorem 3.1

Let $(X_{t})_{t\in[0,T]}$ be the unique solution of the SDE (1.1) and assume that $X\in L^{2}(\varOmega\times[0,T])$. Then $X_{t}$ possesses the chaos expansion (1.3) and the coefficients $x_{\alpha}(t)$ satisfy the system of ordinary differential equations

[TABLE]

Here $b_{\alpha}(t,X_{t})$ (resp. $\sigma_{\alpha}(t,X_{t})$) denotes the $\alpha$-coefficient of the Wiener chaos expansion (1.3) associated with the random variable $b(t,X_{t})$ (resp. $\sigma(t,X_{t})$), and the multi-index $\alpha^{-}(j)$ is defined by (2.10).

Proof.

Using the SDE (1.1) and applying the formula (3.1) we obtain the identity

[TABLE]

Applying the formula (3.1) once again for the random variable $b(s,X_{s})$ we immediately deduce that

[TABLE]

On the other hand, observing the identity $\delta(1_{[0,t]}\sigma(\cdot,X_{\cdot}))=\int_{0}^{t}\sigma(s,X_{s})dW_{s}$, we get by the duality formula (2.11) and (2.9) that

[TABLE]

Putting things together we obtain the identity

[TABLE]

Consequently, the assertion follows after taking the derivative with respect to $t$. ∎

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.

Example 3.2**.**

(i) (Scaled Brownian motion with drift) We start with the toy example of a scaled Brownian motion with drift:

[TABLE]

In this case we obviously have that $x_{\alpha}(t)=0$ for any $\alpha\in\mathcal{I}$ with $|\alpha|\geq 2$. Applying formula (3.2) we obtain the representation

[TABLE]

which is a well known Karhunen–Loéve expansion of the scaled Brownian motion.

(ii) (Geometric Brownian motion) Let us consider the geometric Brownian motion defined via the SDE

[TABLE]

In this setting the propagator system (3.2) translates to

[TABLE]

This system of ordinary differential equations can be solved recursively. For $\alpha=0$ we have $x^{\prime}_{0}(t)=bx_{0}(t)$ and hence $x_{0}(t)=x_{0}\exp(bx)$. If $\alpha$ is the $j$th canonical unit vector in $\mathcal{I}$ (and hence $|\alpha|=1$) we obtain the differential equation

[TABLE]

Hence, $x_{\alpha}(t)=x_{0}\sigma\exp(bx)\int_{0}^{t}e_{j}(s)ds$. Following this recursion we obtain the general formula

[TABLE]

for any $\alpha\in\mathcal{I}$ with $|\alpha|=p$.

For a general specification of the drift coefficient $b$ and diffusion coefficient $\sigma$ 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

[TABLE]

The approximation of $X_{t}$ is now defined via (1.5):

[TABLE]

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

[TABLE]

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\geq 1}$ of $L^{2}([0,T])$ whereas we are interested in the asymptotic analysis for a general basis $(e_{i})_{i\geq 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 $\mathbh{E}[(X_{t}^{p,k}-X_{t})^{2}]$.

Theorem 3.3

Let $(X_{t})_{t\in[0,T]}$ be the solution of the SDE (1.1). Suppose that the coefficient functions $b$ and $\sigma$ satisfy the Lipschitz and linear growth conditions

[TABLE]

Moreover, assume that $b,\sigma\in C^{1,\infty}([0,T]\times\mathbh{R})$, where $C^{1,\infty}([0,T]\times\mathbh{R})$ denotes the space of functions $f:[0,T]\times\mathbh{R}\to\mathbh{R}$ that are once continuously differentiable in the first component and infinitely differentiable in the second component, such that

[TABLE]

for $t\in[0,T]$, any $l=0,1$, and $m\geq 0$. Then it holds that

[TABLE]

where $C=C(t,K)$ is a positive constant and the function $E_{l}(t)$ is defined by

[TABLE]

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)dW_{s}$, we readily deduce that

[TABLE]

Hence, the upper bound on the right-hand side of (3.6) indeed converges to [math] 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\geq 1}$, while the error associated with truncation of basis strongly depends on $(e_{i})_{i\geq 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\leq 4$).

Example 3.4**.**

Let us explicitly compute the last terms of the upper bound (3.6) for two prominent bases of $L^{2}([0,1])$.

(i) (Trigonometric basis) Consider the orthonormal basis $(e_{i})_{i\geq 1}$ given by

[TABLE]

In this setting we obtain that

[TABLE]

for any $j\geq 1$. Consequently, we deduce that

[TABLE]

for some $C>0$.

(ii) (Haar basis) The Haar basis is a collection of functions

[TABLE]

defined as follows:

[TABLE]

In this case we deduce the following representation for $E_{j,n}(t)$:

[TABLE]

Since $\max_{t\in[0,1]}E_{j,n}^{2}(t)=2^{-(n+1)}$ and $E_{j,n}(t)\neq 0$ only for $t\in[2^{-n+1}(j-1),2^{-n+1}j]$, we finally obtain that

[TABLE]

Note that the exponential asymptotic rate for a fixed $p$th order Wiener chaos is equivalent to $\mathcal{O}(k^{-1})$, if we choose $k=2^{n}$ basis elements $(e_{i})$. 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.

4 Proofs

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}$.

Proposition 4.1

Under conditions of Theorem 3.3 we obtain the estimate

[TABLE]

Proof.

We show the assertion of Proposition 4.1 by induction over $n$. For $n=0$ the result is a well-known consequence of the Lipschitz and linear growth conditions (3.5) on the functions $b$ and $\sigma$ (see e.g. [15, Theorem 2.9, page 289] which relates the second moment of $X_{t}$ to the second moment of the initial condition $X_{0}=x_{0}$). For $n=1$ we known from the formula (2.13) that the process $(D_{s}X_{t})_{s\in[0,T]}$ satisfies the SDE

[TABLE]

for $s\leq t$, and $D_{s}X_{t}=0$ for $s>t$. Next, we introduce the two-dimensional stochastic process $(Y^{(1)}_{s;t})_{t\geq s}$ via

[TABLE]

Note that the process $(Y^{(1)}_{s;t})_{t\geq s}$ satisfies a two-dimensional SDE with initial condition given by $Y^{(1)}_{s;s}$. In particular, setting $Y^{(1)}_{s;t}=(Y^{(1.1)}_{s;t},Y^{(1.2)}_{s;t})$, we have the representation

[TABLE]

where $\overline{b}(u,y)=\frac{\partial}{\partial x}b(u,y_{1})y_{2}$ and $\overline{\sigma}(u,y)=\frac{\partial}{\partial x}\sigma(u,y_{1})y_{2}$ for $y=(y_{1},y_{2})$. Applying again the estimates of [15, Theorem 2.9, page 289] to the diffusion process $(Y^{(1)}_{s;t})_{t\geq s}$, we conclude the estimate

[TABLE]

Now, we will perform the induction step. Notice that $D_{s}X_{t}\in\mathbh{D}^{1,\infty}$, so its Malliavin derivative satisfies again an SDE according to [24, Theorem 2.2.2]. We define the stochastic process

[TABLE]

Assume that the assertion of Theorem 3.3 holds for all $k=1,\ldots,n-1$. In order to compute the estimate, we have to consider the initial values of the SDE system satisfied by $Y_{s_{1},\ldots,s_{n};\,t}^{(n)}$. The formula below can be found in [24, Theorem 2.2.2]. Before we state it, we need some additional notation. The stochastic process D^{m}X_{t}=\{D^{m}_{s_{1},\ldots,s_{m}}X_{t}\bigm{|}(s_{1},\ldots,s_{m})\in[0,T]\} depends on the $m$ time points $s_{1},\ldots,s_{m}$. For any subset $J=\{j_{1}<\cdots<j_{\eta}\}$ of $\{1,\ldots,m\}$ with $|J|=\eta\leq m$ elements, denote $s(J)=(s_{j_{1}},\ldots,s_{j_{\eta}})$. Further on, we define

[TABLE]

and

[TABLE]

where the sums run over the set $\mathfrak{P}^{m}$ of all partitions $J_{1}\cup\cdots\cup J_{\nu}$ of $\{1,\ldots,m\}$, where $J_{1},\ldots,J_{\nu}$ are disjoint sets. We determine $\mathfrak{z}(t)=\sigma(t,X_{t})$ as well. With these notations at hand, we find by (2.13) and induction that the $n$th order Malliavin derivative $D^{n}_{s_{1},\ldots,s_{n}}X_{t}$ satisfies the linear SDE

[TABLE]

for $\hat{s}:=\max\{s_{1},\ldots,s_{n}\}\leq t$ and $D^{n}_{s_{1},\ldots,s_{n}}X_{t}=0$ else. Hence, its initial value is given by

[TABLE]

where

[TABLE]

Finally, we apply [15, Theorem 2.9, page 289] as for the case $n=1$ and taking into account that the assertion holds for $k=1,\ldots,n-1$:

[TABLE]

This completes the proof of Proposition 4.1. ∎

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

[TABLE]

where $d(\alpha):=\max\{i\geq 1:~{}\alpha_{i}>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.,

[TABLE]

with $\mathbf{t}^{n}=(t_{1},\ldots,t_{n})$ and symmetric kernel functions $\xi_{n}(\mathbf{t}^{n};t)$ being defined by

[TABLE]

for all $n\in\mathbh{N}_{0}$. We conclude that

[TABLE]

where we abbreviate

[TABLE]

Therefore, we deduce via (4) that

[TABLE]

Finally, using (4.8) and applying Proposition 4.1 we obtain

[TABLE]

The treatment of the term $A_{2}(p,k)$ is more involved. Recalling the definition of $d(\alpha)=\max\{i\geq 1:~{}\alpha_{i}>0\}$, we introduce the characteristic set $(i_{1},\ldots,i_{n})$ of $\alpha\in\mathcal{I}$ for $i_{1}\leq i_{2}\leq\cdots\leq i_{n}$ and $|\alpha|=n$. It is defined as follows: $i_{1}$ is the index of the first non-zero component of $\alpha$. If $\alpha_{i_{1}}=1$ then $i_{2}$ is the index of the second non-zero component of $\alpha$; otherwise $i_{1}=i_{2}=\cdots=i_{\alpha_{i_{1}}}$ and $i_{\alpha_{i_{1}}+1}$ is the second non-zero component of $\alpha$. The same operation is repeated for the index $i_{\alpha_{i_{1}}+1}$ and further non-zero components of $\alpha$. 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,\ldots)\in\mathcal{I}$. Here the characteristic set of $\alpha$ is given by

[TABLE]

For any $\alpha\in\mathcal{I}$ with $|\alpha|=n$ and a basis $(e_{i})_{i\geq 1}$ of $\mathbh{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

[TABLE]

where $(i_{1},\ldots,i_{n})$ is the characteristic set of $\alpha$ and the sum runs over all permutations $\pi$ within the permutation group $\mathfrak{P}^{n}$ of $\{1,\ldots,n\}$. From (3.1) we know that $x_{\alpha}(t)=\mathbh{E}[X_{t}\varPsi^{\alpha}]$. On the other hand, we have the representation

[TABLE]

Now, for any $\alpha\in\mathcal{I}$ with $|\alpha|=n$, we obtain the identity

[TABLE]

by (2.11). Since

[TABLE]

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

[TABLE]

with $t_{0}:=0$ and $t_{n+1}:=t$, by changing the order of integration. Then for any $i_{n}=l\geq 1$ we integrate by parts to deduce

[TABLE]

where

[TABLE]

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\leq j-1$ and $s_{i}=t_{i+1}$ for all $i>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,\ldots,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},\ldots,s_{n-1})$ and $\mathbf{s}^{n-1,n}=(s_{1},\ldots,s_{n-1},s_{n})$. Then

[TABLE]

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

[TABLE]

we obtain from (4.12)

[TABLE]

since $|{\alpha^{-}(i_{|\alpha|}})|=|\alpha|-1$ and $\alpha!\geq\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) and the identity (4.11), we conclude that

[TABLE]

Thus, we get the estimate

[TABLE]

Now, using the Cauchy–Schwarz inequality we deduce that

[TABLE]

Recall the identity $\xi_{n}(\cdot;t)=(n!)^{-1}\mathbh{E}[D^{n}X_{t}]$. Using similar arguments as in the proof of Proposition 4.1, we obtain the inequality

[TABLE]

Since $|s_{j}-s_{j-1}|\leq t$, we finally conclude that

[TABLE]

Consequently, we deduce

[TABLE]

Putting things together we conclude

[TABLE]

Combining (4) and (4.15) we obtain the result of Theorem 3.3.

5 Numerical results and sparse truncation

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,

[TABLE]

on the horizon $t\in[0,1]$ with $\mu=\sigma=x_{0}=1$. In this setting we have the identity

[TABLE]

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).

5.1 Absolute errors of variances

Figure 1 shows the absolute errors between the variances of the analytical solution and different polynomial chaos expansions. One observes the fast decay of the approximation error for fixed $k$ and increasing $p$, and vice versa. The main difference between the observed errors is that the maximum of the trigonometric expansion is always at the end of the time horizon at $t=1$, while the Haar basis leads to almost zero error on $k$ equidistant time points, including $t=0$ and $t=1$.

Table 1 in the Appendix lists approximation errors, numbers of coefficients, and computational times for several choices of $p$ and $k$. The exponential increase in the number of coefficients with increasing $p$ and $k$ leads to an exponential increase in runtime, as expected. To overcome this issue, we propose a heuristic to reduce the computational burden.

5.2 Sparse truncation

The information contained in the coefficient functions $x_{\alpha}(\cdot)$ decays with increasing order $p$ of the basis polynomials $\varPsi^{\alpha}$ and the decaying rate of the Gaussian expansion, i.e., the index $k$ of the basis polynomials $e_{k}$ used for constructing $\varPsi^{\alpha}$ [22].

Hence, we can define sparse index sets for reducing the number of multi-indices $\alpha$ from the full truncated index set $\mathcal{I}_{p,k}$ without losing too much information, i.e., accuracy in a numerical simulation.

Definition 1** (Sparse truncation of first order).**

Let $p$ be the maximum order of the index $\alpha\in\mathcal{I}_{p,k}$. Then the first order sparse index $\mathbh{r}=(r_{1},\ldots,r_{k})$ satisfies $p=r_{1}\geq r_{2}\geq\cdots\geq r_{k}$ and we define the first order sparse index set

[TABLE]

Example 5.1**.**

Let $k=5$ and $p=3$. Then a possible choice of the sparse index is $\mathbh{r}=(3,2,2,1,1)$. For constructing the first order polynomials all five first order basis polynomials can be used. The second order polynomials are comprised by all possible combinations of first order basis polynomials depending on $e_{1},\ldots,e_{5}$ and the second order basis polynomials of $e_{1},e_{2},e_{3}$. Analogously, the third order polynomials are constructed. By using this first order sparse index set $\mathcal{I}_{p,k}^{\mathbh{r}}$ the number of coefficient functions $x_{\alpha}(\cdot)$ appearing within the propagator system (3.2) can be reduced drastically without impairing the solution much. The full index set $\mathcal{I}_{p,k}$ consists of $\frac{(k+p)!}{k!p!}=56$ terms, whereas this first order sparse index set includes $42$ terms.

An even more rigorous reduction of the number of coefficients included in the propagator system can be achieved by using a second order sparse index, i.e., a series of first order sparse indices $(\mathbh{r}^{j})_{j=0,\ldots,p}$ that depend on the actual order of the polynomials $\varPsi^{\alpha}$ with $j=|\alpha|$, [22]. This approach allows to exclude crossing products of random variables $W(e_{i})$ from the construction of higher order basis polynomials $\varPsi^{\alpha}$ that add only negligible information to the system.

Definition 2** (Sparse truncation of second order).**

Let $p$ be the maximum order of the index $\alpha\in\mathcal{I}_{p,k}$. Then the second order sparse index $(\mathbh{r})=(\mathbh{r}^{j})_{j\leq p}$ is a series of first order sparse indices $\mathbh{r}^{j}=(r_{1}^{j},\ldots,r_{k}^{j})$ satisfying $j=r_{1}^{j}\geq r_{2}^{j}\geq\cdots\geq r_{k}^{j}$ for all $j=|\alpha|\leq p$ and we define the second order sparse index set

[TABLE]

Example 5.2**.**

Considering again the setting of the previous example with $k=5$ and $p=3$, one possible choice of a second order sparse index is given via $\mathbh{r}^{1}=(1,1,1,1,1)$, $\mathbh{r}^{2}=(2,2,2,1,0)$, and $\mathbh{r}^{3}=(3,2,0,0,0)$. In constructing basis polynomials of order $|\alpha|=3$ we can use all combinations of basis polynomials depending on the first two random variables $e_{1}$ and $e_{2}$ up to orders $3$ and $2$, respectively. Thus, these are $\sqrt{6}H_{3}(e_{1})$, $\sqrt{2}H_{2}(e_{1})H_{1}(e_{2})$, and $\sqrt{2}H_{1}(e_{1})H_{2}(e_{2})$, compare [22].

Table 2 in the Appendix lists the first and second order sparse indices that were used for the numerical study. The results in Table 1 show that at the price of a slightly higher error, and of course the loss of a guaranteed upper bound as in (3.6), the computational times could be reduced by several orders of magnitude.

Appendix

{funding}

This project has received funding from the \gsponsor[id=GS1,sponsor-id=100010663]European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No \gnumber[refid=GS1]647573), from \gsponsor[id=GS4,sponsor-id=501100001659]Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – \gnumber[refid=GS4]314838170, GRK 2297 MathCoRe, from the project “Ambit fields: probabilistic properties and statistical inference” funded by \gsponsor[id=GS7,sponsor-id=100008398]Villum Fonden, and from CREATES funded by the \gsponsor[id=GS8,sponsor-id=501100001732]Danish National Research Foundation.