自1973年Black和Scholes给出了著名的期权定价公式，金融衍生品的定价问题就成为金融数学的核心研究内容.Black-Scholes(BS)模型将衍生证券的价格表示为标的资产价格的函数.BS模型简单且易于计算，但股价波动率为常数的假设，与实际市场不一致，波动率的微笑曲线验证了这一点.因此，很多学者致力于改进BS模型，随机波动率模型(SV)就是其中一类，即将股价波动率设为另一个随机过程.随机波动率的概念最先由Hull等^[1]提出，之后由Stein等^[2]以及Heston^[3]对其研究工作做了进一步拓展，提出了不同的随机波动率模型.

期权是金融市场中常见的衍生产品，它是一类金融合同，给予了投资者在未来某一时间以约定好的价格买进或者卖出某一标的资产的权利.

随着金融衍生品市场的发展，金融产品种类越来越多样化，结构越来越复杂.随机波动率模型下，大部分金融产品的价格解析解无法求解，需要借助数值方法来求解.在众多的金融衍生品种，期权的计算理论更为基础和根本，很多金融产品的计算研究都以欧式期权为参考^[4].本文主要研究两类随机波动率模型下，基于精确模拟的抽样算法的期权定价及敏感性估计的加速模拟理论.本文以欧式期权的计算为例，文中算法理论可以进一步应用到其他更为复杂的金融产品的计算研究中.

1 蒙特卡罗模拟方法

蒙特卡罗方法以概率论和统计理论为基础，利用随机数来解决实际问题.蒙特卡罗方法在解决实际问题时，先建立一个随机模型，使模型的变量等于所求问题的解，然后通过对变量进行模拟以计算模拟变量的统计特征，最后给出模型的近似解，也就得到了原问题的近似解，其模拟精度可以用估计值的标准误差来衡量.大数定律保证了随着模拟次数的增加，估计值将趋近于正确的值.中心极限定理提供了对于有限次数模拟的误差界限的估计.

设随机变量X，其函数g(X)的数学期望为

$ G = E\left( g \right) = \int {g\left( x \right){\rm{d}}F\left( x \right)} $

(1)

式中:F(x)为随机变量的分布函数.根据随机变量X的分布对其进行抽样模拟，产生X的简单子样，用相应统计量g(x₁), g(x₂), …, g(x_n)的算术平均值

$ {G_n} = \frac{1}{n}\sum\limits_{i = 1}^n {g\left( {{x_i}} \right)} $

(2)

作为G的近似估计值.G_n依概率收敛于所要求的值G，即对∀ε＞0，有

$ \mathop {\lim }\limits_{n \to \infty } P\left( {\left| {{G_n} - G} \right| < \varepsilon } \right) = 1 $

(3)

且有如下近似等式：

$ P\left( {\left| {{G_n} - G} \right| < \frac{\sigma }{{\sqrt n }}{\lambda _\alpha }} \right) \approx \frac{1}{{\sqrt {2{\rm{\pi }}} }}\int_{ - \lambda }^{ + \lambda } {{{\rm{e}}^{ - \frac{{{t^2}}}{2}}}{\rm{d}}t} = 1 - \alpha $

(4)

式中:σ称为模拟量G的标准差；α称为显著水平；λ_α是与显著水平α对应的正态差.

不等式

$ \left| {{G_n} - G} \right| < \frac{\sigma }{{\sqrt n }}{\lambda _\alpha } $

(5)

近似地以概率1-α成立.

如果σ≥0，则Monte Carlo方法的误差ε为

$ \varepsilon = \frac{\sigma }{{\sqrt n }}{\lambda _\alpha } $

(6)

由式(6) 可以看出，σ和n决定了蒙特卡罗方法的误差ε.该误差形式 $\sigma /\sqrt n $ 是蒙特卡罗方法的主要特点.在不改变σ的情况下，要想提高一位数字的精确度，必须要将n提高100倍，极大地增加了计算工作量.然而从另一角度来看，如果使σ减小10倍，就可以减少100倍的工作量.

为提高蒙特卡罗方法的模拟效率，很多学者提出了各种不同的方差减小技术^[5]，如控制变量技术^[6]、重要抽样技术、条件蒙特卡罗、分层抽样算法等，本文主要采用条件蒙特卡罗加速技术的思想.

条件蒙特卡罗主要原理是条件方差公式.计算随机变量V的期望，对一个任意的随机变量H，E[V|H]表示在给定条件H下V的条件期望，这是一个关于H的随机变量，双期望公式为

$ E\left[ V \right] = E\left[ {E\left[ {V\left| H \right.} \right]} \right] $

(7)

条件方差公式

$ {\mathop{\rm var}} \left( V \right) = {\mathop{\rm var}} \left( {E\left[ {V\left| H \right.} \right]} \right) + E\left[ {{\mathop{\rm var}} \left( {V\left| H \right.} \right)} \right] $

(8)

因为E[var(V|H)]≥0, 有

$ {\mathop{\rm var}} \left( V \right) \ge {\mathop{\rm var}} \left( {E\left[ {V\left| H \right.} \right]} \right) $

(9)

可见随机变量E[V|H]的方差比原来随机变量V的方差小，从而在已知精确值E[V|H]的情形下达到了方差减小的目的.

使用条件蒙特卡罗方法可以提高在随机波动率模型下的模拟效率.以波动率为条件对价格进行模拟，使得问题的维度降低，从而提高精度.该方法能很好地适用于对于确定波动率下有解析解的路径依赖的期权的计算.

2 随机波动率模型及精确模拟算法

本节讨论两类随机波动率模型，分别是无跳的随机波动率模型(SV)和带跳的随机波动率模型(SVCJ)及其精确模拟算法. SV模型描述的是连续时间下标的资产价格与瞬时方差率的动态变化过程，但计算机无法直接模拟出连续时间的动态变化，需要先将连续时间离散化为离散时间.欧拉离散方法是普遍使用的离散方法，它可以近似估计出离散时间点股价和方差过程的路径，但欧拉离散常会带来偏差(bias)，有时候达到不能忽略的程度^{[5, 7]}.所以，在研究蒙特卡罗的加速模拟技术的同时，需要研究标的过程的精确模拟算法.

为进一步提高精确抽样的效率，本文在Broadie和Kaya^[7-8]的精确模拟方法的基础上，研究了舍取抽样技术^[9]及其参数的确定问题，讨论了无跳SV和双跳SVCJ两类随机波动率模型的精确模拟算法，考察了几种方差减小技术在欧式期权定价中的应用.

