高阶导数的振动和自动微分

对于每个蒙特卡罗路径，j=1。。Mo计算{Xnk}k=1:n-1，un-1，σn-根据（2.3），（2.6）。0.50.60.70.80.91.11.21.31.41.50 0.2 0.4 0.6 0.8 1路径的第一部分（Euler方案的第n-1项）最后一项图1。V ibrato分解的模拟路径方案。o计算V（un-1） o计算un-1.θi和σn-1.θiby（2.11）、（2.13）和（2.12）o复制MZtimes最后一个时间步，即对于mZ∈ （1，…，MZ）–计算V（un-1+σn-1.√hZ（mZ））和V（un-1.- σn-1.√赫兹（mZ））3。在（2.18）中，通过对所有MZresults进行平均，然后乘以uθi和σθi然后在M条路径上求平均值。备注1。对于简单的案例，如欧洲期权的敏感性，asmall MZsu ffices；这是因为t是外环中关于M的另一个平均值。备注2。对于欧式期权，也可以使用Black-Scholes公式计算（2.15）2.4中的预期值。第二种竞争。假设X、b和σ依赖于两个参数（θ、θ）∈ Θ。计算二阶导数有两种方法。要么在使用引理2.1的同时区分振子（2.15），要么在二阶导数上应用振子。2.4.1。振动微分的二阶导数。让我们对第二个参数θj进行微分（2.15）：θiθjE[V（XT）]=E√Huθiθj·EhV（u+σ√hZ）σ-TZi公司+uθi·θjEhV（u+σ√hZ）σ-TZi公司u=un-1（θ）σ=σn-1（θ）+2h∑θiθj：超高压（u+σ√hZ）σ-T（ZZT- 一） σ-1i+∑θi：θjEhV（u+σ√hZ）σ-T（ZZT- 一） σ-1iu=un-1（θ）σ=σn-1（θ）（2.19）衍生品可以进一步扩张；例如，在一维情况下，在一个繁琐的代数之后，我们得到了：定理2.6。（二阶振动微分）θE[V（XT）]=E“uθEV（u+σ√hZ）Zσ√H+uθEV（u+σ√hZ）Z- 1σh+σθEV（u+σ√hZ）Z- 5Z+2σh+σθEV（u+σ√hZ）Z- 1σ√H+ 2.uθσθEV（u+σ√hZ）Z- 3Zσh（2.20）2.4.2。二阶振子的二阶导数。

同样的振动策略也可以直接应用于二阶导数。与之前一样，衍生品被转换为XT的PDF p：θiθjE[V（XT）]=ZRdV（x）p（x）Pθiθjp（x）dx=ZRdV（x）[ln p公司θiθj+ ln p公司θi ln p公司θj]p（x）dx=EV（x）ln p公司θiθj+ ln p公司θi ln p公司θj（2.21）然后θθE[V（(R)XnT（θ，θ））]=^1θθ（u，σ）=uθuθ^1u（u，σ）+σθσθ^1σ（u，σ）+uθθ^1u（u，σ）+σθθ^1σ（u，σ）+uθσθ+σθuθ^1uσ（u，σ）。我们需要计算两个新术语θθun-1（θ，θ）和θθσn-1（θ，θ）。它需要计算切向过程Yt的θiof的一阶导数，我们称之为Y（2）t（θ，θ）。然后（2.13）是不同的，一个基本但繁琐的计算涉及以下命题：命题2.7。引理2.11中定义的θi相切过程Y（i）tde具有一个由dy（ij）t=hb′θiθj（θ，θ，Xt）+b′θi，x（θ，θ，Xt）Y（j）t+b′θj，x（θ，θ，Xt）Y（i）t+b′x（θ，θ，Xt）Y（i）tY（j）t+b′x（θY（ij）tidt+hσ′θiθj（θ，θ，Xt）+σ′θi，x（θ，θ，Xt）Y（j）t+σ′θj，x（θ，θ，Xt）Y（i）t+σ′′x（θ，θ，Xt）Y（i）tY（j）t+σ′x（θ，θ，Xt）Y（ij）tidWt。最后，在单一情况下θ=θ=θ，这给出了建议2.8。（二阶颤音）θE[V（XT）]=E“uθEV（u+σ√hZ）Zσ√H+uθEV（u+σ√hZ）Z- 1σh+σθ超高压（u+σ√hZ）Z- 5Z+2σh+σθEV（u+σ√hZ）Z- 1σ√H+ 2.uθσθEV（u+σ√hZ）Z- 3Zσh（2.22）备注3。它相当于命题2.6，因此相当于颤音的直接差别。2.5。高阶颤音。借助Fa\'a di Bruno公式，可以将振子AD方法推广到振子相对于参数θ的更高阶微分，并将其推广到具有向量变元的复合函数，如Mishkov[32]所示。2.6。对偶变换、正则性和方差。

