第九章-期权定价的有限差分方法

1第九章期权定价的有限差分方法在本章中，我们将给出几个简单的例子来说明基于偏微分方程（PDE）框架的期权定价方法。具体的方法的是利用第五章中讲述的有限差分方法来解决Black-scholes偏微分方程。在9.1节中，我们会回顾衍生品定价的数值解法以及指出如何利用适当的边界条件来模拟一个特定的期权。在9.2节中我们将会应用简单的显式（差分）方法来求解一个简单的欧式期权。正如你已熟知的那样，这种方法只能解出一些可以从金融角度来解释的不稳定的数值解。在9.3节中我们将可以看到使用完全的隐式方法可以解决这种不稳定问题。在9.4节中我们将介绍Crank-Nicolson方法在障碍期权定价中的应用，它可以看做是一种显式与完全隐式方法的混合。最后，在9.5节中，我们会看到迭代松弛方法可以用于解决使用全隐式方法来解决美式期权定价时由于存在提前执行的可能性而导致的自由边界问题。9.1使用有限差分法解BS方程在2.6.2节中，我们给出了一个标的资产在时间t的价格为)(tS的期权，该期权的价格是一个函数),(tSf，且),(tSf满足偏微分方程（9.1）通过不同的边界条件可以让这个方程刻画不同的期权的特征。在某些地方可能因为假设的改变或者对路径依赖的改变而导致方程式的具体形式改变，但是此处仅仅作为一个起点，帮助读者了解如何应用基于有限差分方法来解决期权定价的问题。正如我们在第五章中遇到的情况那样，要用有限差分方法来解偏微分方程，在此处我们必须建立资产价格和时间的离散网格。设T是2期权的到期日，而Smax是一个足够大的资产价格，在我们所考虑的时间范围内，)(tS的数值不能超过Smax。设定Smax是因为偏微分方程的区域关于资产价格是无边界的。但是为了达到计算的目的，必须要求它是有界的。Smax相当于+∞。网格通过点(S,t)取得，其中(S,t)满足0S，S，S2，……,maxSSM,0t,t，t22，……,TtN。本章中使用网格符号为,我们回顾一下（9.1）方程式的几种不同解法：向前差分向后差分中心（或对称）差分对于第二个差分式子，有至于究竟采用哪种方法进行离散化，我们将在后面的实际操作过程中对显式和隐式的方法作出详细的阐述说明。另一个值得我们注意的问题是如何设置边界条件。对于执行价格3为K的欧式看涨期权，其终值条件为对于执行价格为K的欧式看跌期权，其终值条件为但是当我们涉及到资产价格的边界问题时，这个问题就变得复杂了。因为在数学上要求必须在一个有界区域里来解这个方程，而资产价格的区域是无界的。我们可以通过一些例子来说明这个问题。例9.1我们首先考虑一支简单的欧式看跌期权。当资产的价格)(tS非常大时，看跌期权几乎是一文不值，因为我们几乎可以肯定的是其执行价低于现价：Smax的值必须足够大才能使得假设的边界条件是合理的。当资产的价格0)(tS时，鉴于资产的价格符合几何布朗运动的动力学模型，我们可以认为资产的价格会保持为0；因此到期收益为K，在时间t贴现值为其边界条件可以写为：例9.2我们可以从9.1推导出处理普通欧式看涨期权的方法。不妨设资产的价格0)(tS,那么在任何时刻t，我们都可以知道该期权价值为0：4而对于一个足够大的资产价格)(tS,我们可以肯定的是该期权一定为价内期权，即具有大于零的价值的期权，那么期满行权我们可以得到收益为KTS)(。将终值K贴现到t时刻，并且考虑标的资产的价格是一个函数)(tS，那么我们根据无套利定价原理可以得到其边界条件如下：具体条件可以写为：对于资产价格S较大时，需要一个替代的边界条件，但这时候必须要求期权的Δ为1；在这种情况下，我们所提出的边界条件一般针对未知函数的倒数，而非函数本身。这便是数学物理方程中所谓的Neumann边界条件。我们不会采取这种做法，因为这将使得数值解法复杂化了。当我们面对障碍式期权时，事情可能就变得简单了。对敲出期权来说，例如一个向下触及失效看跌期权，在标的资产价格触碰到障碍时期权作废，价值为0。向上触及失效看涨期权的情况也是类似的，它们的优点都在于我们计算必须涉及到的边界即域是自然有界的（显然存在的）。美式期权的定价就显得更加复杂难以处理了，因为这涉及到了提前执行边界的问题。我们必须考虑什么水平的资产价格和在什么时间是期权的最佳实施点。因此，在我们解决问题的过程中必然产生一个自由边界。这就需要为不同的奇异期权设定不同的边界条件。如何计算出正确的边界条件以及使用数值方法逼近正确的结果，这依赖于对具体期权的选择。5利用显式法求解欧式期权由于初次尝试解方程（9.1），我们不妨考虑一个简单的欧式看跌期权。我们分别用中心差分对S进行求导，用向后差分对时间进行求导来逼近求解。这并不是唯一的方法，但是无论采取什么样的方法都必须保证边界条件是符合特定期权模型的。例9.1中的边界条件差分后的结果可以用以下的方程式来表示：（9.2）需要特别指出说明的是，由于存在终点条件的假设，因此方程必须适当地用回推法求解。不妨设（9.2）式中有j=N，由于存在给定的终点条件，我们有一个未知量,1iNf，可以认为它是关于三个已知量的函数。对其他的每个时间层（i）作同样的考虑。整理方程，我们可以得到：（9.3）此时有6functionprice=EuPutExpl(S0,K,r,T,sigma,Smax,dS,dt)%建立网格，并且在必要时调整增量M=round(Smax/dS);dS=Smax/M;N=round(T/dt);dt=T/N;matval=zeros(M+1,N+1);vetS=linspace(0,Smax,M+1);veti=0:M;vetj=0:N;%建立边界条件matval(:,N+1)=max(K-vetS,0);matval(1,:)=K*exp(-r*dt*(N-vetj));matval(M+1,:)=0;%建立三对角矩阵a=0.5*dt*(sigma^2*veti-r).*veti;b=1-dt*(sigma^2*veti.^2+r);c=0.5*dt*(sigma^2*veti+r).*veti;%求解方程forj=N:-1:1fori=2:Mmatval(i,j)=a(i)*matval(i-1,j+1)+b(i)*matval(i,j+1)+...7c(i)*matval(i+1,j+1);endend%返回价格，有可能在网格外线性插值生成price=interp1(vetS,matval(:,1),S0);图9.1简易欧式期权显式差分方法定价的MATLAB代码定价不是直接利用MATLAB来完成的。图9.1中是定价所使用的MATLAB的程序代码，另外计算还需要给出Smax的值以及上文提到的两个离散步骤。唯一一点需要注意的是，在数学符号中幂指数可以从0开始，但是在MATLAB的矩阵计算中，幂指数是从1开始的。另外，如果资产的初始价格不在网格上，那么我们只能够寻找两个相邻点进行插值。一般我们采取比较粗略的线性插值。但是如果有需要计算分析期权价格的敏感性（在实际案例中经常需要这样做）则需要进行更加复杂的插值计算。[c,p]=blsprice(50,50,0.1,5/12,0.4);PP=4.076EuPutExpl(50,50,0.1,5/12,0.4,100,2,5/1200)ans=4.0669[c,p]=blsprice(50,50,0.1,5/12,0.3);PP=2.8446EuPutExpl(50,50,0.1,5/12,0.3,100,2,5/1200)8ans=2.8288我们可以看到，数值方法能够提供比较精确的结果，为了提高精度，我们也可以采用更加精细的网格进行计算。EuPutExpl(50,50,0.1,5/12,0.3,100,1.5,5/1200)ans=2.8597EuPutExpl1(50,50,0.1,5/12,0.3,100,1,5/1200)ans=-2.8271e+022在此处我们看到在第五章中分析过的不稳定数值解的例子。有一种解决方案就是使用隐式方法进行求解。另外一种方法则是进行稳定性分析，然后在从离散化过程中导出边界。此处我们不会再使用第二种方案，这与第五章中对热传导方程进行的简易变换处理方法十分相似。接着，本章在下一节中将从金融的角度来阐释不稳定性，这暗示我们需要通过改变变量来改变方程式。9.2显式法不稳定性的金融学解释图9.2显式法（a）与隐式法（b）的求解BS方程的直观图示在显式法中，我们通过),(ttSSf,),(ttSf和),(ttSSf三个值得到),(tSf。这看起来与我们在7.4这部分（见图9.2（a））介绍过的三叉树方法有点类似。这样方便我们从显式法推导出新的方,ijf,1ijf1,1ijf（a）,1ijf1,1ijf,ijf（b）1,1ijf1,1ijf9法。如下【参考文献1中第18章】，假设在点(i,j)和点(i,j+1)对S求一阶导数和二阶导数是相等的：另外一种方法可以取得同样的效果，用1,jif代替9.2式右侧的jif,。这种方法引入的误差是有界的，并且这种方法的网格细化至趋于零。改变过后的替换后的差分方程变成了：也可写成（当i=1,2,…,M-1且j=1,2,…,N-1时）：同时有：这个方法仍然是显式的而且仍然受到不稳定性的影响。系数ˆia，ˆib以及ˆic本身也对此作出了很好的解释。回顾之前我们使用过的二叉树和三叉树方法，我们是在风险中性的假设下将后继节点对期权的预期价值进行贴现得到期权的价格。事实上，上述系数，包括)1/(1tr可以看做是长度为t的时间间隔的贴现系数。此外我们还有：10这也说明了上述系数等于概率乘以一个折现系数。他们是否是风险中性的概率呢？我们应当首先测试资产价格在时间间隔t内增长的值：这正是我们所期望的风险中性世界，对于其方差的增加，有：此外，当t极小时有：这也是服从风险中性假设的几何布朗运动的。因此我们可以认为排除一些小问题后，显式法是可以和三叉树方法取得同样的效果的。即“概率”d和0可能是负数。细心的读者可能会发现一个循环，在第五章中我们学习了一个线性组合的系数的稳定性条件。在热传导方程中，我们必须确定这个组合是一个凸组合，即系数是正的并且相加等于1。就像一个概率函数。有一种解决这个问题的方法，就是【参考文献1】中通过改变变量来实现的。通过Z=lnS来代换重写B-S方程，我们得到了实现稳定性的简单条件。然而，变量的改变对某些奇异期权来说却也未必是好事。在下一节中我们将使用隐式法来避免数值不稳定的问题。9.3使用全隐式法对欧式期权定价为了克服显式法带来的不稳定性问题，我们采用全隐式法。首对时间t采用向前差分求其偏导数。我们得到一个网格方程：也可写成（当i=1,2,…,M-1且j=1,2,…,N-1时）：11此时有：此时我们有需要使用3个未知量去计算出一个未知量（见图9.2b）。首先要注意的是，在每一个时间层，我们有M-1个方程式和N-1个未知量。边界条件决定两个未知量，终点条件也给出了最后时间层的值。正如显式法的案例那样，我们必须用回推法，解决当j=N-1,…,1情况下的所有线性方程组。在时间层j，得到的方程组如下：值得我们注意的是，上述的矩阵是个三对角矩阵，并且相对于时间层i来说它是一个不变的常量。因此我们可以使用LU因式分解来解决问题1。所有工作的MATLAB代码在图9.3中给出。1由于该矩阵具有稀疏结构，编写特定的代码来解决这类线性系统更加容易。此处我们只使用现成的MATLAB功能。12functionprice=EuPutImpl(S0,K,r,T,sigma,Smax,dS,dt)%建立网格并在必要时调整增量M=round(Smax/dS);dS=Smax/M;N=round(T/dt);dt=T/N;matval=zeros(M+1,N+1);vetS=linspace(0,Smax,M+1)';veti=0:M;vetj=0:N;%建立边界条件matval(