重组二项式期望值的并行计算

对于一个进程，不需要将工作负载预先分配给各个进程；在tead中，每个进程都可以确定其要处理的2条路径的份额。这可以使用分配给每个流程的唯一排名。每一个进程都会计算出一个局部期望值，该值位于s充足空间的分区上，最终的期望值是通过对所有进程求和来计算的。这种求和是在消息传递接口（MPI）框架中通过reduce操作完成的，reduce操作以高效的方式协调进程之间的通信[7]。利用进程的唯一秩来并行化算法是并行计算中的一个常见主题。例如，Swarztrauber和Sweet[12]使用数据的二进制表示，将计算映射到泊松方程的并行直接解中的过程。这里，我们使用伯努利路径的前导位将路径映射到进程上。假设有M个并行进程的秩M=0，M-1.请注意，在MPI框架中，traditionallystart的排名为0。等级为m的过程称为“过程m”。我们假设M≤ N，M是2的幂。设r=log（M），以便进程的秩M可以用r位二进制表示M=zr写入-1r级-1+···+z+z，其中∈ B、 进程m被分配了所有路径x和pre fix（zr-1.z、 z）；这套2N-rpaths isdenotedBNm={x∈ BN:x=zr-1.xr公司-1=z，xr=z}。请注意，集合BNM构成BN的分区。图3.1显示了从M到BNm的映射图。对于每个伯努利路径x∈ BNm，进程m计算遍历路径p（x）和支付值VN（x）的概率。payoffs Vmon p进程m的本地期望值计算为vm=e-qTXx∈BNmp（x）VN（x）。

（3.1）最后，将局部期望值相加，得出最终结果tv（S）=M-1Xm=0Vm；（3.2）这是通过MPI red uce操作实现的，以获得进程0上的结果。（3.2）的计算需要2N-按照串行计算的要求，在2个并行进程上执行r步，而不是在单个进程上执行2步。如果我们愿意放弃完美的负载平衡，该方法可以扩展到M不完全是2的幂的情况。例如，我们可以考虑一个分区BN=BN∪ · · · ∪ bNK对于某些K>>M。P进程0可以处理bnmf，对于M=0，M，2M。，p进程1可以处理BNM，对于m=1、m+1、2M+1。，等等注意，路径划分的思想适用于每个节点上有两个以上分支的树。这类树的例子包括三项式树[13]。图3.1：进程Bernoulli路径映射4。Monte Carlo估计和方差缩减回想一下，集合BN中的Bernoulli路径数随N呈指数增长。当NBE变大时，即使有相当多的处理器，也不可能精确计算期望值（2.6）。蒙特卡罗（MC）估计提供了一种在不枚举整个样本空间的情况下近似复杂期望值的方法。在本节中，我们讨论一种MC方法，该方法使用第3节中的分区方案来近似产生的M个并行进程。MTH流程负责从BNm中绘制fr，形式=0，M-因此，我们有效地列举了每条路径的第一个r=logM步骤，并通过蒙特卡罗绘制其余步骤。与使用相同抽奖次数的基本MC估计值相比，这可以减少方差。定义θ=E[V（X）]=Xx∈BNV（x）p（x），其中V（x）中的suffix N被删除以便于标记。（2.6）中的选项值可以写为V=e-qTθ。

给定θ的估计值θ，则V的估计值为ΓV=e-qT^θ，其方差为Var（^V）=e-2qTVar（^θ），方差估计值isdVar（^V）=e-2qTdVar（^θ）。因此，我们将重点估计本节剩余部分的θ。让x，xRbe R独立且相同分布（i.i.d.）的伯努利路径采样自BN。然后θ的MC估计量由θ=RRXi=1V（xi）给出，其方差为var（θ）=RVar[V（X）]，（4.1），可通过DVAR（θ）=RRXi=1从MC图中估计五（xi）-^θ.在第3节中，我们将Bernoulli路的sp ace BN划分为BN，BNM。让Dm表示事件[X∈ BNm]发生概率为P（Dm），对于m=0，M- 1、Fu rthermore，考虑分区X=（Z，Y），其中Z∈ 品牌Y∈ BN公司-r、 我们现在可以把θ写成θ=M-1Xm=0E[V（X）| Dm]P（Dm），（4.2），其中P（Dm）=Xy∈BN公司-rP（Z=zm，Y=Y）=P（Z=zm），zm是m的二进制表示，对应于mthprocess的秩。设θ（m）=E[V（X）| Dm]，设X（m），x（m）Rmbe来自BNM上路径分布的i.i.d.样本，m=0，M- 1、假设PM-1m=0Rm=R，因此使用的样本量与基本MCestimator中的一样。估计量θ（m）s=RmRmXi=1V（x（m）i）是方差为mvar[V（x）| Dm]的θ（m）的无偏估计量。将（4.2）中的θ（m）替换为θ（m），得到分块MC估计量^θs=m-1Xm=0^θ（m）sP（Dm）。（4.3）继Rubinstein和Kroese【14】之后，我们选择与P（Dm）成比例的样本大小Rm，作为每m的Rm=R·P（Dm）。通过这种选择，并且忽略R·P（Dm）可能不是精确整数，分割MC估计量的方差可以写为var（θs）=m-1Xm=0Var（^θ（m））[P（Dm）]=RM-1Xm=0Var[V（X）| Dm]P（Dm）。

