格式更好的html版本
蒙特卡洛方法介绍
根据风险中性定价公式（1.18），在时间T回报为x的证券的价值为 \(e^{-rT}\mathbb{E}^R[x] \tag{5.1}\) 用蒙特卡洛方法估计，是指用模拟方法得出随机变量x的样本值，将样本值平均得出数学期望估计值。 蒙特卡洛方法有两个主要缺陷：
- 很难（尽管不是不可能）用来计算提前执行的合约。 即无法在t时刻就整合所有未来路径的信息判断当前是否行权。
- 如果用计算所花时间衡量，这种方法是低效率的（尽管这种方法在计算以多个资产为标的衍生证券价值时，比其他方法的速度要快） 要使得估计标准误足够小，需要足够大的样本量，大量的计算要耗费很多时间。蒙特卡洛方法的困难之处在于如何减少所需样本量。
二叉树模型介绍
对于重复组合形成的树，当N增加时，计算时间以N的线性函数增加；而对非重复组合形成的树，计算时间以N的指数函数增加。遗憾的是，对路径依赖期权来说，这种计算时间的节省是无法实现的，因为一个树（不管是重复组合还是非重复组合的树）中的路径数是$2^N$。
欧式看涨期权价值的数学期望等于： \(\sum_{i=0}^N \frac{N!}{i!(N-i)!}p^i(1-p)^{N-i}\max(u^i d^{N-i} S-K,0) \tag{5.4}\) 将这个数学期望乘以$e^{-rT}$得出期权价值。
利用蒙特卡洛方法对欧式期权进行定价时，相当于用M-点分布来近似S(T)的分布，每个点赋以相等的概率。 二叉树模型中，选取的是S(T)的特定点的集合，并采用（5.4）的方法赋以它们概率，以此来近似S(T)的分布。 当所选取点的个数增加时，蒙特卡洛方法和二叉树方法得出的近似分布都收敛到S(T)的连续时间分布。但是，由于选取了特殊的点并赋以特定的概率，二叉树方法能够采用较少的点而得出同样精度的近似分布，即对于给定的精度，用蒙特卡洛方法需要M次模拟，而二叉树模型中的期数N 要比蒙特卡洛方法中的模拟次数M小得多。因此，二叉树方法的计算速度要快很多。此外，在美式期权定价中，与蒙特卡洛方法相比，二叉树方法更具有优越性。 在路径以来期权的定价中，蒙特卡洛方法更快。在以多个资产为标的的期权定价中，蒙特卡洛方法计算速度要快一些。
美式期权的二叉树模型
最优执行策略是指，当期权内在价值超过期权未来价值期望的贴现时执行期权。 因此，沿二叉树进行反向递归时，检验每个节点是否为最优执行点，如果是，则用期权内在价值代替期望贴现值。 对于美式看跌期权： \(P(i)=\max(K-u^i d^{n-i}S, e^{-r\Delta t}pP(i+1)+e^{-r\Delta t}(1-p)P(i)) \tag{5.9}\)
二叉树模型的参数
考虑长度为T年的N期二叉树模型。每期长度为$\Delta t=T/N$。对连续时间模型，在离散时间区间$\Delta t$上有 \(\Delta \log S =v\Delta t+\sigma \Delta B\) 其中$v=r-q-\sigma^2/2$，B为风险中性概率测度下的布朗运动。在风险中性概率测度下，连续模型中$\Delta \log S$的均值和方差分别为 \(\mathbb{E}^R[\Delta \log S]=v\Delta t , \text{var}^R[\Delta \log S]=\sigma^2 \Delta t\) 因此， \(\frac{\mathbb{E}^R[\Delta \log S]}{\Delta t}=v , \frac{\text{var}^R[\Delta \log S]}{\Delta t}=\sigma^2\) 在二叉树模型中， \(\frac{\mathbb{E}^R[\Delta \log S]}{\Delta t}=\frac{p\log u+(1-p)\log d}{\Delta t}\\ \frac{\text{var}^R[\Delta \log S]}{\Delta t}=\frac{p(1-p)(\log u-\log d)^2}{\Delta t}\) 固定时间总长度T不变，令$N\rightarrow \infty$（等价于令$\Delta t \rightarrow 0$），二叉树模型在某种意义下收敛到连续时间模型的充分条件为 \(\frac{p\log u +(1-p)\log d}{\Delta t}\rightarrow v\\ \frac{p(1-p)(\log u -\log d)^2}{\Delta t}\rightarrow \sigma^2\) 最常用的模型是CRR，即取d=1/u，并且 \(u = \mathrm{e}^{\sigma \sqrt{\Delta t}} \tag{5.10a}\)
\(p = \frac{\mathrm{e}^{(r-q) \Delta t} - d}{u - d} \tag{5.10b}\) 另一个很有名的模型是取p=1/2，并且 \(u = \exp\left(\left(r - q - \frac{1}{2}\sigma^{2}\right)\Delta t + \sigma\sqrt{\Delta t}\right) \tag{5.11a}\)
\[d = \exp\left(\left(r - q - \frac{1}{2}\sigma^{2}\right)\Delta t - \sigma\sqrt{\Delta t}\right) \tag{5.11b}\]Trigeorgis提出，p、u和d的选择应该使二叉树模型中$\Delta \log S$的均值和方差与连续模型中的均值和方差完全相等，即 \(\frac{p\log u+(1-p)\log d}{\Delta t}=\nu \tag{1}\)
\(\frac{p(1-p)(\log u-\log d)^{2}}{\Delta t}=\sigma^{2} \tag{2}\) 两个方程有三个未知数，取d=1/u，可得 \(\log u = \sqrt{\sigma^{2}\Delta t + \nu^{2}(\Delta t)^{2}} \tag{5.13a}\)
\[p = \frac{1}{2} + \frac{\nu\Delta t}{2\log u} \tag{5.13b}\]二叉树模型中的希腊字母
对任何定价模型来说，要估计希腊字母的值，只需用不同的参数值运行两次定价程序，然后用两个定价之差除以参数值的差。 从理论上来说，参数变化量越小，估计出来的希腊字母越精确，但在实际计算中，计算机四舍五入带来的误差会限制采用过小的参数变化量。 \(\delta=\frac{C_u-C_d}{S_u-S_d} \tag{5.14}\)
\(\Gamma=\frac{\delta_u-\delta_d}{(S_u-S_d)/2} \tag{5.15}\) 在二叉树模型中，可以通过更为有效的方法对两个最重要的希腊字母$\delta$和$\gamma$进行估计，而不必采用多次运行定价程序的方法。 如在（5.14）和（5.15）中取$S_u=u^2S$和$S_d=d^2S$
蒙特卡洛方法中的希腊字母I：差值比
为了尽量减少计算误差和计算次数，在计算不同参数下的衍生证券价值时，可以采用相同的随机抽样。
\(Var(X-Y)=Var(X)+Var(Y)-2Cov(X,Y)\) 采用相同的抽样，可以让Cov(X,Y)为0,减少方差 在BS假设下，只需要产生出S(T)就可以采用$[S_u(0)/S(0)]\times S(T)$计算出$S_u(T)$，$S_d(T)$同理，然后采用与计算衍生证券价值标准误相同的方法可以计算希腊字母值的标准误。 资产价格对数过程： \(\begin{aligned} \log S_d(T) &= \log S_d + \left(r - q - \frac{1}{2}\sigma^2\right) T + \sigma B(T) \\ \log S(T) &= \log S + \left(r - q - \frac{1}{2}\sigma^2\right) T + \sigma B(T) \\ \log S_u(T) &= \log S_u + \left(r - q - \frac{1}{2}\sigma^2\right) T + \sigma B(T) \end{aligned}\)
推导结果： \(\begin{aligned} S_d(T) &= \left(\frac{S_d}{S}\right) S(T) \\ S_u(T) &= \left(\frac{S_u}{S}\right) S(T) \end{aligned}\) 因此，在BS假设下，只需要模拟S(T)的值，然后按照上面的公式相乘得出$S_d(T)$和$S_u(T)$ 以上方法适合对于无路径依赖的衍生证券。 对于更一般的欧式衍生证券：
- $\delta$ 差值$(C_u(T)-C_d(T))的样本均值乘以$e^{-rT}/(S_u-S_d)$。 标准误= \(\frac{\mathrm{e}^{-rT}}{S_u - S_d} \sqrt{ \frac{1}{M(M-1)} \left( \sum_{i=1}^M \left[ C_{ui}(T) - C_{di}(T) \right]^2 - M \left[ \overline{C_u(T)} - \overline{C_d(T)} \right]^2 \right) }\) 标准误为$\frac{s}{\sqrt{M}}$，而$s^2=\frac{1}{M-1}\sum_{i=1}^M[x_i-\bar{x}]^2=\frac{1}{M-1}(\sum_{i=1}^M x_i^2-M\bar{X}^2)$
-
$\gamma$ 伽玛的蒙特卡洛估计也是一个样本均值： \(\Gamma = \frac{2}{(S_u - S)(S_u - S_d)} C_u - \frac{2}{(S_u - S)(S - S_d)} C + \frac{2}{(S - S_d)(S_u - S_d)} C_d \tag{5.16}\) 通常情况下取$S_u=(1+\alpha)S,S_d=(1-\alpha)S$，$\alpha$为某个常数。此时，(5.16)为 \(\Gamma=\frac{C_u-2C+C_d}{\alpha^2 S^2} \tag{5.17}\) 而伽玛的标准误为 $$ \frac{\mathrm{e}^{-rT}}{\alpha^2 S^2} \sqrt{\frac{1}{M(M - 1)}} \times \sqrt{\sum_{i=1}^M \left[ C_{ui}(T) - 2C_i(T) + C_{di}(T) \right]^2 - M \left[ \overline{C_u}(T) - 2\overline{C}(T) + \overline{C_d}(T) \right]^2}
$$
蒙特卡洛方法中的希腊字母II：按路径估计
参数估计是“有偏的”，是指估计量期望值不等于参数真值。注意，如果蒙特卡洛估计值是有偏的，即使模拟次数很多、标准误很接近0,蒙特卡洛方法给出的结论也是不正确的。 该节阐述了在期权定价中，如何使用蒙特卡洛方法计算希腊值“Delta”，并重点分析了传统“有限差分法”存在的偏差问题，以及更优的“路径法”估计原理。
核心主题：蒙特卡洛方法计算 Delta 的两种方式及其偏差
整个讨论的背景是：我们无法用公式直接计算期权价格对标的资产价格（S）的敏感度（Delta）时，会采用蒙特卡洛模拟，通过生成大量资产价格的未来路径来估算。关键问题在于如何从模拟的路径中计算出准确的 Delta。
第一部分：传统方法——“有限差分法”及其问题
1. 方法描述
这种方法很直观，类似于数值求导：
- 设定一个初始标的价格
S。 - 对原始价格
S进行一个微小的扰动，得到两个新价格：S_u（向上扰动）和S_d（向下扰动）。 - 分别基于
S,S_u,S_d进行蒙特卡洛模拟，计算对应的期权到期价值C(T),C_u(T),C_d(T)。 - Delta 的估计值就是折扣后的样本均值： \(\frac{C_{u}(T-C_{d}(T)}{S_{u}-S_{d}} \tag{5.18}\)
2. 偏差的来源
文本通过分析期权到期时的三种状态，揭示了这种估计方法存在偏差，即估计值的期望不等于真实的 Delta。
- 情况一：深度虚值
- 条件： $S_u(T) ≤ K$ （即即使价格向上扰动，到期时期权仍无价值）
- 结果： $C_u(T) = C_d(T) = 0$，因此比值 (5.18) 为 0。
- 情况二：深度实值
- 条件： $S_d(T) ≥ K$ （即即使价格向下扰动，到期时期权仍为实值）
- 结果： $C_u(T) = S_u(T) - K$，$C_d(T) = S_d(T) - K$。此时比值 (5.18) 恰好等于 S(T)/S。
- 情况三：临界状态（偏差根源）
- 条件： $S_u(T) > K > S_d(T)$ （即未扰动时，期权可能在虚值和实值之间）
- 结果：只有 $C_u(T)$ 有价值，$C_d(T) = 0$。此时比值 (5.18) 是一个复杂的表达式，且小于
S(T)/S。
核心问题：真实的 Delta 在数学上等于对以下导数求期望： \(\frac{\partial{\max(0,S(T-K))}}}\) 这个导数非常简单：
- 若 $S(T) > K$，导数为 $S(T)/S$。
- 若 $S(T) < K$，导数为
0。
而“有限差分法”的估计量 (5.18) 在临界状态（情况三）下的行为与真实的导数不一致，导致其期望值（公式5.22）与真实Delta的期望值（公式5.21）不同。这就产生了偏差。
偏差的严重后果：即使进行无限次模拟，将标准误差降为零，这个估计值也会系统地偏离真实值，给出错误的答案。
第二部分：解决方案——“路径法”估计
既然知道了问题所在，一个更优的方案被提出：路径法。
1. 方法描述
这种方法直接、高效，完全避免了有限差分法的问题。
- 直接基于初始价格
S生成样本路径。 - 对于每一条路径，如果期权到期时是实值的（即
S(T) > K），则直接计算该路径对 Delta 的贡献为S(T)/S；如果是虚值的，贡献为0。 - Delta 的估计值就是所有这些贡献值的折扣样本平均。
- 用数学公式表示，就是直接估计真实的 Delta（公式5.21）：
$e^{-rT}E^{R}\left[\frac{S(T}{S}\mathbf{1}_{S(T)>K}\right])$
其中
𝟙是指示函数，当S(T) > K时为1，否则为0。
2. 路径法的优势
- 无偏性：因为它直接估计了真实的 Delta 表达式，所以是一个无偏估计。
- 效率高：只需要对原始路径进行模拟，无需像有限差分法那样对每条路径做两次扰动模拟，计算成本更低。
- 标准误差可计算：其结果是一个简单的样本平均值，因此其标准误差可以像常规蒙特卡洛估计一样方便地计算。
总结
| 特性 | 有限差分法 | 路径法 |
|---|---|---|
| 核心思想 | 通过扰动初始价格，模拟“上/下”两条路径，用差分近似导数。 | 直接对每条路径计算其支付函数关于标的价格的解析导数。 |
| 是否存在偏差 | 是，在期权价值临界点处行为与真实导数不符。 | 否，直接估计真实导数的期望。 |
| 计算效率 | 较低，每份期权需模拟三次（原始、上、下）。 | 较高，每份期权只需模拟一次。 |
| 准确性 | 即使模拟次数无限，结果也可能不正确。 | 随着模拟次数增加，结果收敛于真实值。 |
因此，在可以使用路径法的情形下（即支付函数关于标的价格的导数存在且可求），它通常是计算希腊值的更优选择。这种方法由 Broadie 和 Glasserman 提出，并可以扩展到其他希腊值（如 Gamma, Vega等）和其他更复杂的模型。
实现
欧式看涨期权的蒙特卡洛定价
期权价格$e^{-rT}\mathbb{E}^R[x]$ 差值比方法计算的德尔塔$e^{-rT}\frac{\bar{C_u(T)-C_d(T)}}{S_u-S_d}$ 按路径计算的德尔塔$e^{-rT}\mathbb{E}[\frac{S(T)}{S}x]$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
use rand::Rng;
use rand_distr::StandardNormal;
pub fn European_Call_MC(S:f64,K:f64,r:f64,sigma:f64,q:f64,T:f64,M:usize)->(f64,f64,f64){
let LogS0=S.ln();
let drift=(r-q-0.5*sigma.powi(2))*T;
let SigSqrtT=sigma*T.sqrt();
let UpChange:f64=(1.01_f64).ln();
let DownChange:f64=(0.99_f64).ln();
let mut SumCall=0.0;
let mut SumCallChange=0.0;
let mut SumPathwise=0.0;
let mut rng=rand::rng();
for i in 0..M{
let rand_increment:f64=rng.sample(rand_distr::StandardNormal);
let LogS:f64=LogS0+drift+SigSqrtT*rand_increment; //log S(T)的模拟值
let CallV:f64=(LogS.exp()-K).max(0.0); //期权价值
SumCall+=CallV; //期权价值求和
let LogSu=LogS+UpChange; //log S_u(T)的模拟值
let CallVu:f64=(LogSu.exp()-K).max(0.0); //期权价值
let LogSd=LogS+DownChange; //log S_d(T)的模拟值
let CallVd:f64=(LogSd.exp()-K).max(0.0); //期权价值
SumCallChange+=CallVu-CallVd; //期权价值之差
if LogS.exp()>K{
SumPathwise+=LogS.exp()/S; //路径求和
}
}
let CallV:f64=(-r*T).exp()*SumCall/M as f64;
let Delta1:f64=(-r*T).exp()*SumCallChange/(M as f64*0.02*S);
let Delta2:f64=(-r*T).exp()*SumPathwise/M as f64;
(CallV,Delta1,Delta2)
}
GARCH模型的蒙特卡洛定价
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
pub fn European_Call_GARCH_MC(S:f64,K:f64,r:f64,sigma0:f64,q:f64,T:f64,N:usize,kappa:f64,theta:f64,lambda:f64,M:usize)->(f64,f64){
//输入参数：
//S:股票初始价格
//K:执行价格
//sigma0:初始波动率
//q:红利支付率
//T:到期时间
//N:时间区间数
//kappa、theta、lambda:GARCH参数
//M:模拟次数
let deltaT:f64=T/N as f64;
let sqrtT:f64=deltaT.sqrt();
let a=kappa*theta;
let b:f64=(1.0-kappa)*lambda;
let c:f64=(1.0-kappa)*(1.0-lambda);
let LogS0=S.ln();
let mut SumCall=0.0;
let mut SumCallSq=0.0;
let mut rng=rand::rng();
for i in 0..M{
let mut LogS=LogS0;
let mut Sigma=sigma0;
for j in 0..N{
let rand_increment:f64=rng.sample(rand_distr::StandardNormal);
let y:f64=Sigma*rand_increment;
LogS+=(r-q-0.5*Sigma.powi(2))*deltaT+sqrtT*y;
Sigma=(a+b*y.powi(2)+c*Sigma.powi(2)).sqrt(); //更新波动率
}
let CallV:f64=(LogS.exp()-K).max(0.0);
SumCall+=CallV;
SumCallSq+=CallV.powi(2);
}
let CallV=(-r*T).exp()*SumCall/M as f64;
let StdError=(-r*T).exp()*((SumCallSq-SumCall.powi(2)/M as f64)/(M*(M-1)) as f64).sqrt();
(CallV,StdError)
}
其中标准误的计算公式是： \(\sqrt{\frac{1}{N(N-1)}(\sum_{i=1}^{N}x_i^2-N\bar{x}^2)}\)
欧式期权的二叉树定价
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
pub fn Europea_Call_Binomial(S:f64,K:f64,r:f64,sigma:f64,q:f64,T:f64,N:usize)->f64{
let dt=T/N as f64; //时间区间长度
let u=(sigma*dt.sqrt()).exp(); //上升步长
let d=1.0/u; //下降步长
let pu=(((r-q)*dt).exp()-d)/(u-d); //上升概率
let pd=1.0-pu; //下降概率
let u2=u*u;
//计算底部节点（即每次都下降达到的节点）处的股票价格、到达该节点的概率以及公式(5.4)中的第一项
let mut S=S*d.powi(N as i32);
let mut prob=pd.powi(N as i32);
let mut CallV=prob*(S-K).max(0.0);
for i in 1..=N{
S*=u2;
prob*=(pu/pd)*(N-i+1) as f64/i as f64;
CallV+=prob*(S-K).max(0.0);
}
(-r*T).exp()*CallV
}
上面代码中关于i次上升的概率和i-1次上升的概率之比为： \(\frac{p^i(1-p)^{N-i}N!/i!(N-i)!}{p^{i-1}(1-p)^{N-i+1}N!/(i-1)!(N-i+1)!}=\frac{(N-i+1)p}{(1-p)i}\)
美式期权的二叉树定价
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
pub fn American_Put_Binomial(S0:f64,K:f64,r:f64,sigma:f64,q:f64,T:f64,N:usize)->f64{
//输入参数为
//S0:初始股票价格
//K:执行价格
//r:无风险利率
//sigma:波动率
//q:红利支付率
//T:到期时间
//N:时间区间个数
let dt=T/N as f64; //时间区间长度
let u=(sigma*dt.sqrt()).exp(); //上升步长
let d=1.0/u; //下降步长
let pu=(((r-q)*dt).exp()-d)/(u-d); //上升概率
let pd=1.0-pu; //下降概率
let discount=(-r*dt).exp();
let dpu=discount*pu; //贴现上升概率
let dpd:f64=discount*pd; //贴现下跌概率
let u2=u*u;
let mut S=S0*d.powi(N as i32); //底部节点的股票价格
let mut PutV=vec![0.0;N+1];
PutV[0]=(K-S).max(0.0); //底部节点看跌期权价值
for i in 1..=N{
S*=u2;
PutV[i]=(K-S).max(0.0);
}
for i in (0..N).rev(){ //沿时间倒回到0时
let mut S_cur=S0*d.powi(i as i32); //底部节点的股票价格
PutV[0]=(K-S_cur).max(dpd*PutV[0]+dpu*PutV[1]);
for j in 1..=i{
S_cur*=u2;
PutV[j]=(K-S_cur).max(dpd*PutV[j]+dpu*PutV[j+1]);
}
}
PutV[0]
}
德尔塔和伽玛的二叉树估计
```rust pub fn American_Put_Binomial_Delta_Gamma(S0:f64,K:f64,r:f64,sigma:f64,q:f64,T:f64,N:usize)->(f64,f64,f64){ //输入参数为 //S0:初始股票价格 //K:执行价格 //r:无风险利率 //sigma:波动率 //q:红利支付率 //T:到期时间 //N:时间区间个数 //返回(put value,delta,gamma) let dt=T/N as f64; let NewN=N+2; let u=(sigmadt.sqrt()).exp(); let d=1.0/u; let pu=(((r-q)dt).exp()-d)/(u-d); let pd=1.0-pu; let dpu=(-rdt).exp()pu; let dpd=(-rdt).exp()pd; let u2=u*u;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
let mut S=S0*d.powi(NewN as i32);
let mut PutV=vec![0.0;NewN+1];
let mut S=S0*d.powi(NewN as i32);
PutV[0]=(K-S).max(0.0);
for i in (2..=NewN-1).rev(){
let mut S_cur=S0*d.powi(i as i32);
PutV[0]=(K-S_cur).max(dpd*PutV[0]+dpu*PutV[1]);
for j in 1..=i{
S_cur*=u2;
PutV[j]=(K-S_cur).max(dpd*PutV[j]+dpu*PutV[j+1]);
}
}
let Su=S0*u2; //高股票价格
let Sd=S0/u2; //低股票价格
let DeltaU=(PutV[2]-PutV[1])/(Su-S0); //上中点德尔塔
let DeltaD=(PutV[1]-PutV[0])/(S0-Sd); //下中点德尔塔
let distance=S0*(u*u-d*d); //Su和Sd的距离
let delta=(PutV[2]-PutV[0])/distance;
let gamma=2.0*(DeltaU-DeltaD)/distance;
(PutV[1],delta,gamma)
}
```
