在第6.4节中，我们提供了数字证据，表明协方差核的选择通常对产生的模拟器^fk只有很小的影响。推断超参数θ和过程方差σ的经典方法是使用基于第3.1节所述分布的似然函数，通过极大似然估计或惩罚极大似然估计优化边际似然。任何一种情况都会导致非线性优化问题。在连续设置中，超参数在每个阶段都得到理想的拟合，以便随着更多数据的到来，改善似然函数。然而，这在计算上是不切实际的，因为评估可能性需要O（| Dk |）。替代方法是在算法中的某些点上重新设置超参数，例如在10%、20%、，90%（或采用非线性计划，如2%、4%、8%，…）预算的一半已经用完了。后一种方案的优点是，当超参数不太确定时，以及当参数更便宜时，可以更早地进行拟合。上述逻辑还暗示了对保持| Dk |小的算法的偏好。从这个意义上说，像ST-GP这样的更有针对性的算法比SV-GP更可取，SV-GP本质上是在许多场景中传播内部模拟。相关开销涉及预测^fkon Z，即评估mk（z1:N）和sk（z1:N），这是所有采集功能所需要的。这带来了O（| Dk |·| Z |）效应，再次显示了大型设计Dk的成本。它还带来了一个技巧，即通过远离尾部的筛选方案，将Z减少到一个andidate集Zcandkby。筛选的一种方法是重新使用等式（19）和（23）中的权重Wk（zn）来评估情景的相关性。

具体而言，在下面的案例研究中，我们确定候选集Zcandkat stage k viaZcandk={zm∈ Z:WVaRk（zm）/（NXn=1WVaRk（zn））>10-3} ，对于VaR和Zcandk={zn∈ Z:WTVaRk（zn）>10-3} 用于TVaR。因此，在评估采集功能时，仅考虑外部场景（方法ST-GP、SE-GP、SV-GP）。如图所示，在第一轮案例研究的典型运行中，该筛选将预测集从n=10000个场景缩小到第一轮后的| Zcand |=528（657用于TVaR），从而将在Zcand上预测GP代理的开销减少了一倍≈ 此外，候选集| Zcandk |随着k的增长而收缩，因为GP更好地了解尾部场景。在下面的示例中，我们使用了一个固定的阶段预算rk公司≡ r修饰版本可以很容易地包括在各轮中更改批次大小，例如，在GP具有高度全局不确定性的前一轮中的较小批次。选择批量大小还决定了总轮数K，影响GP回归和预测开销，以及UR采集函数的计算成本。一般来说，我们发现（参见第6.1节），算法只关注少量场景，这意味着我们的顺序标准会重复选择相同的zn。因此，大批量生产r=0.1N，与较小的r=0.01N，减少了选择zk+1的开销。回到GP模拟器本身，我们发现（12）中的一个简单趋势函数u（·）对于去趋势是有效的。虽然u（·）的选择对预测^fk几乎没有影响，但它有助于非线性优化程序设置超参数，并且在顺序阶段的每一个开始都是有益的。注意，趋势在一定程度上改变了响应的空间依赖性，从而改变了θjhyperparameters。

参数平均函数也可以通过具有各自系数的基函数指定，因此u（z）=β+Pjβjhj（z）。在这种情况下，β与其他超参数同时估计，这种方法称为通用克里格法（Universal Kriging）[33]。根据我们的经验，在大多数金融环境中，u（·）可以很容易地被确定为投资组合的内在价值（如我们的第一个案例研究），也可以被视为常数β（作为GP的一部分，见第二个案例研究）。5.4与其他方法的比较为了对第4节中定义的拟议算法进行基准测试，我们比较了几种备选方案。我们主要关注其他基于回归的方法，并专注于量化（i）自适应预算分配的作用；（ii）与更简单的1、2或3阶段方法相比的顺序方法。表1.1总结了这些基准以及本节前面讨论的程序。LB：用作下限的完美信息方法。我们假设尾部/分位数是先验已知的，并希望最小化^R的均方误差。对于VaR，这对应于最小化qα处的后验方差，这是通过将整个分位数N分配到真正的分位数场景来实现的。在这种情况下，没有GP替代项，估计方差是蒙特卡罗平均误差，即τ（qα）/N。对于TVaRα，我们将预算N均匀地分配到真尾{zn:f（zn）≤ f（αN）}，并对结果进行GP。2、A3-GP：Liu和Staum的自适应三阶段算法【27】。第一阶段模拟类似于我们方法的太空填充试点场景。第二阶段通过筛选的候选集Zcand来分配uniformlyacross。然后，第三阶段求解R0to最小化方差，如SV-GP和引理1中的^Rlike。我们遵循他们的建议，第2阶段的设计位置为0.02N，第3阶段的预算为N=0.7N。3.

U2-GP：两阶段进近。在第5.2节所述的空间填充第一阶段之后，剩余预算将统一分配到最低的2αN方案中。这是rn=N/（2αN）{m（zn）的自适应分配的最简单版本≤m（2αN）}。注意，U2-GP可以看作是SR-GP的一个版本，K=2轮，L=1，U=2αN。与此方法相比，多阶段程序的增益是量化的。4、U1-GP：一种将N/| Z |均匀分配给每个Z的一级方法∈ Z、 将GP替代项^f设置为最终设计，并使用m1:n估计^R。这是最粗糙的比较器，它没有自适应分配，但仍采用空间平滑。我们发现，对于大型场景集，将aGP配置到整个输出集合在计算上是不可行的，N 1000，所以我们改为{zn:yn≤ y（2000）}，这产生了与T GP相似的开销时间，并且与无GP平滑相比仍有显著改善。我们还考虑了两种非GP方法来研究空间平滑的重要性。第一种被称为U1-SA，与U1-GP相同，但使用样本均值“ynand empiricalvariations^τ（zn）”代替基于GP的后验均值。这是“普通的”嵌套模拟方法。第二种是排序和选择（R&S）方法，Broadie等人[8]介绍了该方法，目的是将f1:nw与给定的L级进行比较。根据我们的符号，在k阶段，该BR SAalgorithm分配Rk对方案zk+1的模拟，该方案使HBRk（zn）最小化。=rnk·|(R)ynk-L |/τ（zn）。上述采集函数将场景样本平均值与L进行比较，通过各自样本不确定性τ（zn）/Rnkw进行归一化，其类似于HECIin（25），但没有模拟器。出于我们的目的，我们将L替换为^RHDk。

