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

因此，我们对方程Φ（v；u（xT），σ（xT））=U或v=Φ的解进行初始猜测-1（U；u（xT），σ（xT））。然后，我们将牛顿法应用于函数Fh，N（v；x，y，xT）以生成一个序列，该序列近似于方程（17）的解：vk+1=vk-Fh，N（vk；x，y，xT）F′h，N（vk；x，y，xT）。我们对循环进行固定次数的迭代，通常为4或5次，如果牛顿法生成的序列未能在预定的迭代次数内收敛，则使用对分搜索法。我们总结了从给定的ytxt=XT的条件分布生成样本的方法。算法3。该算法分别从给定的ytxt=X（1）T、···、X（L）T的条件分布生成L个样本Y（1）T、··、Y（L）T。步骤（1）根据（18）确定l、h和N，步骤（2）在u=-inh，n=1，···，NStep（3）对n=1，··，n和l=1，··，LStep（4）生成IID均匀随机数U（1），··，U（l）步骤（5），每个l=1，··，l，s olve Fh，n（y（l）T；x、 y，x（l）T）=U（l）对于y（l）T5数值结果在本节中，我们将我们的精确模拟方法与其他现有方法进行数值比较：WMSV模型的E uler离散化方案和赫斯顿模型的Broadie-Kaya方法。5.1精确采样与欧拉模式之间的比较随机微分方程组（2）和（3）不接受封闭形式的解。在这种情况下，模拟模型的一种方法是对时间间隔进行离散，并在这些离散的时间网格上模拟另一个与模型近似的过程。欧拉离散是一种基本的离散格式。设0=t<t<··<tN=t是时间间隔[0，t]划分为N个相等的子间隔，即。，t=ti-钛-1=T/N，i=1，·N。

我们通过设置^Xt=x和^Xti来区分Wishart过程（2）=^Xti-1+ (δΣ∑+H^Xti-1+^Xti-1小时)t+q^Xti-1.Wti∑+Wti∑Wtiq^Xti-1.+,哪里Wti=Wti- Wti-1.这里，A+表示对称矩阵A的正部分：我们设置A+=Odiag（λ+，··，λ+d）O如果A=Odiag（λ，··，λd）O正交矩阵。为了得到^Xtipositive半定义，我们在每个时间网格上取正部分。原木价格过程（3）的分解为^Yti=^Yti-1+R-tr[^Xti-1]t+trhq^Xti-1.WtiR+ 兹提皮德- RRi、 在哪里Zti=Zti- Zti-1.欧拉离散格式易于理解和实现。但它有严重的缺点：从Euler方案中提取的样本分布不同于真实分布，可能需要非常精确的离散化，以将偏差减小到可接受的程度。这些将通过数值结果加以说明。我们将从分布、收敛性和性能方面比较精确s模拟方法和Euler离散格式。本文中的模拟实验是在一台装有英特尔Core2 Quad 3.00 GHz处理器和3.25 GBRAM的台式PC上进行的，运行Windows XP Professional。我们在R2009a版本中使用了MATLAB。为了证明精确模拟法和欧拉离散法的经验分布的不同，我们生成了10个样本，并使用这些样本估计了YT的密度函数。为了估计密度函数，我们使用MATLAB函数KSDensity，该函数使用核密度估计计算密度估计。通过对特征函数（4）进行数值反演，得到了理论密度函数。

图1显示了参数集的Yt密度函数：δ=1.1，r=0，y=0，T=1.0，-2.5-2.-1.5-1.-0.50 0.5 100.20.40.60.811.21.41.61.8YT理论计算-0.2-0.15-0.1-0.05 0.05 0.1 0.15 0.2 0.25 0.30.80.911.11.21.31.41.51.61.71.8YT理论AlexactEuler-0.8-0.75-0.7-0.65-0.6-0.55-0.5-0.45-0.4-0.35-0.30.10.150.20.250.30.350.40.450.50.550.6YT理论Alexacteuler0。2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.700.20.40.60.811.21.4YT理论Alexacteulerfigure 1：原木价格的密度估计：左上面板显示了密度估计的整体形状，其他三个面板放大了左上面板的不同部分。x=0.0298 0.01190.0119 0.0108, H=-1.2479-0.8985-0.0820-1.1433,Σ =0.3417 0.34930.1848 0.3090, R=-0.2243-0.1244-0.2545-0.7230.除δ外，该模型的参数集取自Da Fonseca和Grasselli[13]。Eu ler方法的时间步数设置为25。从图1中，我们观察到，从Euler方法得出的样本分布明显不同于真实分布。密度估计的误差源于欧拉格式的离散化。更糟糕的是，没有合适的方法来测量偏差误差。相比之下，精确模拟的误差主要来自方差，并且可以很容易地通过样本方差进行测量。为了说明这种差异，我们给出了两种模拟方法得到的欧式看涨期权价格估计。除δ=3.2外，模型参数与上述相同。我们考虑了一种货币看涨期权，即履约价格K设置为等于他们的票据，DAX指数的校准δ为0.5776，这违反了δ>d的假设- 1.事实上，在这种情况下，Wishart过程X不再是Cuchiero等人意义上的有效过程[12]。方法。ofMC估计数。

