GPU上求解倒向随机微分方程的多步法

摘要

倒向随机微分方程（BSDE）是定价和套期保值的重要工具。低计算时间的高精度定价对于最小化金钱损失变得很有趣。因此，我们探索了在期权定价中并行化高阶多步方案的机会。在多步骤方案中，每个空间网格点的计算是独立的，这一事实促使我们选择使用CUDA的大规模并行GPU计算。在我们的研究中，我们确定了性能瓶颈，并应用适当的优化技术来减少均匀空间域中的计算时间。运行时实验表明，对于在单个GPU NVIDIA GeForce 1070 Ti上的并行实现，可以实现乐观的加速。

1介绍

倒向随机微分方程（BSDEs）由于其关键特性之一，即提供非线性抛物型偏微分方程（PDEs）解的概率表示，在物理和金融等领域得到了广泛应用。我们考虑一个（解耦）正倒向随机微分方程（FBSDE），其形式如下：

$$\textstyle\开始{cases}dX{t}=a（t，X{t{）\，dt+b（t，X{t}）\结束{cases}$$

(1)

哪里\（X_{t}\），\（a \ in \ mathbb{R}^{n}\），b条是一个\（n次d次）矩阵，\（W_{t}\）是一个d日-维布朗运动，\（f（t，X_{t}，y_{tneneneep，z_{t{}）：[0，t]\times\mathbb{R}^{n}\times\ mathbb}R}^}m}\times \ mathbb{R}^{m\ times d}\to\mathbb是驾驶员功能和ξ是终端条件。对于\（a=0）和\（b=1），即\（X{t}=W{t}\），获得表格的标准BSDE

$$\textstyle\begin{cases}-dy{t}=f（t，y_{t}，z_{t{）\，dt-z{t}\，dW{t}，\\y_{t}=\xi=g（W{t}），\end{casesneneneep$$

(2)

哪里\（y_{t}\in\mathbb{R}^{m}\）和\（f（t，y_{t}，z_{t{）：[0，t]\times\mathbb{R}^{m}\times\ mathbb}R}^}m\ times d}\to\mathbb{R}^{m{）.的解的存在性和唯一性(1)已在中得到证明[20]. 如果\（dX_｛t｝\）在(1)定义为资产价格建模的几何布朗运动，以及克是欧洲期权的收益函数，那么(1)与Black-Scholes PDE相关，请参阅[14]. 这就是说\（y_{0}\）给出了期权价格和\（z{t}/b\）代表hegding投资组合，代表期权价格的敏感性\（y{t}\）关于资产价格\（X_{t}\）这在定价理论中也称为Δ-套期保值。

通常，很难找到BSDE的分析解。最近，人们提出了许多数值方法。对于空间网格上的数值方法，我们参考，例如[29–32]. 基于傅里叶方法的BSDE方法在[24，25]. 对于概率方法，研究了蒙特卡罗方法，例如[10，16]和中基于树的方法[6，26]. 还有其他许多人，看[1，2，4，8，9，17]。

高阶数值格式通常需要更多的计算工作量，因此需要有效的并行实现。针对金融领域的定价问题，已经开发了一些基于图形处理单元（GPU）计算的加速策略，但很少有基于BSDE的方法。这些作品可以在[7，11，22]，其中加速策略应用于收敛阶不高于2的数值方法。

为了获得更高的收敛速度，在中提出了时空网格上的第一个多步格式[32]其中，通过及时离散BSDE得到的被积函数通过使用拉格朗日插值多项式进行近似。在[13]，我们成功地在GPU上并行化了该多步方案，并通过Black-Scholes BSDE显示了期权定价的计算时间增益。然而，由于龙格现象，多步方案仅在3个多时间层内稳定。为了获得更好的稳定性和允许更多的时间级，在[27]用样条代替拉格朗日插值多项式。原则上，在该多步方案中可以选择任意多个多级时间电平，通常时间电平越多，精度越高。然而，使用更多的时间级别也需要更多的计算成本。因此，在这项工作中，我们研究了多步方案中的大规模并行GPU计算[27]，使该方案在实践中更加有用。例如，低计算时间的高精度可以最大限度地减少财务问题中的货币损失。可以在中找到一个应用程序，该应用程序显示了金融模型定价和风险管理的有效方法的重要性[三]。

本文的提醒如下。在Sect。 2，我们介绍了多步骤方法[27]用于BSDE的数值解。章节三介绍了数值方法的算法框架，以及由于其均匀域的特殊结构，初步减少计算时间的潜力。在Sect。 4，我们描述了进一步减少计算时间的GPU加速策略。章节5显示了包括金融应用在内的数值结果。最后，我们在第节中给出了结论。 6。

2多步骤方案

在本节中，我们将介绍多步方案[27]定位为在GPU上并行。

2.1前期工作

让\（（\Omega，\mathcal{F}，\mathbb{P}，\{mathcal{F}（F）_{t} \}_{0\ le t\le t}）\）是一个完整的过滤概率空间。在这个空间中，一个标准d日-维布朗运动\（W_｛t｝\）定义，以便过滤\（数学{F}（F）_{t} }_{0\le-t\le-t}\）是天然过滤\（W_{t}\）。我们定义\（|\cdot|\）作为欧氏空间中的标准欧氏范数\（\mathbb{R}^{m}\）或\（\mathbb{R}^{m\乘以d}\）和\（L^{2}=L^{2}_{\mathcal{F}}（0，T；\mathbb{R}^{d}）\）全套\（\mathcal{F}（F）_{t} \）-值为的自适应平方可积过程\（\mathbb{R}^{d}\）此外，让\（\mathcal{F}（F）_{s} ^{t，x}\）对于\（t\le s\le t\）成为σ-布朗运动产生的场\（\{x+W_{r} -西_{t} ，t\ le r \ le s \｝\）从时空点开始\（（t，x）\）。我们定义\（E_｛s｝^｛t，x｝[x]\）作为随机变量的条件期望X（X）在过滤下\（\mathcal{F}（F）_{s} ^{t，x}\），即。\（E_{s}^{t，x}[x]=E[x\vert\mathcal{F}（F）_{s} ^{t，x}]\）。

一对过程\（（（y_{t}，z_{t{）：[0，t]\times\Omega\to\mathbb{R}^{m}\times\ mathbb}R}^}m\ timesd}\）是BSDE的解决方案(1)如果是的话\（\mathcal{F}（F）_{t} \）-自适应，平方可积，并且满足(1)在某种意义上

$$y_｛t｝=\xi+\int_｛t｝^｛t｝f（s，X_｛s｝，y_｛s｝，z_｛s｝）\，ds-\int_｛t｝^｛t｝z_｛s｝\，dW_｛s｝，\ quad t\ in[0，t）$$

（3）

哪里\（f（t，X_{t}，y_{tneneneep，z_{t{}）：[0，t]\times\mathbb{R}^{m}\times\ mathbb}R}^}n}\times \ mathbb{R}^{m\ times d}\to\mathbb是\（\mathcal{F}（F）_{t} \）-右边的第三项是Itó型积分。此解在正则条件下存在[20]。