此外，我们按照[8]的建议，根据收集到的总方差，通过加权|τ（zn）估计局部噪声方差，名称描述参数完全信息下的SLB蒙特卡罗TVaR均匀分配到真实尾部场景SST GP目标MSE准则（22）εk=s（^RHDk）SE-GP预期水平改进准则（25）SV-GP批次方差最小化（第4.3节）SR-GP-1{zn:m（L）k上的均匀≤ m（zn）≤ m（U）k}L=50，U=50（VaR），L=1，U=50（TVaR）SR-GP-2 L=26，U=75（VaR），L=1，U=75（TVaR）A3-GP 3阶段，从[27]阶段3预算的70%；注册护士≡ 10U2-GP两级SR-GP，L=1、U=100和K=2U1-GP均匀分配GP根据“y1:NU1-SA均匀分配使用样本平均数”y1:Nas输出br SA R＆s算法，来自【8】表1：第5.4节所述的比较方法。所有方法（U1-GP除外）使用算法1中描述的相同初始化过程，Ninit=0.01N第一阶段场景。选择顺序方法的预算参数得出K=100。对应于一个非常简单的同调平滑器：？τk（zn）=rnkrnk+？r^τk（zn）+？rrnk+？r？τk，其中？τk=NPNn=1^τk（zn）是所有经验方差的平均值。这里，r是smoothingparameter，我们取下面的[8]中的r=5。第6.3节讨论了与这些非GP基准的比较。评论请注意，按照上述“配方”，可以将任何引入的采集函数与基于样本平均数的方法相结合，产生ST-SA或SV-SA方法。如上所述，没有模拟器的一个限制是无法在没有足够内部模拟的场景下正确预测τ（z）。从另一个方向来看，我们回顾了最小二乘蒙特卡罗方法（如Bauer等人[3]），该方法使用空间平滑，但使用非顺序分配。

我们不与后者进行比较，因为我们的重点是实验设计，而不是纯粹的回归/平滑步骤。6案例研究：Black-Scholes期权组合从一个二维示例开始，其中f（z）可以精确计算。使用已知f在2-dw中工作可以方便地可视化算法，并提供精确的误差计算。考虑一个投资组合，其价值由两项风险资产Zt驱动≡ （St，St）具有几何布朗运动动力学：dSt=βStdt+σStdW（1）t，dSt=βStdt+σStdW（2）t。在风险中性度量Q下，Wiare相关布朗运动与dhW（1），W（2）it=ρdt，β是恒定利率。表2总结了我们的模型参数。投资组合包括看涨期权：长100 K=40-沙地上的罢工要求50 K=85-删除S。我们的风险期限为T=1年，目标VaRα和TVaRα风险度量为α=0.005。通过风险中性定价，投资组合的价值为T，从z开始≡ （z，z）∈ R+isf（z）。=EQh100e-β（T-T）装货单- 40+- 50e-β（T-T）装货单- 85+（ST，ST）=zi。（33）资产头寸初始价格SiStrike KiMaturity TiVolatilityσ为100 50 40 2 25%S-50 80 85 3 35%相关性ρ=0.3利率β=0.04表2：股票和S上Black Scholes投资组合的二维案例研究参数。从Black Scholes公式获得的f等高线图如图3所示。观察到f（z）在左上角最为负，而τ（z）在z和z两个方向都在增加，因此，在不确定性较大且场景更稀少的右上角可能会花费更多的注意力。图中还显示了场景集Z，我们从（ST，ST）的双变量对数正态分布中生成了一个固定样本。注意，这里Z是使用Q分布生成的，因此[0，T]和[T，Ti]上因子的动力学是相同的，即。

物理和风险中性措施是一致的。这只是为了更简单地介绍案例研究；根据我们的经验，Z的作用仅次于其他考虑因素。图3中的红线显示了真实（相对于所示Z）分位数损失f（αN）=f（qα）=-4052.02，指示方法应针对的区域。图3：二维Black-Scholes期权案例研究的真实投资组合值f（z），以及相应的模拟标准偏差τ（z）。红色表示真实分位数损失f（50）=-4052.02.点云表示场景集Z={zn，n=1，…，10000}。在本案例研究中，通过模拟S的对数正态值，S条件为（S，S）=（z，z），并将其插入Payoffy=100e得到输出yn，i-β（T-T）装货单- 40+- 50e-β（T-T）装货单- 85+, （34）式中，T=1、T=2和T=3。我们继续使用全局实验参数N=10，N=10，比较表1中列出的所有方法，从而得出αN=50。注意，sinceN=N，U1-GP为每个场景分配一个内部模拟。在初始化过程中，使用的完全顺序方案为N=0.01N=100，rn=10（即N=N- rNinit=9000），K=100级，因此在所有其他轮次中，rk=N/K=90。GP超参数在阶段k=10、20、…、，对于每个顺序方法为100，对于A3 GP和U2-GP，每个阶段后为100。因为我们有风险RVaR的确切值f（50）和尾部VaR RTVaR=Pn=1f（n），VaR或TVaR的给定最终估计值的偏差和平方误差（SE）都是简单的=^RK- R、 和SE（^RK）=^RK- R.对于（12）中GP的平均函数u（·），我们取T处投资组合的内在价值，即u（z，z）=100e-0.04（z- 40)+- 50e-2·0.04（z- 85)+.

（35）将上述选择与更简单的常数平均函数u（z）=β和其他GP超参数的βfittedalong进行比较；尽管后者确实增加了GP的不确定性sk（zn），但在方案的最终性能方面没有产生统计上的显著差异。6.1比较算法通过执行100次宏复制m=1，…，对算法进行评估，100（即使用一组新的采样yn，i’s重复整个过程），对于固定集合Z的VaR和TVaR，获得估计值^R[1]K，^R[100]K.这会产生各种算法估计量的真实抽样分布，控制内部模拟的内在可变性。这些输出用于计算偏差、方差和SE，结果通过几个表格和图表进行说明。第6.2节提供了顺序方法的详细可视化。表3报告了（i）SD（^R[1:100]K），这是100次宏观复制的经验标准偏差，（ii）与^RK相关的平均估计GP后验不确定度：s.=Pm=1s（^R[m]K），（iii）相对于地面真值的均方根误差，RM SE=qPm=1（^R[m]K- R） 与（iv）average | DK |一起，是算法选择的不同外部场景的数量。注意，对于LB和U1-GP以外的方法，| DK|-100（其中100=Ninitis试验阶段0方案的数量）是初始化后选择的位置数量。图4和图5显示了^R[m]Kand s（^R[m]K）的结果分布箱线图，其中水平线是从Black-Scholes分析计算中获得的真实风险度量R。对于VaR，我们重申，我们的估计值预计会收敛到Harrell-Davis估计值RVaR，HDR在下面的讨论中被用作基本事实。

