一类倒向随机微分方程多步格式的最优收敛阶估计

Sharp Convergence Rates of a Class of Multistep Formulae for Backward Stochastic Differential Equations

摘要: 自解的存在唯一性问题被解决后，倒向随机微分方程逐步被应用于众多研究领域，例如随机最优控制、偏微分方程、金融数学、风险度量、非线性期望等。本文借助等距节点插值型求积公式的误差余项研究一类倒向随机微分方程多步格式的收敛阶问题。通过在空间层采用Chebyshev网格并结合重心Lagrange插值，提高了算法的数值精度。数值实验结果表明新提出的收敛阶是最优的。

Abstract: Since the well-posedness of its solution was established, the backward stochastic differential equa-tion has been applied in many research fields, such as stochastic optimal control, partial differential equations, financial mathematics, risk measurement, nonlinear expectation and so on. This paper studies the convergence rate of a class of multistep formulae for backward stochastic differential equations with the help of remainder of the quadrature rule over the uniform grid. Due to the ap-plication of barycentric Lagrange interpolation with Chebyshev grids, computational accuracy of the multistep formula is greatly improved. Numerical results indicate that the proposed convergence order is sharp.

1. 引言

1973年，Bismut [1] 首次引入了线性倒向随机微分方程的概念，并研究了该类方程解的存在唯一性。基于如下的BSDE，