让我们考虑半线性PDE

$$\frac{\partialu}{\particalt}+\sum_{i=1}^{n} 一个_{i} （t，x）\分形{\部分u}{\部分x{i}}+\分形{1}{2}\sum{i，j=1}^{n}\sum{k=1}^{d} b条_{i，k}（t，x）b{j，k}（t，x）\分数{\部分^{2}u}{\部分x_{i} x个_{j} }+f\bigl（t，x，u，b（t，x）^{\top}D_{x}u\bigr）=0$$

(4)

与终端条件\（u（T，x）=g（x）\）.通过直接应用Itó’s引理，可以得到以下定理。

定理1

（非线性Feynman-Kac定理）

让 \（在C^{1,2}中为u\） 令人满意的(4)假设存在一个常数 C类 这样的话 \（|b（t，x）^｛\top｝D_｛x｝u（t，x）|\leq C（1+|x|）\） 对于每个 \（（t，x）\in[0，t]\times\mathbb{R}^{m}\），然后

$$y_{t}=u（t，X_{tneneneep），\qquad z_{t{=b（t，X_{tneneneei）^{top}D_{X}u（t、X_{t}）$$

(5)

是唯一的解决方案(1).

我们注意到[19]证明了由具有所有订单矩的Lévy过程驱动的BSDEs解的存在性和唯一性，该解可用于Lév y市场的定价。一般非马尔可夫BSDE的非线性Feynman-Kac公式已在[21]，其中考虑路径相关的拟线性抛物偏微分方程。此外，还研究了完全非线性偏微分方程的Feynman-Kac表示，例如[23]。

2.2稳定的半离散格式

让N个是一个正整数，并且\（增量t=t/N）均匀划分时间间隔的步长\（[0，T]\）:\（0=t_{0}<t_{1}<cdots<t_{N-1}<t_{N}=t\），其中\（t_{n}=t_{0}+n\增量t\），\（n=0,1，点，n）。

让k个和\（K_{y}\）是两个正整数，这样\（1\le-k\le-k_{y}\le-N\），分别表示时间层和插值点的数量。BSDE(2)可以表示为

$$y_{t_{n}}=y_{t_n+k}}+\int_{t_n}}^{t_n+k{}f（s，y_{s}，z_{s{）\，ds-\int_}n}}^{t_n+k}}z_{s}\，dW{s}$$

(6)

采用条件期望\（E_{t_{n}}^{x}[\cdot]\）在(6)为了获得解的适应性并使用三次样条多项式逼近被积函数年进程读取（请参阅附录)

$$y_{t_{n}}=E_{t_{n}}^{x}[y_{t_{n+k}}]+\sum_{j=0}^{k_{y} -1个}\biggl[a{j}^{y}\Delta t+\frac{b_{j}^{y{}\Delta t^{2}+\frac{c_{j{}^{y}\ Delta t^3}{3}+\ frac{d_{j}#y}\ Delta t ^{4}}{4}\biggr]+R_{y}^{n}$$

哪里\（R_｛y｝^｛n｝\）是插值误差。对于z（z）过程，使用我和\（K_{z}\）而不是k个和\（K_{y}\），将两边乘以\（增量W_{t_{n+l}}）在(6)取条件期望\（E_{t_{n}}^{x}[\cdot]\），使用三次样条插值的参考方程读取（参见附录)

$$\开始{对齐}0={}&l\增量tE_{t_{n}}^{x}[z_{t_n+l}}]+\sum_{j=0}^{K_{z} -1个}\biggl[a{j}^{z{1}}\增量t+\frac{b_{j}^{z_{1}{\增量t^{2}}{2}+\frac{c_{j{}^{z{1}}\增量t{3}{3}+\frac{d_{j}^{z{1{}}\Delta t^4}{4}\biggr]\\&{}-\sum_{j=0}^{K_{z} -1个}\biggl[a{j}^{z{2}}\增量t+\frac{b_{j}^{z_2}}\三角形t^{2}{2}+\frac{c_{j{}^{z{2}}\增量t^{3}{3}+\frac{d_{j}*^{z_2}}\三角t^4}}{4}\biggr]+R{z}^n}，\end{对齐}$$

哪里\（R_｛z｝^｛n｝=R_｛z_｛1｝｝^｛n｝+R_｛z_｛2｝｝^｛n｝\）是插值误差。在[27]，作者证明了当

$$\begin{aligned}&k=1，\ldot，k_{y}，\quad\text{with}k_{y}=1，2，3，\ldots，N，\\&l=1，\ quad\text{with{k_{z}=1、2、3，\ltots，N.\end{alinged}$$

也就是说，该算法允许任意多个时间级别\（K_{y}\）和\（K_{z}\）.利用三次样条插值条件计算未知系数，我们得到

$$\开始{对齐}y_{t_{n}}={}&E_{t_{n}}^{x}[y_{t_{n+K_{y}}]+\增量tK_y}\和_{j=0}^{K_y{}}\gamma_{K_{y}，j}^{K_y}}E_{t_n}}^x}\bigl[f（t_{n+j}，y_{t}{n+j}}，z{t_{n+j}}）\bigr]+R{y}^{n}，\\z{t_n}}={}&\Biggl（E_{t_n{}}^x}[z{t_}n+1}}]+\sum{j=1}^{K{z}}\gamma_{K{z}，j}^1}E_{t{n}^{x}}\bigl[f（t{n+j}，y{t{n+j}}{n}}{\增量t}，\结束{对齐}$$

(7)

具有\（\gamma_{K_{y}，j}^{K_}y}）和\（\gamma_{K_{z}，j}^{1}\）表示三次样条插值的计算系数（表1和表2给出最多6个时间级别的值）。如所示[27]本地错误(7)由提供

$$\bigl\vert R_{y}^{n}\bigr\vert=\mathcal{O}\bigl$$

前提是（f）和克足够光滑。在(7)我们需要除以Δ吨找到的价值z（z）过程。因此，为了平衡时间截断错误，可以设置\（K{z}=K{y}+1）。

表1系数\（{\gamma_{K_{y}，j}^{K__y}}{j=0}^{K_y}）直到\（K_{y}=6\）

表2系数\（{\gamma_{K_{z}，j}^{1}{j=0}^{K_}z}}）直到\（K_{z}=6\）

