Deterministic Quadrature Formulas for SDEs Based on Simplified Weak Itô–Taylor Steps

Abstract

We study the problem of approximating the expected value ${\mathbb E}f(X(1))$ of a function f of the solution X of a d-dimensional system of stochastic differential equations (SDE) at time point 1 based on finitely many evaluations of the coefficients of the SDE, the integrand f and their derivatives. We present a deterministic algorithm, which produces a quadrature rule by iteratively applying a Markov transition based on the distribution of a simplified weak Itô–Taylor step together with strategies to reduce the diameter and the size of the support of a discrete measure. We essentially assume that the coefficients of the SDE are s-times continuously differentiable and the diffusion coefficient satisfies a uniform non-degeneracy condition and that the integrand f is r-times continuously differentiable. In the case $r \le (\lfloor s/2 \rfloor - 1) \cdot 2d/(d + 2)$, we almost achieve an error of order $\min (r, s)/d$ in terms of the computational cost, which is optimal in a worst-case sense.

Appendix

In this section, we present the sequential support point elimination procedure R of Davis [2] that constitutes the main ingredient of the support cardinality reduction operator $R_{\tau }$ constructed in Sect. 4.4, and we discuss its computational cost. Furthermore, we provide results on smoothness properties of the semigroup of the linear operators $P_t:F^1\rightarrow F^1$ associated with equations $(x_0,a,b)\in {\mathcal C}_{K,\lambda }^{s,s}$, see (22).

1.1 A Sequential Support Point Elimination Procedure

Fix $k\in {\mathbb N}$ and k functions

$$\begin{aligned} h_1,\ldots ,h_k:{\mathbb R}^d\rightarrow {\mathbb R}. \end{aligned}$$

Let $N\in {\mathbb N}$ with $N>k$ and consider a measure

$$\begin{aligned} \mu = \sum _{i=1}^{N} w_i\cdot \delta _{x_i} \in {\mathcal M}({\mathbb R}^d) \end{aligned}$$

with $|{\text {supp}}(\mu )|=N$. Choose any $u\in {\mathbb R}^{k+1}{\setminus } \{0\}$ that satisfies

$$\begin{aligned} \sum _{i=1}^{k+1} u_i\cdot h_j(x_i) = 0,\quad j=1,\ldots ,k, \end{aligned}$$

(39)

and $\{i:u_i< 0 \}\ne \emptyset$, and put

$$\begin{aligned} \alpha = \min _{i: \,u_i <0}\left( -\frac{w_i}{u_i}\right) . \end{aligned}$$

Let