（4.4）相应的方差估计量isdVar（^θs）=RM-1Xm=0“RmRmXi=1V（x（m）i）-^θm#P（Dm）=RM-1Xm=0RmXi=1V（x（m）i）-^θm.为了验证^θ可以使^θ的方差减少，总方差定律给出了var[V（X）]=EDVar[V（X）| D]+VarDE[V（X）| D]（4.5）=M-1Xm=0Var[V（X）| Dm]P（Dm）+VarDE[V（X）| D]=R Var（^θs）+VarDE[V（X）| D]，（4.6），其中最后一个等式为f from（4.4）。用（4.1）中的Var（^θ）替换（4.6）中的左侧，并将两侧除以R，我们得到Var（^θ）=Var（^θs）+VarDE[V（X）| D]R.（4.7）注意，如果V（X）| D不依赖于伯努利路径的第一个R阶，则VarDE[V（X）| D]=0∈ {0，N}，即当进程数M∈ {1，2N}。当M=1时，分区MC方法与基本MC方法相同，当M=N时，它与（2.6）中的精确期望值相同。由于假设支付V（X）依赖于整个路径X，当0<r<N时，（4.7）右侧的第二项大于0，因此，（4.3）中的分割MC估计通常会严格减少方差。当V（X）在整个Dm中是不均匀的，并且在每个Dm中是均匀的时，还原将更加明显。备注1。分区MC方法的一个有趣变化是重用来自BN的相同ofR绘图样本-ron所有流程。再次考虑分区X=（Z，Y），并假设Zand Y是独立的。设x（m）i=（zm，yi），其中y，Y和z（m）=（z（m）r的分布中提取的Y.i.d-1.z（m），z（m））是m的二进制表示，对应于mthprocess的rankof。然后，θ（m）的估计值由θ（m）=RRXi=1V（x（m）i）（4.8）给出，θin（4.2）可以由θ=m无偏估计-1Xm=0￠θ（m）P（Dm）=RRXi=1M-1Xm=0V（x（m）i）P（Dm）=RRXi=1EX | Y[V（Z，yi）]。

（4.9）该估计值的方差为Var（|θ）=RVarYEX|Y[V（X）]，MC样本的方差估计值为Var（|θ）=RRXi=1EX | Y[V（Z，yi）]-θ（m）=RRXi=1M-1Xm=0V（x（m）i）P（Dm）-θ（m）！。我们将∧θ称为共享样本MC估计量。现在，V（X）的方差可以写为Var[V（X）]=VarYE[V（X）| Y]+EYVar[V（X）| Y]=R Var（|θ）+EYVar[V（X）| Y]。（4.10）用（4.1）中的Var（^θ）替换（4.10）中的左侧，并将两侧除以我们得到的Var（^θ）=Var（|θ）+EYVar[V（X）| Y]R.（4.11），当且仅当V（X）| Y不依赖于伯努利路径的第一步，因为我们假设V（X）依赖于整个路径X，Var（|θ）严格小于Var（^θ）。以下结果总结了^θ、^θs和^θ方差之间的关系。定理4.1。假设m=0，…，则Rm=R，M- 1，0<r<N，V（Zk，Y）和V（Zl，Y）与所有k，l呈正相关∈ {0，…，M- 1} 。然后Var（^θs）≤ Var（|θ）≤ Var（^θ）。证据我们已经证明了Var（|θ）≤ Var（^θ）。现在我们有Var（^θs）=RPM-1m=0Var[V（X）| Dm][P（Dm）]，从（4.9）开始，Var（|θ）=RVarY“M-1Xm=0V（Zm，Y）P（Dm）#=RM-1Xm=0变化[V（Zm，Y）][P（Dm）]+RXXk6=lCovV（Zk，Y），V（Zl，Y）P（Dk）P（Dl）=Var（^θs）+RXXk6=lCovV（Zk，Y），V（Zl，Y）P（Dk）P（Dl）≥ Var（^θs）。期权定价的应用我们采用第3节所述的方法，使用R 3.2.2和Julia0.4.6编程环境对看跌期权进行估值。计算是在带有computenodes的分布式集群上运行的，每个集群都有两个Intel E5-2650v2 Ivy Bridge（2.6 GHz，20 MB缓存）处理器和8个coresper节点，每个节点总共有16个核。所有节点都有64 GB的主内存，并通过四数据速率的带内互连进行连接。