d维情形的稳定半离散格式如下：我们表示\（（y^{n}，z^{n{）\）近似于\（（y_｛t_｛n｝｝，z_｛t_｛n｝｝）\），给定随机变量\（（y^｛N-i｝，z^｛N-i｝）\），\（i=0,1，点，K-1）具有\（K=\max\{K_{y}，K_{z}\}）.然后\（（y^{n}，z^{n{）\）可以为找到\（n=n-K，点，0）这样的话

$$开始{对齐}&y^{n}=E_{t_{n}}^{x}\bigl[y^{n+K_{y}}\bigr]+\Delta tK_{y}\sum_j=0}^{K_y}}\gamma_{K_y}，j}^{K_y{}}E_{t_n}}x}\bigl[\bigl（t_{n+j}，y^{n+j{，z^{n+j}\bigr）\bigr]，\\&z^{n}=\Biggl（E_{t_{n}}^{x}\bigl[z^{n+1}\biger]+\sum_{j=1}^{K{z}}\gamma_{K{z}，j}^{1}E_{t{n}{x}\bigl[f\bigl（t_{n+1}，y^{n+j}，z^{n+j}\较大）\三角洲W_{t_{n+j}}^{\top}\bigr]-\sum_{j=1}^{K_{z}\gamma_{K_}z}，j}^{1}E_{t_{n}}^}x}\bigl[z^{n+j}\biger]\Biggr）\Big/\gamma_{K{z}，0}^{1}\end{aligned}$$

(8)

哪里\（y^{n}=（y^}n，\波浪线{m}}）{\波浪线}\times1}\），\（z^{n}=（z^}n，\瓦尔德{m}，\瓦尔特{d}）{\瓦尔德}m}\乘以d}\），\（Delta W_{t_{n+j}}^{top}=（W_{t_{n+j}}^}{\tilde{d}}）_{\tilde{d}\times1}-（W_}t_{n}}^\\tilde{d\times1}），\（tilde{m}=1,2，\ldots，m\）和\（tilde{d}=1,2，\ldots，d\）在下文中，我们仅给出了误差分析的结果，以供参考[27]和[32]。

引理1

中局部截断误差的局部估计(7)满足

$$\bigl\vert R_{y}^{n}\bigr\vert\le C\Delta t^{min\{K_{y{+2,5\}}，\qquad\bigl\overt R_{z}^{n}\biger\vert\leC\Deltat^{min \{K_{z}+2,5\}}$$

哪里 \（C>0\） 是一个常数，取决于 T型，（f），克 和的导数 （f） 和 克。

定理2

假设初始值满足

$$\textstyle\begin{cases}\max_{N-K_{y}+1\len\leN}E[\verty_{t_N}}-y^{N}\vert]=\mathcal{O}（\Delta t^{K_y}+1}），\\max_{N-K_{y}+1 \len\le N}E[\vert y_{t_{N}-y_{N{}\vert]=\mathcal{0}（\ Delta t_{4}），\结束{cases}$$

哪里 \（K_{y}=1,2,3\） 对于第一个方程和 \（K_{y}>3\） 对于第二个。对于足够小的时间步长Δ吨 可以看出

$$\sup_{0\len\leN}E\bigl[\bigl\vert y_{t_{n}}-y^{n}\bigr\vert\bigr]\leC\Delta t^{min\{K_{y}+1，4\}}$$

哪里 \（C>0\） 是一个常数，取决于 T型，（f），克 和的导数 （f） 和 克。

定理3

假设初始值满足

$$\textstyle\begin{cases}\max_{N-K_{z}+1\len\leN}E[\vertz_{t_N}}-z^{N}\vert]=\mathcal{O}（\Delta t^{K_{z}}），\\max_{N-K_{z}+1 \len\le N}E[\vert z_{t{N}}-z^}\vert]=\mathcal{0}（\ Delta t|3}），\结束{案例}$$

哪里 \（K_｛z｝=1,2,3\） 对于第一个方程和 \（K_{z}>3\） 对于第二个，定理中初值的条件 2已实现。对于足够小的时间步长Δ吨 可以看出

$$\sup_{0\len\leN}E\bigl[\bigl\vert z_{t_{n}}-z^{n}\bigr\vert\bigr]\leC\Delta t^{min\{K_{y}+1，K_{z}，3\}}$$

哪里 \（C>0\） 是一个常数，取决于 T型，（f），克 和的导数 （f） 和 克。

备注1

如果（f）不依赖于过程z（z），的最大收敛阶年过程是4和3z（z）过程；如果（f）取决于过程z（z），的最大收敛阶年和z（z）过程为3。

2.3全离散方案

让Δx个表示均匀分区中的步长d日-尺寸实轴，即。

$$\mathbb{R}^{\tilde{d}}=\Bigl\{x_{i}^{\ tilde{d\}}\big\vertx_{i}^{\\tilde{d}}\in\mathbb{R}，i\in\mathbb{Z}，x_{i}^\\tilde}{d}-x_{i}^{\tilde{d}}，\lim_{i\to\pm\infty}x_{i}^{\ tilde{d}}=\pm\infty\Bigr\}$$

哪里