在本示例中，相对于真实分位数的差值约为20，f（αN）=-4052.02，右驾=-4032.21.我们从风险价值开始讨论。统一分配（U1-GP）结果未绘制，因为它们远远超出轴限制；表格显示其RMSE比LB大60倍以上。所有其他方法都比传统的嵌套蒙特卡罗方法有明显的改进。按照U1-GP的顺序进行比较，U2-GP和A3-GP大致显示出改善，这要归功于Var0.005TVaR0.005SD（^RHDK）s RMSE | DK | SD（^RK）s RMSE | DK | LB 44.35 40.42 46.77 1 47.29 53.48 47.42 1ST-GP 50.57 48.55 50.59 121.52 59.12 55.17 61.46 118.27SE-GP 50.48 50.71 74.03 116.12 93.36 87.83 95.70 111.79SV-GP 56.50 48 48.48 28 60.53 305.03 55.78 54.76 56.65 163.08SR-GP-1 63.27 61.90 69.74 112.43 61.48 55.34 61.86 165.27SR-GP-2 50.45 49.8250.52 180.97 61.66 61.56 62.13 193.46A3-GP 61.07 54.18 60.83 292.83 63.57 59.92 63.18 297.44U2-GP 68.76 55 55.91 68.47 194.55 64.92 67.77 64.87 194.64U1-GP 695.33 560.52 2965.05 10909.17 700.43 3003.07 10表3：对于二维Black-Scholes案例研究，样本标准差（SD）超过100个宏重复，平均GP后验标准差s，以及每种方法的RMSE RK，以及平均最终设计尺寸。方法说明见表1。图4：二维Black-Scholes案例研究的VaR估计。左箱线图显示了最终估值的分布情况；右侧是相应GP标准偏差s（^RVaRK）的分布。结果基于每种方法的100次宏复制。增加级数，与U2-GP相比，后者的RMSE减少约20%。按照100个阶段的顺序方案，我们又获得了15-20%的改进–将ST-GP与A3-GP进行比较。可以进行其他一些比较。

在两种SUR类型的策略中，SE-GP相对于ST-GP表现不佳，这表明后一种获取函数更适合分位数学习，可能是因为它能够解释估计的^RHDK中的不确定性。对于基于秩的方法，SR-GP-2表现良好（与ST-GP竞争）表明，在本案例研究中，L=26，U=75的选择对VaR0.005效果良好，而SR-GP-1表现不佳表明其对经验分位数的关注范围狭窄。由于过于激进，SR-GP-1可能会完全错过一些构成RHD的场景，从而导致更多的偏差和分散交叉运行。总体而言，就RMSE而言，表现最好的是ST-GP、SV-GP和SR-GP-2。图5：二维Black-Scholes案例研究的TVaR估计。左箱线图显示了最终估值的分布；右侧是相应的GP标准偏差s（^RTVaRK）。结果基于每种方法的100次宏复制。箱线图我们看到，所有方法（除了预期无偏的LB）都有偏高，即低于报告资本要求。值得注意的是，SE-GP和SR-GP-1的偏差相对较高，SD（^RK）也较大。请注意，如果没有空间平滑，分位数估计值会偏低，因为经验分位数相对于地面真值更极端。通过空间平均，模拟器以另一种方式拉动偏差。表3还比较了s（^R[m]K）和SD（^RK）。前者是^R标准误差的内部仿真器BaseDestinate，而后者是在宏复制中观察到的采样标准偏差。没有办法先验地计算SD（^RK），我们希望S（^RK）可以作为代理。这种准确报告标准错误的能力是仿真的一个重要部分。

当然，s（^R[m]K）本身是随机的，因此在表中我们报告了其平均值，图4显示了其自身的抽样分布。因此，目标是使s（^R[m]K）稳定的Acrosruns，平均值接近SD（^RK）。我们确实观察到，所有方法都相当接近SD（^R[1:100]K），最大的差异出现在SV-GP和A3-GP中。这两种方法都是为了最小化s（^Rk），这可能会导致这些值偏低。从箱线图中，我们注意到SV-GP在整个运行中具有最稳定的s（^RK），而SR-GP-1、A3-GP和U2 GP的标准误差估计值最不可靠。对于有风险的尾部价值，结果几乎相同。我们再次看到SE GP的表现不佳（RMSE比LB大2.01倍），这证实了各自的收购功能不足。还请注意，在TVaR环境中，更具攻击性的SR-GP-1优于SSR-GP-2。这说明SR-GP需要仔细调整L和U值。通过图6中的扇形图提供了顺序方法的最终视觉比较，图6说明了随着k轮进度，估计风险度量的演变。这是允许“在线”使用算法的顺序程序的主要优点之一。因此，用户可以监视^rk，例如判断收敛性、自适应停止模拟或报告临时估计。扇形图显示r[m]kas的宏观复制上的分位数（以m为单位）是k的函数。大多数方法最初偏向高，随着k的增加和采样方差的减小，偏向消失。SV-GP以最快的速度减少偏差。值得注意的是，SR-GP-1并不收敛于GP SE-GP SV-GP SR-GP-1 SR-GP-2VaRαTVaRα图6：对于二维Black-Scholes案例研究，扇形图描述了随着预算的消耗，^RHD、VaRk（顶部）和^RTVaRK（底部）的演变。