2.1 SV模型下期权的加速模拟算法 2.1.1 SV模型介绍

Heston Model是由Steven Heston于1993年提出的描述随机波动率下标的资产动态过程的数学模型^[3].在风险中性测度下，它假设标的资产S_t与瞬时方差率V_t满足以下随机微分方程模型：

$ {\rm{d}}{S_t} = r{S_t}{\rm{d}}t + \sqrt {{V_t}} {S_t}\left[ {\rho {\rm{d}}W_t^{\left( 1 \right)} + \sqrt {1 - {\rho ^2}} {\rm{d}}W_t^{\left( 2 \right)}} \right] $

(1)

$ {\rm{d}}{V_t} = \kappa \left( {\theta - {V_t}} \right){\rm{d}}t + {\sigma _v}\sqrt {{V_t}} {\rm{d}}W_t^{\left( 1 \right)} $

(2)

方程(1) 刻画了标的资产的动态变化，其中S_t表示t时刻的标的资产价格；r表示风险中性漂移率； $\sqrt {{V_t}} $ 表示波动率；方程(2) 刻画了瞬时方差率V_t的动态变化，其中θ表示长期均方差；κ表示方差回归速度；σ_v是这个方差过程对应的波动率；W_t⁽¹⁾与W_t⁽²⁾是两个独立的布朗运动; ρ表示回报过程与波动率过程间的相关系数.

当时间u＜s＜t，标的资产价格S_u与瞬时方差率V_u已知时，方程(1)、(2) 可采用下列欧拉形式离散:

$ \begin{array}{l} {S_t} = {S_u}\exp \left[ {r\left( {t - u} \right) - \frac{1}{2}\int_u^t {{V_s}{\rm{d}}s} + \rho \int_u^t {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right] \cdot \\ \;\;\;\;\;\;\;\exp \left[ {\sqrt {1 - {\rho ^2}} \int_u^t {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 2 \right)}} } \right] \end{array} $

(3)

$ {V_t} = {V_u} + \kappa \theta \left( {t - u} \right) - \kappa \int_u^t {{V_s}{\rm{d}}s} + {\sigma _v}\int_u^t {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} $

(4)

2.1.2 SV模型与精确模拟算法

Broadie和Kaya提出的精确模拟法^[7-8]是一种无偏差方法，旨在将原随机过程根据它的精确分布模拟抽样，以消除将连续时间离散化带来的偏差.划分时间0=t₀＜…＜t_m=T, 0≤i＜j≤m，考虑时间区间 $\left\lfloor {{t_i},{t_j}} \right\rfloor $ .SV模型下精确模拟的具体步骤为：

步骤1 根据V_ti，从V_tj的分布中生成抽样.

1985年Cox提出，当V_ti(t_i＜t_j)已知时，V_tj符合非中心卡方分布，可表示为

$ \begin{array}{l} {V_{{t_j}}} = \frac{{\sigma _v^2\left( {1 - \exp \left( { - \kappa \left( {{t_j} - {t_i}} \right)} \right)} \right)}}{{4\kappa }} \cdot \\ \;\;\;\;\;\;\;\chi _d^2\left( {\frac{{4\kappa \exp \left( { - \kappa \left( {{t_j} - {t_i}} \right)} \right)}}{{\sigma _v^2\left( {1 - \exp \left( { - \kappa \left( {{t_j} - {t_i}} \right)} \right)} \right)}}{V_u}} \right) \end{array} $

(5)

其中χ′_d²(λ)表示自由度为d，非中心参数为λ的非中心卡方分布随机变量，其抽样过程可参考文献[1-8].

步骤2 根据V_ti和V_tj，从 $\int {_{{t_i}}^{{t_j}}} $ V_sds分布中抽样.

当抽样出V_tj后，可以根据V_ti和V_tj抽样出 $\int {_{{t_i}}^{{t_j}}} $ V_sds，Broadie和Kaya通过拉普拉斯变换得到 $\int {_{{t_i}}^{{t_j}}} $ V_sds的条件特征函数Φ(a)^[7].

$ \begin{array}{l} \mathit{\Phi }\left( a \right) = E\left[ {\exp \left( {{\rm{i}}a\int_{{t_i}}^{{t_j}} {{V_s}{\rm{d}}s} } \right)\left| {{V_{{t_i}}},{V_{{t_j}}}} \right.} \right] = \frac{{\gamma \left( a \right)\exp \left( { - 0.5\left( {\gamma \left( a \right) - \kappa } \right)\left( {{t_j} - {t_i}} \right)} \right)\left( {1 - \exp \left( { - \kappa \left( {{t_j} - {t_i}} \right)} \right)} \right)}}{{\kappa \left( {1 - \exp \left( { - \gamma \left( a \right)\left( {{t_j} - {t_i}} \right)} \right)} \right)}} \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\exp \left\{ {\frac{{{V_{{t_i}}} + {V_{{t_j}}}}}{{\sigma _v^2}}\left[ {\frac{{\kappa \left( {1 + \exp \left( { - \kappa \left( {{t_j} - {t_i}} \right)} \right)} \right)}}{{1 - \exp \left( { - \kappa \left( {{t_j} - {t_i}} \right)} \right)}} - \frac{{\gamma \left( a \right)\left( {1 + \exp \left( { - \gamma \left( a \right)\left( {{t_j} - {t_i}} \right)} \right)} \right)}}{{1 - \exp \left( { - \gamma \left( a \right)\left( {{t_j} - {t_i}} \right)} \right)}}} \right]} \right\} \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\frac{{{{\rm{I}}_{0.5d - 1}}\left[ {\sqrt {{V_{{t_i}}}{V_{{t_j}}}} \frac{{4\gamma \left( a \right)\exp \left( { - 0.5\gamma \left( a \right)\left( {{t_j} - {t_i}} \right)} \right)}}{{\sigma _v^2\left( {1 - \exp \left( { - \gamma \left( a \right)\left( {{t_j} - {t_i}} \right)} \right)} \right)}}} \right]}}{{{{\rm{I}}_{0.5d - 1}}\left[ {\sqrt {{V_u}{V_t}} \frac{{4\kappa \exp \left( { - 0.5\kappa \left( {t - u} \right)} \right)}}{{\sigma _v^2\left( {1 - \exp \left( { - \kappa \left( {t - u} \right)} \right)} \right)}}} \right]}} \end{array} $