$\{\begin{array}{l}-\text{d}{y}_t=f\left(t,y_t,z_t\right)\text{d}t-z_t\text{d}{W}_t\\ y_T=\epsilon \end{array}$

Pardoux和Peng [2] 得到了非线性倒向随机微分方程解的存在唯一性。这一重要研究成果奠定了正倒向随机微分方程组的理论基础。从那时起，BSDE得到了许多研究者的广泛研究。在之后的研究中，通过非线性Feynman-Kac公式，Ma，Protter和Yong [3] 提出了研究正倒向随机微分方程组的四步法，该方法证明了在任意时间区间上的解的存在唯一性。1991年，Peng [4] 得到了正倒向随机微分方程(FBSDEs)和偏微分方程之间的直接关系，然后在 [5] 中，他还找到了随机控制问题的最大值原理。

虽然BSDE的解析解往往难以获得，但BSDE的近似数值解相对容易计算，因此在实际应用中，求解近似数值解成为一种非常理想的方法。求解BSDE的数值方法主要有两类。一种是基于BSDE与其相应的抛物型偏微分方程之间的关系，Ma [3] 等利用这种关系，首次提出了一种典型的研究FBSDEs的四步方法(而不是数值方法)。基于在四步法中提出的类似思想，后续又有学者提出了一些数值求解的算法。另一种方法是直接基于BSDE设计的，当函数f不依赖于变量 $Z_t$ 时，Chevance [6] 提出了一种有效的时空离散化方案。Bender [7] 等引入了一种求解BSDEs的正演方案。Peng [8] 提出了一种在合理假设下收敛的迭代线性近似算法。Memin等在 [9] 中提出了一些求解BSDEs的随机游走方法，并证明了这些方法的收敛性。在某些较弱的正则性假设下，Zhang提出了一种改进的求解BSDEs的数值方案，并在 [10] 中研究了其收敛速度，然后Cvitanic等在 [11] 中提出了一些求解FBSDEs的Monte Carlo方法。Zhao等在 [12] 中提出了一种高精度求解BSDEs的数值方法，Zhao等在 [13] 中研究了 $\theta$ 格式的误差估计。2014年，Zhao等在 [14] 中从BSDEs中推导出两个参考ODE，然后逼近参考ODE中的条件期望和导数，得到了带导数的数值格式。同年，Chassagneux在 [15] 中，提出了BSDEs的线性多步方法。Chassagneux等在 [16] 中，提出了BSDEs的Runge-Kutta方法。2017年，Tang等在 [17] 中提出了BSDEs中的延迟修正方法，其关键特点之一是迭代使用低阶数值方法，从而得到高阶数值格式。2019年，Zhang等在 [18] 中导出了一类单步多导数方法。通过几个数值例子说明了该方法的计算效率和高阶精度。为了说明该方法的优点，并与 $\theta$ 方法进行了比较。2021年，Zhou等在 [19] 中结合Feynman-Kac公式配合Itô-Taylor展开，在数值求解 $z_t$ 时使用有限差分，提出了一类BSDEs的高阶一步格式。

目前国内外学者一直致力于对BSDEs数值解的研究，想要寻求具有较高收敛阶的数值方法。2010年，Zhao等在 [20] 中利用时空网格提出了一种稳定的多步方案。Zhao等证明了所提出的半离散格式的收敛性可以是高阶的，而精确的阶数取决于所使用的时间网格点。本文拟在Zhao的基础上采用Chebyshev网格进行空间离散，使用重心Lagrange插值 [21] 近似非网格点的插值，最后，我们在证明过程中考虑了Newton-Cotes系数的代数精度，并给出了多步格式收敛阶的证明。

2. 多步半离散格式

令 $\left\{\Omega ,\mathcal{F},P,{\left\{{\mathcal{F}}_t\right\}}_{0\le t\le T}\right\}$ 是一个完整的、滤波的概率空间，在这上面定义了一个标准的d维布朗运动 $W_t$。如果一个过程 $\left(y_t,z_t\right):\left[0,T\right]\times \Omega \to {ℝ}^m\times {ℝ}^{m\times d}$ 它是 $\left\{{\mathcal{F}}_t\right\}$ -适应，平方可积而且满足(1)，

$y_t=\xi +{\displaystyle {\int }_t^{T}f\left(s,y_s,z_s\right)\text{d}s}-{\displaystyle {\int }_t^T{z}_s\text{d}{W}_s},t\in \left[0,T\right),$ (1)

我们称这个过程为(1)的 $L^2$ -适应解。其中 $f:\left[0,T\right]\times \Omega \times {ℝ}^m\times {ℝ}^{m\times d}\to {ℝ}^m$ 是关于每一个 $\left(y_t,z_t\right)$ 适应的随机过程，右边的第三项是一个Itô型积分。

Pardoux和Peng证明了(1)的解的存在唯一性，我们假设初值 $y_t$ 形式为 $\phi \left(W_t\right)$，然后 的解 $\left(y_t,z_t\right)$ 能够表示为如下形式，

$y_t=u\left(t,W_t\right)z_t=\nabla u\left(t,W_t\right)\forall t\in \left[0,T\right).$

其中， $\nabla u$ 表示 $u\left(t,x\right)$ 相对于空间变量x的梯度，而 $u\left(t,x\right)$ 是以下抛物型偏微分方程的解：

$\frac{\partial u}{\partial t}+\frac{1}{2}{\displaystyle {\sum }_{i=1}^d\frac{{\partial }^{2}u}{\partial x_i^2}}+f\left(t,u,\nabla u\right)=0.$

考虑情况 $m=d=1$，对方程(1)的区间 $\left[0,T\right]$ 进行均匀网格划分，N为一个正常数，有 $0=t_0<t_1<\cdots <t_N=T$， $t_i=t_0+i\Delta t\left(i=0,1,\cdots ,N\right)$ 和时间步长 $\Delta t$，对于两个给定的正常数k和 $K_y$ 且满足 $1\ll k\ll K_y\ll N$，考虑

$y_{t_n}=y_{t_{n+k}}+{\displaystyle {\int }_{t_n}^{t_{n+k}}f\left(s,y_s,z_s\right)\text{d}s}-{\displaystyle {\int }_{t_n}^{t_{n+k}}{z}_s\text{d}{W}_s},$ (2)

在(2)的等式两端同取期望 ${\mathbb{E}}_{t_n}^x[\cdot ]$，我们得到

$y_{t_n}={\mathbb{E}}_{t_n}^x\left[y_{t_{n+k}}\right]+{\displaystyle {\int }_{t_n}^{t_{n+k}}{\mathbb{E}}_{t_n}^x}\left[f\left(s,y_s,z_s\right)\right]\text{d}s.$ (3)

在等式(3)两边同乘以 $\Delta W_S$，再取条件数学期望，根据Itô等距公式，

$0={\mathbb{E}}_{t_n}^x\left[y_{t_{n+l}}\Delta W_{t_{n+l}}\right]+{\displaystyle {\int }_{t_n}^{t_{n+l}}{\mathbb{E}}_{t_n}^x}\left[f\left(s,y_s,z_s\right)\Delta W_S\right]\text{d}s-{\displaystyle {\int }_{t_n}^{t_{n+l}}{\mathbb{E}}_{t_n}^x}\left[z_s\right]\text{d}s.$

注意到 ${\mathbb{E}}_{t_n}^x\left[f\left(s,y_s,z_s\right)\right]$ 是关于s的确定性函数，因此我们选择用Lagrange插值基于插值点 $\left(t_{n+i},{\mathbb{E}}_{t_n}^x\left[f\left(t_{n+i},y_{t_{n+i}},z_{t_{n+i}}\right)\right]\right),i=0,1,\cdots ,K_y$ 去近似这个确定函数，我们有

${\displaystyle {\int }_{t_n}^{t_{n+k}}{\mathbb{E}}_{t_n}^x\left[f\left(s,y_s,z_s\right)\right]\text{d}s}={\displaystyle {\int }_{t_n}^{t_{n+k}}{P}_{K_y}^{t_n,x}\left(s\right)\text{d}s}+R_y^n,$

其中

$R_y^n={\displaystyle {\int }_{t_n}^{t_{n+k}}\left\{{\mathbb{E}}_{t_n}^x\left[f\left(s,y_s,z_s\right)\right]-P_{K_y}^{t_n,x}\left(s\right)\right\}\text{d}s}.$

当时间为等距划分时，我们有

${\displaystyle {\int }_{t_n}^{t_{n+k}}{P}_{K_y}^{t_n,x}\left(s\right)\text{d}s}=k\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^k{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y_{t_{n+i}},z_{t_{n+i}}\right)\right],$ (4)