在100个算法运行时，红线对应于^R[m]kover的0.1、0.25、0.5、0.75和0.9分量分位数。在一些跑步中，强调了探索的重要性。对于TVaR，我们观察到，经过几轮之后，偏差（^Rk）通常已经很低，这表明估计水平集比估计分位数更容易。我们还看到，ST GPAR最初很快就成功地学习了TVaR，但随后其性能停滞不前；k=100SV时，GP可获得更好的偏差和更低的SD（Rk）。SE-GP在TVaR方面的表现较差，特别是在SD（RK）非常高的方面，证实了前面的讨论。检查批量大小的影响r、 我们重新运行了ST-GPr=0.001N（即每轮9个新的内部模拟），观察到的s和RMSE分别为48.37和59.22。这些值在统计上与相同的ST-GP方案无法区分r=0.01N。换句话说，批量大小为r=0.001N（即K=1000发）与r=0.01Ndoes在本实验中没有提供显著的改善。我们发现，较小的批次最终会使用一个过程，该过程连续多次选择相同的场景，因此它的行为与较大的批次相似。6.2比较顺序方案为了补充图1，图1给出了一个fixeddataset中不同方案针对的方案的快照，图7显示了代表运行结束时ST-GP、SE-GP和SV-GP产生的总体预算分配r1:NK。对于VaR0.005，ST-GP和SE-GP在50级左右的效应浓度=αN时表现相似。作为参考，初始化后的Nbudgetafter的93%、73%和99.1%分配给了{z:f（25）≤ f（锌）≤ f（75）}分别用于ST-GP、SE-GP和SV-GP，最大复制量maxNRNK为3790、1260和990。

这突出了高适应性，在单个场景上花费了高达30%的预算（与在分位数场景上花费100%的LB基准相比）。对于TVaR，95%、100%和98.1%的Nbudget用于{z:f（1）≤ f（锌）≤ f（50）}分别用于ST-GP、SE-GP和SV-GP，最大分配为1890、4600和531。特别是，ST-GP和SV-GP成功地从左尾的几乎所有场景中识别和采样（或多或少具有可比性RNK），而SE-GP导致了命中和未命中设计，因为f（1:50）的真实位置是相关的。尽管所有SE-GP样本都处于真正的尾部，但其设计往往非常集中（只添加了大约十几个非试点场景），创建了小的空间集群，每个集群中探索一个场景。VaR Tvar图7：二维Black-Scholes案例研究的序贯方法的最后阶段后，复制计数f1:n的rnKversus真实等级。y-轴位于对数刻度上。所有三个方案共享完全相同的初始化阶段，可以通过场景（尤其是在每个图的右侧）看到r=10。否则，对于ST-GP和SE-GP，所有rnK都是r=90。对于SV gp，RNK上没有可低至1的约束。从不同的角度来看，表3列出了平均总设计尺寸| DK |。我们观察到，ST-GP和SE-GP只使用了大约120个场景（回想一下，初始化期间已经选择了Ninit=100个场景），而SV-GP使用了大约300个场景。这是因为ConstructionalSV-GP跨多个场景定位RNKA，这也可以在图7中观察到。特别是，SV-GP rnK=1刚刚“探索”了许多场景（参见ZSVin第4.3节的讨论），但并未真正使用。

由此产生的更大的设计规模转化为更大的计算开销。值得注意的是，SE-GP比ST-GP更集中，似乎无法有效地探索Z。图8说明了在k=1、20、100时，使用SV-GP方案对VAR和TVaR的sk（z）演变。随着k的增加，后验不确定度sk（zn）减小；此外，我们还看到了内部模拟的目标位置，对于VaR，在分位数qα（真秩50）附近sk（zn）最低，而在左尾（真秩≤ 50）用于TVaR。两个变量0.005TVaR0.005之间有明显区别。图8：使用SV-GP，在k=1、20、100阶段的后验GP标准偏差sk（zn）进行一次运行，根据f的真实秩排序（z1:N）。y-轴位于对数刻度上。较低的sk（zn）表明在zn的空间邻域中，内部模拟的密度较高。面板，带VaR-0.005关注等级50多于等级1-10，反之亦然，T V aR0.005。图9中也观察到了这种影响：VaRα在qα处的方差和偏差最小，TVaRα在整个左尾处的方差和偏差最小。6.3空间建模VAR TVAR的收益图9：使用ST-GP，K=100后单次运行的输出。两个面板都显示了真实的f（z），以及点估计值mK（z）和'ynK，以及95%可信/置信区间，所有这些都是根据f的真实等级（z1:N）排序的。在每个“ynKis”周围有一个95%的误差条，使用hetGP估计值τ（zn）/rnkforth模拟方差，而每个mK（zn）都有基于sK（zn）的误差条。图9显示了典型ST-GP运行K=100轮结束时GP替代物的空间特征。回想一下，我们的hetGP模型平滑了样本平均值和样本方差τK（zn）。这有双重效果，即提高学习f（zn）的准确性和降低该估计的不确定性。

实际上，在图中我们可以看到，相对于经验“ynK”，mK（zn）非常接近f（zn），而且，GP后验方差sK（zn）远低于经验^τK（zn）。我们强调，该图将底层的2-D空间简化为一维实现，因此相邻秩的z在空间上可能相距很远，因此具有非常不同的（mK（z），sK（z））值。特别是，误差条较大的场景通常在空间上远离真实的分位数轮廓（TVaR情况下为左尾），因此算法发现它不值得采样，并保持rnKlow。对于RNK高的情况，两个误差条更接近。两个面板的右侧都明显存在较大的不确定性和偏差，因为算法明确针对左尾，因此在其他地方具有较高的准确性。VaR的最左端也观察到了这一特征。为了更好地量化学习R的改善，由于图9中暗示的跨场景信息的空间借用，我们首先比较了两种非自适应算法：标准嵌套抽样（称为U1-SA）与U1-GP，前者通过GP替代物使其平滑。这种比较体现了纯粹的“空间”效应，而不会因自适应设计而产生任何混淆。表4报告了各种模拟预算的结果。我们可以看到，使用空间模型可以在改善偏差和RMSE方面将速度提高约5倍。例如，N=500000的U1-SA的RMSE为152.49，与模拟次数仅为五分之一的U1-GP的theRMSE相当（N=100000的为149.09）。对于真正大的预算，由于样本平均值开始产生f（z）的合理估计，差距缩小了。

