随机波动模型中期权定价的高阶ADI格式

（22）我们希望在紧凑的模具上找到uyy和uyy的二阶精确近似值，以找到uyy和uyy的四阶精确表达式。重新安排（20），我们得到：uyy=（vy）2β- κ（vy）α（θ）- vy+g.通过重复应用链式规则，uyy（xi，yj）和uyy-yy（xi，yj）的二阶精确近似式为：uyy（xi，yj）=（2（vyj）αkvyj- 2（vyj）αθαk+2（vyj）αkvyj（vyj）2βvyjδy0ui，j+（2（vyj）αkvyj- 2（vyj）2ββv- 2（vyj）αθkyj（vyj）2βvyjδyui，j+（vyj）2βδy0gi，j+O(y） ，（23）uyy-yy（xi，yj）=2（2（vyj）αkvyj- 2（vyj）2ββv- 2（vyj）αθkyj（vyj）4βvyj-4β（vyj）2βyj！δy0gi，j+（vyj）2βvyj（2（vyj）ααkv+4（vyj）αkv-2（vyj）ααθkyj+2（vyj）αkv）-2β（2（vyj）αkvyj- 2（vyj）αθαk+2（vyj）αkvyj（vyj）2βvyj-2（vyj）ααkvyj- 2（vyj）αθαk+2（vyj）αkvyj（vyj）2βvyj+（2（vyj）αkvyj- 2（vyj）2ββv- 2（vyj）αθkyj）（2（vyj）ααkvyj- 2（vyj）αθαk+2（vyj）αkvyj（vyj）4βvyj！δy0ui，j+2（vyj）αkvyj- 2（vyj）αθαk+2（vyj）αkvyj（vyj）2βvyj+（vyj）2βvyj（2（vyj）αkvyj+4（vyj）αkvyj- 4（vyj）2ββvyj- 2（vyj）αθαk- 2（vyj）αθk）-2β（2（vyj）αkvyj- 2（vyj）2ββv- 2（vyj）αθkyj（vyj）2βvyj-2（vyj）αkvyj- 2（vyj）2ββv- 2（vyj）αθkyj（vyj）2βvyj+（2（vyj）αkvyj- 2（vyj）2ββv- 2（vyj）αθkyj（vyj）4βvyj！δyui，j+（vyj）2βδygi，j+O(y） 。（24）其中δy0和δyde表示标准的二阶中心差分算子。分别将（23）和（24）代入（21）和（22）中，得到了紧集s tencil上uy（xi，yj）和uyy（xi，yj）的四阶精确近似值（此处未给出）。通过将这些四阶精确近似代入（20）并将u项和g项分别分离到左侧和右侧，我们得到了一个线性系统，它可以用矩阵表示为m:Ayu=bygw，其中u=（u2,2，u2,3，…，uN）-1米-1） ，g=（g2,2,g2,3，…，gN）-1米-1).

我们不在y方向施加任何边界条件，而是使用相同的格式离散边界网格点，并通过外推处理产生的鬼影点；关于边界条件的详细信息见第7节。系数矩阵Ay和BY是具有以下结构的块三对角矩阵：Ay=A1,1yA1,2y0 0A2,1yA2,2yA2,3y0 0。。。。。。。。。0安-3，N-4yAN-3，N-3yAN-3，N-2y0 AN-2，N-3yAN-2，N-2y,由=B1,1yB1,2y0 0B2,1yB2,2yB2,3y0 0。。。。。。。。。0亿-3，N-4yBN-3，N-3yBN-3，N-20亿日元-2，N-3yBN-2，N-2y,其中，每个Aj，jy=diag[ai，j]和Bj，jy=diag[bi，j]都是对角矩阵，这些对角线上的值如下所示：ai，j±1=2y（vyj）2β-12（vyj）2βvyj- 2（vyj）2ακvyj+2（vyj）2β+αακvyj- 2（vyj）2β+αβκvyj+4（vyj）2αθkvyj+2（vyj）4ββv- 2（vyj）2β+αθακvyj+2（vyj）2β+αθβκvyj+2（vyj）2β+ακvyj- 2（vyj）2αθκyj+（vyj）4ββv±-（vyj）ακvyj+（vyj）αθκ2vY-24yβ（vyj）- 2（vyj）2ακvyjy+（vyj）2β（vyj）αακvyjY- 4（vyj）2β+αβκvyjy+4（vyj）2αθακvyjY- 2（vyj）2ακvyjY- （vyj）2β+αθακvy+4（vyj）2β+αθαβκvy+（vyj）2β+ακvyjY- 4（vyj）2β+αβκvyjY- 2（vyj）2αθακyjy+2（vyj）2αθκvyjy+（vyj）2β+αθακvY!, （25）ai，j=6（vyj）2β+2- 2（vyj）2αkvyj+2（vyj）2β+αkvyj+2（vyj）4ββv- 2（vyj）2β+αβkvyj+4（vyj）2αθkvyj- 2（vyj）2β+αθαkvyj+（vyj）4ββv+2（vyj）2β+αθβkvyj+2（vyj）2β+αkvyj- 2（vyj）2αθkyj- 6（vyj）4β+2）, （26）bi，j±1=±-2（vyj）αkvyjY- 4（vyj）2ββvyjy+2（vyj）αθkvyjy24（vyj）2β+2y+，bi，j=。(27)5. 显式步骤的高阶方案ADI方案（8）的第一步和第四步仅对之前的近似值进行操作，以显式计算更新的近似值。这些步骤中的不同运算符采用（5）右侧的形式。对于混合导数项，不可能利用微分算子的结构在紧凑的计算模板上获得四阶近似。

因此，为了在（8）的显式步骤中保持方案的四阶精度，每个微分算子F、Fand和Fare中的导数近似于经典的四阶中心微分算子，这些算子在空间域中操作在一个较大的5×5模板上。对于F（u）=vyuxx- （vy）- r） ux，我们有以下方案：hvyUx+R-vyUxii，j=2r- vyj24十、-vyj24十、ui，j-2+8vyj- 16r24x+16vyj24十、ui，j-1.-30vyj24徐，j+16r- 8 vyj24x+16vyj24十、ui，j+1+vyj- 2r24十、-vyj24十、ui，j+2+O(x） 。对于F（u）=（vy）2βuyy+κ（vy）α（θ-vy）vuy，我们有：h（vy）2βUy+κ（vy）α（θ）- vy）vUyii，j=-κ（vyj）α（θ）- vyj）12vY-（vyj）2β24Y用户界面-2，j+8κ（vyj）α（θ）- vyj）12vy+16（vyj）2β24Y用户界面-1，j-30（vyj）2β24是的，j+-8κ（vyj）α（θ）- vyj）12vy+（vyj）2β24Yui+1，j+κ（vyj）α（θ）- vyj）12vY-（vyj）2β24Yui+2，j+O(y） 。最后，对于混合d导数项F=ρ（vy）β+uxy，使用以下计算模板：hρ（vy）β+U十、yii，j=-64ρ（vyj）β+144十、是的-1，j-1+64ρ（vyj）β+144十、是的-1，j+1+64ρ（vyj）β+144十、yui+1，j-1.-64ρ（vyj）β+144十、yui+1，j+1-ρ（vyj）β+144十、是的-2，j-2+8ρ（vyj）β+144十、是的-2，j-1.-8ρ（vyj）β+144十、是的-2，j+1+ρ（vyj）β+144十、是的-2，j+2+8ρ（vyj）β+144十、是的-1，j-2.-8ρ（vyj）β+144十、是的-1，j+2-8ρ（vyj）β+144十、yui+1，j-2+8ρ（vyj）β+144十、yui+1，j+2+ρ（vyj）β+144十、yui+2，j-2.-8ρ（vyj）β+144十、yui+2，j-1+8ρ（vyj）β+144十、yui+2，j+1-ρ（vyj）β+144十、yui+2，j+2+O(十、y） 。使用这些四阶近似，可以直接计算（8）中的第一步和第四步。ADI格式的每个解的空间边界值由边界条件确定，空间域中的所有剩余点都需要计算模板。对于显式步骤，当我们希望在沿空间域内网格边界的任意点近似微分算子F（u）时，5×5点计算模板超过了空间边界。

例如，如果我们希望评估F（u2,2），我们将要求来自空间域之外的重影点的贡献，如图1中的要点所示u4,1u4,2u4,3u4,4ou3,1u3,2u3,3u3,4ou2,1u2,2u2,3u2,4ou1,1u1,2u1,3u1,4⊙ o o o o图1：示例：使用计算域左下角的5×5点计算模板计算F（u2,2）；计算域外部的重影点（从域内部外推值）用项目符号（o，o,⊙), 边界上的网格点以罗马文字设置。我们从网格点u（xi，yj）推断信息，其中i=1，M- 1，j=1，N- 1.确定这些g主点的值，以评估沿空间域内网格边界的任意点处的微分算子F（u）。为了计算这些重影点的值，我们对三种情况使用以下五点外推公式：x=L边界（o）：ui，0=5ui，1- 10ui，2+10ui，3- 5ui，4+ui，5+O(x） ，y=l边界(o) : u0，j=5u1，j- 10u2，j+10u3，j- 5u4，j+u5，j+O(y） ，x=L，y=L角(⊙) : u0,0=5u1,1- 10u2,2+10u3,3- 5u4,4+u5,5+O(十、y） 。在x=K和y=K边界以及其余三个角s处的外推是类似处理的。6.从给定Un开始求解高阶ADI方案-1，ADI格式（8）包括六个近似步骤，以获得下一时间级别的解。使用第5节中导出的5×5点计算模板，可以显式求解第一个近似值。我们的解的第二个近似值，用Y表示，必须隐式地求解：Y=Y+φt（F（Y）- F（联合国）-1)) <==> F（Y）- 联合国-1) =φt（Y）- Y） 。（28）我们应用第4节中建立的四阶紧致格式来解（28）。

对于mwe获得的矩阵（Y- 联合国-1） =Bxφt（Y）- Y）+ d、 收集左手边未知的Y形项和已知的Y形项，Un形项-1和d在右边我们得到（Bx- φ税）Y=BxY- φ塔逊-1.- φtd。为了解决这个问题，我们将三对角矩阵（Bx）求逆- φ税）。对于ADI格式的第三步，我们进行类似的处理，并使用第4节中介绍的高阶紧致格式来简单地求解。ADI方案的第四步、第五步和第六步分别类似于第一步、第二步和第三步。注意，矩阵（Bx- φ税务）在方案（8）中出现两次，分别出现在第二步和第五步。同样地- φ（泰）在第三和第六步中出现了。因此，使用LU分解，在方案（8）的每个时间步中只需要两个矩阵逆。此外，由于偏微分方程（5）中的系数不依赖于时间，因此矩阵是常数，因此在进行时间迭代之前，不能对它们进行LU分解，以获得高效的算法。将第4节和第5节中介绍的四阶空间离散化与二阶时间分裂（8）相结合，得到了一个高阶ADI格式，在spa ce中的一致性阶数为2和4。7.在x=Land x=Kwe imposeu（L，yj，τk）=1的情况下，狄里克莱条件的边界条件- erτ+L，j=1，2，N、 k=1，2，u（K，yj，τK）=0，j=1，2，N、 k=1，2。使用齐次Neumann条件（6d）和（6e），它们在极限y内正确→ ∞还有y→ 在（有限）边界y=L>0和y=Kwould处分别为0，将导致沿这些边界的误差。

因此，我们不在这两个边界处施加任何边界条件，而是使用内部的计算模板离散偏导数。边界上的未知值是通过从内部的值外推来设置的。这引入了一个数值误差，需要考虑的是，外推的阶数应足够高，以免影响精度的总体阶数。我们参考Gustafsson[18]来讨论逼近阶对全局收敛速度的影响。我们使用以下公式：uki，1=5uki，2- 10uki，3+10uki，4- 5uki，5+uki，6+O(y） ，uki，N=5uki，N-1.- 10uki，N-北卡罗来纳州2+10uki-3.- 5uki，N-4+uki，N-5+O(y） .8。数值实验在本节中，我们报告数值实验的结果。我们估计了高阶ADI格式的数值收敛阶，然后进行了附加实验以验证其稳定性。8.1. 数值收敛我们对高阶ADI格式的收敛阶进行了数值研究。由于期权定价问题的初始条件，即支付函数V（S，σ，T）在S=E时是非光滑的，我们通常无法观察到高阶收敛[28]。平滑初始条件的一种简单方法是选择网格，使初始条件的非平滑点不是mes h的点。这种网格的构造总是可以以简单的方式进行。按照这种方法，我们的方案可以直接考虑非光滑的支付，并且我们观察到了高阶数值收敛。或者，可以使用合适的平滑算子来实现类似的效果，参见[28,14]。对于Conv-nie nc e，我们选择大小相等的空间步长h=x=y、 在水平和垂直方向创建e Venlyspace网格。

我们在（8）中设置参数φ=0.5，在所有数值中为101001502002500.050.10.150.20.250102030405060σSV（S，σ，0.5）。图2：欧洲看跌期权价格的数值解，T=0.5参数值执行价格E=100到期时间T=0.5利率r=0.05波动率v=0.1平均反转速度κ=2波动率长期平均值θ=0.1相关系数ρ=-0.5随机网格比γ=0.5随机波动率漂移参数α=0随机波动率扩散参数β=0.5表1：数值实验的默认输入参数-4.-3.5-3.-2.5-2.-1.5-14-13-12-11-10-9-8.-7.-6.-5.-4对数（h）对数（l）∞ 误差）ρ=-0.5毫-阿迪- 顺序：3.89ρ=-0.5秒-订购ADI- 顺序：1.97ρ=-0.2毫-阿迪- 顺序：3.92ρ=-0.2秒-订购ADI- 顺序：1.97ρ=0.1 HO-阿迪- 顺序：3.91ρ=0.1秒-订购ADI- 订单：1.96图3:l∞-在不同相关参数ρ实验值下，高阶ADI格式与标准二阶空间ADI格式的误差比较。图2显示了在T=0.5时，使用表1中的参数得出的欧元期权价格的数值解。我们计算了l范数误差ε和最大范数误差ε∞与细网格上的数值参考解相关的数值解。我们确定抛物线网格比γ=t/hto是抛物线偏微分方程（5）的自然常数。然后，对于一些表示常数的m和C，我们期望这些误差渐近收敛为ε=chm。这意味着ln（ε）=ln（C）+mln（h）。因此，ε对h的双对数图应与斜率为m的直线相切。这为实验确定方案的顺序提供了一种方法。我们期望在空间中观察到近似有序（h）的数值收敛速度。

为了进行比较，我们使用基于（8）的标准二阶ADI方案，结合空间中的二阶中心差分离散，进行了额外的实验。图3显示了l∞-错误与空间步长h=y=x、 我们观察到，数值收敛阶与格式的理论阶非常吻合。在所有情况下，对于给定的网格宽度h，高阶ADI方案优于标准的二阶ADI方案。或者换句话说，为了实现选定的精度水平，我们可以对高阶ADI方案使用同轴网格，而不是标准的二阶方案，后者可以解算更小的线性系统，因此计算效率更高。8.2. 数值稳定性分析在本节中，我们将研究时间步长的选择是否存在任何稳定性限制t用于高阶ADI方案。与标准二阶格式不同，algebraichγ0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.20.20.30.40.50.60.70.80.910.511.522.533.544.5x 10-3hγ0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.20.20.30.40.50.60.70.80.910.511.522.533.544.5x 10-3图4：ρ=0（左）和ρ=-抛物线网格比γ=0.5（右）τ/手网格宽度H高阶紧致格式的数值稳定性分析的复杂性非常高，因为已建立的稳定性概念意味着高阶紧致格式的强大代数问题。因此，文献[14,9,16]中关于高阶紧致格式的稳定性结果很少。

这一点在更高的空间维度中更为明显，因为大多数现有的高阶紧致格式分析稳定性结果的研究仅限于一维分析环境。对于具有混合导数项和常数系数的微分方程（无对流），ADI方法（8）的稳定性分析（在空间[22]中采用标准二阶离散）表明，它是无条件稳定的。[22]中的分析基于对简化线性测试方程稳定性的研究，这意味着假设所有涉及的离散矩阵都是标准的和协同计算的。高阶紧致s模式的离散化矩阵通常不满足这些假设，而且，在目前非常数系数的情况下，情况更复杂。因此，彻底的稳定性分析超出了本文件的范围。相反，为了验证该方案的稳定性，我们进行了额外的数值稳定性测试。我们注意到，在我们的数值实验中，我们始终观察到s表的行为。我们计算抛物线网格比γ=t/手动调整网格宽度h。绘制平面上的相关l-范数误差应允许我们根据γ或高单元雷诺数（大h）发生的振荡来检测稳定性限制。[9,10]中也使用了这种数值稳定性研究方法。我们使用表1中的参数给出了欧洲看跌期权的结果。对于我们的稳定地块，我们使用γ=k/10和k=2，10和空间网格点的降序序列。图4显示了相关参数ρ=0和ρ=-0.5.

我们观察到，抛物线网格比γ对l-误差的影响很小，相对误差不超过5×10-3是两个稳定性图的值。我们可以推断，对于这两种情况，γ似乎都不存在稳定条件。随着h值的增加，也会导致更高的c ell Reynolds数，误差会逐渐增大，数值解不会出现振荡。这些观察结果通过改变抛物线γ=t/h0。2 0.4 0.6 0.8 1.0HO-ADI l-误差3.8871 3.8870 3.8868 3.8866 3.8 864标准ADI l-误差2.4521 2.4519 2.4517 2.4514 2.4 510HO-ADI l∞-错误3.8960 3.8961 3.8961 3.8962 3.8964标准ADI l∞-误差1.9744 1.9744 1.9744 1.9743 1.9742表2：变抛物线网格比γ=t/hmesh比值γ。表2中报告的数值收敛阶表明，高阶格式的数值收敛阶在l范数和l∞-范数为4，与抛物线网格比γ无关。结论通过将空间中的四阶（紧致和非紧致）有限差分格式与Hundsdorfer和Verwer的二阶ADI时间步格式相结合，我们构造了一种新的数值方法来求解随机波动率模型的期权定价问题。使用Heston随机波动率模型和一般参数来近似欧佩恩看跌期权价格的数值实验证实了该方案在空间和时间上的数值收敛性，而抛物线网格比率范围的结果表明无条件稳定。确认TBD确认Leve rhulme信托研究项目“高阶非线性偏微分方程的新型离散化”（RPG-2015-69）的部分支持。