其中

$b_{k_y,i}^k=\frac{1}{k\Delta t}{\displaystyle {\int }_{t_n}^{t_{n+k}}{\displaystyle {\prod }_{\begin{array}{l}j=0\\ j\ne i\end{array}}^{K_y}\left(\frac{s-t_{n+j}}{t_{n+i}-t_{n+j}}\right)\text{d}s}},$

即

$b_{k_y,i}^k=\frac{1}{k}{\displaystyle {\int }_0^t{\displaystyle {\prod }_{\begin{array}{l}j=0\\ j\ne i\end{array}}^{K_y}\left(\frac{s-j}{i-j}\right)\text{d}s}},i=0,1,\cdots ,K_y.$ (5)

通过上面的推导我们得到关于 $y_t$ 的多步半离散格式为

$y_{t_n}={\mathbb{E}}_{t_n}^x\left[y_{t_{n+k}}\right]+k\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^k{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y_{t_{n+i}},z_{t_{n+i}}\right)\right]+R_y^n.$ (6)

类似可以得到关于 $z_t$ 的多步半离散格式为

$0={\mathbb{E}}_{t_n}^x\left[z_{t_{n+l}}\right]+{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^l{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y_{t_{n+i}},z_{t_{n+i}}\right)\Delta W_{t_{n+i}}\right]-{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^l{\mathbb{E}}_{t_n}^x}\left[z_{t_{n+i}}\right]+R_z^n,$ (7)

其中