表4还比较了样本分位数^τ（^RHD）（由U1 SAto计算标准误差）处的平均局部经验方差与^RHDK基于GP的平均标准误差。我们观察到，对于相同的预算，不确定性减少了50%以上，这再次归功于空间借用。U1-SA U1-GP BR-SAN^τ（^RHD）RMSE s RMSE^τ（^RHD）RMSE1·10-6578.12 585.01 2965.05-2·101263.18 3660.92 359.05 1385.44-5·10650.62 1569.38 202.97 344.55 135.60 186.821·10410.72 759.87 146.49 148.62 53.50 61.582·10272.09 384.27 107.35 112.63 34.72 43·165 62.42 163.86 81.26 77.22 19.45 21.581·10112.89 92.44 54.28 66.69 12.98 13.87表4：二维黑色Scholes投资组合案例研究，使用嵌套MC（U1-SA）、Broadie等人[8]的BR-SA算法和统一GP替代物U1-GP，对^rke的后验标准差和RMSE进行100次宏观复制的平均值。U1-SA是基于样本平均数的标准嵌套蒙特卡洛法；U1-GP使用基于单阶段均匀分配的模拟输出GP，BR-SA使用MC样本平均值以及顺序R&S程序，见第5.4节。对于N=10，内部模拟预算为rn=1，不允许对-SA方法的τ（qα）进行估计。从不同的角度，我们可以将BR-SA方法与ST-GP方法进行比较，BR-SA方法对内部模拟进行非均匀分配。虽然BR-SA显著优于U1-SA（例如，N=5·10的BR-SA与预算大10倍的U1-SA相比），但它仍然远远落后于基于代理的顺序方法。BR SA需要10-20倍的预算才能开始匹配上一节中仅使用n=10的方法的RMSE。因此，虽然BR-SA基本上没有回归/顺序设计开销，但增加模拟效果的成本是巨大的。

综上所述，在本例中，空间建模产生了数量级的模拟节省。（我们还注意到，通过其构造，BR-SA首先在每个场景分配两个内部模拟，然后才进入其自适应阶段≤ 2N降低为U1-SA；即使在小型模拟预算下，我们的方案也会产生自适应分配。）6.4 GP仿真器的选择要使用GP仿真器，必须选择协方差核以及将模型拟合到仿真输出的特定方法。为了研究各自的影响，我们比较了两种不同的核族Mat'ern-5/2与高斯核族的使用，参见（14）和（32）以及两种不同的建模模拟方差τ（z）的方法。具体而言，我们比较了使用点估计^τ（zn）来估计τ（zn）（在DiceKriging R包中实现）的随机克里格方法和hetGP的联合响应噪声面方法。后者更复杂，因此带来更多的计算开销。表5提供了使用ST-GP采集功能和上述三种GP模拟器选择时的最终RMSE以及几个阶段的偏差。BiasApproach Kernel RMSE k=1 k=10 k=20 k=50 k=100hetGP Mat'ern-5/2 57.52 166.065 113.925 97.751 37.158 28.066hetGP Gaussian 68.06 91.475 104.427 74.874 52.716 48.249SK Mat'ern-5/2 69.31 914.728 206.267 113.716 69.821 48.670表5：对于二维Black-Scholes投资组合案例研究，不同TGP模型/Kernel系列的平均RMSE。我们还报告了一系列中间阶段k=1、10、20、50、100的平均偏差偏差（^Rk）。所有方法都使用ST-GP规则，并基于100个宏复制。我们观察到，使用经验方差估计^τ（zn）（SK）在初始阶段会产生明显的高偏差。反过来，hetGP模型能够有效地消除偏差。

虽然随着算法的发展，偏差会迅速下降，但SK emulator在k=100阶段仍有50%以上的偏差（和20%以上的RMSE）。因为在实践中，人们不知道偏差何时得到有效缓解，这是使用SK仿真器的一个显著缺点，并强调了正确建模噪声表面τ（·）的重要性。在比较核函数时，高斯核函数与Mat'ern-5/2相比，显示出略高的RMSE和偏差，尽管没有明显的边缘。这与民间传说的观点相吻合，即Mat'ern是kernelfamily的默认选择。7案例研究：随机利率和死亡率下的终身年金为了检查在更复杂的环境中，每种方法的相对性能是否保持不变，我们转向更高维度、更复杂支付函数的案例研究。尽管复杂性增加，但我们仍保持相同的预算N=10，这预计会扩大拟议方法之间的相对差异。本节中的设置考虑到年金受益人签订合同，在S年内开始收取款项，从S年开始每年持续支付，直到个人死亡。（实际上，为最终付款设置了一些折扣。）法规要求对保险公司本合同的T=1年价值的分位数进行分析。投资组合损失的主要驱动因素是利率风险（低利率增加了年金支付的现值）和死亡率风险（延长寿命提高了年金的价值）。这些因素由6-D随机状态捕获：三因素利率模型和三因素随机死亡率模型。我们希望在地平线T上找到该年金的现值≤ S、 付款开始前atS。为此，设T（x）为今天x岁个体的剩余随机寿命。

t处的年费支付基于t（x）>t，其可能性P（t，x）表示asP（t，x）=P r（t（x）>t）=1-t型-1Xu=0P r（u<T（x））≤ u+1）。=1.-t型-1Xu=0q（u，x+u），（36），其中q（u，x+u）是死亡率，即年龄为x+u的个体在u年死亡的概率（或者解释为今天年龄为x+u和x+u+1的个体死亡）。综合特殊的死亡率经验（在（36）中以P r（·）表示），我们关注影响未来死亡率演变的系统因素。因此，我们将q（u，x+u）视为随机死亡因子（Zt）驱动的arandom变量，即q（Zt；t，x）。这意味着，（36）中的生存概率P（Z[T，T]；T，x）取决于水平线和总年份T之间的整个路径Z[T，T]。设βT为时间T的瞬时利率，它也取决于（Zt）的（其他）分量。在T（以T（x）>T为条件，并假设寿命和利率风险是独立的）下，年金的净现值为f（z）=-E“xu-xXt=Se-RtTβuduP（Z[T，T]；T，x）ZT=z，T（x）>T#。（37）对于生存概率，我们选择以下模型，在Cairns等人【9】中称为（M7）：logit q（Zt；t，x）=κt+κt（x- \'\'x）+κt（十）- (R)x）- ^σx+ γt-x、 （38）在（38）中，随机驱动因素是κjt，j=1，2，3以及γt-x；其余为固定参数。具体而言，(R)x是模型适用的平均年龄，^σxis是适用年龄的平均值xof（x- 被解释为年龄效应。该模型允许为x年龄添加一个插件，而突变因素是周期效应κIt，它捕获了历年的死亡率演变，以及γt-X这是队列效应。这些构成了离散时间ARIMA模型，遵循通常的选择，即每个κII是一个带漂移的随机游动。

