Wishart多维随机波动率模型的精确模拟

该算法生成赫斯顿模型的状态变量XT和YTO。步骤（1）根据XT的分布生成样本步骤（2）根据给定XT的I=RTXtdt的条件分布生成样本步骤（3）生成标准正态随机数Z步骤（4）设置YT=y+ρσXT- 十、+ （r）-κθρσ）T+ρκσ-I+p（1）- ρ） 伊兹。对于步骤（1），需要根据XT的分布生成样本。幸运的是，XT的分布在文献中是众所周知的（例如，见Glasserman[20]）。xT定律可以表示为xT=σ（1- E-κT）4κχ′2δ4κe-κTσ（1）- E-κT）x,其中δ=4κθσ，χ′2δ（λ）表示具有δ自由度的非中心卡方随机变量和非中心参数λ。有关非中心卡方随机变量的采样，请参阅[20]。布罗迪和卡娅的傅里叶反演技术思想在第（2）步开始发挥作用。他们表明，给定xT=xtdt时，txtdt的条件特征函数φ（·x，xT）可以写成φ（λ| x，xT）=Ex，yhexpiλRTXtdtXT=xTi=γ（λ）e-(1/2)(γ(λ)-κ） T（1）- E-κT）κ（1）- E-γ（λ）T）×expnx+xTσhκ（1+e）-κT）1- E-κT-γ（λ）（1+e）-γ（λ）T）1- E-γ（λ）TioI0。5δ-1.√xxT4γ（λ）e-0.5γ（λ）Tσ（1）-E-γ（λ）T）I0。5δ-1.√xxT4κe-0.5κTσ（1-E-κT）,式中γ（λ）=√κ- 2σiλ。然后，他们对条件特征函数φ（·| x，xT）进行数值反转，以获得条件分布Px的近似值RTXtdt≤ v|XT=XT≈ Fh，N（v | x，xT）=hvπ+πNXn=1sin（hnv）nRe|（hn | x，xT）, （13） 其中h>0是离散化步长。最后，他们将逆变换方法应用于Fh，N（·| x，xT）来模拟积分RTXTDT，也就是说，他们生成了一个统一的随机数，并对RTXTDT的以下g方程进行了数值求解RTXtdt | x，xT= U.4精确模拟方法在本节中，我们详细介绍了WMSV模型的精确采样算法。

由于过程（X，Y）是一个时间齐次马尔可夫过程，给定初始值（X，Y）=（X，Y）的（XT，YT）的精确模拟方法可以直接推广到给定任意t<t的（XT，YT）的精确模拟方法。因此，我们只考虑给定初始值（X，Y）=（X，Y）时（XT，YT）的模拟。对于WMSV模型，很难实现Broadie和Kaya方法的简单扩展。困难来自波动因子过程X的维数。在Broadie和Kaya方法中，需要从给定XT=XT的TXTDT的条件（单变量）d分布生成样本，这可以通过傅里叶反演技术实现。相反，WMSV模型的综合波动系数rRTXTDT具有d×d矩阵变量分布，这使得几乎不可能通过傅里叶反演技术从分布中生成样本。因此，我们没有遵循Broadie和Kaya的方法，而是通过定理2.3直接实现目标。我们的方法大致包括以下两个步骤。算法2。该算法生成WMSV模型的L对状态变量（XT，YT）。步骤（1）从XT的分布生成L个样本X（1）T，···，X（L）T步骤（2）对于每个L=1，···，L，从Yt的条件分布生成一个样本Y（L）T给定XT=X（L）在以下小节中，我们将详细介绍这两个步骤。4.1从定理2.3所示的XTA分布中取样，XTA具有自由度δ、协方差矩阵V（0，0）和非中心参数V（0，0）的矩阵的非中心Wish art分布-1Ψ(0, 0)xψ（0,0）。具有整数自由度（即δ）的Wishart分布∈ N） 在多元统计分析的文献中有广泛的研究（例如。

参见[21,27,29]）。如果δ是一个大于或等于d的整数，精确模拟的一种方法是将正规随机矩阵平方[8]。设Nt为下列等式的解Dnt=NtHdt+dBt∑，带NN=x，（14），其中B是标准δ×d矩阵布朗运动。我们可以很容易地检查NN满足随机微分方程（2）。方程（14）允许一个封闭形式的解=N+ZtdBs∑e-嘘以斯，0≤ t<∞.请注意，NTA的行是独立的正态随机向量，具有公共协方差矩阵V（0，0）和NTis净H的平均值。因此，对于整数自由度δ，可以通过采样δ×d i.i.d.正态随机变量来实现XT的精确采样。还有一种更复杂、更有效的方法，可以模拟自由度为整数的非中心Wishart分布[21]。据我们所知，Ahdida和Alfonsi[2]直到最近才设计了一种具有非整数值自由度的非中心Wishart分布的精确抽样方案。他们使用了一种分裂的方法，将愿望艺术过程分解为一个极小的生成器。他们的Rexact采样方案要求最多采样d（d- 1） /2 i.i.d.n正常随机变量和d n中央卡方随机变量。在本文的数值实验中，我们使用了Ahdida和Alfonsi[2]的精确抽样方法。4.2从Ytgevent xT的条件分布中抽样本小节专门讨论从Ytgevent=xT的条件分布中抽样的步骤。在我们的方法中，这一步是最技术和最耗时的一步。在给定的条件下，我们采用了反转技术来解释函数。

我们遵循类似的方法，但我们直接反转给定Ytxt=XT的条件特征函数，以避免转换矩阵变量随机变量的特征函数的困难。设φ（·；x，y，xT）为条件特征函数，即φ（λ；x，y，xT）=Ex，yheiλYTXT=xTi，表示λ∈ R.设F（·；x，y，xT）为相应的分布函数：F（v；x，y，xT）=Px，yYT≤ 五、XT=XT.通过Levy的反演公式，可以从特征函数中恢复分布，即-∞ < l<v，我们有F（v；x，y，xT）=F（l；x，y，xT）+2πZ∞-∞ν（λ；x，y，xT）e-iλl- E-iλviλdλ=F（l；x，y，xT）+πZ∞Reh~n（λ；x，y，xT）e-iλl- E-iλviλidλ=F（l；x，y，xT）+πZ∞Imh~n（λ；x，y，xT）（e）-iλl- E-iλv）idλ。注意，我们可以通过取一个小的l，使F（l；x，y，xT）小到可以忽略。因此，积分项在上述表达式中占主导地位。上面的积分可以用数值积分法来近似。我们用梯形法则数值计算分布函数。众所周知，梯形规则适用于振荡积分，因为误差往往会被抵消（见阿巴特和惠特[1]第4节）。整条正实线上的积分可以用以下方法近似。为了简单起见，我们把被积函数写成g（λ）。Z∞g（λ）dλ=g（0+）h+h∞Xn=1g（nh）- ed（h）=g（0+）h+hNXn=1g（nh）- 教育署（h）- et（N），其中ed（h）和et（N）分别是由于连续变量的离散化和有限和的运行而产生的误差。每个误差都可以通过取足够小的h和足够大的N来任意减小，我们在第4.3节中详细讨论了这些误差的控制。人们很容易发现g（0+）=v- l。

因此，我们用以下有限三角级数来近似分布函数：Fh，N（v；x，y，xT）=h（v-l）2π+πNXn=1Imh~n（nh；x，y，xT）n（e-inhl- E-inhv）i=h（v-l）2π+πNXn=1重新ν（nh；x，y，xT）]nsin（nhv）- 辛（nhl）(15)-πNXn=1伊姆河ν（nh；x，y，xT）]ncos（nhv）- cos（nhl）.我们讨论了如何计算条件特征函数φ（λ；x，y，xT）。条件特征函数可通过取u=-式（7）中的iλ。公式（7）包括普通微分方程组（8）的解ψ、φ、ψ和dV。如第2.3节所示，系统（8）允许Heston模型的封闭形式解。一般来说，由于矩阵乘法的非交换性，系统（8）不允许闭式解。但系统（8）可以通过数值方法有效地求解。特别是，我们在实验中使用了MATLAB函数ode45来求解方程。需要注意的是，函数ψ、φ、ψ和V需要在（0,0），（0，-ih），··，（0，-iNh），这样的评估点可以在样本X（1）T，···，X（L）T上均匀地选择。因此，对于所有模拟运行，仅在这些网格点处求解系统（8）一次就足够了，如果模拟运行的数量L相对大于N，则此数值步骤不会造成计算负担。表达式（7）需要计算复数的幂，通常是一个多值函数：（z）ν=（|z|eiArg（z））ν=（|z|ei（Arg（z）+2mπ））ν=|z|eiνArg（z）+2mνπ，m∈ Z、 主要论点-π<Arg（z）≤ π. 如果λ7→ z（λ）是R上不为0的复值连续函数，则存在唯一的连续函数λ7→ arg（z（λ））使得-π<arg（z（0））≤ π.