$\begin{array}{c}{R}_z^n={\displaystyle {\int }_{t_n}^{t_{n+l}}{\mathbb{E}}_{t_n}^x}\left[f\left(s,y_s,z_s\right)\Delta W_S\right]\text{d}s-{\displaystyle {\int }_{t_n}^{t_{n+l}}{\mathbb{E}}_{t_n}^x\left[z_s\right]\text{d}s}+{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^l{\mathbb{E}}_{t_n}^x}\left[z_{t_{n+i}}\right]\\ -l\Delta t{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^l{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y_{t_{n+i}},z_{t_{n+i}}\right)\Delta W_{t_{n+i}}\right].\end{array}$

这里的推导过程我们用到了等式(8)

${\mathbb{E}}_{t_n}^x\left[y_{t_{n+l}}\Delta W_{t_{n+l}}\right]=l\Delta t{\mathbb{E}}_{t_n}^x\left[z_{t_{n+l}}\right].$ (8)

3. 最优收敛阶证明

在本小节中，我们将对多步半离散方案进行误差分析，即生成器函数f独立于随机变量 $z_t$，对于方程，

$y_t=\phi \left(W_t\right)+{\displaystyle {\int }_t^{T}f\left(s,y_s\right)\text{d}s}-{\displaystyle {\int }_t^T{z}_s\text{d}{W}_s}$ (9)

为了简单的分析，我们只考虑 $m=d=1$ 的情况，但下面得到的所有结果也适用于一般的情况。考虑如下的简化形式，

$y_{t_n}={\mathbb{E}}_{t_n}^x\left[y_{t_{n+k_y}}\right]+k_y\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^k{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y_{t_{n+i}}\right)\right]+R_y^n,$

$0={\mathbb{E}}_{t_n}^x\left[z_{t_{n+1}}\right]+{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^l{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y_{t_{n+i}}\right)\Delta W_{t_{n+i}}\right]-{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^l{\mathbb{E}}_{t_n}^x}\left[z_{t_{n+i}}\right]+\frac{1}{\Delta t}{R}_z^n.$ (10)

相应的我们的多步方案为：

$y^n={\mathbb{E}}_{t_n}^x\left[y^{n+K_y}\right]+K_y\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^{K_y}{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y^{n+i}\right)\right],$

$0={\mathbb{E}}_{t_n}^x\left[z^{n+1}\right]+{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^1{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y^{n+i}\right)\Delta {W^{\prime }}_{t_{n+i}}\right]-{\displaystyle {\sum }_{i=0}^{K_z}{b}_{K_z,i}^1{\mathbb{E}}_{t_n}^x}\left[z^{n+i}\right].$ (11)

其中 $n=N-K,\cdots ,0$。

对于后续的证明，这里我们用到了Zhao [20] 中的假设1，假设2，引理1，引理3。对于Zhao [20] 中的引理2，我们考虑Newton-Cotes系数的代数精度，进行了修改如下。

引理2 设 $R_y^n$ 和 $R_z^n$ 为参考方程(10)中定义的局部截断误差。然后在假设2下，我们有以下的局部估计：

$\left|R_y^n\right|\le C{\left(\Delta t\right)}^{k_y+2}$， $k_y$ 为奇数时；

$\left|R_y^n\right|\le C{\left(\Delta t\right)}^{k_y+3}$， $k_y$ 为偶数时；

$\left|R_z^n\right|\le C{\left(\Delta t\right)}^{k_z+2}$。

其中 $C>0$ 是一个一般常数，只依赖于T， $\phi$ 的上界和f及其它们的导数。

上述引理的证明与 [13] 中的引理3.2的证明相似，因此我们在这里省略了它。但考虑到多步格式(11)中y的求解方程系数为Newton-Cotes系数，结合Lagrange插值的局部截断误差和Newton-Cotes系数代数精度我们可以得到以上 $R_y^n$ 的结论。

现在我们在下面的定理中给出了 $y_{t_n}-y^n$ 的误差估计：

