layout: post
title: "《衍生证券教程》第五章 蒙特卡洛来模拟方法和二叉树模型 笔记"
description: ""
author: "Hao"
category: true
tagline: "衍生品"
tags: [finance math , derivative pricing]
根据风险中性定价公式（1.18），在时间T回报为x的证券的价值为
e^{-rT}\mathbb{E}^R[x] \tag{5.1}用蒙特卡洛方法估计，是指用模拟方法得出随机变量x的样本值，将样本值平均得出数学期望估计值。
蒙特卡洛方法有两个主要缺陷：
对于重复组合形成的树，当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
为了尽量减少计算误差和计算次数，在计算不同参数下的衍生证券价值时，可以采用相同的随机抽样。
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)
以上方法适合对于无路径依赖的衍生证券。
对于更一般的欧式衍生证券：
\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) }\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} \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}
参数估计是“有偏的”，是指估计量期望值不等于参数真值。注意，如果蒙特卡洛估计值是有偏的，即使模拟次数很多、标准误很接近0,蒙特卡洛方法给出的结论也是不正确的。
该节阐述了在期权定价中，如何使用蒙特卡洛方法计算希腊值“Delta”，并重点分析了传统“有限差分法”存在的偏差问题，以及更优的“路径法”估计原理。
整个讨论的背景是：我们无法用公式直接计算期权价格对标的资产价格（S）的敏感度（Delta）时，会采用蒙特卡洛模拟，通过生成大量资产价格的未来路径来估算。关键问题在于如何从模拟的路径中计算出准确的 Delta。
这种方法很直观，类似于数值求导：
S。S 进行一个微小的扰动，得到两个新价格：S_u（向上扰动）和 S_d（向下扰动）。S, S_u, S_d 进行蒙特卡洛模拟，计算对应的期权到期价值 C(T), C_u(T), C_d(T)。\frac{C_{u}(T-C_{d}(T)}{S_{u}-S_{d}} \tag{5.18}文本通过分析期权到期时的三种状态，揭示了这种估计方法存在偏差，即估计值的期望不等于真实的 Delta。
情况一：深度虚值
情况二：深度实值
情况三：临界状态（偏差根源）
S(T)/S。核心问题：真实的 Delta 在数学上等于对以下导数求期望：
\frac{\partial{\max(0,S(T-K))}}}这个导数非常简单：
0。而“有限差分法”的估计量 (5.18) 在临界状态（情况三）下的行为与真实的导数不一致，导致其期望值（公式5.22）与真实Delta的期望值（公式5.21）不同。这就产生了偏差。
偏差的严重后果：即使进行无限次模拟，将标准误差降为零，这个估计值也会系统地偏离真实值，给出错误的答案。
既然知道了问题所在，一个更优的方案被提出：路径法。
这种方法直接、高效，完全避免了有限差分法的问题。
S 生成样本路径。S(T) > K），则直接计算该路径对 Delta 的贡献为 S(T)/S；如果是虚值的，贡献为 0。𝟙 是指示函数，当 S(T) > K 时为1，否则为0。| 特性 | 有限差分法 | 路径法 |
|---|---|---|
| 核心思想 | 通过扰动初始价格，模拟“上/下”两条路径，用差分近似导数。 | 直接对每条路径计算其支付函数关于标的价格的解析导数。 |
| 是否存在偏差 | 是，在期权价值临界点处行为与真实导数不符。 | 否，直接估计真实导数的期望。 |
| 计算效率 | 较低，每份期权需模拟三次（原始、上、下）。 | 较高，每份期权只需模拟一次。 |
| 准确性 | 即使模拟次数无限，结果也可能不正确。 | 随着模拟次数增加，结果收敛于真实值。 |
因此，在可以使用路径法的情形下（即支付函数关于标的价格的导数存在且可求），它通常是计算希腊值的更优选择。这种方法由 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]
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)
}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)
} 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)
} 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)} 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
} 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}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]
}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]
}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=(sigma*dt.sqrt()).exp();
let d=1.0/u;
let pu=(((r-q)*dt).exp()-d)/(u-d);
let pd=1.0-pu;
let dpu=(-r*dt).exp()*pu;
let dpd=(-r*dt).exp()*pd;
let u2=u*u;
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)
}
```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=(sigma*dt.sqrt()).exp();
let d=1.0/u;
let pu=(((r-q)*dt).exp()-d)/(u-d);
let pd=1.0-pu;
let dpu=(-r*dt).exp()*pu;
let dpd=(-r*dt).exp()*pd;
let u2=u*u;
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)
}
```