误差估计时间步长模拟运行（秒）精确N/A50000 0.192415 0.1430×10-212.14100000 0.191451 0.1009×10-224.28500000 0.192143 0.4535 × 10-3121.41000000 0.191513 0.3202 × 10-3242.8欧拉50000 0.195202*0.1456 × 10-2179.4100000 0.194829*0.1034 × 10-2358.7500000 0.193827*0.4603 × 10-31793.51000000 0.194197*0.3259 × 10-33587.050000 0.194016 0.1461 × 10-2358.4100000 0.193264 0.1027 × 10-2716.8500000 0.193008*0.4570 × 10-33584.21000000 0.193073*0.3234 × 10-37168.5表1：δ=3.2的看涨期权价格的蒙特卡罗估计：理论期权价格为0.191575。星号数字是指理论价格位于95%置信区间之外的数字。S=eY=1。我们通过将Carr和Madan方法[9]应用于拉普拉斯变换公式（4）来计算期权的理论价格。表1显示了期权价格的蒙特卡洛估值。Euler方法给出的结果不可靠：理论价格仅在Euler方法建立的两个（八个）95%置信区间内。此外，欧拉方法的计算时间比精确模拟方法要长得多。随着δ的减小，欧拉方法的精度变差。对于较小的δ，离散Wishart过程更频繁地穿过s y mmetric正微分方程的圆锥体边界。每次它通过边界时，它的n负部分都会被截断，使其成为正方，这种截断可能会导致严重的偏差误差。为了说明这种情况，wesetδ=1.1，所有其他参数的设置同上。表2显示，对于δ=1.1，乌勒离散化方法很难给出可靠的估计。

理论上的代价永远不在于欧拉方法建立的95%置信区间，而在于精确模拟的95%置信区间。5.2与Broadie和Kaya方法的比较为了将我们的方法应用于Heston模型，我们使用了显式公式（12）而不是（7）。在（12）中，我们不需要求解普通微分方程，我们可以使用内置的MATLAB函数，而不是Koev和Edelman[26]的算法来计算矩阵参数的超几何函数。我们在第4.3节中指出，在我们的方法中，需要仔细选择网格大小h和截断数N。任何基于傅里叶反演的模拟方法都是如此，例如Broadie和Kaya的方法[7]。在Broadie和Kaya的论文中，他们没有给出如何选择h和N的任何指南，他们认为适当的h和N可以通过反复试验找到。但很难找到合适的h和N，这将通过数值实验加以说明。MATLAB函数besseli使用Amos[3]的算法计算第一类Iν（·）的修改贝塞尔函数。方法。ofMC估计数标准误差估计时间步长模拟运行（秒）精确N/A50000 0.114129 0.7381×10-315.34100000 0.113630 0.5215 ×10-330.68500000 0.112916 0.2320 ×10-3153.401000000 0.113085 0.1642 × 10-3306.80Euler50000 0.117013*0.7654 ×10-3180.33100000 0.117020*0.5441 ×10-3360.67500000 0.116698*0.2428 ×10-31803.341000000 0.116940*0.1718 ×10-33606.6750000 0.116940*0.7541 ×10-3360.54100000 0.115963*0.5321 ×10-3721.09500000 0.115333*0.2383 ×10-33605.451000000 0.115363*0.1684 ×10-37210.90表2：δ=1.1的看涨期权价格的蒙特卡罗估计：理论期权价格为0.113000。