定理1 设 $y_{t_n}$ 和 $y^n$ 分别为BSDE和多步格式的解。在假设2成立的基础上，并且初始值满足 ${\max }_{N-K_y<j\le N}\mathbb{E}\left[\left|y_{t_j}-y^j\right|\right]=O\left({\left(\Delta t\right)}^{K_y+1}\right)$。然后对于足够小的时间步长 $\Delta t$，我们有

$\underset{0\le n\le N}{\sup }\mathbb{E}\left[\left|y_{t_n}-y^n\right|\right]\le C{\left(\Delta t\right)}^{K_y+1}$， $k_y$ 为奇数时；

$\underset{0\le n\le N}{\sup }\mathbb{E}\left[\left|y_{t_n}-y^n\right|\right]\le C{\left(\Delta t\right)}^{K_y+2}$， $k_y$ 为偶数时；

其中 $C>0$ 是一个一般常数，依赖于T， $\phi$ 的上界和f及其它们的导数。

证明：令 $e_y^n=y_{t_n}-y^n$ 对于 $n=N,N-1,\cdots ,0$。从(10)和(11)我们得到

$e_y^n={\mathbb{E}}_{t_n}^x\left[e_y^{n+K_y}\right]+K_y\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^{K_y}{\mathbb{E}}_{t_n}^x}\left[f\left(t_{n+i},y_{t_{n+i}}\right)-f\left(t_{n+i},y^{n+i}\right)\right]+R_y^n.$ (12)

然后我们有

$\left|e_y^n\right|\le {\mathbb{E}}_{t_n}^x\left[\left|e_y^{n+K_y}\right|\right]+L{K}_y\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^{K_y}{\mathbb{E}}_{t_n}^x}\left[\left|e_y^{n+i}\right|\right]+\left|R_y^n\right|,$ (13)

其中L是函数 $f\left(t,y_t\right)$ 相对于 $y_t$ 的Lipschitz常数。对不等式(13)两边同取数学期望 $\mathbb{E}[\cdot ]$，我们得到

$\mathbb{E}\left[\left|e_y^n\right|\right]\le \mathbb{E}\left[\left|e_y^{n+K_y}\right|\right]+L{K}_y\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^{K_y}\mathbb{E}}\left[\left|e_y^{n+i}\right|\right]+\mathbb{E}\left[\left|R_y^n\right|\right].$ (14)

令 $N_k=\left[\frac{N-n}{K_y}\right]$。对于满足 $1\le s\le N_k$ 的整数s，我们也有同样的估计

$\mathbb{E}\left[\left|e_y^{n+\left(s-1\right)K_y}\right|\right]\le \mathbb{E}\left[\left|e_y^{n+s{K}_y}\right|\right]+L{K}_y\Delta t{\displaystyle {\sum }_{i=0}^{K_y}{b}_{K_y,i}^{K_y}\mathbb{E}}\left[\left|e_y^{n+\left(s-1\right)K_y+i}\right|\right]+\mathbb{E}\left[\left|R_y^{n+\left(s-1\right)K_y}\right|\right].$ (15)

然后我们将上面不等式(15)按照 $s=1,2,\cdots ,N_k$ 相加可以得到

$\mathbb{E}\left[\left|e_y^n\right|\right]\le \mathbb{E}\left[\left|e_y^{n+N_k{K}_y}\right|\right]+2L{K}_y\Delta t{\displaystyle {\sum }_{i=0}^{N_k{K}_y}\mathbb{E}}\left[\left|e_y^{n+i}\right|\right]+{\displaystyle {\sum }_{i=0}^{N_k-1}\mathbb{E}}\left[\left|R_y^{n+i{K}_y}\right|\right],$ (16)

整理有

$\left(1-2L{K}_y\Delta t\right)\mathbb{E}\left[\left|e_y^n\right|\right]\le \mathbb{E}\left[\left|e_y^{n+N_k{K}_y}\right|\right]+2L{K}_y\Delta t{\displaystyle {\sum }_{i=n+1}^{N_k{K}_y}\mathbb{E}}\left[\left|e_y^i\right|\right]+{\displaystyle {\sum }_{i=0}^{N_k-1}\mathbb{E}}\left[\left|R_y^{n+i{K}_y}\right|\right].$ (17)