$$\mathbb{R}^{\tilde{d}}=\mathbb{R}^{1}\times\mathbb}R}^}2}\times\cdots\times\mathbb{R{^{d}\quad\text{和}\qua2\tilde{d}=1，2，\ldots，d$$

让\（x{\mathbf{i}}=（x^{1}_{i{1}}，x^{2}_{i_{2}}，\ldot，x^{d}_{i{d}}）\）对于\（\mathbf{i}=（i{1}，i{2}，\ldots，i{d}）\in\mathbb{Z}^{d}\）.我们表示\（年）^{无}_{\mathbf{i}}，z^{无}_{\mathbf{i}}）\）近似于\（（y_｛t_｛n｝，x_｛\mathbf｛i｝｝｝，z_｛t_｛n｝，x_｛\mathbf｛i｝｝）\），给定随机变量\（年^{北}_{\mathbf{i}}，z^{北}_{\mathbf{i}}）\），\（l=0,1，点，K-1）具有\（K=\max\{K_{y}，K_{z}\}）.然后\（年）^{无}_{\mathbf{i}}，z^{无}_{\mathbf{i}}）\）可以为找到\（n=n-K，点，0）这样的话

$$\开始{对齐}（&y）^{无}_{\mathbf{i}}=\hat{电子}_{t_{n}}^{x_{mathbf{i}}\bigl[\hat{y}^{n+K_{y}}\bigr]+\Delta tK_{y}\sum_{j=1}^{K_y}}\gamma_{K_y}，j}^{K_{y{}}\hat{电子}_{t{n}}^{x{mathbf{i}}\bigl[f\bigl（t{n+j}，what{y}^n+j{，what{z}^{n+j}\bigr）\bigr]+\Delta tK{y}\gamma_{K{y{，0}^{K{y}}f\bigle（t}n}，y^{无}_{\mathbf{i}}，z^{无}_{\mathbf{i}\bigr），\\&z^{无}_{\mathbf{i}}=\Biggl（\hat｛E｝_{t{n}}^{x{mathbf{i}}\bigl[\hat{z}^{n+1}\bigr]+\sum_{j=1}^{K{z}}\gamma_{K{z}，j}^{1}\hat{电子}_{t{n}}^{x{mathbf{i}}\bigl[f\bigl（t{n+j}，what{y}^{n+j}，hat{z}^{n+j}\bigr）\Delta W{t{n+j}}^{top}\biger]-\sum_{j=1}^{K{z}}\gamma_{K{z}，j}^{1}{电子}_{t_{n}}^{x_{mathbf{i}}\bigl[\hat{z}^{n+j}\bigr]\Biggr）\big/\gamma_{K{z}，0}^{1}，\end{aligned}$$

(9)

哪里\（帽子{电子}_{t{n}}^{x{i}}[\cdot]\）用于表示条件期望的近似值。每个空间网格点的计算\（x{\mathbf{i}}\）在(9)每个时间层都是独立的\（t{n}\）因此，并行化策略与空间离散化完全相关，这将在以下章节中讨论。

条件期望中的函数包括d日-布朗运动的维数概率密度函数，可以选择例如高斯-埃尔米特求积规则，以仅用几个空间点实现高精度。条件期望可以通过巧妙的表格插值得到足够精确的近似值

$$\帽子{电子}_{t_{n}}^{x_{mathbf{i}}\bigl[\hat{y}^{n+k}\bigr]=\frac{1}{\pi^{\frac}{d}{2}}\sum_{\Lambda=1}^{L}\omega_{\Lambda}\hat}（x_{\mathbf}i}+\sqrt{2k\Delta t}a_{\Lambda}）$$

(10)

哪里\（{y}^{n+k}\）是在空间点处插值\（（（x_{\mathbf{i}}+\sqrt{2k\Delta t}a{\Lambda}）\）基于\（y^{n+k}\）值，\（（\omega{\Lambda}，a{\Lambeda}）\）对于\（\Lambda=（\Lambda_{1}，\Lambda _{2}，\ ldots，\ Lambda_{d}）\）是Hermite次多项式的权重和根我（请参见[12]),\（\omega{\Lambda}=\prod_{\tilde{d}=1}^{d}\omega_{\Lambda{\tilde{d}}\），\（a{\Lambda}=（a{Lambda{1}}、a{\Lambda}}，\ldots，a{\Lambda{d}}）和\（\sum_{\Lambda=1}^{L}=\sum_{\Lambda{1}=1，\ldots，\Lambda{d}=1}^}L，\ldot，L}\）以类似的方式，可以在(9).

三算法框架

根据第节中所示的数值算法。 2，数值求解BSDE的整个过程可以分为三个步骤。

• 1 构建时空离散域。

  我们划分时间段\（[0，T]\）进入之内N个时间步长使用\（增量t=t/N）即。，\（N+1）时间层和空间域\（\mathbb{R}^{d}\）使用步长Δx个（另请参阅第2.3). 我们使用截断域$[-16, 16 ]$和$[-8, 8 ]$对于网格空间，前者用于一维示例，后者用于二维示例。请注意，也可以使用较大的域，但近似值不会得到改进。也就是说，这些截断的区域对于我们的数值实验来说已经足够了。此外，为了平衡时间和空间方向上的误差，我们调整Δx个和Δ吨这样他们就满足了平等\（（\Δx）^｛r｝=（\Δt）^｛q+1｝\），其中\（q=\min\{K_{y}+1，K_{z}\}\）和对表示计算条件期望时用于生成非网格点的插值方法的全局误差。为了更好地说明，我们展示了以下领域的可视化\（K{y}=K{z}=1\），\（d=1）在图中1，包括非栅格点。

  图1

  

  时空域

• 2 计算 K（K） 初始解 \（\boldsymbol{K=\max\{K_{y}，K_{z}\}}\）。

  通常，只给出终端值，一个需要计算另一个\（K-1）初始值。为了获得这些初始值，我们首先\（K=1）并选择一个非常小的时间步长Δ吨。

• 三。 计算数值解 \（\粗体符号{（y_{0}^{0}，z_{0{^{0{）}\） 反向使用方程式(9)。

  请注意年这个过程是通过Picard迭代隐式完成的。

在步骤1中，一些生成的非网格点可能位于截断域之外（参见图1). 对于这些点，我们将边界上的值作为近似值（常数外推）。注意，也可以对外部的那些点进行外推，但需要更长的计算时间。如中所述(10)，我们使用插值来近似下列值年和z（z）在非网格点。插值时计算时间的缺点是在插值过程中找到新点的位置。一种自然搜索算法是循环遍历所有网格点，并找出该点所属的间隔。在最坏的情况下\（\mathcal{O}（M^{d}）\）需要工作。幸运的是，高斯-海米特求积的结构为非网格点创造了对称性。回想一下，每个新点都是作为\（X_{\lambda_{\tilde{d}}}=X_{i{\tilde{d}{+\sqrt{2k\Delta t}a{\lampda_{\ tilde{d}}}）这意味着\（\运算符名称｛int｝（\ frac｛X_｛\lambda_｛\tilde｛d｝｝-X_｛\main｝｝｛\Delta X｝）\）对于\（[x{\min}，x{\max}]\中的x{i{\tilde{d}}）和\（M-\运算符名称｛int｝（\frac｛X_｛\lambda_｛\tilde｛d｝｝｝-X_｛\main｝｝｛\Delta X｝）\）对于\（[x{\max}，x{\min}]\中的x{i{\tilde{d}}）给出了网格间隔的左边界\（X_{\lambda_{\tilde{d}}}）属于，具有\（\operatorname{int}（x）\）给出的整数部分x个。因此，算法中的这一步可以在\（\mathcal{O}（d）\）也就是说，没有for-loop。这种好处来自空间领域的一致性。这大大减少了总计算时间，这将在数值实验中得到证明。

在步骤2中，我们不考虑2K（K）(K（K）对于年和K（K）对于z（z）)每个新计算的插值，但只有2个。假设我们在时间层\（t_{n-K}\）.计算年和z（z）这个时间层上的值，需要计算条件期望K（K）时间层。在我们的数值实验中，我们考虑了1维和2维示例，对于具有大内存的设备（GPU），高维问题可以并行化。三次样条插值用于查找一维情况下的非网格值，而双三次插值用于查找二维情况下的双三次值。例如年过程是\（A_{y}\ in \mathbb{R}^{K\次（4^{d}\次M^{d{）}\），存储所有系数。当我们在时间层时\（t_{n-K-1}\），只考虑与先前计算值相对应的样条插值。然后，矩阵的列\（A_{y}\）向右移动+1，以删除最后一列，并在第一列中输入当前计算的系数。新的\（A_{y}\）用于当前步骤。遵循相同的程序，直到\（t_{0}\）这也减少了算法的工作量。

在最后一步（步骤3）中，我们考虑将条件期望的计算合并到(9)这对于减少并行实现中的计算时间非常重要。这将在下一节中详细介绍。

4并行算法

在本节中，我们将介绍前面章节中介绍的多步方案的并行化策略。基本上，我们直接将问题并行化，同时关注最佳的计算统一设备体系结构（CUDA）执行模型，即创建阵列以便对齐和合并访问，减少对全局内存的冗余访问，在需要时使用寄存器等。。

算法的第一步和第二步在主机中实现。第三步，与操作员一起执行\（\mathcal{O}（d）\）插值搜索过程中的工作以及该过程中系数移位的思想在GPU设备中得到了充分的实现。从调用(9)需要以下步骤来向后计算每个未知时间层上的近似值：

• 非网格点的生成 \（\boldsymbol{X_{\Lambda}=X_{\mathbf{i}}+\sqrt{2k\Delta t}a_{\Lambda}}\）。

  在均匀空间域中，非网格点只需要生成一次。为此，创建一个内核，每个线程在其中生成\（L^{d}\）每个空间方向的点（高斯端点）。这个内核在总算法时间中增加了微不足道的计算时间。

• 数值的计算 \（\hat{\boldsymbol{y}}\） 和 \（\hat{\boldsymbol{z}}\） 在非网格点。

  这是算法中最耗时的部分，它与为给定的插值函数找到相应的插值函数有关年和z（z）点，以便在非网格点中插值它们的值。对于一维情况，我们考虑了三次样条插值。自(9)涉及两个线性系统的解，双共轭梯度稳定(稳定双共轭梯度法) [28]由于矩阵是三对角的，所以使用了迭代方法。为此，我们考虑CUDA基本线性代数子程序(cuBLAS公司)和CUDA Sparse(cuSPARSE公司)图书馆[5]. 对于向量的内积、二次范数和加法，我们使用cuBLAS公司库。对于矩阵向量乘法，我们使用cuSPARSE公司由于系统矩阵的结构，使用压缩稀疏行格式的库。此外，我们创建了一个核来基于求解的系统计算样条系数。最后，创建了一个核来应用该算子进行插值搜索，以查找非网格点的值。请注意，每个线程都被指定查找\（m+m\乘以d\）值(米对于年和\（m\乘以d\）对于z（z）). 对于二维示例，我们考虑了双三次插值。我们需要为每个点计算16个系数。基于双三次插值思想，我们需要一阶导数和混合导数。使用四阶精度的有限差分格式（中心为内部点，向前和向后为边界点）进行近似。因此，在每个线程计算这些值的位置创建了一个内核。此外，为了找到16个系数，需要对每个点应用矩阵向量乘法。因此，每个线程使用另一个内核执行矩阵-向量乘法。最后，创建一个内核来应用插值搜索过程的操作符，以查找非网格点上的值，每个线程在这些点上进行计算\（m+m\乘以d\）值。

• 条件期望的计算。

  如上所述，我们合并了条件期望的计算。对于右手边的第一个条件期望(9)，我们创建一个内核，其中每个线程通过使用(10). 对于其他内核，我们将它们的计算（三个条件期望）合并到一个内核中，即\（帽子{电子}_{t{n}}^{x{mathbf{i}}[\hat{z}^{n+j}]\），\（帽子{电子}_{t{n}}^{x{\mathbf{i}}[f（t{n+j}，\hat{y}^{n+j}，\ hat{z}^{n+j}）]\）和\（帽子{电子}_{t{n}}^{x{\mathbf{i}}[f（t{n+j}，\hat{y}^{n+j}，\ hat{z}^{n+j}）\Delta W{t{n+j}}]\），用于\（j=1,2，\点，K\）。这减少了从全局内存多次访问数据的次数，并为线程提供了更多的工作，使其不会处于空闲状态。注意，一个线程计算\（2倍m\倍d+m\）值。

• 计算 z（z） 值。

  中的第二个方程(9)使用并计算每个线程\（m\乘以d\）值。

• 计算 年 值。

  中的第一个方程式(9)使用，每个线程计算米值，使用Picard迭代过程。

在下一节中，我们将给出一些数值示例，以展示GPU计算的加速计算及其在期权定价中的应用。

5数值结果

我们使用CUDA C编程实现了并行算法。将并行计算时间（仅一次运行）与CPU上的串行计算时间进行比较。此外，还计算了加速比。中央处理器（CPU）是Intel（R）Core（TM）i5-4670 3.40Ghz，具有4个内核，其中只有一个内核用于串行计算。GPU是NVIDIA GeForce 1070 Ti，总内存为8GB GDDR5。

我们选择Hermite多项式的次数\（长=32\）对于一维示例和\（L=8）对于二维的。对于Picard互动，我们选择\（p=30\）使用这些值，求积误差和迭代误差非常小，不会影响收敛速度。

示例1

考虑非线性BSDE[32]。

$$\textstyle\begin{cases}-dy{t}=（-y_{t}^{3}+\frac{5}{2}y_{t}^{2}-\frac{3}{2{y_{t}）\，dt-z{t}\，dW{t}，\\y_{t}=\frac{exp（W{t}+t）}{\exp（W{t}+t）+1}。\结束{cases}$$

分析解决方案为

$$\textstyle\begin{cases}y_{t}=\frac{\exp（W{t}+t）}{\ exp（W{t}+t）+1}，\\z{t}=\frac}\exp。\结束{cases}$$

精确解\（T=1\）是\（（y{0}，z{0}）=（\压裂{1}{2}，\压裂{1'{4}）\）。在表中三，我们说明了在插值非网格值时，对于高N个。请注意，计算时间以秒为单位。我们看到，如果没有操作符来查找非网格点在空间中的位置，计算的串行时间会非常长。

表3初步结果\（N=256），\（K{y}=K{z}=3\）例如1

在表中4，我们使用每个块256个线程显示结果\（K=K_｛y｝=K_｛z｝\），\（t_{0}=0\）和\（T=1\）如前所述，通过使用更高的步骤方案可以获得更好的准确性。对于\（K>3\），由于定理，误差适度减少2.时间层的值越高N个可以为GPU分配更多的工作，从而提高应用程序的速度。对于\（K>3\），因为空间数据点的数量来自M（M）都是一样的。这些在图中显示\（\log_{10}（|y_{0,0}-y_{0}^{0}|)\)，\（\log_{10}（|z_{0,0}-z_{0}^{0}|)\)并在以下方面加速N个对于\（K=1，\ldots，6\）在图中2.最高加速比（26×）是针对具有\（K=6\）和\（N=1024）。

图2

示例结果图1

表4实施例的结果1

示例2

考虑非线性BSDE[32]

$$\textstyle\begin{cases}-dy{t}=\frac{1}{2}（\exp（t^{2}）-4ty_{t} -3个\exp（t^{2} -年_{t} \exp（-t^{2}））+z_{t}^2}\exp。\结束{cases}$$

解析解为

$$\textstyle\boot｛cases｝y_｛t｝=\ln（\sin（W_｛t｝）+3）\exp（t^｛2｝），\\z_｛t｝=\exp（t^｛2｝）\frac｛\cos（W_｛t｝）｝｛\sin（W_｛t｝）+3｝。\结束{cases}$$

精确解\（T=1\）是\（（y{0}，z{0}）=（ln（3），frac{1}{3}）。每个块使用256个线程的结果\（K=K{y}=K{z}\），\（t_{0}=0\）和\（T=1\）如表所示5。我们得出的结论与前一示例中的结论相同，只是精确度降低了。这是由于定理中给出的收敛顺序三，最多减少到3，因为驱动程序功能取决于z（z）过程。此外，由于驱动程序功能更加复杂（即访问的数据越多，使用的功能单元也越多，等等），因此与前一示例相比，我们获得了更高的加速比。加速比为35倍。我们显示\（\log_{10}（|y_{0,0}-y_{0}^{0}|)\)，\（\log_{10}（|z_{0,0}-z_{0}^{0}|)\)并在以下方面加速N个对于\（K=1，\ldots，6\）如图三。

图3

示例结果图2

表5示例结果2

为了优化应用程序，我们使用了以下迭代方法：

• 1 将探查器应用于应用程序以收集信息

• 2 确定应用程序热点

• 三。 确定性能抑制因素

• 4 优化代码

• 5 重复上述步骤，直到达到所需性能

我们考虑了这个案件\（N=1024）和\（K{y}=K{z}=6\）为了使优化过程更加清晰，我们在每次优化迭代后显示了详细信息，如下所示。

• 在优化过程的第一次迭代中，我们使用nvprof公司（NVIDIA命令行分析器）。结果如表所示6（a） ●●●●。应用程序热点是\（nrm2\_\mathit{kernel}\）内核，它计算稳定双共轭梯度法算法。这已经进行了优化。因此，为了克服这个瓶颈，我们使用了点内核\（\tathit{dot}\_\mathit{kernel}\）.计算时间从8.04秒减少到0.86秒。第一次迭代后的新加速比由35倍变为57倍。

  表6示例迭代优化过程结果2

• 在第二次迭代中，应用程序的下一个瓶颈是为进程计算非网格值的内核年和z（z）(\（\mathit{sp}\_\mathit{inter}\_ \mathrm{non}\_\mathit{grid}\_d\_\mathit{no}\_\ mathit}for}\）)之后每一个时间层向后。内核的性能受到算术和内存操作延迟的限制。因此，我们考虑了循环交换和循环展开技术。这减少了相应内核和与之相关的其他内核的计算时间，如表所示6（b） 。通过这一点，我们减少了\（\mathit{sp}\_\mathit{inter}\_\ mathit}non}\_\fathit{grid}\_d\_\mathit{no}\_\tathit}for}\）从2.48秒减少到1.16秒。默认情况下，我们已将以下项的计算时间从2.28秒减少至1.46秒\（\mathit{calc}\_f\_\mathit{和}\_c\_\mathit{exp}\）（门派第三条核心4)因为我们需要改变非网格点的存储和访问方式，并减少\（\mathit{calc}\_c\_\mathit{exp}\d\）（计算条件期望）从0.22秒到0.04秒。新的加速比为69倍。从表中可以观察到6（c） 同样，应用程序瓶颈是同一个内核。因此，不值得进一步优化应用程序。

• 最后，为了提高并行度，我们将块维度从256个线程减少到128个。最后的加速比是70倍。

示例3

考虑一下Black-Scholes FBSDE[14]

$$\begin{aligned}\textstyle\begin}{cases}dS_{t}=\mu{t}S_{t{，dt+\sigma{t}S_{t｝，dW{t}，\quad S_{0}=x，\-dy{t}=-（r{t}y{t{t}+（\mu{t}-r{t}+\delta（t，S_{t}））\frac{z{t}{t}}{t}}）\，dt+z{t}\，dW{t}，\\y_{t}=（S_{T} -克 )^{+}. \结束｛cases｝\displaystyle\end｛aligned｝$$

对于常数参数（即。\（r_｛t｝=r\），\（\mu_{t}=\mu\），\（\sigma{t}=\sigma），\（\delta{t}=\delta\）)，看涨期权的解析解为

$$\开始{对齐}\textstyle\开始{cases}y_{t}=V（t，S_{t}）=S_{t{\exp（-\delta（t-t））N（d_{1}）-K\exp a，\\d_{1/2}=\frac{ln（\frac{S_{t}}{K}）+（r\pm\frac{\sigma^{2}}{2}）（t-t）}{\sigma\sqrt{t-t}}，\end{cases}\displaystyle\end{aligned}$$

哪里\（N（\cdot）\）是累积标准正态分布函数。在这个例子中，我们考虑\（T=0.33\），\（K=S_{0}=100\），\（r=0.03\），\（\mu=0.05\），\（δ=0.04\）和\（σ=0.2）（取自[32])用解决方案\（（y_{0}，z_{0}）\doteq（4.367110.0950）\）。请注意，终端条件对于z（z）过程。因此，对于接近执行价格的离散点K（K）（也称为货币区域）z（z）进程将在下一个时间层上导致较大错误。为了克服这个非光滑性问题，我们考虑平滑初始条件，参照克雷斯的方法[15]. 对于示例的前面部分三，我们使用了解析解

$$S_{t}=S_{0}\exp\biggl（\biggl（\mu-\frac{\sigma^{2}}{2}\biggr）t+\sigma W{t}\bigr）$$

为了确保统一的股价域，我们切换到对数股价域\（X_{t}=\ ln（S_{t}）\）。在表中7我们展示了使用这种转换的重要性。

表7初步结果\（N=256\），\（K{y}=K{z}=3\）例如三

每个块使用256个线程的结果\（K=K{y}=K{z}\），\（t_{0}=0\）和\（T=0.33\）如表所示8与前面的BSDE示例一样，对于考虑的最大步数，可以获得最高的精度K（K）以及时间层的数量N个，即6步方案\（N=256），其中我们也有12倍的最高加速。我们绘制了\（\log_｛10｝（|y_{0,0}-y_{0}^{0}|)\)，\（\log_{10}（|z_{0,0}-z_{0}^{0}|)\)并在以下方面加速N个对于\（K=1，\ldots，6\）在图中4。