星号数字是指理论价格位于95%置信区间之外的数字。本节旨在展示（13）中h和N的各种选择如何影响Broadie和Kaya方法的准确性和性能，并将其性能与我们的方法进行比较。为此，我们使用我们的方法生成了10个Heston模型样本，我们还使用Broadie和Kaya的方法为h和N的各种组合生成了相同数量的样本。我们使用它们获得了欧洲看涨期权价格的蒙特卡罗估计。模型参数设置为tox=0.010201，y=log（100），κ=6.21，θ=0.019，σ=0.61，ρ=-0.7，r=0.0319，T=1.0，K=100。这组p参数取自杜菲等人[16]：它于1993年11月2日被校准为标准普尔500的期权价格，并在布罗迪和卡亚[7]中使用。表3显示了模拟结果。从表中可以看出，h=32和N=25似乎是实验输出中h和N的最佳选择。但是，即使选择了h和N，Broadie和Kaya的方法也比我们的方法慢。通过这个实验，我们发现Broadie-Kaya方法需要很多实验才能找到h和N的最佳选择。更糟糕的是，h和N的最佳选择对模型参数的变化非常敏感。为了证明其方法的灵敏度，我们给出了除T=0.25外其他参数相同的另一个实验结果。表4显示，对于h和N的相同情况，真正的期权价格决不取决于其方法的置信区间。相比之下，我们根据（18）自动选择h和N，使我们的方法能够适应参数变化。6.结论性意见我们提出了一种精确模拟WISHART多维随机波动率模型（WMSV）资产价格和波动率因子的方法。

我们的模拟方法受Broadie和Kaya[7]的傅里叶反演技术的启发，并基于我们方法中的条件分析，h和N由（18）方法自动选择，h N Nh MC估计标准误差时间（秒）KK N/A N/A N/A 6.8125 7.4302×10-359.77BK100 800 6.8276*7.4377 × 10-3280.4200 1600 6.8059 7.4206 ×10-3591.1400 3200 6.7993 7.4208 ×10-31172.050 800 6.8162 7.4299 ×10-3123.4100 1600 6.8025 7.4220 ×10-3297.0200 3200 6.8211*7.4257 × 10-3584.825 800 6.8098 7.4302 ×10-368.050 1600 6.8011 7.4144 × 10-3132.6100 3200 6.8191 7.4286 ×10-3293.613 832 6.8346*7.4252 × 10-341.1625 1600 6.8222*7.4173 × 10-372.5250 3200 6.8184 7.4137 × 10-3130.8表3：T=1.0时买入期权价格的蒙特卡罗估计：理论期权价格为6.8061。星号数字是指理论概率位于95%置信区间之外的数字。KK指的是我们的方法，BK指的是布罗迪和卡娅的方法。方法HN Nh MC估计标准误差时间（秒）KK N/A N/A N/A 2.6720 2.9214×10-378.18BK100 800 2.7883*2.8438 × 10-3279.5200 1600 2.7597*2.8631 × 10-3608.0400 3200 2.7463*2.8705 × 10-31290.250 800 2.7853*2.8435 × 10-3124.0100 1600 2.7587*2.8664 × 10-3306.6200 3200 2.7483*2.8796 × 10-3647.725 800 2.7802*2.8408 × 10-367.350 1600 2.7610*2.8678 × 10-3138.7100 3200 2.7488*2.8798 × 10-3327.913 832 2.7845*2.8433 × 10-340.6325 1600 2.7651*2.8695 × 10-375.1950 3200 2.7445*2.8736 × 10-3150.9表4：T=0.25时买入期权价格的蒙特卡罗估计：理论期权价格为2.6709。带星号的数字是理论价格位于95%置信区间之外的数字。KK指的是我们的方法，BK指的是布罗迪和卡娅的方法。给定波动水平的资产价格的拉普拉斯变换。