这里我们令 $D=\frac{1}{1-2L{K}_y\Delta t}$， $N_1=\frac{2L{K}_y}{1-2L{K}_y\Delta t}$，以及 $M_0={\max }_{N-K_y<j\le N}\mathbb{E}\left[\left|e_y^j\right|\right]+{\displaystyle {\sum }_0^{N_k-1}\mathbb{E}}\left[\left|R_y^{n+i{K}_y}\right|\right]$。

对于足够小的时间步长 $\Delta t$，D和 $N_1$ 明显是正数，并由一个正的常数为界。然后通过引理3和上式，得到下面的不等式：

$\mathbb{E}\left[\left|e_y^n\right|\right]\le {\text{e}}^{N_{1}T}\left(D{M}_0+\Delta t{K}_y{N}_1{M}_0\right)={\text{e}}^{N_{1}T}\left(D+\Delta t{K}_y{N}_1\right)M_0.$

通过引理2，我们可以知道

$\mathbb{E}\left[\left|R_y^{n+i{K}_y}\right|\right]\le C{\left(\Delta t\right)}^{k_y+2}$， $k_y$ 为奇数时；

$\mathbb{E}\left[\left|R_y^{n+i{K}_y}\right|\right]\le C{\left(\Delta t\right)}^{k_y+3}$， $k_y$ 为偶数时，

即

${\displaystyle {\sum }_{i=0}^{N_k-1}\mathbb{E}\left[\left|R_y^{n+i{K}_y}\right|\right]}\le N\mathbb{E}\left[\left|R_y^{n+i{K}_y}\right|\right]\le CN{\left(\Delta t\right)}^{K_y+2}=C{\left(\Delta t\right)}^{K_y+1}$， $k_y$ 为奇数时；

${\displaystyle {\sum }_{i=0}^{N_k-1}\mathbb{E}\left[\left|R_y^{n+i{K}_y}\right|\right]}\le N\mathbb{E}\left[\left|R_y^{n+i{K}_y}\right|\right]\le CN{\left(\Delta t\right)}^{K_y+3}=C{\left(\Delta t\right)}^{K_y+2}$， $k_y$ 为偶数时。

在定理的初始条件 ${\max }_{N-K_y<j\le N}\mathbb{E}\left[\left|y_{t_j}-y^j\right|\right]=O\left({\left(\Delta t\right)}^{K_y+1}\right)$ 下，我们可以得到结论

$M_0\le C{\left(\Delta t\right)}^{K_y+1}$， $k_y$ 为奇数时；

$M_0\le C{\left(\Delta t\right)}^{K_y+2}$， $k_y$ 为偶数时。

综上，得证。

我们在下面的定理中给出 $z_{t_n}-z^n$ 的误差估计。

定理2 设 $z_{t_n}$ 和 $z^n$ 分别为BSDE和多步格式的解。在假设2成立的基础上，并且初始值满足 ${\max }_{N-K_z<j\le N}\mathbb{E}\left[\left|z_{t_j}-z^j\right|\right]=O\left({\left(\Delta t\right)}^{K_z}\right)$。然后对于足够小的时间步长 $\Delta t$，我们有

$\underset{0\le n\le N}{\sup }\mathbb{E}\left[\left|z_{t_n}-z^n\right|\right]\le C{\left(\Delta t\right)}^{\min \left(K_y+1,K_z\right)},$

其中 $C>0$ 是一个一般常数，依赖于T， $\phi$ 的上界和f及其它们的导数。

对于 $z_{t_n}-z^n$ 的误差估计证明见Zhao [20] 中定理2的证明。

4. 数值实验