$$\begin{aligned} \widetilde{w}_i = {\left\{ \begin{array}{ll} w_i + \alpha \cdot u_i, &{}\quad \text {if } i\in \{1,\ldots ,k+1\},\\ w_i, &{}\quad \text {if }i\in \{k+2,\ldots ,N\}, \end{array}\right. } \end{aligned}$$

and define

$$\begin{aligned} \widetilde{\mu }= \sum _{i=1}^{N} \widetilde{w}_i\cdot \delta _{x_i}. \end{aligned}$$

Clearly, $\widetilde{w}_1,\ldots , \widetilde{w}_{N} \ge 0$ and $\widetilde{w}_i=0$ for some $i\in \{1,\ldots ,k+1\}$. Hence $\widetilde{\mu }\in {\mathcal M}({\mathbb R}^d)$ and

$$\begin{aligned} {\text {supp}}(\widetilde{\mu })\subset {\text {supp}}(\mu ),\, |{\text {supp}}(\widetilde{\mu })|\le N-1. \end{aligned}$$

Furthermore, by (39),

$$\begin{aligned} \sum _{i=1}^{N} \widetilde{w}_i\cdot h_j(x_i) = \sum _{i=1}^{N} w_i\cdot h_j(x_i),\quad j=1,\ldots ,k. \end{aligned}$$

Iteratively repeating this procedure, we obtain a measure $R(\mu )\in {\mathcal M}({\mathbb R}^d)$ with

$$\begin{aligned} {\text {supp}}(R(\mu ))\subset {\text {supp}}(\mu ),\, |{\text {supp}}(R(\mu ))|\le k \end{aligned}$$

and

$$\begin{aligned} \int _{{\mathbb R}^d} h_j\, \mathrm{d}\mu = \int _{{\mathbb R}^d} h_j \, \mathrm{d}R(\mu ),\quad j=1,\ldots ,k. \end{aligned}$$

Clearly, at most $N-k$ steps are needed to obtain $R(\mu )$ and in each step the number of basic computational operations needed to compute the support points and the corresponding weights of the actual measure $\widetilde{\mu }$ is bounded by $c\cdot k^3$. Therefore, the total number $\text {op}(R(\mu ))$ of basic computational operations needed to compute $R(\mu )$ satisfies

$$\begin{aligned} \text {op}(R(\mu )) \le c\cdot |{\text {supp}}(\mu )|\cdot k^3. \end{aligned}$$

Let $\gamma \in {\mathbb N}$. The particular choice

$$\begin{aligned} k= \left( {\begin{array}{c}2\gamma +1 +d\\ d\end{array}}\right) ,\,\, \{h_1,\ldots ,h_k\} = \{p_\alpha :\alpha \in {\mathbb N}_0^d, |\alpha |_1\le 2\gamma +1\} \end{aligned}$$

(40)

with the monomials

$$\begin{aligned} p_\alpha (x) =x^\alpha ,\quad x\in {\mathbb R}^d, \end{aligned}$$

yields the operator R employed in Sect. 4.4.

Remark 7

In the case of (40) with $d=1$, the system of linear equations (39) is given by

$$\begin{aligned} Au=0, \end{aligned}$$

(41)

where A is a $(2\gamma +2)\times (2\gamma +3)$ Vandermonde matrix with entries $A_{j,i}=x_i^{j-1}$. Consider the $(2\gamma +3)\times (2\gamma +3)$ Vandermonde matrix $\bar{A} = (x_i^{j-1})_{1\le j,i\le 2\gamma +3}$. Then $\bar{A}$ is invertible and the last column of $\bar{A}^{-1}$ solves the Eq. (41). Using an explicit representation of $\bar{A}^{-1}$, see, e.g.,  [15, Section 14(c)], one obtains

$$\begin{aligned} \text {ker}(A)=\left\{ \left( z\cdot \prod _{k\not =i}(x_k-x_i)^{-1} \right) _{i=1, \ldots , 2\gamma +3}:z\in {\mathbb R}\right\} \end{aligned}$$

for the null space $\text {ker}(A)$ of A.

In the case of (40) with $d>1$, explicit formulas for a solution u of Eq. (39) seem to be unknown up to now, such that suitable numerical methods have to be employed.

1.2 Smoothness Properties of the Associated Semigroup

Let $s,r\in {\mathbb N}$, $\beta \ge 0$ and $K\ge 1\ge \lambda >0$. We fix an equation

$$\begin{aligned} (x_0,a,b)\in {\mathcal C}_{K,\lambda }^{s,s} \end{aligned}$$

(42)

and we consider the semigroup of linear operators $(P_t)_{t\in [0, \infty )}$ associated with a and b, see (22).

Lemma 8

For every $f\in {F}^{r,\beta }$ and $t\in (0,1]$ we have

$$\begin{aligned} P_tf\in C^{s+1}({\mathbb R}^d) \end{aligned}$$

and for all $y\in {\mathbb R}^d$ and all $\alpha \in {\mathbb N}_0^d$ with $0\le |\alpha |_1\le s+1$,

$$\begin{aligned} |(P_tf)^{(\alpha )}(y)|\le c\cdot \Vert f\Vert _{r, \beta }\cdot (1+|y|^\beta )\cdot \frac{1}{t^{(|\alpha |_1-\min (|\alpha |_1,r))/2}}, \end{aligned}$$

(43)

where $c= c(d,m,s,r,\beta ,K, \lambda )$.

Proof

In the sequel, unspecified positive constants c may only depend on the parameters $d,m,s,r,\beta ,K$ and $\lambda$.

By (42) the distribution of $X^y(t)$ has a Lebesgue density for all $y\in {\mathbb R}^d$ and $t\in (0,1]$. More precisely, put $\sigma = bb^T$, let

$$\begin{aligned} V = \{(v,y,t,z)\mid 0\le v< t \le 1,\,y,z\in {\mathbb R}^d\} \end{aligned}$$

and consider the partial differential equation

$$\begin{aligned} \frac{1}{2}\sum _{i,j=1}^d \sigma _{i,j}(y)\cdot \frac{\partial ^2 u}{\partial y_i \partial y_j}(v,y)+\sum _{i=1}^d a_{i}(y)\cdot \frac{\partial u}{\partial y_i}(v,y)+\frac{\partial u}{\partial v}(v,y) = 0 \end{aligned}$$

(44)

for $(v,y)\in [0,1)\times {\mathbb R}^d$. By (42) the Eq. (44) has a unique fundamental solution

$$\begin{aligned} G:V\rightarrow [0,\infty ) \end{aligned}$$

and we have

$$\begin{aligned} {\mathbb P}_{X^y(t)}(A) = \int _A G(0,y,t,z)\mathrm{d}z \end{aligned}$$

(45)

for every Borel set $A\subset {\mathbb R}^d$ and all $y\in {\mathbb R}^d$ and $t\in (0,1]$. See [5, Theorem 6.5.4].

Moreover, by [4, Chapter 9,Theorem 7] we have $G(0,\cdot ,t,z)\in C^{s+1}({\mathbb R}^d)$ for all $t\in (0,1]$ and $z\in {\mathbb R}^d$ with

$$\begin{aligned} \Bigl |\frac{\partial ^{\alpha }G}{\partial y^\alpha }(0,y,t,z)\Bigr |\le \frac{c}{ t^{(|\alpha |_1+d)/2}} \cdot \exp \left( -c\cdot \frac{|z-y|^2}{t}\right) \end{aligned}$$

(46)

for all $\alpha \in {\mathbb N}_0^d$ with $0\le |\alpha |_1 \le s+1$ and $y\in {\mathbb R}^d$.

Using (45) and (46), we obtain $P_t f\in C^{s+1}({\mathbb R}^d)$ with

$$\begin{aligned} (P_tf)^{(\alpha )}(y)=\int _{{\mathbb R}^d} \frac{\partial ^{\alpha }G}{\partial y^\alpha }(0,y,t,z)\cdot f(z)\mathrm{d}z \end{aligned}$$

(47)

for all $f\in {F}^{r,\beta }$, $t\in (0,1]$, $y\in {\mathbb R}^d$ and $\alpha \in {\mathbb N}_0^d$ with $0\le |\alpha |_1 \le s+1$.

Fix f, t, y of the latter type and $\alpha \in {\mathbb N}_0^d$ with $1\le |\alpha |_1 \le s+1$. Put

$$\begin{aligned} q=\min (|\alpha |_1,r) \end{aligned}$$

as well as

$$\begin{aligned} I =\int _{{\mathbb R}^d} \left( f(z) -\sum _{\rho \in {\mathbb N}_0^d:0\le |\,\rho |_1\le q-1}\frac{f^{(\rho )}(y)}{\rho !} \cdot (z-y)^{\rho } \right) \cdot \frac{\partial ^{\alpha }G}{\partial y^\alpha }(0,y,t,z)\mathrm{d}z \end{aligned}$$

and

$$\begin{aligned} I_\rho =\frac{f^{(\rho )}(y)}{\rho !} \cdot \int _{{\mathbb R}^d} (z-y)^{\rho } \cdot \frac{\partial ^{\alpha }G}{\partial y^\alpha }(0,y,t,z)\mathrm{d}z \end{aligned}$$

for $\rho \in {\mathbb N}_0^d$ with $0\le |\rho |_1\le q-1$.

Using integration by parts as well as (46) and (47) we obtain

$$\begin{aligned} I_\rho= & {} (-1)^{|\rho |_1}\cdot f^{(\rho )}(y) \cdot \int _{{\mathbb R}^d} \frac{\partial ^{\alpha -\rho }G}{\partial y^{\alpha -\rho }}(0,y,t,z)\mathrm{d}z\\= & {} (-1)^{|\rho |_1}\cdot f^{(\rho )}(y) \cdot (P_t 1)^{(\alpha -\rho )}(y)=0 \end{aligned}$$

for all $\rho \in {\mathbb N}_0^d$ with $0\le |\rho |_1\le q-1$. Hence

$$\begin{aligned} (P_tf)^{(\alpha )}(y) = I. \end{aligned}$$

Clearly,

$$\begin{aligned} I = \sum _{\rho \in {\mathbb N}_0^d:|\rho |_1=q}\int _{{\mathbb R}^d} f^{(\rho )}(\theta _{z})\cdot (z-y)^{q}\cdot \frac{\partial ^{\alpha }G}{\partial y^\alpha }(0,y,t,z)\mathrm{d}z, \end{aligned}$$

where $\theta _{z} \in {\mathbb R}^d$ satisfies $|\theta _{z}-y|\le |z-y|$. Since $f\in F^{r,\beta }$ we get, similarly to (30),

$$\begin{aligned} |f^{(\rho )}(\theta _{z})|\le c\cdot \Vert f\Vert _{r, \beta }\cdot (1+|y|^\beta +|z-y|^\beta ) \end{aligned}$$

for all $\rho \in {\mathbb N}_0^d$ with $|\rho |_1=q$ and all $z\in {\mathbb R}^d$. Using (46) we therefore obtain

$$\begin{aligned} |I|&\le \frac{c}{ t^{(|\alpha |_1+d)/2}} \cdot \Vert f\Vert _{r, \beta }\cdot \int _{{\mathbb R}^d} (1+|y|^\beta +|z-y|^\beta )\cdot |z-y|^{q}\cdot \exp \left( -c\cdot \frac{|z-y|^2}{t}\right) \mathrm{d}z\\&\le \frac{c}{ t^{(|\alpha |_1+d-q)/2}} \cdot \Vert f\Vert _{r, \beta }\cdot \int _{{\mathbb R}^d} (1+|y|^\beta +|z-y|^\beta )\cdot \exp \left( -c\cdot \frac{|z-y|^2}{t}\right) \mathrm{d}z\\&\le \frac{c}{ t^{(|\alpha |_1-q)/2}} \cdot \Vert f\Vert _{r, \beta }\cdot (1+|y|^\beta ), \end{aligned}$$

which completes the proof of (43) in the case $1\le |\alpha |_1\le s+1$. The case $\alpha =0$ is treated analogously.$\square$