通常，γ·作为ARMA模型。为了再现性，我们使用R软件包StMoMo【37】，其中包含用于拟合和模拟的工具（38）。也就是说，我们利用了数据中包含的英格兰和威尔士（E＆W）死亡率，以及之前凯恩斯等人研究的死亡率。由于P（Z[T，T]；T，x）的表达式（36）具有在（u，x+u）年龄/年评估的方程式（38），因此它导致了固定队列效应γx，与ARMA过程的历史值相匹配。（βt）的利率动态来自【13】，通过随机波动率ζ和随机漂移αtdβt=（(R)β）的三因素Cox-Ingersolross模型定义- αt）dt+pβtζtdWβt，（39）dαt=（(R)α- αt）dt+√αtζtdWαt，dζt=（(R)ζ- ζt）dt+pζtИdWζt，其中Wβ、Wα和Wζ是独立的标准布朗运动。总的来说，这意味着一个6-D马尔可夫状态过程Zt=（βt，αt，ζt，κt，κt，κt）。因此，获得收益实现Yi（zn，i·）的程序是首先模拟路径（zn，it）t≥t来自分发（Zt）t≥T | ZT=zn，插入（κjt）T的实现≥T、 j=1，2，3到（38）来评估（36），并最终使用这些生存概率以及模拟的（βT）T≥Tto计算（zn，i·）=xu-xXt=Se-RtTβn，iuduP（zn，i[T，T]；T，x）。（40）注意，κ是离散时间的，而利率模型使用连续t。对于后者，我们使用一种简单的带离散化的前向Euler方法t=0.1。总的来说，一个现金流（·）的评估大约需要0.01115秒（即10个整体内部模拟大约需要2分钟），而第一个案例研究的评估需要0.000513秒。

这种计算成本较高的模拟器更能代表现实生活中的模拟引擎，并减少了模拟器安装、选择和预测方面的间接成本影响。7.1结果对于本节的其余部分，我们考虑水平T=1，并分析方程式（37）中的f（Z），即未来一年的终身年金净现值。我们取x=55，年号开始日期为S=10，即个人在x=65岁退休。如前所述，我们模拟附录B中的算法1，以确定N=10的场景集Z。方程（39）中的利率参数为“β=0.04，”“α=0.04，”“ζ=0.02，φ=0.05，E＆W死亡率模型在年龄范围x上进行了拟合∈ [55，89]使用StMoMo包[37]。与之前的案例研究一样，GP超参数在阶段10、20、…、，100，并且在3 GP和U2-GP的每一轮之后。对于平均函数，我们拟合常数u（z）≡ β.在设置（37）下，f（z）没有闭合形式的评估，因此我们通过模拟获得基准RB。该值通过运行SR-GP方法确定，N=2·10，K=200轮（因此rk=1·10；我们也采用非常保守的L=1，U=200）。结果为RVaRB=-16.0498和RTVaRB=-16.3739，估计标准偏差为s（RVaRB）=0.00202，s（RTVaRB）=0.00188。我们重复第6节的方法，对HVAR0.005和TVaR0.005执行100次固定Z的宏复制。

偏差和后验不确定度的箱线图如图10所示，估计值标准偏差、RMSE和最终设计尺寸的数值如表6所示。VaR0.005TVaR0.005SD（^RHDK）s RMSE | DK | SD（^RK）s RMSE | DK | TimeST GP 0.0394 0.0455 0.0403 151.83 0.0461 0.0404 0.0472 147.10 330SE-GP 0.0427 0.0493 0.0459 143.27 0.2853 0.2717 0.2850 101.62 295SV-GP 0.0382 0.0406 0.0380 497.81 0.0408 0.038 93 0.0407 254.35 403SR-GP-10.0467 0.0470 0.0485 135.03 0.0430 0.0402 0.0430 184.73 219SR-GP-2 0.0391 0.0437 0.0434 217.36 0.0447 0.0431 0.0450 224.96 198A3-GP0.0434 0.0497 0.0436 298.15 0.0464 0.0490 0.0461 298.55 115U2-GP 0.0598 0.0598 0.0596 194.03 0.0684 0.0601 0.0689 195.81 112U1-GP 0.5020 0.4156 0.5853 100.4940 0.4705 0.6709 10表6：基于100次宏观复制的6天人寿年金案例研究结果。我们报告了^R[1:100]K的样本标准偏差、平均GP后验标准偏差s和^RK的RMSE，以及每种方法的平均最终设计尺寸。最后一列以秒为单位报告运行时间。方法说明见表1。在相对绩效方面，结果与第一个案例研究略有不同。对于VAR0.005，SV-GP再次具有最低的RMSE值，尽管ST-GP仅差6%。有趣的是，A3-GP的表现和SR-GP-2一样好。由于A3-GP分配ramong 200个位置，这表明SR-GP-2的带宽L=26，U=75对于这种更复杂的设置来说可能太紧了，特别是因为这一次它只略优于激进的SR-GP-1。与之前一样，U1-GPFailes（RMSE比SV-GP高出15倍以上），虽然U2-GP有很大的改进，但与其他方法相比，它的RMSE仍然高得多。图10中的方框图显示，ST-GP、SV-GP和SR-GP-2提供了最稳定的不确定度量化（s（^RHDK）的紧密分布和接近实际D（^RK））。