为数值求解 $\left(y_t,z_t\right)$，我们也需要空间离散，这里我们选择Chebyshev网格点。为此，我们在时间层 $t_n$ 上引进 $R^d$ 空间上的一个剖分 $D_h^n=\left\{X_i|X_i\in R^d\right\}$。用 $h^n>0$ 表示剖分 $D_h^n$ 中网格点的密度：

$h^n=\underset{x\in {ℝ}^q}{\max }\underset{x_i\in D_h^n}{\min }\left|x-x_i\right|=\underset{x\in {ℝ}^q}{\max }dist\left(x,D_h^n\right).$

同理，推导全离散下 $\left(y_t,z_t\right)$ 的格式为

$y_i^n={\hat{\mathbb{E}}}_{t_n}^{xi}\left[{\hat{y}}^{n+K_y}\right]+K_y\Delta t{\displaystyle {\sum }_{j=0}^{K_y}{b}_{K_y,j}^{K_y}{\hat{\mathbb{E}}}_{t_n}^{xi}}\left[f\left(t_{n+j},{\hat{y}}^{n+j},{\hat{z}}^{n+j}\right)\right]+K_y\Delta t{b}_{K_y,0}^{K_y}f\left(t_n,y_i^n,z_i^n\right),$

$0={\hat{\mathbb{E}}}_{t_n}^{xi}\left[{\hat{z}}^{n+1}\right]+{\displaystyle {\sum }_{j=0}^{K_z}{b}_{K_z,j}^1{\hat{\mathbb{E}}}_{t_n}^{xi}}\left[f\left(t_{n+j},{\hat{y}}^{n+j},{\hat{z}}^{n+j}\right)\Delta {W^{\prime }}_{t_{n+j}}\right]-{\displaystyle {\sum }_{j=0}^{K_z}{b}_{K_z,j}^1{\hat{\mathbb{E}}}_{t_n}^{xi}}\left[{\hat{z}}^{n+j}\right]-b_{K_z,0}^1{z}_i^n.$

期望的近似采用Gauss-Hermite求积公式，有如下等式

${\hat{\mathbb{E}}}_{t_n}^{x_i}\left[{\hat{y}}^{n+k}\right]=\frac{1}{{\pi }^{\frac{d}{2}}}{\displaystyle {\sum }_{j=1}^L{\omega }_j{\hat{y}}^{n+k}}\left(x_i+\sqrt{2k\Delta t}{a}_j\right).$

我们进行了一些数值实验，用来说明方案的准确性和有效性。为了方便起见，我们将在时间网格上使用均匀划分，在空间离散上使用Chebyshev非均匀网格进行空间离散化。Chebyshev非均匀网格可以有效地避免Runge现象使空间离散插值近似的更加准确。在我们的数值实验中，我们设定 $T=1$。误差 $e_y^n=y_{t_n}-y^n$ 和 $e_z^n=z_{t_n}-z^n$ 通常来自三个方面：时间离散化，空间插值和用Gauss-Hermite求积公式对条件期望进行近似。我们的主要目标是测试对于时间离散的收敛速度，所以需要尽量平衡其它误差。我们设定Gauss-Hermite正交点的数量为 $L=10$ 用来尽可能减小Gauss-Hermite求积公式带来的误差影响。我们采用重心Lagrange插值来尽可能减小空间离散带来的误差。实验中，我们选用高阶的单步方法或者直接选取精确值作为启动值，利用重心Lagrange插值来获得启动值对应的空间层的非网格点的近似值，应用到多步全离散格式中，循环以上步骤，从而获得数值解 $\left(y^0,z^0\right)$。通过对数值误差进行线性最小二乘拟合，得到了对时间步长 $\Delta t$ 的收敛速度(CR)。

例1考虑如下的线性BSDE