Open MPI 1.8.5[15]被用作MPI框架的底层实现。R是一种统计计算环境，有助于高级数据分析，提供图形功能和解释性高级编程语言[16]。除了核心R环境中可用的统计、计算和p编程特性之外，用户社区提供的许多包还提供了其他功能。Rmpi【17】和pbdMPI【18】包可用于从R编写MPI程序。本节中显示的结果基于Rmpi，但pbdMPI在我们的经验中表现类似。packageRcpp[19]有助于将C++代码集成到R程序中，从而以增加编程负担为代价大幅提高性能。我们还没有在我们的实现中探索Rcpp，但注意到它的潜在用途。Julia是最近开发的一种编程语言，在科学计算、数据分析和高性能计算领域越来越受欢迎。它是一种编译语言，使用低级虚拟机实时技术[21]生成编译到机器级的源代码的优化版本。Julia pr在核心环境中以及通过用户社区提供的包提供了许多计算和统计功能。我们在Julia中使用了MPI包来运行MPI程序。在Julia中，也可以通过CxxWrap和Cpp等包与C++集成，但我们尚未对其使用进行探讨。我们的实现在MPI包中使用了动态Julia代码。因为ejulia被编译成机器级代码，所以用Julia编写的程序应该比用R编写的等效程序性能更好。

本节后面的绩效结果证实了我们的直觉。清单1显示了并行化方法的Julia实现的一个片段。由于我们的R和Julia实现的结构相似，所以我们没有显示类似的R代码列表。在第1行中，我们加载MPI包。由于我们的实现遵循SPMD范式，因此在所有进程上运行相同的代码。第5行请求运行代码的进程的级别，第6行请求MPI通信中的进程总数。正如第11行的while循环所示，每个进程都在2N上工作-totalNBernoulli路径的路由。请注意，第13行中的完整伯努利路径的构造是通过在代码运行到当前（N- r） -维伯努利路径。第14行对calcpath prob的函数调用计算使用第13行构造的伯努利路径的概率。第15行的callto calc payoff函数计算期权支付；它们的函数定义没有显示出来，因为它们独立于并行化方法。最后，在第21行，将各个p进程的期望值相加，以获得进程0.1导入MPI2的最终答案。3 MPI。在i t（）中，4 comm=MPI。CO MM\\u W ORLD5 id=MPI。Comm\\u ra nk（Comm）6 M=MPI。通信大小（co m m）7。8 r=log2（M）9 l\\u n=转换（Int64，2^（n-r））10。11 wh ile i<l\\u n12 node=i13 path=cat（2，in teg er\\ba se\\b（id，2，r），i nte ge r\\u b base \\u b（node，2，N-r））14 p\\u vt=ca lc\\u pat h\\u pro b（path，prob）15 vt=calc\\u payof f（s，K，u，d，opt\\u type，pat h）16 v+=p\\u v t*vt17 i+=118 end19。20 v=exp（-q*T）*v21减少d\\U v=MPI。还原e（v，MPI。