与第一个案例研究一样，SE-GP在^RK及其报告的标准误差中会导致更多的跨期离散。TVaR结果相似；我们评论了一些明显的差异。SV-GP仍然具有最低的RMSE和最低的SD（^RK）。另一方面，SE-GP在其采集功能和不稳定s（^RK）方面也存在同样的问题。SR-GP方法的表现次之（尽管在s（^RK）不稳定的情况下，其连续确定性定量的可靠性明显较低），A3-GP和STGP紧随其后。我们将SV-GP的相对“胜利”归因于其更具差异性的分配，即更大的DK，这表明在更复杂的例子中估计TVaR需要更广泛的撒网，而不是仅针对尾部。总的来说，本案例研究证实了基于GP的顺序方法的性能，即使在高维环境中也是如此。估计的^RKGP后验不确定度s（^RK）VaR0.005TVaR0.005Figure 10：对于终身年金案例研究，最终^rkeptimates（左）和s（^RK）（右）的箱线图显示了每种方法的100次宏观复制。第一行：风险价值；最下面一行：TVaR。左侧的范围虚线是风险度量的参考值。8结论我们首次全面检查了在一组固定的外部场景下使用高斯过程模拟器和嵌套模拟来估计VaRα和Tvarα的顺序方法。无论是假设水平（VaRα）已知【30、4、15】和/或无噪声输出【26】，使用非顺序算法【27】，还是不使用空间模型进行平滑【8】，仿真文献中现有的方法都无法满足该应用。我们记录了相对于普通嵌套蒙特卡罗的两个数量级的节约，由于自适应分配，节约了大约10倍，由于空间建模，节约了5倍。

这些收益对于小型模拟预算来说是最强的N≈ N如果现有方法基本上失败或不适用。作为另一个优势，使用GP模拟器还可以得到本质上无偏的风险估计器（与有偏的非光滑版本相比），并提供标准误差的可靠估计。我们还利用了HETGP中新引入的heteroskedastic GP框架【6】。当重复计数rk较低时，噪声面τ（·）的改进估计会减少^f中的偏差，并产生一个更稳健的^RkThrough估计量。所有全序贯方法都比简单的两阶段法和三阶段法有显著的改进【27】。我们将GP仿真与顺序设计相结合的主要缺点是回归/预测开销，但通过批处理和仔细选择候选集Zcandk，这一开销显著减少。此外，在调用Y（·）代价高昂的更现实的示例中，计算采集函数、拟合和预测所花费的总时间变得微不足道。我们确定了两种性能最好的方法。SV-GP是[27]的完全顺序适应，在两个案例研究中，在RMSE方面表现最好。SV-GP还导致输出的后验不确定度sK（^RK）的变化最小，使该值成为真实采样标准偏差的可靠代理。ST-GP的表现也差不多，事实上在前几轮中表现最好，参见图6的扇形图。SV-GP的实现涉及到一个单独的辅助优化问题，并带来了更多的开销和编码复杂性。相比之下，ST-GP的标准很简单，并且还将仿真器模型大小保持在最小（Dk是GP开销的关键决定因素）。

一个更简单/更快的方案是基于银行的SR-GP，该方案也表现良好，但其结果在很大程度上取决于土地U的选择，而土地U的选择又因案例研究和VaRα与TVaRα的不同而不同。这是SV-GP和ST-GP的另一个优点，它们不需要用户指定的参数。SE-GP的混合性能提供了一个警示，即设计良好的采集功能对整体性能很重要。我们得出的结论是，当预算N非常有限时，ST-GP是一个很好的选择，而对于更大的预算，像SV-GP这样的方案是最好的。有多种方法可以进一步发展。通过微调各前瞻变量优化期间使用的权重w，可以进一步优化表现最佳的SV-GP。例如，【27】中的相关三阶段方法使用通过模拟获得的经验权重估计，以确定TVaR估计器的权重。修改优化标准以避免rnk≥ 1将消除场景探测并大幅缩小SV-GP的设计尺寸| Dk |。或者，自适应选择RK将减少仿真开销。相反，可以扩展ST-GP或SE-GP方法，在ZK+1和rk，即将配料内部化到设计施工中（而目前，相应的采集功能在技术上是一步先行）。也欢迎对^RK偏差的各种收敛速度进行理论分析，尽管这对已知阈值L的等高线查找这一更简单的问题来说已经是一个极大的挑战。初始化阶段也需要进一步分析。

回想一下，我们使用的固定预算为r=0.1N，通过空间填充方案确定的九个代表性场景。原则上R应尽可能小，让算法本身分配模拟的剩余部分，但在实践中，GP超参数的合理起点至关重要。更好的规则/启发式选择R需要考虑。此外，在Z的几何极值点的意义上，可能需要进一步关注各自的极值区域，而不是给定Z的空间，例如，参见【11】中统计深度函数的使用。最后但并非最不重要的一点是，可以对我们框架下的GP模拟器进行各种调整。回想一下，代理作出了许多假设，例如高斯观测噪声和协方差结构中的空间平稳性。在实践中可能会违反这些规定；如果模型规格不准确是一个问题，那么可以考虑在GP生态系统GP GLM（例如，t分布噪声）、树状图和本地GPs中构成一个活跃研究领域的众多替代方案。参考文献【1】Bruce Ankenman、Barry L Nelson和Jeremy Staum，《模拟元建模的随机克里格法》，运筹学，58（2010），第371-382页。[2] 巴塞尔委员会，《交易账簿的基本审查：经修订的市场风险框架》，咨询文件，2013年10月版。[3] Daniel Bauer、Andreas Reuss和Daniela Singer，关于基于嵌套模拟的solvencycapital需求计算，ASTIN公告：IAA杂志，42（2012），第453-499页。[4] Julien Bect、David Ginsbourger、Ling Li、Victor Picheny和EmmanuelVazquez，《故障概率估计的计算机实验顺序设计》，统计与计算，22（2012），pp。

773–793.[5] Mickael Binois和Robert Gramacy，《hetGP：复制下的异方差高斯过程建模与设计》，2017年。CRAN包v.1.0.0。[6] Mickael Binois、Robert B Gramacy和Michael Ludkovski，《大型模拟实验的实用异方差高斯过程建模》，技术报告，arXiv预印本XIV:1611.059022016。[7] Kurt M Bretthauer、Anthony Ross和Bala Shetty，《分层抽样中优化分配的非线性整数规划》，欧洲运筹学杂志，116（1999），第667-680页。[8] Mark Broadie、Du Yiping和Ciamac C Moallemi，《通过nestedsequential simulation进行有效风险评估》，《管理科学》，57（2011），第1172-1194页。[9] Andrew JG Cairns、David Blake、Kevin Dowd、Guy D Coughlan、David Epstein和Marwa Khalaf Allah，《死亡率密度预测：六种随机死亡率模型的分析》，《保险：数学与经济学》，48（2011），第355-367页。[10] Tony F Chan、Gene Howard Golub和Randall J LeVeque，《更新公式和成对算法以计算样本方差》，COMPSTAT 1982年第五届研讨会，1982年，图卢兹，斯普林格，1982年，第30-41页。[11] Matthieu Chauvigny、Laurent Devineau、St'ephane Loisel和V'eroniqueMaume Deschamps，《多因素快速远程但非极端分位数：对偿付能力II和企业风险管理的应用》，《欧洲精算杂志》，1（2011），第131-157页。[12] Chun Hung Chen、Donghai He、Michael Fu和Loo Hay Lee，《选择最佳子集的有效模拟预算分配》，通知《计算杂志》，20（2008），第579-595页。[13] Lin Chen，《随机均值和随机波动率：利率期限结构的三因素模型及其在衍生品定价和风险管理中的应用》，载《利率动力学、衍生品定价和风险管理》，第卷。