其中γ(a)= $\sqrt {{\kappa ^2} - 2{\sigma _v}{\rm{i}}{a}} $ ，I_v(x)是第一类修正贝塞尔函数.

假设在已知V_ti和V_tj的情况下，随机变量V(t_i, t_j)和 $\int {_{{t_i}}^{{t_j}}} $ V_sds同分布，那么利用傅里叶展开可以由条件特征函数Φ(a)求解出V(t_i, t_j)的累积分布函数F(x)和密度函数f(x)^[7].

$ \begin{array}{l} F\left( x \right) \equiv P\left( {V\left( {u,t} \right) \le x} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;\;\frac{2}{{\rm{\pi }}}\int_0^{ + \infty } {\frac{{\sin \left( {ux} \right)}}{u}{\mathop{\rm Re}\nolimits} \left[ {\mathit{\Phi }\left( a \right)} \right]{\rm{d}}a} \end{array} $

$ f\left( x \right) \equiv \frac{2}{{\rm{\pi }}}\int_0^{ + \infty } {\cos \left( {ux} \right){\mathop{\rm Re}\nolimits} \left[ {\mathit{\Phi }\left( a \right)} \right]{\rm{d}}a} $

分布函数F(x)和密度函数f(x)形式复杂，采取直接抽样技术不容易操作，本文讨论舍取抽样技术^[9]及该技术中参数的选取问题，以提高抽样效率.

注意密度函数f(x)是关于0到无穷区域的积分，在Matlab中，它无法直接求解，可以利用梯形公式化为有限区间上的积分.为了有效使用梯形公式，先使用Matlab绘制出被积函数p(x)≡ $\frac{2}{{\rm{\pi }}}$ ·cos(ux)Re[Φ(a)]与变量a的图像，其中S₀=100;r=0.05;V₀=0.015;κ=2;θ=0.01;σ=0.05;ρ=－0.75，见图 1.

由图 1可以看出，选取的积分区间与T有关.此外κ，V₀等初始量也会影响积分区间的选择.因此对于不同的初始值要考虑选取合适的积分区域.在计算机模拟时，可以利用梯形公式计算有限区间上这个积分的值近似得到密度函数，但需要特别注意对积分区间的划分.

由梯形公式得到 $\int {_{{t_i}}^{{t_j}}} $ V_sds的密度函数f(x)后，使用舍取抽样技术抽取符合 $\int {_{{t_i}}^{{t_j}}} $ V_sds分布的样本点，步骤如下：

(1) 选择合适的密度函数g(x)，使对于所有常数c＞1，f(x)＜cg(x)成立.

(2) 从g(x)中抽样出x，并从均匀分布U~(0, 1) 中抽样出u.

(3) 检查u＜ $\frac{{f\left( x \right)}}{{cg\left( x \right)}}$ 是否成立.如果不等式成立，接受x作为 $\int {_{{t_i}}^{{t_j}}} $ V_sds的抽样；如果不等式不成立，拒绝x并重复以上过程.

为选取合适的g(x)，先绘制出 $\int {_{{t_i}}^{{t_j}}} $ V_sds密度函数图像，见图 2，其中S₀=100;T=1.00;r=0.05;V₀=0.015;κ=2;θ=0.01;σ=0.05;ρ=－0.75.

由图 2可知, $\int {_{{t_i}}^{{t_j}}} $ V_sds的密度函数f(x)的图像接近Gamma分布，因此可以考虑选取g(x)是符合Gamma分布的密度函数^[9]，对应参数可以由 $\int {_{{t_i}}^{{t_j}}} $ V_sds的特征函数Φ(a)求解，即为提高舍取抽样效率，匹配f(x)与g(x)的一阶矩与二阶矩得到Gamma分布的参数如下：

$ g\left( x \right) = \frac{{{\beta ^\alpha }}}{{\mathit{\Gamma }\left( \alpha \right)}}{x^{\alpha - 1}}\exp \left( { - \beta x} \right), $

$ E\left[ X \right] = \frac{{\mathit{\Phi '}\left( 0 \right)}}{i}, $

$ E\left[ {{X^2}} \right] = - \mathit{\Phi ''}\left( 0 \right), $

$ \alpha = - \frac{{{{\left[ {\mathit{\Phi '}\left( 0 \right)} \right]}^2}}}{{{{\left[ {\mathit{\Phi '}\left( 0 \right)} \right]}^2} - \mathit{\Phi ''}\left( 0 \right)}}, $

$ \beta = \frac{{\mathit{\Phi '}\left( 0 \right)}}{{\left[ {{{\left[ {\mathit{\Phi '}\left( 0 \right)} \right]}^2} - \mathit{\Phi ''}\left( 0 \right)} \right]i}} $

选取c=1.1时，f(x), cg(x)的图像见图 3，其中S₀=100;T=1.00;r=0.05;V₀=0.015;κ=2;θ=0.01;σ=0.05;ρ=－0.75.

所以c=1.1满足舍取抽样法的条件，可以用它抽样出对应的 $\int {_{{t_i}}^{{t_j}}} $ V_sds的样本点.

步骤3 根据V_ti, V_tj， $\int {_{{t_i}}^{{t_j}}} $ V_sds计算 $\int {_{{t_i}}^{{t_j}}} \sqrt {{V_s}} $ dW_s⁽¹⁾的值.

当有V_ti，V_tj和 $\int {_{{t_i}}^{{t_j}}} $ V_sds的样本后，可以从以下等式抽样出 $\int {_{{t_i}}^{{t_j}}} \sqrt {{V_s}} $ dW_s⁽¹⁾：

$ \begin{array}{l} \int_{{t_i}}^{{t_j}} {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} = \frac{1}{{{\sigma _v}}}\left( {{V_{{t_j}}} - {V_{{t_i}}} - \kappa \theta \left( {{t_j} - {t_i}} \right) + } \right.\\ \;\;\;\;\;\;\;\left. {\kappa \int_{{t_i}}^{{t_j}} {{V_s}{\rm{d}}s} } \right) \end{array} $

步骤4 以 $\int {_{{t_i}}^{{t_j}}} $ V_sds， $\int {_{{t_i}}^{{t_j}}} \sqrt {{V_s}} $ dW_s⁽¹⁾为条件求出S_ti.

定义 $\left\lfloor {{t_i},{t_j}} \right\rfloor $ 间的平均方差为