图4

示例结果图三

表8示例结果三

与并行多步方案相比[32]英寸[13]，数值结果更准确，因为我们可以考虑\（K>3\）此外，我们优化了为Black-Scholes BSDE创建的内核\（N=256）和一个高K（K），即\（K=K{y}=K{z}=6\）优化迭代过程与示例中相同2.最终加速比为31×。请注意，此加速适用于256个时间层。在示例中2，我们针对1024个时间层进行了优化。高精度\（\mathcal{O}（10^{-12}）\）具有较高的值K（K）计算时间比6.47秒要少。

显然，在GPU上实现高维BSDE的并行化更有趣。

示例4

我们从二维BSDE开始[30]

$$\textstyle\begin{cases}-dy{t}=（y_{t}-z_{t} A类）\，dt-z{t}\，dW{t}，\\y_{t}=\sin（MW{t}+t），\end{cases}$$

哪里\（W{t}=（W{t}^{1}，W{t{^{2}）^{\top}\），\（z{t}=（z{t}^{1}，z{t{^{2}）\），\（A=（\frac｛1｝｛2｝，\frac｛1｝｛2｝）^｛\top｝\）和\（M=（1,1）\）。

解析解为

$$\textstyle\begin{cases}y_{t}=\sin（MW{t}+t），\\z{t}=（\cos（MW_at}+t），\cos。\结束{cases}$$

精确解\（T=1\）是\（（y{0}，（z{0}^{1}，z{0{2}））=（0，（1,1））\）。每个块使用256个线程的结果\（K=K{y}=K{z}\），\（t_{0}=0\）和\（T=1\）如表所示9。请注意\（|z）_{0,0}-z_{0}^{0}|\)在二维情况下由\（压裂{1}{2}（|z{0,0}^{1} -z_{0}^{0,1}|+|z{0,0}^{2} -z（-z）_{0}^{0,2}| )\).图\（\log_{10}（|y_{0,0}-y_{0}^{0}|)\)，\（\log_{10}（|z_{0,0}-z_{0}^{0}|)\)并在以下方面加速N个对于\（K=1，\ldots，6\）如图所示5。在减少误差方面，我们得出的结论与之前相同，但加速除外。在二维中不清楚后者在增长时的幅度增加K（K），因为我们考虑了最大时间层数\（N=64）。最高的加速比约为59倍，需要大约3 GB的内存。