通过追踪主变元Arg（z（λ））并在不连续点处加或减2π，可以很容易地构造这样一个连续。我们可以用∧的幂| t | n | n | n | n | n | n | n | n | n | n | n | n | n | n | n。很明显λ7→ detV（0，-iλ）是一个不能达到0的连续函数，我们可以应用上述观察值来计算（detV（0，-iλ）δ/2。MATLAB函数ode45基于Dormed and Prince[15]的显式R unge-Kutta（4,5）公式。条件特征函数（7）涉及矩阵参数的超几何函数。矩阵参数的超几何函数由以下幂函数定义f（b；y）=∞Xk=0X |ι|=k（b）ιCι（y）k！。（16） |ι和（b）ι的定义见A.1节。分区多项式Cι（y）不是矩阵y的多项式，而是其特征值的多项式。所以，我们有时用下面的方式写超几何函数和区域多项式：Cι（y）=Cι（α，·α，αd），F（b；y）=F（b；α，·α，αd），其中α，·α，y.Koev和Edelman[26]的α-dare特征值利用分区多项式和广义超几何系数的组合性质，开发了计算（16）的截断版本的分析算法：mF（b；α，·α，αd）=mXk=0X|||=k（b）ιCι（α，·αd）k！。由于级数（16）的分母增长速度快于阶乘阶，级数收敛得很快，截断的形式给出了一个很好的近似值。在第4.3节中，我们提供了详细的截断错误分析。Koev和E delman算法的MATLAB实现可在作者的主页[25]中找到。但是rou tin e只得到真正的输入参数。因此，我们修改了适用于复杂输入参数的代码。为了对原木价格YT进行采样，我们使用了逆变换方法d。