$ \bar \sigma _j^2 = \frac{{\left( {1 - {\rho ^2}} \right)\int_{{t_i}}^{{t_j}} {{V_s}{\rm{d}}s} }}{{{t_j} - {t_i}}} $

定义辅助变量

$ {\xi _j} = \exp \left( { - \frac{{{\rho ^2}}}{2}\int_{{t_i}}^{{t_j}} {{V_s}{\rm{d}}s} + \rho \int_{{t_i}}^{{t_j}} {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right) $

已知方差路径与t_i时标的资产价格S_ti，可以得到这时S_tj的表达式为

$ {S_{{t_j}}} = {S_{{t_i}}}{\xi _j}\exp \left[ {\left( {r - \frac{{\bar \sigma _j^2}}{2}} \right)\left( {{t_j} - {t_i}} \right) + {{\bar \sigma }_j}\sqrt {{t_j} - {t_i}} Z} \right] $

由此就可以模拟出SV模型下标的资产过程.

2.1.3 SV模型下欧式期权定价的条件蒙特卡罗算法

考虑同样到期日为T，执行价格为K，初始标的资产价格为S₀，波动率为常数σ的欧式看涨期权，它的价格可以用Black-Scholes公式表示，记为B_S(S₀, σ).

在SV模型精确模拟方法的前3步可以抽样出 $\int \begin{array}{l} T\\ 0 \end{array} $ V_sds， $\int \begin{array}{l} T\\ 0 \end{array} \sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}$ ，Wiilard^[10]提出以它们为条件时，S_T符合对数正态分布.

定义[0, T]间的平均波动率为

$ \bar \sigma = \sqrt {\frac{{\left( {1 - {\rho ^2}} \right)\int_0^T {{V_s}{\rm{d}}s} }}{T}} $

定义辅助变量

$ \xi = \exp \left( { - \frac{{{\rho ^2}}}{2}\int_0^T {{V_s}{\rm{d}}s} + \rho \int_0^T {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right) $

这时S_T满足

$ {S_T} = {S_0}\xi \exp \left[ {\left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T + \bar \sigma \sqrt T Z} \right] $

那么在 $\int \begin{array}{l} T\\ 0 \end{array} $ V_sds和 $\int \begin{array}{l} T\\ 0 \end{array} \sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}$ 为条件时，看涨期权价格可以写为B_S(S₀ξ, σ).

利用条件蒙特卡罗可知SV模型下欧式期权的价格为

$ \begin{array}{l} C = E\left[ {\exp \left( { - rT} \right){{\left( {S\left( T \right) - K} \right)}^ + }} \right] = \\ \;\;\;\;\;\;E\left[ {E\left[ {\exp \left( { - rT} \right){{\left( {S\left( T \right) - K} \right)}^ + }\left| {\int_0^T {{V_s}{\rm{d}}s} } \right.,} \right.} \right.\\ \;\;\;\;\;\;\left. {\left. {\int_0^T {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right]} \right] = \\ \;\;\;\;\;\;E\left[ {{S_0}\xi N\left( {{d_1}} \right) - \exp \left( { - rT} \right)KN\left( {{d_2}} \right);\bar \sigma } \right] = \\ \;\;\;\;\;\;E\left[ {{B_S}\left( {{S_0}\xi ,\bar \sigma } \right);\bar \sigma } \right] \end{array} $

其中：S(T)表示T时刻股价；K表示敲定价；

$ {d_{1,2}} = \frac{{\ln \frac{{{S_0}\xi }}{K} + \left( {r \pm \frac{1}{2}{{\bar \sigma }^2}} \right)T}}{{\bar \sigma \sqrt T }}. $

2.2 SVCJ模型下期权的加速模拟算法 2.2.1 SVCJ模型的精确模拟算法

SVCJ模型是在Heston SV模型的基础上增加了方差的跳跃过程.它假设标的资产S_t与瞬时方差率V_t满足以下动态随机模型：

$ \begin{array}{l} {\rm{d}}{S_t} = \left( {r - \lambda \bar \mu } \right){S_t}{\rm{d}}t + \sqrt {{V_t}} {S_t}\left[ {\rho {\rm{d}}W_t^{\left( 1 \right)} + } \right.\\ \;\;\;\;\;\;\;\;\;\left. {\sqrt {1 - {\rho ^2}} {\rm{d}}W_t^{\left( 2 \right)}} \right] + {S_t}\left( {{J_s} - 1} \right){\rm{d}}{N_t} \end{array} $

(6)

$ {\rm{d}}{V_t} = \kappa \left( {\theta - {V_t}} \right){\rm{d}}t + {\sigma _v}\sqrt {{V_t}} {\rm{d}}W_t^{\left( 1 \right)} + {J_v}{\rm{d}}{N_t} $

(7)

方程(6) 刻画了标的资产的动态变化，其中S_t表示t时刻的标的资产价格；r表示风险中性漂移率； $\sqrt {{V_t}}$ 表示波动率；方程7刻画了瞬时方差率V_t的动态变化，其中θ表示长期均方差，κ表示方差回归速度；σ_v是这个方差过程对应的波动率；W_t⁽¹⁾与W_t⁽²⁾是两个独立的布朗运动过程；ρ表示回报过程与波动率过程间的相关系数；N_t是参数为λ_J的泊松过程，它与布朗运动独立；J_s表示标的资产价格相对跳跃幅度；J_v表示方差的跳跃幅度.特别地，标的资产价格与方差的跳跃是相关的，即它们是同时发生跳跃的，如果在时刻t发生跳跃，那么S_t+=S_t－J_s，V_t+=V_t－+J_v，它们之间有相关系数ρ_J.跳跃幅度J_v符合期望为μ_v的指数分布.当J_v已知时，J_s符合对数正态分布，ln J_s~N(μ_s+ρ_JJ_v, σ_s²)，且参数μ_s和μ相关，满足表达式

$ {\mu _s} = \ln \left[ {\left( {1 + \bar \mu } \right)\left( {1 - {\rho _J}{\mu _v}} \right)} \right] - \frac{1}{2}\sigma _s^2 $

划分时间0=t₀＜…＜t_m=T, 0≤i＜j≤m，考虑时间区间 $\left\lfloor {{t_i},{t_j}} \right\rfloor $ ，具体步骤如下：

步骤1 模拟出区间 $\left\lfloor {{t_i},{t_j}} \right\rfloor $ 间的方差跳跃步数n_j，更新时间为t_i＜t_{jump_1}＜…＜t_{jump_nj}＜t_j.

步骤2 根据V_ti，从V_tj的分布中生成样本V_{tjump_1}, …, V_{tjump_nj}, V_tj.

步骤3 根据V_ti和V_{tjump_1}, …, V_{tjump_nj}, V_tj，生成 $\int \begin{array}{l} {t_{{\rm{jump\_1}}}}\\ {t_i} \end{array} $ V_sds, $\int \begin{array}{l} {t_{{\rm{jump\_2}}}}\\ {t_{{\rm{jump\_1}}}} \end{array} $ V_sds, …, $\int \begin{array}{l} {{t}_{i}}\\ {t_{{\rm{jump\_}}{{n}_j}}} \end{array} $ V_sds.

步骤4 模拟出每次跳跃大小J_{v, jump_1}, …, J_{v, jump_nj}，更新 $\tilde V$ _{tjump_i}=V_{tjump_i}+J_{v, tjump_i}, 1≤i≤n_j.

步骤5 根据V_ti，V_tj， $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} $ V_sds计算 $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} \sqrt {{V_s}} {\rm{d}}{W}_s^{\left( 1 \right)}$ 的值.

步骤6 以 $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} $ V_sds， $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} \sqrt {{V_s}} {\rm{d}}{W}_s^{\left( 1 \right)}$ ，跳跃步数n_j，跳跃大小J_{v, jump_1}, …, J_{v, jump_nj}为条件求出S_ti.