$\{\begin{array}{l}-\text{d}{y}_t=\left(-y_t^3+2.5y_t^2-1.5y_t\right)\text{d}t-z_t\text{d}{W}_t\\ y_T=\frac{\exp \left(W_T+T\right)}{\exp \left(W_T+T\right)+1}\end{array}$

其解析解可以表示为

$\{\begin{array}{l}{y}_t=\frac{\exp \left(W_t+t\right)}{\exp \left(W_t+t\right)+1}\\ z_t=\frac{\exp \left(W_t+t\right)}{{\left(\exp \left(W_t+t\right)+1\right)}^2}\end{array}$

我们使用多步全离散格式来解决例1中的线性BSDE。在表1和表2中，我们分别列出了 $\left|y_0-y^0\right|$ 和 $\left|z_0-z^0\right|$ 的数值误差，以及其收敛速度。

Table 1. Error and convergence rate of | y 0 − y 0 | in Example 1

表1. 例1中 $\left|y_0-y^0\right|$ 的误差和收敛速度

Table 2. Error and convergence rate of | z 0 − z 0 | in Example 1

表2. 例1中 $\left|z_0-z^0\right|$ 的误差和收敛速度

我们设定终端时间为 $T=1$。那么精确解 $\left(y_t,z_t\right)$ 在 $t=0$ 时， $\left(y_0,z_0\right)$ 是 $\left(\frac{1}{2},\frac{1}{4}\right)$。此外，我们将时间层划分设定为 $N=4,8,\cdots ,64$。从表1和表2的数值结果中，我们可以得出结论，多步全离散格式对于求解BSDE是稳定和准确的。此外，实验数据表明，当 $K_y$ 为奇数时，y的收敛阶可以达到 $K_y+1$，当 $K_y$ 为偶数时，y的收敛阶可以达到 $K_y+2$，对于z，收敛阶可以达到 $\min \left(K_y+1,K_z\right)$。

例2考虑如下的线性BSDE

$\{\begin{array}{l}-\text{d}{y}_t=-y_t\left(1-y_t\right)\left(\frac{3}{4}-y_t\right)\text{d}t-z_t\text{d}{W}_t\\ y_T=\frac{1}{\exp \left(-W_T-\frac{T}{4}\right)+1}\end{array}$

其解析解可以表示为

$\{\begin{array}{l}{y}_t=\frac{1}{\exp \left(-W_t-\frac{t}{4}\right)+1}\\ z_t=\frac{\exp \left(-W_t-\frac{t}{4}\right)}{{\left[\exp \left(-W_t-\frac{t}{4}\right)+1\right]}^2}\end{array}$

此例子中，我们设定终端时间为 $T=1$。精确解在 $t=0$ 时为 $\left(\frac{1}{2},\frac{1}{4}\right)$。误差 $\left|y_0-y^0\right|$ 和 $\left|z_0-z^0\right|$ 分别在表3和表4中给出，并且给出了相应的收敛速率。表中数据表明，y的收敛阶与我们定理1的理论推导一致，z的收敛阶符合定理2给出的结论。

Table 3. Error and convergence rate of | y 0 − y 0 | in Example 2

表3. 例2中 $\left|y_0-y^0\right|$ 的误差和收敛速度

Table 4. Error and convergence rate of | z 0 − z 0 | in Example 2

表4. 例2中 $\left|z_0-z^0\right|$ 的误差和收敛速度

5. 结论

本文介绍了一类求解倒向随机微分方程的时空网格上的多步格式。我们在时间上采用Lagrange插值离散积分，在空间上采用Chebyshev网格和重心Lagrange插值，此外，我们还采用了Gauss-Hermite求积公式来求解近似条件数学期望。我们在生成元函数f与随机变量 $z_t$ 无关时，给出了多步半离散格式的最优收敛阶估计。最后，我们通过一些数值实验验证了该多步格式的最优收敛阶：当 $K_y$ 为奇数时，y的收敛阶为 $K_y+1$，当 $K_y$ 为偶数时，y的收敛阶为 $K_y+2$，z的收敛可以达到 $\min \left(K_y+1,K_z\right)$。

基金项目

国家自然科学基金项目(11901133)。

NOTES

^*通讯作者Email: [e-mail]