所提出的模拟方法可用于生成路径无关或轻度路径依赖的无偏PRICE估计，如百慕大期权和具有离散红利的美式看涨期权。通过数值实验验证了精确模拟方法的准确性和速度。数值结果表明，精确模拟方法给出了准确的模拟结果，并且比粗糙的欧拉离散化方法更快。对Heston模型的数值比较表明，我们的方法比Broadie和Kaya的原始模拟方法更具适应性。在Broadie和Kaya[7]的论文中，作者将赫斯顿模型的精确模拟方法扩展到了其他有效跳跃扩散模型。使用与他们相同的应用程序，WMSV模型的准确模拟方法可以扩展到单资产矩阵a f nejump扩散模型[28]，前提是模型中的跳跃具有恒定的强度，并且跳跃大小可以精确采样。这些扩展非常简单，wellin Broadie和Kaya[7]对此进行了解释。附录。1非中心Wishart分布和特殊函数本节旨在回顾非中心Wishart分布和一些特殊函数的定义。有关多元分布和特殊函数的详细讨论，请参阅Muirhead[29]。非中心Wishart分布的概率密度函数包含两个特殊函数：多元伽马函数和矩阵参数的超几何函数。我们从多元伽马函数开始，它是用锥S++自由度对称正定d×d矩阵上的积分定义的。定义A.1。

用Γd（a）表示的多元伽马函数定义为Γd（a）=ZS++dexp(-tr[y]）（det[y]）a-（d+1）/2（dy）、f或Re（a）>（d）- 1).注意，当d=1时，多元伽马函数变成了通常的伽马函数，Γ（a）=Γ（a）。为了定义矩阵参数的超几何函数，我们需要引入分区多项式，它是根据正整数的划分定义的。设k为正整数。k的一个分区ι写成ι=（k，k，···），其中pjkj=kand k≥ K≥ ··· ≥ 0.我们按字典顺序排列分区：设ι=（k，k，···）和ι=（l，l，···）为两个分区，如果kj>ljj，我们写出两部分不相等的第一个索引j。如果ι=（k，···，kd）和ι=（l，···，ld）是两个具有ι>ι的分区，则我们可以得到单项式αk···αkddi，其权重高于单项式αl··αldd。定义A.2。设y是一个d×d复对称矩阵，其特征值为α，··，αd是k的不超过d个部分的划分。对应于ι=（k，···，kd）的y的区域多项式，用Cι（y）表示，是特征值α，··，αd中k次的对称齐次多项式，其中（i）Cι（y）中的最高权重项是αk····················································。（ii）Cι（y）是微分算子的本征函数依偎y=dXj=1αjαj+dXj=1dXi=1i6=jαjαj- αiαj.换句话说，Cι（y）满足对于某些λ，yCι（y）=λCι（y）。（iii）由于ι在k的所有分区上都不同，分区多项式在（tr[y]）k的展开中具有单位系数，即X |ι|=kCι（y）=（tr[y]）k=（α+·+αd）k，其中|ι|表示ι的各部分之和，即|ι|=k+·+kd。定义A.3。

矩阵变元的超几何函数由pfq（a，···，ap；b，··，bq；y）给出=∞Xk=0X |ι|=k（a）ι··（ap）ι（b）ι··（bp）ιCι（y）k！，其中，一个分区的广义超几何系数（a）ι=（k，···，kd）由（a）ι=dYj=1给出A.-（j）- 1)kj和（a）k=a（a+1）·（a+k）- 1） ，（a）=1。这里，函数的参数是复对称d×d矩阵，参数aj，Bj是任意复数。分母参数不允许为零、整数或半整数≤（d）-1） ，否则s中的一些分母将消失。请注意，在d=1的情况下，k只有一个分区不超过d，即（k）。因此，矩阵参数的超几何函数变成了实变量的一般超几何函数。与非中心卡方分布类似，具有整数自由度的非中心Wishart分布通过正态随机向量与自身相乘来定义。考虑δ(∈ N） Rd的独立随机向量，用Z，···，Zδ表示，在奇异协方差矩阵C上具有均值μ，···，μδ和公共N的多变量非线性分布。随机矩阵的分布xw=δXk=1ZkZkis称非中心Wishart分布为自由度δd、协方差矩阵C和非中心参数矩阵Ohm = C-1Pδk=1ukuk、 Ifδ≥ d、 那么W有一个概率密度函数，它的形式是：（det[y]）（δ-D-1） /2dδ/2Γd（δ）（det[C]）δ/2expn-trC-1y+Ohm属于δ;OhmC-1y, Y∈ 当δ是大于d的任何实数时，函数（21）仍然是密度函数-1.这一观察结果使非中心Wishart分布的概念在整数自由度上得以扩展。定义A.4。假设δ>d- 1，C∈ S++d，以及Ohm d×d矩阵是这样的吗Ohm issymmetric正半定义。