前5步模拟和SV模型基本相同，下面详细解释第6步.

已知区间 $\left\lfloor {{t_i},{t_j}} \right\rfloor $ 的 $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} $ V_sds， $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} \sqrt {{V_s}} {\rm{d}}{W}_s^{\left( 1 \right)}$ ，方差跳跃步数n_j，对应跳跃大小J_{v, jump_1}, …, J_{v, jump_nj}后，可抽样出S_ti.

定义 $\left\lfloor {{t_i},{t_j}} \right\rfloor $ 间的平均方差为

$ \bar \sigma _j^2 = \frac{{{n_j}\sigma _s^2 + \left( {1 - {\rho ^2}} \right)\int_{{t_i}}^{{t_j}} {{V_s}{\rm{d}}s} }}{{{t_j} - {t_i}}} $

定义辅助变量

$ \begin{array}{l} {\xi _j} = \exp \left( {\sum\limits_{k = 1}^{{n_j}} {\left( {{\mu _s} + {J_{v,k}}{\rho _J} + \frac{{\sigma _s^2}}{2}} \right) - \lambda \bar \mu \left( {{t_j} - {t_i}} \right) - } } \right.\\ \;\;\;\;\;\;\;\left. {\frac{{{\rho ^2}}}{2}\int_{{t_i}}^{{t_j}} {{V_s}{\rm{d}}s} + \rho \int_{{t_i}}^{{t_j}} {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right) \end{array} $

已知方差路径与t_i时刻标的资产价格S_ti，可以得到t_j时刻标的资产价格S_tj的表达式

$ {S_{{t_j}}} = {S_{{t_i}}}{\xi _j}\exp \left[ {\left( {r - \frac{{\bar \sigma _j^2}}{2}} \right)\left( {{t_j} - {t_i}} \right) + {{\bar \sigma }_j}\sqrt {{t_j} - {t_i}} Z} \right] $

由此就可以模拟出SVCJ模型的标的资产价格路径.

2.2.2 SVCJ模型下期权定价的条件蒙特卡罗算法

在SVCJ模型精确模拟方法中可以抽样出 $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} $ V_sds， $\int \begin{array}{l} {t_j}\\ {t_i} \end{array} \sqrt {{V_s}} {\rm{d}}{W}_s^{\left( 1 \right)}$ ，同SV中的推导，利用条件蒙特卡罗可知SVCJ模型下欧式期权价格为

$ \begin{array}{l} C = E\left[ {\exp \left( { - rT} \right){{\left( {S\left( T \right) - K} \right)}^ + }} \right] = \\ \;\;\;\;\;\;E\left[ {E\left[ {\exp \left( { - rT} \right){{\left( {S\left( T \right) - K} \right)}^ + }\left| {\int_0^T {{V_s}{\rm{d}}s} } \right.,} \right.} \right.\\ \;\;\;\;\;\;\left. {\left. {\int_0^T {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right]} \right] = \\ \;\;\;\;\;\;E\left[ {\sum\limits_{n = 0}^\infty {\frac{{{{\left( {{{\lambda '}_J}T} \right)}^n}}}{{n!}}\exp \left( { - {{\lambda '}_J}T} \right)} \left[ {{S_0}\xi N\left( {{d_{{n_1}}}} \right) - } \right.} \right.\\ \;\;\;\;\;\;\left. {\left. {K\exp \left( { - rT} \right)N\left( {{d_{{n_2}}}} \right)} \right];\bar \sigma } \right] = \\ \;\;\;\;\;\;E\left[ {\sum\limits_{n = 0}^\infty {\frac{{{{\left( {{{\lambda '}_J}T} \right)}^n}}}{{n!}}\exp \left( { - {{\lambda '}_J}T} \right){B_S}\left( {{S_0}\xi ,{\sigma _n},{r_n}} \right);\bar \sigma } } \right] \end{array} $

其中：

$ {d_{{n_1},{n_2}}} = \frac{{\ln \frac{{{S_0}\xi }}{K} + \left( {{r_n} \pm \frac{1}{2}\sigma _n^2} \right)T}}{{{\sigma _n}\sqrt T }}, $

$ \sigma _n^2 = {{\bar \sigma }^2} + \frac{{n\sigma _s^2}}{T}, $

$ {{\lambda '}_J} = {\lambda _J}\left( {1 + \bar \mu } \right), $

$ {r_n} = r - {\lambda _J}\bar \mu + \frac{n}{T}\left( {{\mu _s} + J_k^v{\rho _J} + \frac{{\sigma _s^2}}{2}} \right) $

此时定义[0, T]间的平均波动率为

$ \bar \sigma = \sqrt {\frac{{\left( {1 - {\rho ^2}} \right)\int_0^T {{V_s}{\rm{d}}s} }}{T}} $

定义辅助变量

$ \xi = \exp \left( { - \frac{{{\rho ^2}}}{2}\int_0^T {{V_s}{\rm{d}}s} + \rho \int_0^T {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right) $

3 期权定价的数值模拟结果

下面将分别在精确模拟算法基础上考察原始蒙特卡罗法、对偶变量法和条件蒙特卡罗方法对欧式看涨期权的价格估计的加速效果.同时，SV、SVCJ模型中涉及到初始参数有K、T、r、V₀、κ、σ、θ、ρ、σ_s、μ、λ_J、μ_v、ρ_J，不同的参数也会影响蒙特卡罗法下期权估计值和误差，以下也将参数对期权价格估计和误差估计的影响做出分析.

表 1为SV模型下不同模拟路径数m对蒙特卡罗法、对偶变量法和条件蒙特卡罗中欧式看涨期权价格误差估计效果.其中，S₀=100，K=100，T=1.00，r=0.05，V₀=0.015，σ=0.15，θ=0.02，κ=4，ρ=－0.25.

由表 1可知，随着模拟路径数m的增加，3种方法的价格估计值都将趋于收敛；条件蒙特卡罗方法较其他两种方法有方差减小的效果.

本文继续考察了误差减小效果随相关系数ρ，波动率σ，方差回归速度κ，长期均方差θ和执行时间T等参数变化的影响，见表 2~4，表中R1表示对偶变量法相对原始蒙特卡罗法误差减小倍数; R2表示条件蒙特卡罗相对原始蒙特卡罗法误差减小倍数；R3表示条件蒙特卡罗法相对对偶变量法误差减小倍数.

表 2中，S₀=100, K=100, T=1.00, r=0.05, V₀=0.015, κ=4, σ=0.15, θ=0.02, m=10, 000.由表 2可知，随着相关系数的增加，3种方法的价格估计值逐渐减小；条件蒙特卡罗法较其他两种方法具有较小的方差，达到了方差减小的效果，特别当相关系数趋于0时，条件蒙特卡罗法的方差减小的效果非常明显.

表 3中，S₀=100，K=100，T=1，r=0.05，V₀=0.015，κ=4，θ=0.02，ρ=－0.25，m=10 000.由表 3可知，随着波动率的增加，蒙特卡罗法，对偶变量法，条件蒙特卡罗的价格估计值逐渐增加，对应方差呈逐渐减小趋势; 条件蒙特卡罗的方差减小倍数逐渐增加，较其他两种方法具有较小的方差，有更好的方差减小效果.

表 4中，S₀=100, K=100, T=1.00, r=0.05, V₀=0.015, σ=0.15, θ=0.02, ρ=－0.25, m=10 000.由表 4可知，随着方差回归速度的增加，3种方法的价格估计值逐渐增加，对应方差也呈逐渐增加趋势; 条件蒙特卡罗法的方差减小倍数逐渐减少，其方差始终低于其他两种方法，说明较其他两种方法，条件蒙特卡罗法具有更好的方差减小效果.

4 期权敏感度分析

期权敏感性指期权价格相对于定价参数变动的敏感程度，敏感性分析可以帮助投资者选取适当的期权组合，减小风险.蒙特卡罗法中求解Greeks的方法可总结为3类：有限差分近似、顺向微分估计(pathwise method, PW)和似然比估计(likelihood ration method, LR).其中有限差分近似虽然容易操作，但它会为估计值带来偏差^[5]，本文将使用PW方法和LR方法，考察欧式期权价格随参数变化的影响.PW是基于对收益函数的微分，LR则是基于对标的资产价格密度函数的微分，它们都能求解出Greeks的无偏估计量，LR的应用性更广泛，PW的结果更精确.下面讨论用LR与PW方法计算Δ、γ和ρ.Δ是期权价格关于标的资产初始价格S₀的导数，它表示期权价格与标的资产初始价格之间的关系；γ是期权价格关于标的资产初始价格S₀的二阶导数，它表示Δ与标的资产初始价格之间的变化关系；ρ是期权价格关于无风险利率r的导数，它表示期权价格与无风险利率之间的变化关系.

考虑执行价格K，到期日T的欧式看涨期权，它的价格应为

$ \begin{array}{*{20}{c}} {C = E\left[ {\exp \left( { - rT} \right){{\left( {{S_T} - K} \right)}^ + }} \right],记}\\ {{C_i} = \exp \left( { - rT} \right){{\left( {{S_T} - K} \right)}^ + }} \end{array} $

在精确模拟方法中抽样出 $\int \begin{array}{l} T\\ 0 \end{array} {V_s}{\rm{d}}{s,}\int \begin{array}{l} T\\ 0 \end{array} \sqrt {{V_s}} {\rm{d}}W_s^{(1)}$ 的条件下，标的资产价格可以视为常数波动率下动态过程，可表示为如下形式：

$ {S_T} = {S_0}\xi \exp \left[ {\left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T + \bar \sigma \sqrt T Z} \right],Z \sim N\left( {0,1} \right) $

σ, ξ的具体表达式可参考第3节，与S₀, r无关.

4.1 PW方法

PW方法假设微分与期望具有可交换性，在此条件下可以推导出欧式看涨期权Greeks(Δ，ρ)的解析解。

4.1.1 Δ

已知Δ是期权价格关于标的资产初始价格S₀的导数，利用链式法则将导数写为

$ \frac{{{\rm{d}}{C_i}}}{{{\rm{d}}{S_0}}} = \frac{{{\rm{d}}{C_i}}}{{{\rm{d}}{S_T}}}\frac{{{\rm{d}}{S_T}}}{{{\rm{d}}{S_0}}} $

其中：

$ \frac{{{\rm{d}}{C_i}}}{{{\rm{d}}{S_T}}} = \left\{ \begin{array}{l} \exp \left( { - rT} \right),\;\;\;\;{S_T} > K\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{S_T} \le K \end{array} \right. $

$ \begin{array}{l} \frac{{{\rm{d}}{S_T}}}{{{\rm{d}}{S_0}}} = \xi \exp \left[ {\left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T + \bar \sigma \sqrt T Z} \right] = \frac{{{S_T}}}{{{S_0}}} \Rightarrow \\ \;\;\;\;\;\;\;\Delta :\frac{{{\rm{d}}C}}{{{\rm{d}}{S_0}}} = E\left[ {\frac{{{\rm{d}}{C_i}}}{{{\rm{d}}{S_0}}}} \right] = E\left[ {\exp \left( { - rT} \right)\left\{ {{S_T} > K} \right\}\frac{{{S_T}}}{{{S_0}}}} \right] \end{array} $

4.1.2 ρ

ρ是期权价格关于无风险利率r的导数，它表示期权价格与无风险利率之间的变化关系，可写为

$ \frac{{{\rm{d}}{C_i}}}{{{\rm{d}}r}} = \left\{ \begin{array}{l} \frac{{{\rm{d}}\left[ {\exp \left( { - rT} \right)\left( {{S_T} - K} \right)} \right]}}{{{\rm{d}}r}},\;\;\;\;{S_T} > K\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{S_T} \le K \end{array} \right. $

其中：

$ \begin{array}{l} \frac{{{\rm{d}}\left[ {\exp \left( { - rT} \right)\left( {{S_T} - K} \right)} \right]}}{{{\rm{d}}r}} = - T\exp \left( { - rT} \right){S_T} + \\ \;\;\;\;\;\;\exp \left( { - rT} \right)\frac{{{\rm{d}}{S_T}}}{{{\rm{d}}r}} + KT\exp \left( { - rT} \right) \end{array} $

$ \begin{array}{l} \frac{{{\rm{d}}{S_T}}}{{{\rm{d}}r}} = T{S_0}\xi \exp \left[ {\left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T + \bar \sigma \sqrt T Z} \right] = T{S_T} \Rightarrow \\ \;\;\;\;\;\;\rho :\frac{{{\rm{d}}C}}{{{\rm{d}}r}} = E\left[ {\frac{{{\rm{d}}{C_i}}}{{{\rm{d}}r}}} \right] = E\left[ {\exp \left( { - rT} \right)\left\{ {{S_T} > K} \right\}KT} \right] \end{array} $

故使用PW方法得一条模拟路径中欧式看涨期权Greeks的估计值为

$ \begin{array}{l} \Delta :\exp \left( { - rT} \right)\left[ {{S_T} \ge K} \right]\frac{{{S_T}}}{{{S_0}}}\\ \rho :\exp \left( { - rT} \right)\left[ {{S_T} \ge K} \right]KT \end{array} $

4.2 LR方法

LR方法是基于条件密度函数得到Greeks的估计值，假设S_T的条件密度函数为g_θ(x)，其中θ为概率密度函数的参数，那么欧式看涨期权价格C=E[exp(-rT)(S_T－K)⁺]可以写为

$ C = \int {\exp \left( { - rT} \right){{\left( {x - K} \right)}^ + }{g_\theta }\left( x \right){\rm{d}}x} $

如果微分与积分可交换次序，那么期权价格C关于参数θ的导数可以表示为

$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}C}}{{{\rm{d}}\theta }} = \int {\exp \left( { - rT} \right){{\left( {x - K} \right)}^ + }\frac{{{\rm{d}}{g_\theta }\left( x \right)}}{{{\rm{d}}\theta }}{\rm{d}}x} = }\\ {\int {\exp \left( { - rT} \right){{\left( {x - K} \right)}^ + }\frac{{{{g'}_\theta }\left( x \right)}}{{{g_\theta }\left( x \right)}}{g_\theta }\left( x \right){\rm{d}}x = } }\\ {E\left[ {\exp \left( { - rT} \right){{\left( {{S_T} - K} \right)}^ + }\frac{{{{g'}_\theta }\left( {{S_T}} \right)}}{{{g_\theta }\left( {{S_T}} \right)}}} \right]} \end{array} $

其中exp(-rT)(S_T－K)⁺ $\frac{{g{{'}_\theta }\left( {{S_T}} \right)}}{{{g_\theta }\left( {{S_T}} \right)}}$ 是 $\frac{{{\rm{d}}{C}}}{{{\rm{d}}\theta }}$ 的无偏估计量，称 $\frac{{g{{'}_\theta }\left( {{S_T}} \right)}}{{{g_\theta }\left( {{S_T}} \right)}}$ 为评分函数(score function)，它独立于收益函数exp(-rT)(S_T－K)⁺.

在精确模拟方法中抽样出 $\smallint \begin{array}{*{20}{l}} T\\ 0 \end{array}{V_s}{\rm{d}}s,\smallint \begin{array}{*{20}{l}} T\\ 0 \end{array}\sqrt {{V_s}} \cdot {\rm{d}}W_s^{(1)}$ 后，标的资产价格可以视为常数波动率下动态过程，可表示为如下形式：

$ {S_T} = {S_0}\xi \exp \left[ {\left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T + \bar \sigma \sqrt T Z} \right],Z \sim N\left( {0,1} \right) $

此时S_T符合对数正态分布，它的累积分布函数可以写为

$ \begin{array}{l} P\left( {{S_T} \le x} \right) = P\left( {{S_0}\xi \exp \left[ {\left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T + \bar \sigma \sqrt T Z} \right] \le x} \right) = \\ P\left( {\left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T + \bar \sigma \sqrt T Z \le \ln \left( {\frac{x}{{{S_0}\xi }}} \right)} \right) = \\ P\left( {Z \le \frac{{\ln \left( {\frac{x}{{{S_0}\xi }}} \right) - \left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T}}{{\bar \sigma \sqrt T }}} \right) = \\ \mathit{\Phi = }\left( {\frac{{\ln \left( {\frac{x}{{{S_0}\xi }}} \right) - \left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T}}{{\bar \sigma \sqrt T }}} \right) \end{array} $

对应密度函数为

$ \frac{1}{{x\bar \sigma \sqrt T }}\varphi \left( {\frac{{\ln \left( {\frac{x}{{{S_0}\xi }}} \right) - \left( {r - \frac{{{{\bar \sigma }^2}}}{2}} \right)T}}{{\bar \sigma \sqrt T }}} \right) $

该密度函数可以简写为g(x)= $\frac{1}{{x\bar \sigma \sqrt T }}$ ·φ(d(x))，φ(·)是标准正态密度函数，d(x)的表达式为

$ d\left( x \right) = \frac{{\ln \left( {x/\left( {{S_0}\xi } \right)} \right) - \left( {r - \frac{1}{2}{{\bar \sigma }^2}} \right)T}}{{\bar \sigma \sqrt T }} $

在不同参数下，密度函数的微分化简计算出的评分函数为

$ \begin{array}{l} \frac{{\frac{{\partial g\left( x \right)}}{{\partial {S_0}}}}}{{g\left( x \right)}} = \frac{{\varphi '\left( {d\left( x \right)} \right)\frac{{\partial d\left( x \right)}}{{\partial {S_0}}}}}{{\varphi \left( {d\left( x \right)} \right)}} = \frac{{d\left( x \right)}}{{{S_0}\bar \sigma \sqrt T }}\\ \frac{{\frac{{{\partial ^2}g\left( x \right)}}{{\partial {S_0}}}}}{{g\left( x \right)}} = \frac{{\partial \frac{{\left[ {\varphi \left( {d\left( x \right)} \right)\frac{{d\left( x \right)}}{{{S_0}}}} \right]}}{{\partial {S_0}}}}}{{\varphi \left( {d\left( x \right)} \right)}} = \frac{{{d^2} - d\bar \sigma \sqrt T - 1}}{{{S_0}{{\bar \sigma }^2}T}}\\ \frac{{\frac{{\partial g\left( x \right)}}{{\partial r}}}}{{g\left( x \right)}} = \frac{{\varphi '\left( {d\left( x \right)} \right)\frac{{\partial d\left( x \right)}}{{\partial r}}}}{{\varphi \left( {d\left( x \right)} \right)}} = - \frac{{d\left( x \right)T}}{{\bar \sigma \sqrt T }} \end{array} $

对应敏感性指数LR估计值为

$ \Delta :\exp \left( { - rT} \right){\left( {{S_T} - K} \right)^ + }\left( {\frac{d}{{{S_0}\bar \sigma \sqrt T }}} \right) $

$ \mathit{\Gamma }:\exp \left( { - rT} \right){\left( {{S_T} - K} \right)^ + }\left( {\frac{{{d^2} - d\bar \sigma \sqrt T - 1}}{{{S_0}{{\bar \sigma }^2}T}}} \right) $

$ \rho :\exp \left( { - rT} \right){\left( {{S_T} - K} \right)^ + }\left( { - T + \frac{{d\sqrt T }}{{\bar \sigma }}} \right) $

PW方法无法直接计算出Gamma的值，可以结合LR方法求出欧式看涨期权Gamma的混合估计值，LRPW表示先使用PW方法估计一阶导数，再使用LR方法估计二阶导数，PWLR意义类似.

LR-PW

$ \begin{array}{l} {\rm{LR}} - {\rm{PW}}\\ \;\;\;\;\;\mathit{\Gamma }:\exp \left( { - rT} \right)\left[ {{S_T} \ge K} \right]K\left( {\frac{d}{{S_0^2\bar \sigma \sqrt T }}} \right) \end{array} $

PW-LR

$ \begin{array}{l} {\rm{PW}} - {\rm{LR}}\\ \;\;\;\;\;\mathit{\Gamma }:\exp \left( { - rT} \right)\left[ {{S_T} \ge K} \right]\frac{{{S_T}}}{{S_0^2}}\left( {\frac{d}{{\bar \sigma \sqrt T }} - 1} \right) \end{array} $

5 条件蒙特卡罗方法计算Greeks

本节使用蒙特卡罗加速技术来计算Greeks的值.在SV模型下根据条件蒙特卡罗的思想，欧式看涨期权的解析解为

$ \begin{array}{l} C = E\left[ {{B_{\rm{S}}}\left( {{S_0}\xi ,\bar \sigma } \right);\bar \sigma } \right] = E\left[ {{S_0}\xi N\left( {{d_1}} \right) - } \right.\\ \;\;\;\;\;\;\left. {\exp \left( { - rT} \right)KN\left( {{d_2}} \right);\bar \sigma } \right] \end{array} $

其中：

$ {d_{1,2}} = - \frac{{\ln \frac{{{S_0}\xi }}{K} + \left( {r \pm \frac{1}{2}{{\bar \sigma }^2}} \right)T}}{{\bar \sigma \sqrt T }}, $

$ \bar \sigma = \sqrt {\frac{{\left( {1 - {\rho ^2}} \right)\int_0^T {{V_s}{\rm{d}}s} }}{T}} , $

$ \xi = \exp \left( { - \frac{{{\rho ^2}}}{2}\int_0^T {{V_s}{\rm{d}}s} + \rho \int_0^T {\sqrt {{V_s}} {\rm{d}}W_s^{\left( 1 \right)}} } \right) $

它的敏感性参数可以直接求解微分得到，即

$ \Delta :\frac{{\partial C}}{{\partial {S_0}}} = E\left[ {\frac{{\partial {B_{\rm{S}}}\left( {{S_0}\xi ,\bar \sigma } \right)}}{{\partial {S_0}}};\bar \sigma } \right] = E\left[ {\xi N\left( {{d_1}} \right);\bar \sigma } \right] $

$ \begin{array}{l} \mathit{\Gamma }: = \frac{{{\partial ^2}C}}{{\partial S_0^2}} = E\left[ {\frac{{{\partial ^2}{B_{\rm{S}}}\left( {S_0^2\xi ,\bar \sigma } \right)}}{{\partial {S_0}}};\bar \sigma } \right] = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;E\left[ {\xi N'\left( {{d_1}} \right)\frac{1}{{{S_0}\bar \sigma \sqrt T }};\bar \sigma } \right] \end{array} $

$ \begin{array}{l} \rho :\frac{{\partial C}}{{\partial r}} = E\left[ {\frac{{\partial {B_{\rm{S}}}\left( {{S_0}\xi ,\bar \sigma } \right)}}{{\partial r}};\bar \sigma } \right] = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;E\left[ {KT\exp \left( { - rT} \right)N\left( {{d_2}} \right);\bar \sigma } \right] \end{array} $

表 5所示为3种方法中期权价格估计与Greeks分析.

由表 5可知，在SV和SVCJ模型下，蒙特卡罗法、对偶变量法和条件蒙特卡罗法对期权价格估计方差和Greeks方差中，条件蒙特卡罗法始终有较小方差，较其他两种方法，条件蒙特卡罗法有方差减小效果.

6 结论

本文在SV和SVCJ两种随机波动率模型下，基于标的资产过程精确模拟算法，考察了欧式看涨期权定价和敏感性Greeks的加速模拟理论计算问题，讨论了不同情形下条件蒙特卡罗法和对偶变量法两种方差减小技术的加速效果.

从模拟结果看，不同起始参数对条件蒙特卡罗的方差减小效果有一定影响，但条件蒙特卡罗法对期权价格估计和Greeks估计的方差始终小于其他两种方法的方差，说明精确模拟方法和条件蒙特卡罗法结合对期权的定价估计和Greeks估计有很好的结果.相比常用的欧拉离散与蒙特卡罗法，它能够得到无偏且更小方差的估计值.

本文算法可以很方便地解决其他更加复杂的产品的计算问题，如障碍期权的定价和敏感性估计问题、篮子期权的计算问题等.