在本节中，为了简单起见，我们假设d=q=1。从m振动开始φ（u，σ）=E[f（u+σ√假设f Lipschitz与Lipschitz系数[f]Lip连续，我们得到^1u（u，σ）=Ef（u+σ√hZ）Zσ√H= Ef（u+σ√赫兹）- f（u- σZ√h）Z2σ√H. （2.23）因此，方差满足VaRf（u+σ√赫兹）- f（u- σ√赫兹）Z2σ√H≤ E“f（u+σ√赫兹）- f（u- σ√赫兹）Z2σ√H#≤ [f] LipE“（2σ√hZ）4σhZ#=[f]LipE[Z]=3[f]Lip。当E[Z]=0时，我们也有^1u（u，σ）=Ef（u+σ√赫兹）- f（u）Zσ√H. （2.24）然后，Varf（u+σ√赫兹）- f（u）Zσ√H≤ E“f（u+σ√赫兹）- f（u）Zσ√H#≤σh[f]LipEh（σ√hZ）Zi=[f]LipE[Z]=3[f]lip备注4。公式（2.23）和（2.24）的方差相等，但后者的计算成本较低。如果f是可微的且f′具有多项式增长，那么我们也有^1u（u，σ）=E[f′（u+σ√hZ）]。（2.25）因此，Varhf′（u+σ√hZ）i≤ Ef′（u+σ√赫兹）≤kf′k∞.备注5。设f]Lip表示f的Lipschitz常数。如果f′有界，我们有[f]Lip=kf′k∞那么（2.25）中的表达式的方差小于（2.23）和（2.24）。假设f′与Lipschitz系数[f′]Lip是Lipschitz连续的。我们可以提高（2.25）的效率，因为evarhf′（u+σ√hZ）i=Varhf′（u+σ√赫兹）- f′（u）i≤ Ef′（u+σ√赫兹）- f′（u）≤ [f′]LiphσE[Z]≤ [f′]Liphσ备注6。假设f（x）=1{x≤K} 显然，我们无法区分预期和之前所见的方差估计。2.6.1。指示器功能。假设f（x）=1{x≤K} 。

简化假设K≤ u，我们有f（u+σ√赫兹）- f（u- σ√赫兹）=新西兰≤K-uσ√ho公司- 1nZ≥u-Kσ√ho公司= 1nZ/∈香港-uσ√h、 u-Kσ√hio，因此f（u+σ√赫兹）- f（u- σ√赫兹）Zσ√H=σ√h | Z | 1nZ/∈香港-uσ√h、 u-Kσ√你好。对于方差，我们有VaRf（u+σ√赫兹）- f（u- σ√赫兹）Zσ√H≤ E“f（u+σ√赫兹）- f（u- σ√赫兹）Zσ√H#.作者：Cauchy Schwarz我们可以写“f（u+σ√赫兹）- f（u- σ√赫兹）Zσ√H#=2σhEZf（u+σ√赫兹）- f（u- σ√赫兹）=2σhEZnZ公司/∈香港-uσ√h、 u-Kσ√hio公司≤2σhE[Z]PZ/∈K- uσ√h、 u- Kσ√H≤√2σh2P级Z≥u- Kσ√H.然后√2σh2P级Z≥u- Kσ√H=√2σhZ+∞u-Kσ√他-udu√2π！。现在 a>0，P（Z≥ （a）≤E-aa公司√2π，所以当a→ +∞,风险值f（u+σ√赫兹）- f（u- σ√赫兹）Zσ√H≤σhre-（u-K）4σh（2π）qu-Kσ√H≤（2π）σhre-（u-K）4σh√u- K-→σ→0如果u6=K，则为0+∞ 否则事实上，这种估计可以用不可微f得到，这证明了可控震源技术的威力。振动加自动微分（VAD）的二阶导数。导致公式（2.22）的差异可通过AD自动得出；然后，一个人只需编写一个实现命题2.19公式的计算机程序，并对计算机程序应用自动微分。我们在此重新命名AD.3.1的基础。自动区分。考虑一个函数z=f（u）在C或C++中由double f（double u）{…}实现为了找到z′u的近似值，可以调用Cdouble dxdu=（f（u+du）-f（u））/dubecausez′u=f′（u）=f（u+du）- f（u）du+O（| du |）。一个好的精度应该通过选择小的来实现。然而，算术截断限制了精度，并表明适当选择du并不容易，因为在某个阈值上，有限差分公式的精度会因几乎为零与几乎为零的比率而退化。如Squire等人图2所述。精度（dzdu的对数-对数图-cos（1.）|使用正向有限差分公式计算，以评估u=1时的sin′（u）。图3：。

与图2相同，但具有使用复杂增量的细化效果；这两种测试都是用Maple-14进行的[38]，一种简单的补救方法是使用复虚增量，因为ref（u+idu）- f（u）idu=参考（u+idu）idu=f′（u）- Ref′′（u+iθdu）duleads to f′（u）=Re[f（u+idu）/（idu）]，其中分子不再是两个ms差值的结果。事实上，测试表明，当您→ 0（图3）。因此，可以选择du=10-8使用aO（10）渲染最后一个术语-16） 准确度，从而获得基本准确的结果。使用此公式的成本是f（）的两次评估，编程需要重新定义标准模板库inC++.3.2的所有双精度函数。AD处于直接模式。概念上更好的想法是基于这样一个事实，即计算机程序的每一行都是可微分的，除非在分支语句的切换点，如if和sqrt函数的零点等。用dx表示变量x的微分，A*b的微分是da*b+A*db，sin（x）的微分是cos（x）dx，等等。通过操作符重载，可以将此代数构建到一个C++类中，这里称为ddouble：类ddouble{public:double val[2]；ddouble（double a=0，double b=0）{val[1]=b；val[0]=a；}ddouble运算符=（const ddouble&a）{val[1]=a.val[1]；val[0]=a.val[0]；返回*this；}ddouble运算符-（const ddouble&a，const ddouble&b）{返回ddouble（a.val[0]-b.val[0]，a.val[1]-b.val[1]）；}ddouble运算符*（const ddouble&a，const ddouble&b）{返回ddouble（a.val[0]*b.val[0]，a.val[1]*b.val[0]+a.val[0]*b.val[1]）；}；因此，所有D双变量都有一个2数组数据：val[0]包含变量的值，val[1]包含其差值的值。

请注意，DdoubleassConstructor默认情况下会将零签名为val[1]。要了解它是如何工作的，请考虑图4中的C++示例，它调用函数f（u，ud）=（u-ud）对于u=2和ud=0.1。图5显示了相同的程序，其中double已更改为ddouble，而u的初始化意味着其差值等于1。打印语句现在显示fw的微分，如果所有参数的微分都被线性化为0，则fw也是其相对于u的导数，u的微分du=1。用alldouble f（double u，double u\\d）{double z=u-u\\d；return z*（u-u\\d）；}编写类doubleint main（）{double u=2，u\\u d=0.1；cout<<f（u，u\\u d）<<endl；return 0；}图4：。一个小型C++程序，用于计算（u- ud）在u=2时，ud=0.1。ddouble f（ddouble u，ddouble u\\u d）{ddouble z=u-u\\d；返回z*（u-u\\d）；}int main（）{ddouble u=ddouble（2，1.），u\\u d=0.1；cout<<f（u，u\\d）。val[1]<<结束；返回0；}图5：。同一程序现在计算DDU（u- ud）在u=2时，ud=0.1。函数和常见的算术运算符有点乏味，但并不困难。可以从www.ann下载一个示例。朱西厄。fr/Pironeau。该方法易于推广到高阶导数。例如，对于二阶导数，a.val[4]将存储a，其与第一和第二参数da、da和第二微分da的差值，其中两个参数可以是相同的。a*b的第二个区别是a*db+da*db+da*分贝+分贝*da等等。请注意，如果main（）中的firstline替换为ddouble u=2，则也可以使用相同的程序计算dfdudc。，ud=D双倍（0.1,1）；。然而，如果两个导数都需要DFDU、dfdudare，那么要么程序必须两次运行，要么必须修改ddouble类以处理偏导数。

在任何一种情况下，计算n个偏导数的成本大约是原始程序的n倍；反向模式没有这种数字复杂性，如果在directmode中使用带有特征的表达模板，则必须使用n>5的反向模式，并且n>5的其他w ise【35】。3.3。反向模式下的AD。考虑F′θ，其中（u，θ）→ F（u，θ）∈ 兰德u∈ Rd和θ∈ 注册护士。假设u是一个姿态良好的d线性系统的解Au=Bθ+c。应用于实现SF的c++程序的直接微分模式将求解线性系统n次，代价至少是dn操作。变分法的数学解从w′θdθ=(θF）dθ+(uF）du，Adu=Bdθ，然后引入s p∈ RDATP溶液=(uF）并置写入(uF）du=（ATp）Tdu=pTBdθ=> F′θdθ=(θF+pTB）dθ。p的线性系统只需求解一次，即至少执行O（d）运算。因此，由于线性系统通常是最昂贵的操作，当n较大时，第二种方法最为有利。仅由赋值组成的C程序可视为变量的三角形线性系统。循环可以展开并视为赋值和测试等。然后，通过上述方法，程序的ithline乘以piand，p从最后一行开始计算；但最大的困难是在计算p时，变量值的簿记。例如，对于f=u+ud的导数，对应于u g由{u=2*ud+4；u=3*u+ud；}驱动的ud，第二行中的u与第三行中的u不同，程序应重写为u1=2*ud+4；u=3*u1+ud；。那么p的系统是p2=1；p1=3*p2；导数是2*p1+p2+1=8。在这项研究中，我们使用了R.J.Hogan在inHogan【25】中描述的librar y adept 1.0。

该库的优点在于，反向模式的编程与上面介绍的直接模式非常相似；所有可区分的变量都被声明为D双变量，与之相关的变量在初始化时被指示，如上所述。3.4。不可微函数。在金融领域，无差别无处不在。例如，K中的二阶导数（x-K） +在x=Kas时不存在，但r的二阶导数∞f（x）（x）-K） +dx是f（K）。分布理论扩展了导数的概念：重边函数H（x）=1{x≥0}表示导数在z e roδ（x）处的狄拉克质量。通过使用δa（x）定义的函数δa（x）近似0处的Dir ac质量，可以将自动微分扩展到一定程度，以处理这种困难=√aπe-xa。现在，假设f在x=z处是不连续的，而在其他地方是不连续的；thenf（x）=f+（x）H（x- z） +f-（x） （1）- H（x- z） ）hencef′z（x）=（f+）′z（x）H（x- z） +（f-)′z（x）（1- H（x- z） （）- （f+（z）- F-（z） ）δ（x- z） 除非添加最后一项，否则二阶灵敏度的计算将不正确。如果在AD库中，斜坡函数x+定义为xH（x），其导数为H（x），如果H定义为其导数等于δaA，并且如果在计算金融资产的程序中，其写入的是t（x- K） +=r安培（x- K） ，则AD库计算的K中的二阶导数将是δa（x- K） 。更重要的是，它也会计算∞f（x）（x）- K） +dx≈NNXi=1f（ξi）δa（ξi- K） 式中，ξ是积分的N个求积点或程序员用来近似积分的蒙特卡罗点。然而，这种伎俩并不能解决所有问题，必须谨慎行事；例如，写下（x- K） +=（x- K） H（x- K） 不会产生正确的结果。此外，精度对a的值相当敏感。备注7。

请注意，有限差分（FD）并没有受到这个问题的困扰，这意味着具有复杂增量的FD对于一阶灵敏度来说是一种不错的方法。对于二阶灵敏度，“非常小超过非常小”的问题仍然存在。VAD和Black Scholes M模型。在本节中，我们实现并测试VAD，并给出一个描述该方法实现的概念性算法（自动完成）。我们关注的指标取决于SDE的解决方案，而不是SDE本身的解决方案。让我们以Black-Scholes模型中的标准欧式看涨期权为例。4.1。VAD的概念算法。1、基本资产及其相切过程Y的时间步长h=tn的一般模拟路径=十、θ相对于参数θ，k=0，N- 2：\'Xnk+1=\'Xnk+rh\'Xnk+\'Xnkσ√hZk+1，\'Xn=X，\'Yn=十、θ？Ynk+1=？Ynk+rh？Ynk+θ（rh）(R)Xnk+(R)Ynkσ√h类+θσ√H(R)XnkZk+1。（4.1）2。对于每个模拟路径（a），生成MZlast时间步长（\'XT=\'Xnn）\'Xnn=\'Xnn-1（1+相对湿度+σ√hZn）。（4.2）（b）通过可控震源，使用抗感应技术计算θ的一阶导数（公式（2.19），σ（Xt）等于Xtσ）VT公司θ=un-1.θ（VT+- VT公司-)Zn？Xnn-1σ√h类+σn-1.θ（VT+- 2VTo+VT-)锌- 1英寸Xnn-1σ√h、 （4.3）带VT±，o=（\'XT±，o- K） +，（\'XT±=\'Xnn-1+右侧Xnn-1±σ′Xnn-1.√hZn'XTo='Xnn-1+右侧Xnn-1.（4.4）和un-1.θ=(R)Ynn-1（1+右侧）+Xnn-1.θ（rh）σn-1.θ=(R)Ynn-1σ√h+Xnn-1.θ（σ√h） （4.5）如果θ=T或θ=r，我们必须加上θ（e-rT）VT至上述结果。（c） 在执行步骤4.3的计算机程序上应用自动微分方法，以计算在某个θ处与θ相关的二阶导数*.（d） 计算每条路径的平均值，即MZ上的平均值。3、计算结果向量的平均值（在M条模拟路径上）并对其进行贴现。4.2。希腊人。

Delta衡量的是溢价E[V（XT）]相对于现货价格X的变化率。Gamma衡量的是Delta相对于现货价格变化的变化率。Gamma对于投资组合的Delta对冲非常重要。Vanna是溢价相对于σ和X的二阶导数。TheVanna通过观察波动性的变化来衡量Delta的变化率。4.3。数字测试。为了生成Random编号，我们选择了C++STL版本11中提供的标准Mersenne捻线器生成器。我们取MZ=1，即我们只模拟每条路径的最后一个时间步；对于所有测试用例，除了欧罗巴-斯科尔斯模型中的调用契约。然而，对于Black-Scholes模型中的欧式调用，我们使用了带或不带Brownian桥的多时间步Witheueler方案。以下数值试验中考虑的参数为K=100，σ=20%，r=5%，T=1年。风险定价的初始价格从1到200不等。Monte C arlo参数设置为10万个模拟时间，2 5个时间步。4.3.1。初步数值试验。在这里，我们重点讨论了VAD对具有恒常性和漂移的标准欧洲看涨期权合约Gamma的数值精度，其中有一个解析Black-Scholes公式。因为颤音的颤音与颤音+AD（VAD）相似，所以比较两者是没有意义的。回想一下（命题2.6和2.8），它相当于将Vibra ato应用于可控震源，或将自动微分应用于可控震源。