图5

示例结果图4

表9实施例的结果4

由于GPU内存限制（8 GB），我们优化了以下情况\（N=64）和\（K{y}=K{z}=6\）在一维情况下，我们使用不同的核来计算非网格点、条件期望和进程值年和z（z）。在这里，由于内存限制，我们合并了这些内核，并将其命名为\（\mathit{main}\_\mathit{body}\）内核。在第一次迭代中，我们使用无价值专业人士结果见表10。应用程序热点是\（\mathit{main}\_\mathit{body}\）内核。内核的性能受到内存操作的限制，因为双三次样条系数的访问不是最优的。然而，我们无法优化这一部分，否则请求对齐和合并的内存操作会移动脊椎系数，线程访问错误的系数。此外，我们不能增加块维度，因为没有足够的资源。因此，最终加速约为59×。

表10主要内核性能示例4

最后，我们考虑零行权欧洲价差期权，这是一个二维问题。

示例5

零打击欧洲利差期权BSDE为

$$开始{对齐}\textstyle\begin{cases}dS_{t}=\mu S_{t{，dt+\sigma S_{tneneneep，dW{t}，\quad S_{0}=x，\\E[dW{t1}]=\rho d_{t}，\\-dy{t}=-（ry_{t}+z{t}A^{-1}M^{top}），dt+z_{t{}，dW}，\\y_{t}=（S_{t}^{1} -S型_{T} ^{2}）^{+}。\结束{cases}\displaystyle\end{aligned}$$

哪里\（S_{t}=（S_{t}^{1}，S_{t{^{2}）^{\top}\），\（\mu=（\mu{1}，\mu{2}）\），\（\西格玛=（\西格玛{1}，\西格马{2}）\），\（W{t}=（W{t}^{1}，W{t{^{2}）^{\top}\），\（z_｛t｝=（z_｛t｝^｛1｝，z_｛t｝^｛2｝）\），$\mathrm{<trans\; data-src="A">A类</trans>}<trans\; data-src="=">=</trans><trans\; data-src="(">(</trans>\begin{array}{cc}{\mathrm{<trans\; data-src="\sigma ">\sigma </trans>}}_{\mathrm{<trans\; data-src="1">1</trans>}}& \mathrm{<trans\; data-src="0">0</trans>}\\ \mathrm{<trans\; data-src="\rho ">\rho </trans>}{\mathrm{<trans\; data-src="\sigma ">\sigma </trans>}}_{\mathrm{<trans\; data-src="2">2</trans>}}& {\mathrm{<trans\; data-src="\sigma ">\sigma </trans>}}_{\mathrm{<trans\; data-src="2">2</trans>}}\sqrt{\mathrm{<trans\; data-src="1">1</trans>}<trans\; data-src="-">-</trans>{\mathrm{<trans\; data-src="\rho ">\rho </trans>}}^{\mathrm{<trans\; data-src="2">2</trans>}}}\end{array}<trans\; data-src=")">)</trans>$和\（M=（\mu{1}-r，\mu{2}-r）\）.解析解由Margrabe公式给出[18]

$$\开始{对齐}\textstyle\开始{cases}y_{t}=V（t，S_{t}）=S_{t{^1}N（d_{1}）-S{t}^2}N（e_{2}）2}}）\pm\frac{\tilde{\sigma}^{2}}{2}（t-t）}{\ tilde{\sigma}\sqrt{t-t}}，\end{cases}\displaystyle\end{aligned}$$

哪里\（\ tilde{\sigma}=\sqrt{\sigama{1}^{2}+\sigma{2}^{2\sigma{1}\sigma_2}\rho}\）和\（V S_{t}=（部分V}{部分S}{部分t}^{1}}S{t}^}，部分V}{\部分S{t{t}^{2}}S_{t{2}）^{top}）在这个例子中，我们考虑\（T=0.1），\（S_{0}^{1}=S_{0{2}=100\），\（r=0.05\），\（mu{1}=mu{2}=0.1\），\（σ{1}=0.25\），\（σ{2}=0.3\）和\（\rho=0.0\）用解决方案\（（y{0}，（z{0}^{1}，z{0{2}））（15.48076，（14.3510，-12.6779））结果见表11与示例中相同的参数4。最高的加速比约为64倍，这同样需要约3 GB的大内存。The plots of\（\log_{10}（|y_｛0,0｝-y_{0}^{0}|)\)，\（\log_｛10｝（|z_{0,0}-z_{0}^{0}|)\)并在以下方面加速N个对于\（K=1，\ldots，6\）在图中6。

