第九章 期权定价的有限差分方法
在本章中，我们将给出几个简单的例子来说明基于偏微分方程（PDE）框架的期权定价方法。具体的方法的是利用第五章中讲述的有限差分方法来解决Black-scholes偏微分方程。在节中，我们会回顾衍生品定价的数值解法以及指出如何利用适当的边界条件来模拟一个特定的期权。在节中我们将会应用简单的显式（差分）方法来求解一个简单的欧式期权。正如你已熟知的那样，这种方法只能解出一些可以从金融角度来解释的不稳定的数值解。在节中我们将可以看到使用完全的隐式方法可以解决这种不稳定问题。在节中我们将介绍Crank-Nicolson方法在障碍期权定价中的应用，它可以看做是一种显式与完全隐式方法的混合。最后，在节中，我们会看到迭代松弛方法可以用于解决使用全隐式方法来解决美式期权定价时由于存在提前执行的可能性而导致的自由边界问题。
使用有限差分法解BS方程
在节中，我们给出了一个标的资产在时间的价格为的期权，该期权的价格是一个函数，且满足偏微分方程
（）
通过不同的边界条件可以让这个方程刻画不同的期权的特征。在某些地方可能因为假设的改变或者对路径依赖的改变而导致方程式的具体形式改变，但是此处仅仅作为一个起点，帮助读者了解如何应用基于有限差分方法来解决期权定价的问题。
正如我们在第五章中遇到的情况那样，要用有限差分方法来解偏微分方程，在此处我们必须建立资产价格和时间的离散网格。设T是期权的到期日，而Smax是一个足够大的资产价格，在我们所考虑的时间范围内，的数值不能超过Smax。设定Smax是因为偏微分方程的区域关于资产价格是无边界的。但是为了达到计算的目的，必须要求它是有界的。Smax相当于+∞。网格通过点(S,t)取得，其中(S,t)满足
，，，……,,
, ，2，……,。
本章中使用网格符号为,我们回顾一下（）方程式的几种不同解法：
向前差分
向后差分
中心（或对称）差分
对于第二个差分式子，有
至于究竟采用哪种方法进行离散化，我们将在后面的实际操作过程中对显式和隐式的方法作出详细的阐述说明。
另一个值得我们注意的问题是如何设置边界条件。对于执行价格为的欧式看涨期权，其终值条件为
对于执行价格为K的欧式看跌期权，其终值条件为
但是当我们涉及到资产价格的边界问题时，这个问题就变得复杂了。因为在数学上要求必须在一个有界区域里来解这个方程，而资产价格的区域是无界的。我们可以通过一些例子来说明这个问题。
例 我们首先考虑一支简单的欧式看跌期权。当资产的价格非常大时，看跌期权几乎是一文不值，因为我们几乎可以肯定的是其执行价低于现价：
Smax的值必须足够大才能使得假设的边界条件是合理的。当资产的价格时，鉴于资产的价格符合几何布朗运动的动力学模型，我们可以认为资产的价格会保持为0；因此到期收益为K，在时间贴现值为
其边界条件可以写为：
例 我们可以从推导出处理普通欧式看涨期权的方法。不妨设资产的价格,那么在任何时刻，我们都可以知道该期权价值为0：
而对于一个足够大的资产价格,我们可以肯定的是该期权一定为价内期权，即具有大于零的价值的期权，那么期满行权我们可以得到收益为。将终值贴现到时刻，并且考虑标的资产的价格是一个函数，那么我们根据无套利定价原理可以得到其边界条件如下：
具体条件可以写为：
对于资产价格较大时，需要一个替代的边界条件，但这时候必须要求期权的Δ为1；在这种情况下，我们所提出的边界条件一般针对未知函数的倒数，而非函数本身。这便是数学物理方程中所谓的Neumann边界条件。我们不会采取这种做法，因为这将使得数值解法复杂化了。
当我们面对障碍式期权时，事情可能就变得简单了。对敲出期权来说，例如一个向下触及失效看跌期权，在标的资产价格触碰到障碍时期权作废，价值为0。向上触及失效看涨期权的情况也是类似的，它们的优点都在于我们计算必须涉及到的边界即域是自然有界的（显然存在的）。美式期权的定价就显得更加复杂难以处理了，因为这涉及到了提前执行边界的问题。我们必须考虑什么水平的资产价格和在什么时间是期权的最佳实施点。因此，在我们解决问题的过程中必然产生一个自由边界。这就需要为不同的奇异期权设定不同的边界条件。如何计算出正确的边界条件以及使用数值方法逼近正确的结果，这依赖于对具体期权的选择。
利用显式法求解欧式期权
由于初次尝试解方程（），我们不妨考虑一个简单的欧式看跌期权。我们分别用中心差分对S进行求导，用向后差分对时间进行求导来逼近求解。这并不是唯一的方法，但是无论采取什么样的方法都必须保证边界条件是符合特定期权模型的。例中的边界条件差分后的结果可以用以下的方程式来表示：
（）
需要特别指出说明的是，由于存在终点条件的假设，因此方程必须适当地用回推法求解。不妨设（）式中有j=N，由于存在给定的终点条件，我们有一个未知量，可以认为它是关于三个已知量的函数。对其他的每个时间层（i）作同样的考虑。整理方程，我们可以得到：
（）
此时有
function price=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=*dt*(sigma^2*veti-r).*veti;
b=1-dt*(sigma^2*veti.^2+r);
c=*dt*(sigma^2*veti+r).*veti;
%求解方程
for j=N:-1:1
for i=2:M
matval(i,j)=a(i)*matval(i-1,j+1)+b(i)*matval(i,j+1)+...
c(i)*matval(i+1,j+1);
end
end
%返回价格，有可能在网格外线性插值生成
price=interp1(vetS,matval(:,1),S0);
图 简易欧式期权显式差分方法定价的MATLAB代码
定价不是直接利用MATLAB来完成的。图中是定价所使用的MATLAB的程序代码，另外计算还需要给出Smax的值以及上文提到的两个离散步骤。唯一一点需要注意的是，在数学符号中幂指数可以从0开始，但是在MATLAB的矩阵计算中，幂指数是从1开始的。另外，如果资产的初始价格不在网格上，那么我们只能够寻找两个相邻点进行插值。一般我们采取比较粗略的线性插值。但是如果有需要计算分析期权价格的敏感性（在实际案例中经常需要这样做）则需要进行更加复杂的插值计算。
>> [c,p] = blsprice(50,50,,5/12,);
>> P
P =
>> EuPutExpl(50,50,,5/12,0 .4,100,2,5/1200)
ans =
>> [c,p] = blsprice (50,50,,5/12, ) ;
>> P
P =
>> EuPutExpl(50,50,,5/12,,100,2,5/1200)
ans =
我们可以看到，数值方法能够提供比较精确的结果，为了提高精度，我们也可以采用更加精细的网格进行计算。
>> EuPutExpl(50,50,,5/12,,100,,5/1200)
ans =
>> EuPutExpl1(50,50,,5/12,,100,1,5/1200)
ans =+022
在此处我们看到在第五章中分析过的不稳定数值解的例子。有一种解决方案就是使用隐式方法进行求解。另外一种方法则是进行稳定性分析，然后在从离散化过程中导出边界。此处我们不会再使用第二种方案，这与第五章中对热传导方程进行的简易变换处理方法十分相似。接着，本章在下一节中将从金融的角度来阐释不稳定性，这暗示我们需要通过改变变量来改变方程式。
显式法不稳定性的金融学解释
图 显式法（a）与隐式法（b）的求解BS方程的直观图示
在显式法中，我们通过,和三个值得到。这看起来与我们在这部分（见图（a））介绍过的三叉树方法有点类似。这样方便我们从显式法推导出新的方法。如下【参考文献1中第18章】，假设在点(i,j)和点(i,j+1)对S求一阶导数和二阶导数是相等的：
另外一种方法可以取得同样的效果，用代替式右侧的。这种方法引入的误差是有界的，并且这种方法的网格细化至趋于零。
改变过后的替换后的差分方程变成了：
也可写成（当i=1,2,…,M-1且j=1,2,…,N-1时）：
同时有：
这个方法仍然是显式的而且仍然受到不稳定性的影响。系数，以及本身也对此作出了很好的解释。回顾之前我们使用过的二叉树和三叉树方法，我们是在风险中性的假设下将后继节点对期权的预期价值进行贴现得到期权的价格。事实上，上述系数，包括可以看做是长度为的时间间隔的贴现系数。
此外我们还有：
这也说明了上述系数等于概率乘以一个折现系数。他们是否是风险中性的概率呢？我们应当首先测试资产价格在时间间隔内增长的值：
这正是我们所期望的风险中性世界，对于其方差的增加，有：
此外，当极小时有：
这也是服从风险中性假设的几何布朗运动的。因此我们可以认为排除一些小问题后，显式法是可以和三叉树方法取得同样的效果的。即“概率” 和可能是负数。细心的读者可能会发现一个循环，在第五章中我们学习了一个线性组合的系数的稳定性条件。在热传导方程中，我们必须确定这个组合是一个凸组合，即系数是正的并且相加等于1。就像一个概率函数。
有一种解决这个问题的方法，就是【参考文献1】中通过改变变量来实现的。通过Z=lnS来代换重写B-S方程，我们得到了实现稳定性的简单条件。然而，变量的改变对某些奇异期权来说却也未必是好事。在下一节中我们将使用隐式法来避免数值不稳定的问题。
使用全隐式法对欧式期权定价
为了克服显式法带来的不稳定性问题，我们采用全隐式法。首对时间t采用向前差分求其偏导数。我们得到一个网格方程：
也可写成（当i=1,2,…,M-1且j=1,2,…,N-1时）：
此时有：
此时我们有需要使用3个未知量去计算出一个未知量（见图）。首先要注意的是，在每一个时间层，我们有M-1个方程式和N-1个未知量。边界条件决定两个未知量，终点条件也给出了最后时间层的值。正如显式法的案例那样，我们必须用回推法，解决当j=N-1,…,1情况下的所有线性方程组。在时间层j，得到的方程组如下：
值得我们注意的是，上述的矩阵是个三对角矩阵，并且相对于时间层i来说它是一个不变的常量。因此我们可以使用LU因式分解来解决问题 。所有工作的MATLAB代码在图中给出。
function price = 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(:,N+1) = max(K-vetS,0);
matval(1,:) = K*exp(-r*dt*(N-vetj));
matval(M+1,:) = 0;
% 建立三对角矩阵系数
a = *(r*dt*veti-sigma^2*dt*(veti.^2));
b = 1+sigma^2*dt*(veti.^2)+r*dt;
c = *(r*dt*veti+sigma^2*dt*(veti.^2));
coeff = diag(a(3:M),-1) + diag(b(2:M)) + diag(c(2:M-1),1);
[L,U] = lu(coeff);
% 求解线性方程组
aux = zeros(M-1,1);
for j=N:-1:1
aux(1) = - a(2) * matval(1,j); % other term from BC is zero
matval(2:M,j) = U \ (L \ (matval(2:M,j+1) + aux));
end
% 返回价格，如果初始价格在网格外则可能进行性线性插值
price = interp1(vetS, matval(:,1), S0);
图 全隐式方法求解简易欧式期权的MATLAB代码
>> [c,p] = blsprice(50,50,,5/12,);
>> P
P =
>> EuPutlmpl(50,50,,5/12,,100,,5/2400)
ans =
这个计算结果是非常精确的，这是因为使用了排除不稳定因素存在的更加精细的网格的结果。另外一种方法就是通过Crank-Nicolson方法来提高精度，这将在后面的章节中做进一步的介绍。
使用Crank-Nicolson方法对障碍式期权定价
Crank-Nicolson方法已经在节中提到过，它通过结合使用显式法和隐式法来提高数值解的精确度。在B-S模型的求解中应用此方法可以得到：
上述方程可以改写为：
此时有：
我们首先考虑在中介绍过的向下触及失效期权，假设该障碍期权是处于连续监测下的。在这种情况下我们需要考虑≤S≤Smax，边界条件是：
综合考虑上述边界条件，我们可以将方程（）改写为如下矩阵形式：
同时，有：
相应的MATLAB代码在图中列出。
这个结果可以同中通过定价分析公式得到的结果进行比较。
>> DOPut (50,50,,5/12,,40)
ans =
>> DOPutCK(50,50,,5/12,,40,100,,1/1200)
ans =
障碍期权还有很多种，可以在【参考文献9】中找到更多的关于障碍式期权定价的偏微分方程应用。
function price = DOPutCK(S0,K,r,T,sigma,Sb,Smax,dS,dt)
% 建立网格并且在必要的时候调整增量
M = round((Smax-Sb)/dS);
dS = (Smax-Sb)/M;
N = round(T/dt);
dt = T/N;
matval = zeros(M+1,N+1);
vetS = linspace(Sb,Smax,M+1)';
veti = vetS / dS;
vetj = 0:N;
% 建立边界条件
matval(:,N+1) = max(K-vetS,0);
matval(1,:) = 0;
matval(M+1,:) = 0;
% 建立系数矩阵
alpha = *dt*( sigma^2*(veti.^2) - r*veti );
beta = -dt**( sigma^2*(veti.^2) + r );
gamma = *dt*( sigma^2*(veti.^2) + r*veti );
M1 = -diag(alpha(3:M),-1) + diag(1-beta(2:M)) - diag(gamma(2:M-1),1);
[L,U] = lu(M1);
M2 = diag(alpha(3:M),-1) + diag(1+beta(2:M)) + diag(gamma(2:M-1),1);
% 求解线性方程组
for j=N:-1:1
matval(2:M,j) = U \ (L \ (M2*matval(2:M,j+1)));
end
% 返回价格，这个结果可能因为初始资产价格在网格外而有线性插值产生。
price = interp1(vetS, matval(:,1), S0);
图 Crank-Nicolson方法对向下触及失效期权定价的MATLAB代码。
美式期权定价
有限差分方法在对简单的欧式期权进行定价时的确是十分方便的，但是这在美式期权定价中却显得并不是很实用。我们将这个方法应用于美式期权定价的时候，精确的公式就不大适用了。这是因为美式期权定价时必须考虑存在提前行权的可能性。根据无套利原理，在(S,t)空间中任一点的期权的价格不能低于其内在价值（即立刻行权的所得收益）。对于一个一般的美式期权，这意味着：
从一个实际操作来看，在显式法中要充分考虑这种情况是很困难的。但是我们可以做一些修改以便能够简单的应用节中提到的方法。在计算得到后，我们可以检测提前行权的概率。并且可以像我们处理二叉树模型那样设：
由于存在着不稳定性问题，我们可能更加倾向于使用隐式法。但是这样处理后又导致了另一个复杂的问题，即根据上述关系要求必须是已知的，而在隐式法中却并不能已知要解决这个问题，我们可以采取一个迭代的方法来求解这个线性方程组，而不是基于LU分解直接求解。在节中我们介绍了Gauss-Seidel松弛迭代算法。为了方便我们先回顾一下这个方法。给出一个线性方程组，如下：
我们将使用下面的迭代方法，初始点为:
其中k是迭代计数器，而ω是松弛变量。该迭代将会一直进行到满足收敛条件，如：
其中ε是容限参数。
现在我们将尝试使用Crank-Nicolson方法对美式期权进行定价。我们要解决的是和（）式大致相同的一个线性方程组，但是这里的边界条件有一点不同，价值为0的期权不用再考虑是否有障碍。我们要用回推法求解的线性方程组如下：
其中右边有：
式子中额外增加的一项是为了将看跌期权独有的边界条件带入模型方程。在使用松弛方法的时候应当考虑矩阵的三对角性质，并且这能够适应美式期权提前执行的特征。设为内在价值，i=1，2，…，M-1，同时。对于每个时间层j，我们进行如下迭代：
当计算从一个时间层进行到下一个时间层的时候，将上一个时间层的计算结果作为下一个时间层迭代计算的初始向量可能是更加合理的。完成上述步骤需要的代码在图中将列出。使用该代码有点麻烦，原因在于MATLAB的向量幂运算是从1开始的，而现实却要求从0开始。这样我们便不能建立一个包含所有的向量矩阵，同时系数矩阵也无法储存。上述迭代最好是直接通过、和计算得出。
该代码的实际效果可能要比MATLAB金融工具箱中附带的使用二叉树方法对美式期权定价的binprice函数（见节）要好。
>> tic, [pr, opt] = binprice(50,50 ,,5/12,1/1200, ,0) ; , toc
Elapsed time is seconds.
>> opt(1,l)
ans =
>> tic ,AmPutCK( 50,50,,5/12,,100,1,1/600,,),toc
ans =
Elapsed time is seconds.
>> tic,AmPutCK(50,50,,5/12,,100,1,1/600,,),toc
ans =
Elapsed time is seconds.
>> tic,AmPutCK(50,50,,5/12,,100,1,1/600,,),toc
ans =
Elapsed time is seconds.
>> tic,AmPutCK(50,50,,5/12,,100,1,1/1200,,),toc
ans =
Elapsed time is seconds.
>> tic,AmPutCK(50,50,,5/12,,100,1,1/100,,),toc
ans =
Elapsed time is seconds.
function price = AmPutCK(S0,K,r,T,sigma,Smax,dS,dt,omega,tol)
M = round(Smax/dS); dS = Smax/M; % 建立网格
N = round(T/dt); dt = T/N;
oldval = zeros(M-1,1); %Gauss-Seidel更新向量
newval = zeros(M-1,1);
vetS = linspace(0,Smax,M+1)';
veti = 0:M; vetj = 0:N;
% 建立边界条件
payoff = max(K-vetS(2:M),0);
pastval = payoff; % values for the last layer
boundval = K*exp(-r*dt*(N-vetj)); % 边界值
% 建立系数矩阵和式子右边矩阵
alpha = *dt*( sigma^2*(veti.^2) - r*veti );
beta = -dt**( sigma^2*(veti.^2) + r );
gamma = *dt*( sigma^2*(veti.^2) + r*veti );
M2 = diag(alpha(3:M),-1) + diag(1+beta(2:M)) + diag(gamma(2:M-1),1);
% 使用SOR方法求解线性方程组
aux = zeros(M-1,1);
for j=N:-1:1
aux(1) = alpha(2) * (boundval(1,j) + boundval(1,j+1));
% 建立右端矩阵并进行初始化
rhs = M2*pastval(:) + aux;
oldval = pastval;
error = realmax;
while tol < error
newval(1) = max ( payoff(1), ...
oldval(1) + omega/(1-beta(2)) * (...
rhs(1) - (1-beta(2))*oldval(1) + gamma(2)*oldval(2)));
for k=2:M-2
newval(k) = max ( payoff(k), ...
oldval(k) + omega/(1-beta(k+1)) * (...
rhs(k) + alpha(k+1)*newval(k-1) - ...
(1-beta(k+1))*oldval(k) + gamma(k+1)*oldval(k+1)));
end
newval(M-1) = max( payoff(M-1),...
oldval(M-1) + omega/(1-beta(M)) * (...
rhs(M-1) + alpha(M)*newval(M-2) - ...
(1-beta(M))*oldval(M-1)));
error = norm(newval - oldval);
oldval = newval;
end
pastval = newval;
end
newval = [boundval(1) ; newval ; 0]; % 加入缺少的值
% 返回价格，这个价格可能因为初始资产价格在网格外而由线性插值生成。
price = interp1(vetS, newval, S0);
图 Crank-Nicolson方法求解美式看跌期权MATLAB代码
通过以上例子我们可以看到，松弛变量ω在迭代方法中对结果的收敛性有十分显著的影响。就运算速度而言，上述有限差分方法（CK）似乎要比二叉树方法快，但是我们必须要注意一些问题。当我们将其与隐式法比较时，两者是不相上下的。而且CPU对MATLAB代码的解码方式也可能影响到计算效率。不管怎样，一个足够大而且精确的网格要比二叉树的点能够提供更加有效的敏感度估计（这在求解B-S方程中会用到）。此外，在处理一些奇异期权的定价问题的时候，上述提到的有限差分方法可能是个更好的选择。
扩展阅读
·许多金融工程应用偏微分方程的例子可以在【参考文献6】和【参考文献7】中找到，其中也包括了应用有限差分方法。另外也可以在【参考文献2】中找到其它许多有用的例子。
·我们对B-S方程直接采用了有限差分方法求解，同时，一个变量的微小调整可能有助于分析结果的稳定性。这可以参考【参考文献3】中的某些例子。在这个参考文献中你还可以找到一个有限差分方法的改进方案，这是一个比简单有限差分方法更加精确的方法。
·【参考文献7】和【参考文献8】是专门介绍金融工程中的有限差分方法的书籍。
·如果你对有限元方法感兴趣，也可以阅读【参考文献5】。
参考文献
【1】 . Hull. Options, Futures, and Other Derivatives (5th ed.). Prentice Hall, Upper Saddle River, NJ, 2003.
【2】 . Kwok. Mathematical Models of Financial Derivatives. Springer Verlag, Berlin, 1998.
【3】. R. Seydel. Tools f o r Computational Finance. Springer-Verlag, Berlin,2002.
【4】. D. Tavella and C. Randall. Pricing Financial Instruments: The Finite Difference Method. Wiley, New York, 2000.
【5】 J. Topper. Financial Engineering with Finite Elements. Wiley, New York,2005.
【6】. P. Wilmott. Derivatives: The Theory and Practice of Financial Engineering. Wiley, Chichester, West Sussex, England, 1999.
【7】 P. Wilmott. Quantitative Finance (vols. I and II). Wiley, Chichester,West Sussex, England, 2000.
【8】 Y.-I. Zhu, X. Wu, and I.-L. Chern. Derivative Securities and Difference Methods. Springer, New York, 2004.
【9】. R. Zvan, . Vetzal, and . Forsyth. PDE Methods for Pricing Barrier Options. Journal of Economic Dynamics and Control, 24: 1563-1590,2000.
由于该矩阵具有稀疏结构，编写特定的代码来解决这类线性系统更加容易。此处我们只使用现成的MATLAB功能。
PAGE
PAGE 24