《经济学与数学系统讲稿》第435页，柏林斯普林格，1996年，第1-36页。[14] Xi Chen、Barry L Nelson和Kyung Kuk Kim，《条件价值风险及其敏感性的随机克里格法》，2012年冬季模拟会议论文集，IEEE，2012年，第1-12页。[15] Cl'ement Chevalier、Julien Bect、David Ginsbourger、Emmanuel Vazquez、VictorPicheny和Yann Richet，基于快速并行克里格法的逐步不确定性减少，并应用于偏移集的识别，Technometrics，56（2014），第455-465页。[16] Cl'ement Chevalier，Victor Picheny和David Ginsbourger，Kriginv：基于kriging的批量顺序反演策略的高效且用户友好的实现，《计算统计与数据分析》，71（2014），第1021-1034页。[17] Marcus C Christiansen和Andreas Niemeyer，《偿付能力II中偿付能力资本要求的基本定义》，ASTIN公告：IAA杂志，44（2014），第501–533页。[18] Alexander Forrester、Andras Sobester和Andy Keane，《通过代理建模进行工程设计：实用指南》，John Wiley&Sons，2008年。[19] Peter I Frazier，《区分区排序和选择的全序列淘汰程序，正确选择概率的严格限制》，《运筹学》，62（2014），第926-942页。[20] Gene H Golub和Charles F Van Loan，《矩阵计算》，第3卷，JHU出版社，2012年。[21]Michael B Gordy和Sandeep Juneja，《投资组合风险度量中的嵌套模拟》，《管理科学》，56（2010），第1833-1848页。Robert B Gramacy和Herbert KH Lee，《高斯过程和极限线性模型》，计算统计与数据分析，53（2008），第123-136页。[23]Frank E Harrell和CE Davis，一种新的无分布分位数估计量，Biometrika，69（1982），pp。

635–640.【24】Bogumil Kami'nski，《随机克里格元模型的更新方法》，欧洲运筹学杂志，247（2015），第859–866页。[25]Joseph Hyun Tae Kim和Mary R Hardy，《量化和纠正估计风险度量中的偏差》，Astin Bulletin，37（2007），第365–386页。[26]T Labopin Richard和V Picheny，《估计黑箱函数百分位数的序贯实验设计》，中国统计局，即将出版（2016）。arXiv预印本XIV:1605.05524。[27]Ming Liu和Jeremy Staum，《有效嵌套模拟预期短缺的随机克里格》，风险杂志，12（2010），第3-27页。【28】Amandine Marrel、Bertrand Iooss、Franc,ois Van Dorpe和Elena Volkova，《利用高斯过程建模复杂计算机代码的有效方法》，计算统计与数据分析，52（2008），第4731-4744页。[29]Jeremy Oakley，《估算不确定计算机代码输出的百分位数》，《皇家统计学会期刊：C辑（应用统计学）》，53（2004），第83-93页。[30]Victor Picheny、David Ginsbourger、Olivier Roustant、Raphael T Haftka和Nam Ho Kim，《精确逼近目标区域的自适应实验设计》，机械设计杂志，132（2010），第071008页。卡尔·爱德华·拉斯穆森和克里斯托弗·威廉姆斯，《机器学习的高斯过程》，麻省理工学院出版社，2006年。[32]James Risk和Michael Ludkovski，《长寿风险产品定价和享乐的统计模拟器》，《保险：数学与经济学》，68（2016），第45-60页。【33】Olivier Roustant、David Ginsbourger和Yves Deville，《DiceKriging，DiceOptim：通过基于kriging的元建模和优化分析计算机实验的两个R包》，统计软件杂志，51（2012），pp。

1–55.[34]Thomas J Santner、Brian J Williams和William I Notz，《计算机实验的设计与分析》，斯普林格科学与商业媒体，2013年。[35]Michael E Sfakianakis和Dimitris G Verginis，《非参数估计量的新家族》，《统计学中的通信：模拟与计算》，37（2008），第337-345页。[36]Simon J Shecker和James Stephen Marron，《核分位数估计》，美国统计协会杂志，85（1990），第410-416页。[37]Andres M.Villegas、Pietro Millossovich和Vladimir K.Kaishev，StMoMo：随机死亡率建模的AnR包，2017年。R软件包版本0.4.0。方差最小化计算引理1。随着rnk+1的增加，我们有以下近似值：（C+candk+1）-1.≈（C+（k）-1+（C+（k）-1(k- candk+1）（C+（k）-1.（41）证明。回想一下candk+1是对角线，输入τ（zn）/（rnk+r0nk）和kis对角线，条目为τ（zn）/rnk。我们重写candk+1=k+Bk+1(-一） Bk+1，其中Bk+1是元素bnn的对角矩阵τ（zn）rrnk公司-rnk+r0nk我是identitymatrix。通过Woodbury矩阵反演公式[20]，我们得到（C+candk+1）-1=（C+k+Bk+1(-一） 黑色+1）-1=（C+（k）-1.-（C+（k）-1Bk+1黑色+1（C+（k）-1Bk+1- 我-1Bk+1（C+（k）-1.（42）当rnk+1较大时，Bk+1和khave项目相对较小，因此Bk+1（C+（k）-1Bk+1≈因此，（42）中的中间矩阵逆为≈ (-一） 和（C+candk+1）-1.≈ （C+（k）-1+（C+（k）-1Bk+1Bk+1（C+（k）-1.