图6

示例结果图5

表11示例结果5

请注意，加速比示例中的要快4因为我们有额外的正向SDE，更多的工作是从GPU线程执行的。但我们无法进一步优化，因为主要限制是GPU内存。在表中12我们展示了主要内核的性能。我们的结果表明，多步方案[27]GPU计算在二维情况下也表现得很好。当考虑单个GPU时，由于内存限制，增加维度变得困难。然而，这种方法可以应用于GPU集群，在该集群中实现更高加速的可能性是巨大的。这突出了在多步骤方案中实现集群GPU计算的事实[27]即使对于高维问题，也可以在较短的计算时间内提供非常精确的结果。

表12主要内核性能示例5

6结论

在这项工作中，我们将在[27]用于在GPU上求解BSDE。首先，我们分析了算法并提出了减少计算时间的方法。最重要的缩减工作是找到插值位置的优化操作。这对于减少计算时间至关重要。为了进一步加速，我们研究了如何在发现性能瓶颈并应用优化技术后优化应用程序。我们的数值结果表明，多步方案非常适合大规模并行GPU计算，并且对于期权定价及其风险管理等实时应用非常有效。我们的并行化策略也适用于其他多步方案，并使这些方案在实践中更加有用。使用稀疏网格和集群GPU计算解决高维问题是我们正在进行的工作的任务。

