块三角矩阵指数的增量计算

直到终止isO（Pnk=0kb）=o（nb）的操作数的总体复杂度界，它与仅对Gn应用缩放和平方的复杂度界相匹配∈ R（n+1）b×（n+1）b，也就是O（（nb））。总之，算法5计算Fnby所需的操作数与仅将sa me缩放和平方设置应用于计算exp（Gn）大致相同，而算法5以增量方式显示所有指数exp（G），exp（Gn）在迭代过程中，满足我们在引言中概述的需求。2.4自适应缩放在算法4和5中，我们假设缩放功率s作为输入参数给出，并且通过计算exp（G），…，它是固定的，exp（Gn）。这与缩放和平方方法中通常使用的方法不同，参见第2.1节。一方面，s必须足够大，以使rk，m（2-sGl）≈exp（2-sGl），对于0≤ l≤ n、 另一方面，如果s选择得太大，则rk，m（2）的估值-sGl）可能因尺寸过大而变得不准确。因此，如果s是固定的，并且范数kGlk随着l的增加而增长，正如人们通常所期望的那样，则无法保证所有l的精确近似。算法6 exp（G），exp（G）。使用自适应缩放输入：Pad'e近似参数k，m，范数界θ。输出：F≈ exp（G），F≈ exp（G）。1： s← max{0，log（kGk）}2：计算融合缩放和平方，存储算法43的中间产物：对于l=1，2。do4：如果kGlk>θ，则5：重新划分Gn=^Gn-（10）6中的las：使用^Gn重新启动算法-l、 7：结束if8：从Fl计算Fl-1使用算法49：如果满足终止标准，则10：返回11：结束if12：结束最短的校准和平方设计，因此根据输入矩阵的范数选择s【11，7，9】。

例如，在[9]中描述的Higham算法中，它是满足k2的最小整数-sGlk公司≤ θ≈ 5.37。。。。（9） 为了将我们的增量评估技术与这种缩放和平方设计相结合，因此必须在评估过程中动态选择缩放功率。假设s满足步骤1中的标准（9）- 1，但不是在步骤l。然后，我们暗示放弃算法4中的所有累积数据结构，增加以匹配Gl的界限（9），并使用重新分区输入矩阵重新启动算法5=G0,0···G0，lG0，l+1···G0，n。。。。。。。。。。。。Gl、lGl、l+1··Gl、nGl+1、l+1··Gl+1、n。。。。。。Gn，n=^G0,0^G0,1···G0，n-l^G1,1···G1，n-l、 。。。。。。^Gn-l、 n个-l|{z}=：^Gn-l、 （10）算法6总结了该过程。结果表明，这种重新启动过程引起的计算超负荷是适度的。在第2.3节复杂性讨论中引入的符号中，通过Higham的缩放和平方方法计算exp（Gn）的操作数为O（log（kGnk）（nb））。由于算法6中最多有log（kGnk）res tartsin，因此增量计算所有指数exp（G），…，的运算总数，exp（Gn）可以由O（log（kGnk）（nb））中的函数限定。我们在第4节中评估了算法6的实际性能。在我们的期权定价应用中，结果表明，矩阵的范数不会显著增长（见第3.2节和第3.3节），即使比例因子固定，也可以计算出所有矩阵指数的相当精确的近似值（见第4.2节）。3多项式模型中的期权定价本节的主要目的是解释某些期权定价技术如何要求对块三角矩阵的矩阵指数进行顺序计算。

描述必须更加简短；更多详情，请参考教科书【5】。因为我们在初始时间t=0时进行评估，所以在时间τ>0时某个期权的价格有助于计算公式的表达式-rτE[f（Xτ）]，（11），其中（X）0≤t型≤τ是一个d维随机过程，用于模拟金融资产在时间间隔内的价格l[0，τ]，f:Rd→ R是所谓的支付函数，R代表一个固定利率。在下面，我们考虑由SDE描述的随机过程，其形式为dxt=b（Xt）dt+∑（Xt）dWt，（12），其中W表示d维布朗运动，b:Rd7→ Rd，∑：Rd7→ Rd×d.3.1多项式扩散模型在过去几年中，多项式扩散模型已成为金融应用（包括期权定价）中的通用工具。在下文中，我们提供了一个简短的总结，并参考菲利波维奇和拉尔森的论文[6]，了解数学基础。对于多项式微分，假设（12）中的向量B和矩阵a的系数：=∑∑TsatisfyAij∈ Pol（Rd），bi∈ i的Pol（Rd），j=1，d、 （13）这里，Poln（Rd）表示总次数最多为n的d元多项式集，即Poln（Rd）：=nX0≤|k级|≤nαkxk | x∈ Rd，αk∈ Ro，其中我们使用多索引符号：k=（k，…，kd）∈ Nd，| k |：=k+···+KD和xk：=xk。xkdd。在下文中，Pol（Rd）表示Rd上所有多变量多项式的集合。与A和b相关，我们定义了部分微分算子G byGf=Tr（Af） +英国电信f、 （14）表示（12）的所谓生成器，请参见[12]。可以直接验证（13）意味着Poln（Rd）在G下对于任何n是不变的∈ N、 即GPoln（Rd） 波兰（Rd）。（15） 备注3.1。在许多应用中，人们对位于状态空间E上的（12）的解感兴趣 Rdto输入，例如，非负性。

