第7章--蒙特卡洛模拟金融衍生产品定价.doc

第7章 蒙特卡洛模拟金融衍生产品定价 本章介绍蒙特卡洛模拟期权定价的内容，要求读者掌握随机数生成方式，了解蒙特卡洛定价就是模拟风险中性测度下标的资产的运动过程，学会蒙特卡洛方法模拟欧式期权定价，掌握提高模拟精度的常用方法。7.1 随机模拟基本原理 1977年，菲力埔伯耶勒(Phelim Boyle)提出了模拟方法求解金融资产定价问题，其想法是假设资产价格分布是随机波动，如果知道了这个波动过程，就可以通过随机模拟不同的路径，每做完一次模拟，就产生了一个最终资产价值，再进行若干次这样的过程，那么所得到的结果就是一个最终的资产价值分布，从这个分布中我们可以得到期望的资产价格。7.1.1 随机数生成函数1均匀分布随机数生成函数MATLAB中的unidrnd函数可以生成1到N的均匀分布随机数。 调用方式R=unidrnd(N);R=unidrnd(N,m);R=unidrn(N,m,n);其中，N所要生成的随机数个数，m确定输出随机矩阵R的行数,n确定输出随机矩阵R的列数2生成服从连续均匀分布的随机数如果需要生成服从连续分布的随机数，则需调用unifrnd函数，其调用格式为调用方式1Runifrnd(A，B)生成位于A、B之间的一个随机数。调用方式2Runifrnd(A,B,m)生成位于A、B之间的随机数。m=m1,m2表示行数列数。调用方式3Runifrnd(A,B,m,n),m,n分别表示行数、列数unifrnd(1,2,5,6),unifrnd(1,2,5,6)712 生成正态分布随机数 调用方式R=normrnd(mu,sigma)R=normrnd(mu,sigma,m)R=normrnd(mu,sigma,m,n)713 特定分布随机数发生器 MATLAB中有统一格式的随机数发生器，函数名称为random，可生成许多服从不同分布的随机数。y=random(name,A1,A2,A3,m,n)表8.1 生成特定分布的随机数函数参数表beta分布：beta,二项分布：bino,卡方:chi2,指数分布：exp,F-分布：f, Gamma:gamLognormal:logn, uniform:unif;Poisson:poiss,T:t;Normal-norm;Noncentral F -ncf, Noncentral-nct714 蒙特卡洛模拟方差削减技术蒙持卡洛模拟精度与模拟次数密切相关，模拟次数越高其精度越高，但是次数增加又会增加计算量。实践证明明减少模拟方差可以提高稳定性，减少模拟次数。有很多种方法可以减小方差，如对偶变量技术、控制变量技术、分层抽样、矩匹配、条件蒙特卡洛模拟等，但最简单并且应用最为广泛的是对偶变量技术与控制变量技术。 对偶变量技术就是先随机抽样得到一组数据，然后以此为基础构造出另一组对偶变量。下面以正态分布为例介绍对偶变量技术。首先从正态分布变量中随机抽取N个样本值，分别为，由此可以得到个模拟值，那么衍生证券蒙特卡洛估计值为 以为基础，构造对偶随机数，是与相互对偶的随机数，由正态分布的性质知，也是服从正态分布，由对偶随机数生成的估计值为对和取平均得到新的估计如果随机抽样的样本模拟得到的估计值比较小，那么与之对偶的随机抽样样本得到的估计值可能会偏大，二者的平均值就可能会接近真实值。如果，那么从上面的不等式可以看出，利用对偶技术可以增加估计稳定性，提高了估计精确度。7.1.5 随机模拟控制变量技术 摔制变量技术就是将与所估计的未知变量密切相关的另一个已知量的真实值和估计值之间的差异作为控制量，以提高估计精度。在定价实践中，将这两种衍生证券用相同的随机抽样样本和时间间隔，实施同样的蒙特卡洛模拟过程，能够得到两个模拟估计值，以第二种衍生证券真实值与估计值之间的差异作为控制变量，最后得到第一种衍生证券的蒙特卡洛估计值。 假定是需要估汁的第一种衍生证券的价值，是价值容易估计的第二种衍生证券的价值，第一种证券与第二种证券相似，而和分别是第一种衍生证券和第二种衍生证券在同样的随机抽样样本的蒙特卡洛估计值，那么利用控制变量技术得到第一种衍生证券的价格估计值为这里就是控制变量，它实际上是第一种衍生证券的蒙特卡洛模拟的估计误差，且上述方程的方差之间的关系为 如果，一定有因此，当两种衍生证券的协方差很大时，或者当两种衍生证券的价格高度相关时，上述关系是成立的，两种衍生证券的正相关性越强，估计效率越理想。然的从实际应用的角度看，这种控制变量技术的应用十分有限，因此，下面是更一般的控制变量技术，其控制变量的形式为方差为 这是关于控制变量系数的二次三项式，下面的目标是能够找到特殊的使方差最小.这时只要取就可以保证方差最小，这种控制变量技术的缺点是需要提前知道协方差的信息，而这一般需要靠经验实现。7.2 蒙持卡洛方法模拟期权定价(编程)考721 蒙特卡洛方法模拟欧式期权定价在期权计算中，我们可以利用风险中性的方法计算期权的价格。风险中性定价形式如下：其中，是期权的价格，是到期日的现金流，是风险中性测度。 如果知道了风险中性测度就可以模拟全路径，也可模拟终端价格，例如计算障碍期权等路径依赖型期权时可以模拟全路径，而欧式期权可模拟终端价格。如果标的资产服从几何布朗运动 那么风险中性定价的关键在于寻找风险中性测度。对于几何布朗运动，可以证明风险中性测度下，标的资产运动过程如下：对于欧式看涨期权，到期日欧式看涨期权现金流如下：其中，是执行价，是无风险利率，是标准差，是正态分布的随机变量。对到期日的现金流用无风险利率贴现，就可以知道期权的价格。【例7-1】假设股票价格服从几何布朗运动，股票现在价格，欧式期权执行价，无风险利率，股票波动的标准差，期权的到期日，试用蒙特卡洛模拟方法计算该期权价格。下面用MATLAB编写一个子程序blsmc进行计算function eucall,varprice,ci=blsmc(s0,K,r,T,sigma,Nu)% 蒙特卡洛方法计算欧式看涨期权的价格% 输入参教% sO:股票价格% K:执行价% r 无风险利率% sigma 股票波动的标准差% Nu 模拟的次数% 输出参数% eucall: 欧式看涨期权的价格% varprice模拟期权价格的方差% ci 95概率保证的期权价格区间randn(seed,O);%定义随机数发生器种子是0，这样可以保证每次模拟的结果相同nuT=(r-0.5*sigma2)*T;sit=sigma*sqrt(T);discpayoff=exp(-r*T)*max0,s0*exp(nuT+sit*randn(Nu,l)-K);%期权到期时的现金流eucall,varprice,ci=normfit(discpayoff)调用子程序可得到欧式看涨期权价格。c,var,ci=blsmc(50,52,0.1,5/12,0.4,1000)c = 5.444522854313244var = 9.136131838483149ci = 4.877583065095581 6.011462643530908从上面的结果可以看到，蒙特卡洛模拟得到的期权价格为5.4445，样本正态拟合的方差为9.1361，95的置信区间为4.8776，6.0115，模拟波动的区间还是很大的。我们用了normfit函数对模拟的结果用正态分布函数进行拟合，这不是必需的，主要是为了考察模拟结果的稳定性，如果不需要考察结果是否稳定，也可直接对模拟的结果求均值。此时可将最后一句改为price=mean(discpayoff)。欧式期权的公式解如下：call,put= blsprice(50,52,0.1,5/12,0.4)call = 5.191081351846794put = 5.068933121521976公式解5.1911和模拟值5.4445二者之间还是存在较大的差距，增加模拟的次数为10000次时结果如下：c,var,ci=blsmc(50,52,0.1,5/12,0.4,10000)c = 5.133786430504636var = 8.933464757068791ci = 4.958672541451405 5.308900319557867模拟结果入5.1338，可以看到期权模拟精度有了显著提高，95%的的置信区间为4.9587，5.3089，置信区间较1000时大大缩小，模拟可靠性增加。下面用对偶方法计算欧式看涨期权的价格。function eucall,varprice,ci=dul_blsmc(s0,K,r,T,sigma,Nu)% 蒙特卡洛方法计算欧式看涨期权的价格(对偶法）% 输入参教% sO:股票价格% K:执行价% r 无风险利率% sigma 股票波动的标准差% Nu 模拟的次数% 输出参数% eucall: 欧式看涨期权的价格% varprice模拟期权价格的方差% ci 95概率保证的期权价格区间randn(seed,O);%定义随机数发生器种子是0，这样可以保证每次模拟的结果相同nuT=(r-0.5*sigma2)*T;sit=sigma*sqrt(T);rand=randn(Nu,1);discpayoff=exp(-r*T)*max(0,s0*exp(nuT+sit*rand)-K);discpayoff1= exp(-r*T)*max(0,s0*exp(nuT+sit*(-1)*rand)-K);%期权到期时的现金流duleucall,dulvarprice,dulci=normfit(discpayoff;discpayoff1)调用子程序可得到欧式看涨期权价格。c,var,ci=dul_blsmc(50,52,0.1,5/12,0.4,10000);结果如下c = 5.209419581125931var = 8.902468445392321ci = 5.086032462943382 5.332806699308480模拟的结果为5.2094，离精确值5.19l1非常接近，说明对偶技术还是非常有效的722 蒙特卡洛方法模拟障碍期权定价障碍期权是特殊形式的期权，例如确定一个障碍值，在期权的存续期内有可能超过该价格，也有可能低于该价格，对于敲出期权而言，如果在期权的存续期内标的资产价格触及障碍值时，期权合同可以提前终止执行；相反地对于敲入期权丽言，如果标的资产价格触及障碍值时，期权合同开始生效。注意障碍值可以低于标的资产现在的价格，也可以高于。如果称为上涨期权，反之称为下跌期权。对于下跌敲出看跌期权，该期权首先是看跌期权，股票价格是，执行价是，买入看跌期权就首先保证以执行价卖掉股票，下跌敲出障碍期权相当于在看跌期权的基础上附加提前终止执行的条款，内容是当股票价格触及障碍值时看跌期权就提前终止执行。因为该期权对于卖方有利，所以其价格应低于看跃期权的价格。下面考虑下跌敲入看跌期权，同样地该期权首先是看跌期权，下跌敲入期权相当于在看跌期权的基础上，附加何时生效的条款，内容是当股票的价格触及障碍值时，看跌期权开始生效，综合地看，标准的看跌期权合同可以拆分为两份产品，分别是下跌敲出看跌期权与下跌敲入看跃期权，用公式表示如下其中，是标准看跌期权价格，与分别表示下跌敲入看跌期权与下跌敲出看跌期权的价格。如果下跌敲出看跌期权提前终止时卖方补偿一些费用给买方，上述公式表示的平价义系就就不再有效。当障碍值确定时，障碍期权存在公式解，其形式如下：其中，是股票价格，是障碍值，是看跌期权执行价，是存续期，是无风险利率，是波动率的标准差，其他参数如下：a=(Sb/S0)(-1+2*r/sigma2); b=(Sb/S0)(1+2*r/sigma2);d1=(log(S0/K)+(r+sigma2/2)*T)/(sigma*sqrt(T);d2=(log(S0/K)+(r-sigma2/2)*T)/(sigma*sqrt(T);d3=(log(S0/Sb)+(r+sigma2/2)*T)/(sigma*sqrt(T);d4=(log(S0/Sb)+(r-sigma2/2)*T)/(sigma*sqrt(T);d5=(log(S0/Sb)-(r-sigma2/2)*T)/(sigma*sqrt(T);d6=(log(S0/Sb)-(r+sigma2/2)*T)/(sigma*sqrt(T);d7=(log(S0*K/Sb2)-(r-sigma2/2)*T)/(sigma*sqrt(T);d8=(log(S0*K/Sb2)-(r+sigma2/2)*T)/(sigma*sqrt(T);利用上面的公式编写下跌敲出障碍期权价格的程序如下 function P=DownOutPut(S0,K,r,T,sigma,Sb)a=(Sb/S0)(-1+2*r/sigma2); b=(Sb/S0)(1+2*r/sigma2);d1=(log(S0/K)+(r+sigma2/2)*T)/(sigma*sqrt(T);d2=(log(S0/K)+(r-sigma2/2)*T)/(sigma*sqrt(T);d3=(log(S0/Sb)+(r+sigma2/2)*T)/(sigma*sqrt(T);d4=(log(S0/Sb)+(r-sigma2/2)*T)/(sigma*sqrt(T);d5=(log(S0/Sb)-(r-sigma2/2)*T)/(sigma*sqrt(T);d6=(log(S0/Sb)-(r+sigma2/2)*T)/(sigma*sqrt(T);d7=(log(S0*K/Sb2)-(r-sigma2/2)*T)/(sigma*sqrt(T);d8=(log(S0*K/Sb2)-(r+sigma2/2)*T)/(sigma*sqrt(T);P=K*exp(-r*T)*(normcdf(d4)-normcdf(d2)-a*(normcdf(d7)-normcdf(d5).-S0*(normcdf(d3)-normcdf(d1)-b*(normcdf(d8)-normcdf(d6);【例7-2】我们考虑一个欧式看跌股票期权。股票的价格为50，看跌期权执行价为50，无风险利率为0.1，时间为5个月，股票年波动率的标准差为0.4。首先公式的解如下：S0=50;K=50;r=0.1;T=5/12;sigma=0.4;call,put=blsprice(S0,K,r,T,sigma)call = 6.116508129330871put = 4.075980984787780 看跌期权价格为4.076. 对于上述看跌期权，进一步地，我们考虑障碍值等于40时下跌敲出期权的价格：P=DownOutPut(S0,K,r,T,sigma,40)P = 0.542418625858436 由于该下跌敲出看跌期权提供的条件过于优厚，买方承受大量风险，作为回报价格较看跌期权便宜许多。 下面用蒙特卡洛方法模拟下跌敲出看跌期权价格，在模拟中我们给出了模拟次数为NRep1,每次模拟时间分为NSteps步离散，障碍值为变量Sb，其现金流如下： 当时，CashFlow=0。我们可以先模拟路径，然后让小于Sb的路径的现金流为0，程序如下：% DOPutMC(s0,k,r,T,sigma,sb,NSteps,NRep1)function P,aux,CI= DOPutMC(s0,k,r,T,sigma,sb,NSteps,NRep1)% 利用蒙特卡洛方法对欧式下跌敲出期权定价% 输入参数% sO:股票价格% k:执行价% r 无风险利率% T:期权存续期 % sigma 股票波动的标准差% Sb 障碍值% NSteps 时间离散数目% NRepl: 模拟次数% 输出参数% P 下跌敲出看跌期权价格% aux 蒙特卡洛模拟的方差% CI 价格的95%-置信区间% 用蒙特卡洛方法模拟风险中性下股价路径dt=T/NSteps;nudt=(r-0.5*sigma2)*dt;sidt=sigma*sqrt(dt);randn(seed,0);rand=randn(NRep1,NSteps);rand1=nudt+sidt*rand;rand2=cumsum(rand1,2);% 沿列方向逐列累加path=s0*exp(rand2);% 利用路径进行定价payoff=zeros(NRep1，1);for i=1:NRep1ax=path(i,:); if min(ax)sb payoff(i)=0; %如果路径中的任意一点价格低于障碍值，现金流为0 else payoff(i)=max(0,k-ax(NSteps); endendP,aux,CI=normfit(exp(-r*T)*payoff);%P为期权价格运行程序：s0=50;k=50;r=0.1;T=5/12;sigma=0.4;sb=40;NSteps=600;NRep1=10000;P,aux,CI=DOPutMC(s0,k,r,T,sigma,sb,NSteps,NRep1)P = 0.592069149986891aux = 1.668695291705626CI = 0.559359363792887 0.624778936180894该期权模拟的价格为0.5594，和前面的公式解0.5424尚存在差距，增加模拟次数，代码如下：s0=50;k=50;r=0.1;T=5/12;sigma=0.4;sb=40;NSteps=600;NRep1=20000;P,aux,CI=DOPutMC(s0,k,r,T,sigma,sb,NSteps,NRep1)P = 0.559719264298951aux = 1.602236168134015CI = 0.537512467987729 0.581926060610173结果和公式解比较接近，如果将存续期改为2/12年，再考察其价格变化，代码如下：P,aux,CI=DOPutMC(s0,k,r,2/12,sigma,sb,60,50000)P = 1.352650569271812aux = 2.366421196429638CI = 1.331907851945229 1.3733932865983957.2.3 蒙特卡洛方法模拟亚式期权定价亚式期权是一种路径依赖型期权，它的收益函数依赖于期权存续期内标的资产的平均价格。平均价格分算术平均和几何平均两种，对于离散算术平均价格定义为 （8.1）其中是离散时间样本点。 离散几何平均价格定义为 （8.2）亚式看涨期权到期现金流为其中，是执行价，是时刻的股价。【例7-3】股票价格为50，亚式看涨期权执行价为50，存续则为5个月，期权到期现金流是每月均价与执行价之差，股票波动率的标准差为0.4，无风险利率为0.1，下面用蒙特卡洛方法计算该亚式期权价格：该期权定价程序为：function P,aux,CI= AsianMC(s0,k,r,T,sigma,NSteps,NRep1)% 利用蒙特卡洛方法对亚式期权定价% 输入参数% sO:股票价格% k:执行价% r 无风险利率% T:期权存续期 % sigma 股票波动的标准差% NSteps 时间离散数目% NRepl:模拟路径数% 输出参数% P 下跌敲出看跌期权价格% CI 蒙特卡洛模拟的方差% NCrossed 价格的95%-置信区间% 用蒙特卡洛方法模拟风险中性下股价路径dt=T/NSteps;nudt=(r-0.5*sigma2)*dt;sidt=sigma*sqrt(dt);randn(seed,0);rand=randn(NRep1,NSteps);rand1=nudt+sidt*rand;rand2=cumsum(rand1,2);% 沿列方向逐列累加path=s0*exp(rand2);% 利用路径进行定价payoff=zeros(NRep1，1);for i=1:NRep1 payoff(i)=max(0,mean(path(i,:)-k);endP,aux,CI=normfit(exp(-r*T)*payoff); 运行程序s0=50;k=50;r=0.1;T=5/12;sigma=0.4;NSteps=5;NRep1=50000;P,aux,CI= AsianMC(s0,k,r,T,sigma,NSteps,NRep1）P = 3.962173517636958aux = 5.966879973484831CI = 3.909871203577914 4.014475831696003这是一个比较粗糙的估计，我们可以用控制变量(control variable)技术提高估计精度。构造显然它和收益函数是相关的，而，这样用控制变量法可以将前面的程序改为% AsianMCCV.mfunction P, CI= AsianMCCV(S0,X,r,T,sigma,NSamples,NRep1,NPilot)% 利用蒙特卡洛方法对亚式期权定价% 输入参数% sO:股票价格% X:执行价% r 无风险利率% T:期权存续期 % sigma 股票波动的标准差% NSample 取样数目% NRepl:模拟路径数% Npilot: 路径数目% 输出参数% P 下跌敲出看跌期权价格% CI 蒙特卡洛模拟的方差% 用蒙特卡洛方法模拟风险中性下股价路径 TryPath=AssetPaths1(S0,r,sigma,T,NSamples,NPilot);%生成路径StockSum=sum(TryPath,2);PP=mean(TryPath(:,2:(NSamples+1),2);TryPayoff=exp(-r*T)*max(0,PP-X);MatCov=cov(StockSum,TryPayoff);c=-MatCov(1,2)/var(StockSum);dt=T/NSamples;ExpSum=S0*(1-exp(NSamples+1)*r*dt)/(1-exp(r*dt);% MCControlVars=zeros(NRep1,1);for i=1:NRep1 StockPath=AssetPaths1(S0,r,sigma,T,NSamples,1); Payoff(i)=exp(-r*T)*max(0,mean(StockPath(2:(NSamples+1)-X); ControlVars(i)=Payoff(i)+c*(sum(StockPath)-ExpSum);end P,aux,CI=normfit(ControlVars);function Spaths=AssetPaths1(S0,mu,sigma,T,NSteps,NRep1)dt=T/NSteps;nudt=(mu-0.5*sigma2)*dt;sidt=sigma*sqrt(dt);Increments=nudt+sidt*randn(NRep1,NSteps);LogPaths=cumsum(log(S0)*ones(NRep1,1),Increments,2);Spaths=exp(LogPaths);下面是控制变量模拟的结果：randn(seed,0);P,CI=AsianMCCV(50,50,0.1,5/12,0.4,5,45000,50000)计算结果为P = 3.980870822887436CI = 3.958134505545162 4.003607140229710根据前面蒙特卡洛模拟的结果有randn(seed,0);P,aux,CI=AsianMC(50,50,0.1,5/12,0.4,5,50000)P = 3.962173517636958aux = 5.966879973484831CI = 3.909871203577914 4.014475831696003可以看出控制变量使得估计的置信区间缩小了。7.2.4 蒙特卡洛模拟经验等价勒测度等价鞅理论是金融衍生品定价的重要方法，衍生品价格就是在等价鞅测度下对衍生品现金流用无风险利率贴现，可以用下面方程表示：其中，表示股票价格，为无风险利率，表示0时刻信息，表示等价敬测度下的期望算子。 Duan和simonate(1998)提出了乘数调整法，该方法可以保证等价鞅性质，当标的资产服从几何布朗运动时，等价鞅下的资产运动满足下列条件：其中，这样模拟的第条路径为这时我们可以构造一个等价鞅测度测度实际上在第条路径的权重上乘以，当样本容量趋于无穷大时有是风险中性测度下的新样本。对于例8-1我们重新编写程序如下%function Emseucall,Var,CI=Emsmc(s0,K,r,T,sigma,Nu)% 蒙特卡洛经验等价鞅方法计算欧式看涨期权的价格% 输入参教% sO:股票价格% K:执行价% r 无风险利率% sigma 股票波动的标准差% Nu 模拟