附录：参考方程的详细推导

我们推导了以下参考方程年和z（z）第节。 2.2。接收的参考方程年，我们需要获得适应性。因此，我们采用条件期望\（E_{t_{n}}^{x}[\cdot]\）在(6)并获得

$$y_{t_{n}}=E_{t_{n}}^{x}[y_{t_{n+k}}]+\int_{t_n}}^}n+k{}E_{t_n{n}{x}\bigl[f（s，y_{s}，z_{s{）\bigr]\，ds$$

(11)

其中使用了Itó积分的鞅性质。近似积分(11)，Teng等人[27]使用三次样条多项式来近似被积函数。基于支撑点\（（（t{n+j}，E_{t{n}}^{x}[f（t{n+j}，y{t{n+j}}，z{t{n+j{}）]），\（j=0，\ldots，K_{y}），我们有

$$\int_{t_{n}}^{t_{n+k}}E_{t_{n}}^{x}\bigl[f（s，y_{s}，z_{s{）\bigr]\，ds=\int_}_n}^{t{n+k}}\波浪线｛S｝_｛K_｛y｝｝^｛t_｛n｝，x｝（s）\，ds+R_｛y｝^｛n｝$$

其中三次样条插值为

$$\波浪线｛S｝_{K_{y}}^{t_{n}，x}（s）=\总和{j=0}^{K_{y} -1个}\波浪线{s}_{K_{y}}^{t_{n}，x，j}（s）$$

哪里

$$\波浪线{s}_{K{y}}^{t{n}，x，j}（s）=a{j}^{y}+b{j}^y}$$

具有

$$s\in[t_{n+j}，t_{n+j+1}]，\quad j=0，\ldots，K_{y} -1。 $$

显然，残差读数

$$R_{y}^{n}=\int_{t_{n}}^{t_{n+k}}\bigl（E_{t_{n}}^{x}\bigle[f（s，y_{s}，z_{s{）\bigr]-\tilde｛S｝_{K_{y}}^{t_{n}，x}\较大）\，ds$$

我们计算

$$\开始{对齐}\int_{t_{n}}^{t_{n+k}}\波浪线｛S｝_{K{y}}^{t{n}，x}（s）\，ds&=\int_{t{n}}^}{t{n+K}}\sum_{j=0}^{K_{y} -1个}\波浪线{s}_{K_{y}}^{t_{n}，x，j}（s）\，ds\\&=sum_{j=0}^{K_{y} -1个}\int_{t_{n}}^{t_{n+k}}\波浪线{s}_{K_{y}}^{t_{n}，x，j}（s）\，ds\\&=sum_{j=0}^{K_{y} -1个}\ int_{t{n+j}}^{t{n+j+1}}\波浪线{s}_{K_{y}}^{t_{n}，x，j}（s）\，ds.\end{aligned}$$

并获得参考方程年如第节所述。 2.2

$$y_{t_{n}}=E_{t_{n}}^{x}[y_{t_{n+k}}]+\sum_{j=0}^{k_{y} -1个}\biggl[a{j}^{y}\Delta t+\frac{b{j}^{y{}\ Delta t^{2}{2}+\frac{c{j}|y}\ Deltat^3}{3}+\frac{d_{j}#y}\ Delta t^4}{4}\biggr]+R{y}^n}$$

接收的参考方程z（z）过程，我们使用我而不是k个在(11)，将两边乘以\（增量W_{t_{n+l}}），取条件期望\（E_{t_{n}}^{x}[\cdot]\）以获得

$$0=E_{t_{n}}^{x}[y_{t_{n+l}}\增量W_{t_n+l}}]+\int_{t_}n}}^{t_n+l}E_{t_n}}^}x}\bigl[f（s，y_{s}，z_{s{）\Delta W_{sneneneep \bigr]\，ds-\int_t{n}^{n+l}}E_{t_{n}}^{x}[z_{s}]\，ds$$

此处使用Itóisometry。再次，使用三次样条插值和

$$开始{对齐}E_{t_{n}}^{x}[y_{t_{n+l}}\Delta W_{t_n+l}]=&\frac{1}{\sqrt{2\pi l\Delta t}}\int_{-\infty}^{\infty}u（t_{n+l{，x+v）\exp\biggl（-\frac}v^{2}{2l\Deltat}\biggr）\，dv，\\=&\frac{l\Delta-t}{\sqrt{2\pi l\Delta t}}\int_{-\infty}^{\infty}\frac}\partialu}{\partial x}（t_{n+l}，x+v）\exp\biggl（-\frac<v^{2}}{2l\Delta t}\biggr）\，dv，\\=&l\ Delta tE_{t_{n}}^{x}[z_{t_{n+l}}]，\end{对齐}$$

我们得到了z（z）过程

$$\开始{对齐}0={}&l\增量tE_{t_{n}}^{x}[z_{t_n+l}}]+\sum_{j=0}^{K_{z} -1个}\biggl[a{j}^{z{1}}\增量t+\frac{b_{j}^{z_{1}{\增量t^{2}}{2}+\frac{c_{j{}^{z{1}}\增量t{3}{3}+\frac{d_{j}^{z{1{}}\Delta t^4}{4}\biggr]\\&{}-\sum_{j=0}^{K_{z} -1个}\比格尔[a｛j｝^｛z_｛2｝｝\三角洲t+\三角洲t^｛z_｛2｝｝\三角洲t^｛2｝｝+\三角洲c^｛j｝^｛z_｛2｝｝\三角洲t^｛3｝｝｝+\三角洲t^｛4｝\比格尔]+R_｛z｝^｛n｝。\结束{对齐}$$