这个问题在文献[6]中得到了广泛的研究，其中（12）o在几种类型的状态空间E上解的存在唯一性 图中显示了大型类A和b的Rdand。现在，让我们确定Poln（Rd）的多项式基础Hn={h，…，Hn}，其中n=dim Poln（Rd）=n+dn, 和w riteHn（x）=（h（x），hN（x））T。设gn表示关于Poln（Rd）的线性算子的H的矩阵表示。根据定义，Gp（x）=Hn（x）TGn~p。对于任何p∈ 坐标向量为p的Poln（Rd）∈ Rn相对于Hn。B y定理3.1在[6]中，相应的多项式矩可以从e[p（Xτ）]=Hn（X）TeτGn~p计算得出。（16）上述设置对应于引言中描述的场景。我们有一个子空间序列pol（Rd） Pol（Rd） 波尔（路） ···  Pol（Rd）和多项式保持性（15）意味着矩阵r e表示的是大小为1，d，1+d, . . . ,n+d- 1n.在本节的其余部分中，我们将介绍两种不同的定价技术，它们需要增量计算形式（16）的多项式矩。3.2 Jacobi模型基于矩的期权定价Jacobi随机波动率模型是多项式微分模型的特例，其特征是SDEdYt=（r- Vt/2）dt+ρpQ（Vt）dW1t+pVt- ρQ（Vt）dW2t，dVt=κ（θ- Vt）dt+σpQ（Vt）dW1t，其中q（v）=（v- vmin）（vmax- 五）(√vmax（最大值）-√vmin），对于某些0≤ vmin<vmax。这里，w1和w2是独立的标准布朗运动，模型参数满足条件κ≥ 0，θ∈ [vmin，vmax]，σ>0，r≥ 0，ρ∈ [-1，1]。Ackerer等人[1]在他们的论文中，在期权定价的背景下使用了该模型，其中资产的价格由St:=EY和VT表示平方随机波动率。

下面，我们简要介绍他们提出的定价技术，并解释其如何涉及多项式矩的增量计算。在具有欧洲索赔贴现支付函数f的Jac obi模型下，初始时间t=0的期权价格（11）可表示为xn≥0fnln，（17）其中{fn，n≥ 0}是f和{ln，n的傅里叶系数≥ 0}是hermitoments。如【1】所述，傅立叶系数可以递归计算。使用（16）计算厄米矩。具体来说，考虑Poln（R）的单项式基础：Hn（y，v）：=（1，y，v，y，yv，v，…，yn，yn-1v，vn）T.（18）Thenln=Hn（Y，V）TeτGn~Hn，（19），其中~Hn包含关于（18）的坐标√nhn公司y- uwσw,用实参数σw，uwand表示第n个斜接多项式hn。在限定数量的条款后截断总和（17），可以获得期权价格的近似值。算法7描述了一种基于求和绝对值选择截断的启发式方法，使用算法5增量计算所需的矩。算法7 Jac-obi随机波动率模型下欧洲所有期权的期权定价输入：模型和支付参数，容差输出：近似期权价格1：n=02：计算l，f；设置价格=lf。3： 而| lnfn |>·Price do4:n=n+15：使用算法4.6计算经验（τGn）：计算埃尔米特矩lnusing（19）。7： 计算[1]中所述的傅立叶系数fn。8： 价格=价格+lnfn；9： 最后，第2节讨论了Gnis的范数估计，它有助于在定标和平方方法中先验地选择S平方参数。在考虑引理3.2的情况下，以下引理为模型提供了这样的估计。设Gnbe（14）中定义的算子G相对于Poln（R）的基（18）的矩阵表示。

定义α：=σ（1+vminvmax+vmax+vmin）2(√vmax（最大值）-√vmin）。然后，以n（r+κ+κθ）为界的Gnis的矩阵1-范数- σα）+n（1+|ρ|α+2σα）。证据雅可比模型中的算符G的形式为gf（y，v）=Tr（A（v）f（y，v））+b（v）f（y，v），其中b（v）=r- v/2κ（θ- 五）, A（v）=vρσQ（v）ρσQ（v）σQ（v）.设置S：=(√vmax（最大值）-√vmin），我们考虑生成器G对基本元素ypvq的作用：Gypvq=yp-2vq+1pp- 1.- yp公司-1vq+1p+qρσS+ yp公司-1vqpr+qρσvmax+vminS- yp公司-1vq-1pqρσvmaxvminS- ypvqqκ+q- 1σS- ypvq公司-2qq- 1σvmaxvminS+ypvq-1季度κθ+q- 1σvmax+Vmin.对于Gn的矩阵1-范数，需要确定（p，q）的值∈ M：={（p，q）∈ N×N | p+q≤ n}，其中Gypvq坐标向量的1-范数变得最大。考虑到相关模型参数的非负性，并将ρ替换为|ρ|，我们得到的上界如下：pp- 1+p+q |ρ|σS+ pr+q |ρ|σvmax+vminS+pq |ρ|σvmaxvminS+qκ+q- 1σS+ qq- 1σvmaxvminS+qκθ+q- 1σvmax+Vmin=pr+qκ（θ+1）+p+2pq |ρ|α+q（q- 1） σα≤n（r+κ+κθ）+n+2pq |ρ|α+n（n- 1） σα。这就完成了证明，注意到M上pq的最大值是以M上的n/4为界的。引理3.2的结果预测，一般来说，Gn的范数增长是相等的。该预测在数值上符合实际相关的参数设置。3.3 Heston模型基于矩的期权定价Heston模型是多项式微分模型的另一种特殊情况，其特点是SDEdYt=（r- Vt/2）dt+ρpVtdW1t+pVtp1- ρdW2t，dVt=κ（θ- Vt）dt+σpVtdW1t，模型参数满足条件κ≥ 0，θ≥ 0，σ>0，r≥ 0，ρ∈[-1，1]。如前所述，资产价格通过St:=eYt建模，而VT表示平方s到仓促波动率。Lasserre et a l【10】开发了一种基于动量和半有限规划（SDP）的通用期权定价技术。

在下文中，我们简要解释了主要步骤以及其中需要增量计算力矩的内容。因此，我们仅限于Hes-ton模型和Europeancall期权的具体情况。考虑支付函数f（y）：=（ey- eK）+对于特定的日志，s trike valueK。设ν（dy）为Yτ-随机变量（Yτ，Vτ）联合分布的边缘分布。将重新严格的度量sν和ν定义为ν=ν|(-∞,K] 和ν=ν|[K，∞). 通过用n项后截断的泰勒级数近似付息函数中的指数，期权价格（11）可以写成ν和ν中的某个线性函数L，即e[f（Yτ）]=L（n，ν，···，νn，ν，··，νn），其中νm表示第i个测度的第m个矩。然后，可以通过解决最优化问题SDPN来计算期权价格的上下限：=最小/最大L（n，ν，···，νn，ν，···，νn）受νj+νj=νj，j=0，···，nν是(-∞, K] ，ν是[K]上的Borel测度，∞).（20） 当通过moment和localizingmatrices将最后两个条件写入（20）时，会出现两个SDP，对应于所谓的截断Stieltjes矩问题。公式（16）在此设置中用于计算矩νj。适当增加松弛顺序可以让我们找到更清晰的界限（这很简单，因为增加n会增加更多约束）。当边界充分闭合时，可以停止旋转。

算法8总结了最终的定价算法。算法8基于SDP和动量松弛的欧式期权定价输入：模型和支付参数，容差输出：近似期权价格1：n=1，差距=12：当差距>do3：计算经验（τGn）使用算法44：计算n阶矩使用（16）5：求解相应的SDPN以获得上下限nd6：差距=|上限- LowerBound | 7：n=n+18：end，而下面的引理将引理3.2的结果扩展到了Heston模型。引理3.3。设Gn是上述关于Poln（R）的基（18）引入的算子G的矩阵表示。然后，Gnis的矩阵1-范数有界于n（r+κ+κθ-σ） +n（1+|ρ|σ+σ）。证据与引理3.2.4的证明类似，我们在Matlab中实现了本文中描述的算法，并将其与Higha m的缩放和平方方法进行了比较，Higha m的缩放和平方方法通常采用13度的径向Pad'e近似，以下称为“expm”。块三角矩阵算法、算法5（固定缩放参数）和算法6（自适应缩放参数）的实现，基于相同的缩放和平方设计，以下称为“incexpm”。所有实验均在标准笔记本电脑上运行（Intel Core i5，2核，256kB/4MBL2/L3缓存）使用单个c计算线程。0 500 1000 1500 2000 2500矩阵大小-2-1自适应ss=6s=12expm0 500 1000 2000 2500矩阵大小-15-14-13-12自适应ss=6s=12图1：随机块三角矩阵的incexpm和expm比较。左：计算前导部分的累积运行时间。右：incexpm w.r的相对误差。t、 扩展。4.1随机块三角矩阵我们首先评估随机生成的块上三角矩阵Gn的运行时间和准确性∈ R2491×2491。

有46个对角块，大小在n 20和80之间变化。矩阵的谱包含在区间内[-80，-0.5]和条件良好的特征基X（κ（X）≈ 10 0）。图1（左）显示了所有前导指数增量计算的挂钟时间。具体而言，如果0≤ l≤ n、 每个数据点显示计算l+1矩阵指数exp（G），exp（G），…，所需的时间vs.dl=b+······+bl，exp（Gl）当使用oexpm时（只需将其分别应用于每个矩阵）；o采用算法6中的自适应缩放策略的incexpm；o固定缩放功率为6的incexpm（expm用于G的缩放）；o具有固定缩放功率12的incexpm（expm用于Gn的缩放）。正如所料，incexpm比单纯地将expm分别应用于每个矩阵要快得多；表1还显示了l=n的总时间。作为参考，我们注意到，Matlab的expm的运行时间仅应用了最终矩阵Gnis13.65s，这非常接近于具有缩放参数se t至12的incexpm的运行时间（有关渐近复杂性的讨论，请参见第2.3节）。事实上，仔细观察incexpm的运行时文件可以发现，更复杂的数据结构所导致的计算开销在平方阶段通过利用块三角矩阵结构得到了很大程度的补偿，而Matlab的sexpm并没有从中自动实现。还值得注意的是，自适应缩放策略的运行时间大约只有使用固定缩放参数6运行算法的运行时间的两倍，尽管其渐进复杂性更差。图1右侧显示了由incexpm获得的近似值的精度。

我们假设expm为参考，并测量这两种近似值之间的相对距离，即kexpm（Gl）- incexpm（Gl）kFkexpm（Gl）kF，表1：在大小为2491的随机块三角形r矩阵上，expm和incexpm获得的运行时间和相对误差。算法时间相对值。er ROREXM 163.60incexpm（自适应）20.01 3.27e-15incexpm（s=6）9.85 2.48e-13incexpm（s=12）13.70 6.17e-140 500 1000 1500 2000矩阵大小-3-2-1自适应ss=7expm0 500 1000 2000矩阵大小-15-14自适应ss=7图2：算法7中Jacobi模型上下文中产生的块上三角矩阵的incexpm和expm的比较。左：计算前导部分的累积运行时间。右：incexpm w.r.t.expm的相对错误或。在每次迭代时，l（图1中将小于机器精度的数量设置为u，用于绘图）。人们注意到，在整个计算序列中，自适应策略的近似值保持在expm附近。对于该策略，观测到的误差下降到u与算法6中的重新启动相关；这一步的近似值与expm的近似值完全相同。即使对于固定的标度参数6和12，获得的近似值也相当准确。4.2期权定价的应用我们现在显示了使用算法7计算期权价格的结果，参数集为sv=0.04，x=0，σw=0.5，uw=0，κ=0.5，θ=0.04，σ=0.15，ρ=-0.5，vmin=0.01，vmax=1，r=0，τ=1/4，k=对数（1.1）。我们使用公差=10-3用于停止算法7。我们探讨了使用不同的算法计算算法7第5行中的矩阵exp单数：具有自适应缩放的incexpm、具有固定缩放参数s=7的incexpm（对应于n=60的引理3.2的上界）和expm。与图1类似，观察到的累积运行时间和错误如图2所示。