清单1：parall el-Bernoulli路径算法的Julia实现。我们以看跌期权为例来说明OUR方法。我们设定了K=10的执行价。资产的当前价格和波动率分别为S=5和σ=0.30。无风险利率为q=6%，到期时间为一年。表5.1（a）和表5.2（a）分别显示了问题大小N=16、20、24、28和32时，R和Julia实现的墙时钟运行时。虽然这两种实现都能很好地随进程数M扩展，但JuliaiImplementation的速度大约是R的10倍。我们在单个进程（M=1）上的R程序f或N=32导致了循环中的过流，该循环计算了自2- 1是可存储在R中的最大整数值。因此，表5.1（a）中记录了该特定字母的运行时间为不适用。如果Tm是M个进程的运行时间，则M的加速比SM和效率Em分别定义为T/Tm和SM/M。如果该程序完全扩展到M个进程，则可获得理想值SM=M和EM=1。这些数字表明了该计划的可扩展性。由于我们的R程序在n=32的情况下没有在单个进程上运行，因此我们认为这种情况下的加速比为2·T/TM，M=2，64和形式=1和M=2时，加速比分别取1和2。表5.1（b）和5.2（b）显示了加速，表5.1（c）和5.2（c）分别显示了我们的R和JuliaiImplements的效率数字。图5.1中的曲线图分别显示了表5.1（b）和表5.1（c）中的加速比和效率数字，图5.2显示了表5.2中的相应曲线图。这些图从视觉上证实了我们的猜测，即对于我们的问题，Julia比R更有效。请注意，对于固定的问题规模，在一定数量的任务之外，加速的优势会降低。

这是因为协调任务的开销开始支配进行有用计算所花费的时间；参见Pacheco【7】f或更多详细信息。这可以从表5.1（a）中的N=16和表5.2（a）中的N=16和N=20中看出。在这两种情况下，由于verysmall的总运行时间主要由MPI reduceoperation消耗的近乎恒定的时间决定，因此s peedup和效率数字明显低于大型问题。表5.1（b）和表5.1（c）中N=32的相对较高的加速比和效率数是将M=2的加速比设置为2，并计算相对于M=2的其余加速比和效率数。由于Julia的运行时间大约比R的运行时间快10倍，如果我们估计R的运行时间为N=32，M=1，并相应地计算加速比和效率数字，我们会注意到数字会下降，并且与其他案例相当。我们实现了第4节中描述的亚式期权和回溯期权的蒙特卡罗估计方法，这两种期权都是路径相关的[2]。在亚洲期权中，期权支付函数中的资产价格统计（到期时间）替换为算术平均值{St:t=1，…，N}。因此，在二叉树模型中，亚洲看跌期权的收益由v（x）=max{K给出- S*, 0}，（5.1），其中S*=NPNt=1St（x），St（x）是贝努利路径x上紧随时间t的资产价值。在回溯操作中，执行价格K或资产价格统计到期时间在支付函数中分别替换为{St}的最大值或最小值。这里我们考虑一个固定的L ookback看跌期权，其收益由V（x）=max{K给出- S*, 0}，其中S*= 最小值{St（x）：t=1，…，N}。

我们使用大小为N的二叉树模型，对第4节中给出的亚洲和固定回望看跌期权的基本MC估计进行了实现，以研究估计是否收敛到准确的预期值（2.6）。我们进一步实现了分区和共享第4节中的样本MC，以研究方差缩减特性。表5.3显示了亚洲看跌期权和da固定回望看跌期权的基本MC估计和相应方差估计，参数K=100、S=20、q=6%、σ=3.0和T=1，使用树大小N=32的二项模型。通过精确计数计算得出的亚洲看跌期权的期权值为82.115，固定回望看跌期权的期权值为93.196。用于MC估计的样本大小从2增加到2，小于路径总数的0.01%。从表5.3可以看出，两个选项的MC估计值收敛到各自的精确值。此外，正如预期的那样，方差估计值随着样本量R的增加而减少。表5.4显示了亚洲和固定回溯看跌期权的分割MC估计值^VS以及相应的方差估计值，使用总样本量R=1024，并在1到64之间改变过程数。请注意，随着过程数量的增加，每个过程的样本量rm减少。表5.4中所示的估计值是1000次重复的平均值。正如预期的那样，表5.3显示，当R=1024时，分区MC估计的方差估计大多小于相应的基本MC估计。表le5.5显示了N=32，R=Rm=1024，m=0，…，亚洲看跌期权的分割样本和共享样本MC估计之间的方差估计比较，M- 1、共享样本表5.1：R实现的不同时间步数的运行时。