我们生成一个统一的随机数U，然后数值求解方程YTFh，N（YT；x，y，xT）=U。（17）该方程可以通过数值方法有效地求解，例如牛顿法，因为函数Fh，是一个严格递增函数。我们将在以下小节中提供详细信息。4.3一些实施问题为了实施我们的方法，我们需要解决一些问题。我们必须在（15）中确定l、h和N的适当值，并在有限系列（16）中确定近似的求和数。我们还需要解决如何解方程（17）。在我们的方法中，仔细选择l、h和N的值是至关重要的，因为它们决定了计算特征函数φ（·；x，y，xT）的点。由于条件分布函数的近似误差，太少的评估点可能会导致蒙特卡罗模拟的大偏差。太多的点可能会使我们的方法太慢，因为特征函数的计算是我们方法中最耗时的步骤。如前一小节所述，近似值（15）涉及三种不同类型的误差：离散化误差ed（h）、截断误差或et（N）和左尾概率F（l；x，y，xT）。我们提出了一种基于给定XT=XT的条件均值和标准差来控制这些误差的方法。回想一下，通过对特征函数φ（·；x，y，xT）：u（xT）=Ex，y[YT | xT=xT]=Im进行微分，可以很容易地找到平均值和标准偏差ν′（0；x，y，xT）,σ（xT）=Ex，y[（YT- u）|XT=XT]=-重新ν′（0；x，y，xT）-Im（ν′（0；x，y，xT））.在Abate和Whitt[1]的第5节中，他们使用泊松求和公式证明离散化误差由两个简单的尾部概率0决定≤ 教育署（h）≤1.-F（2π/h+l；x，y，xT）+ F(-2π/h+v；x、 y，xT）代表v-l<2π/h。

假设我们想要计算F（v；x，y，xT），误差小于。设F（l；x，y，xT）≤ /4和F（u；x，y，xT）≥ 1.- /4. 取h=2πu-l。那么，对于l≤ 五、≤ u，2πh+l=u，和-2πh+v≤ -2πh+u=l，所以| F（l；x，y，xT）-ed（h）|≤ max{F（l；x，y，xT），ed（h）}≤ /2.然后我们转向截断误差et（N）。由于（15）中的和是振荡的，最后一项的绝对值给出了截断误差的估计：|et（N）|≈πImh~n（Nh；x，y，xT）N（e）-iNhl-E-iNhv）我≤2 |||||（Nh；x，y，xT）|π因此，我们可以选择一个大的N，以便||||（Nh；x，y，xT）|≤ π/4，使截断误差约小于/2。这些误差范围和估计值在理论上很有吸引力，但不适合实际使用。例如，我们可以计算φ（Nh；x，y，xT），但这非常耗时。此外，我们需要计算F（·；x，y，xT），这是事先未知的。为了克服这些困难，我们利用了一种从数值实验中获得的启发式思想。给定的YTxt=XT的条件分布与正态分布大致相似，因为波动因子的终值XT已知，且随机性的主要来源是关于布朗运动的随机积分。但它们也表现出一些差异。范数分布相对于均值是对称的，但YT的条件分布是不对称的，并且YT的条件特征函数的衰减速度随着相关参数R的增加而变慢。因此，我们从正态分布所建议的l、h和N开始，然后修改它们以考虑这种影响。

正态分布建议我们采用以下形式的l、h和N：l=Φ-1./4; u（xT），σ（xT）= u（xT）+σ（xT）Φ-1（/4），h=2πΦ-1.1.-/4; u（xT），σ（xT）- Φ-1./4; u（xT），σ（xT）-1= -πσ（xT）Φ-1./4-1，N=lhσ（xT）p-2对数（π/4）m，其中Φ-1（·；u，σ）和Φ-1（·）分别是正态分布N（u，σ）和标准正态分布的反函数。和A. 表示大于a的最小整数。然后我们用平均值u（xT）和相关参数R对其进行修改。因为我们想要选择点u=-inh，n=1，··，n在样本X（1）T，··，X（L）T中均匀分布，我们以保守的方式取L，h和n：l=minl=1，··，lu（X（l）T）+σ（X（l）T）Φ-1(/4),h=-π×minl=1，··，Lc |u（X（l）T）|+σ（X（l）T）Φ-1./4-1，N=maxl=1，··，L1+c√tr[RR] hσ（X（l）T）p-2对数（π/4）,（18） 对于适当的常数c≥ 0和c≥ 特别是，我们在实验中取c=0.1和c=0.5。现在，我们考虑决定在哪里截断超几何级数（16）的问题。整数的划分可以按字典顺序排列，即，如果ι=（k，···，kd）和ι=（l，···，ld）是两个整数的划分，我们将在第一个索引iat中写入ι>如果ki>li，则部分变为不相等。利用这一阶关系，我们给出了超几何函数运行误差的上界。证据见附录A.4。提议4.1。

对于b>（d）- 1） 和一个整数m≥ D- 1，我们有α，α，αd∈ C | F（b；α，··，αd）-mF（b；α，··，αd）|≤（|α|+··+|αd |）m+1（m+1）！（b） ^ι（m+1）e |α|+··+|αd |，（19）其中^ι（m+1）是m+1划分为不超过dparts的分区中的最小分区。请注意，k分成不超过d个部分的分区中，最小的分区^ι（k）=（k，···，kd）是分布最均匀的分区，k+1的最小分区^ι（k+1）是通过将1添加到ki，或将1添加到kjdi提供给kj的kjat而获得的分区-1另一方面。例如，（2,2,2,2,2,1），（2,2,2,2,2）和（3,2,2,2,2）分别是9,10和11分成5部分的最小分区。通过这一观察，我们可以计算超几何系数（b）^ι（k）如下：（b）（1）=b和（b）^ι（k+1）=（b） ^ι（k）（b+k）如果k=····=kd（b）（k）b+kj-（j）- 1)如果k=··=kj-16=kj=···=kd，（20）其中^（k）=（k，···，kd）。我们使用界（19）作为截断级数（16）的标准。请注意，参数B是一个常数δ，满足命题4.1的假设。由于它是一个常数，我们需要在所有模拟运行中只计算一次最小分区的超几何系数，并且它不会增加我们方法的计算负担。因此，上界（19）可以通过较小的额外计算预算来计算。为了解方程（17），我们使用牛顿法。为了增强收敛性，需要仔细选择初始猜测。如前所述，条件分布F（v；x，y，xT）近似为具有平均值u（xT）和方差σ（xT）的正态分布。