显式有限差分法求欧式看涨期权的matlab代码

Question

用显式有限差分法求欧式看涨期权的matlab代码:
具体代码如下，我用blsprice()求black-scholes equation的准确值，但是发现用我代码算出的计算值和准确值不一样，在观察了每一个格点上的准确值和计算值后发现问题出在对边界值的处理上。
因为设定的x最大值（Xplus）是把它当做无穷大来考虑，然后用exp((k+1)x/2+(k+1)^2*t/4)来赋值无穷大时的边界值，但是实际上这个x的取值是有限大的，这样算出来的期权值就比准确值偏大很多，因为blsprice()在对边界值的处理上不会考虑x=xmax时是无穷大的情况，我就陷入困惑不知道该怎么算了。所以是我的计算值计算出错了还是准确值出错了呢？

E=10;r=0.05;vol=0.2;T=0.5;
Xplus=1;Xminus=-1;
N=51;M=21;
k=2*r/(vol*vol)
dtau=0.5*vol^2*T/(M-1)
dx=(Xplus-Xminus)/(N-1)
a=dtau/(dx)^2

u=zeros(M,N);              
%initial condition
for i=2:N-1 
    u(1,i)=max(exp(0.5*(k+1)*(Xminus+(i-1)*dx))-exp(0.5*(k-1)*(Xminus+(i-1)*dx)),0);      
end   
%boundary condition
for  j=1:M                 
     u(j,1)=0;
     u(j,N)=exp(0.5*(k+1)*(Xminus+(N-1)*dx)+0.25*(k+1)^2*(j-1)*dtau);
end        
%calculate other u values
for j=1:M-1                
    for i=2:N-1
        u(j+1,i)=((1-2*a)*u(j,i)+a*u(j,i+1)+a*u(j,i-1)); 
    end
end
%期权价格
vprice=ones(M,N)
for j=1:M
    for i=1:N
        vprice(j,i)=E*u(j,i)*exp(-0.5*(k-1)*(Xminus+(i-1)*dx)-0.25*(k+1)^2*(j-1)*dtau);
    end
end
vprice
%股票价格的向量
Shares=ones(N,1)
for n=0:(N-1)
    Shares(n+1)=E*exp(Xminus+n*dx);
end

%exact value
[call,put]=blsprice(Shares,E,r,T,vol)
%compare
plot(vprice(M,:))
hold
plot(call)

图片为计算值和准确值的比较

Answer 1

你好，同学！因为x最大值（Xplus）是把它当做无穷大来考虑，但是实际操作不可能，所以导致x处给出得边界条件本身就不准确。所以我们可以把边界设远一点，比如1.2倍或者1.5倍的Xplus处，但是程序处理的还是Xminus到Xplus的区间，这样计算的误差就会比较小。

E=10;r=0.05;vol=0.2;T=0.5;
Xplus=1; Xminus=-1;
alpha = 1.5; Xplus = alpha*Xplus; N0 = 51;
N=(N0-1)*alpha+1; M=21;
k=2*r/(vol*vol)
dtau=0.5*vol^2*T/(M-1)
dx=(Xplus-Xminus)/(N-1)
a=dtau/(dx)^2
u=zeros(M,N);              
%initial condition
for i=2:N-1 
    u(1,i)=max(exp(0.5*(k+1)*(Xminus+(i-1)*dx))-exp(0.5*(k-1)*(Xminus+(i-1)*dx)),0);      
end   
%boundary condition
for  j=1:M                 
     u(j,1)=0;
     u(j,N)=exp(0.5*(k+1)*(Xminus+(N-1)*dx)+0.25*(k+1)^2*(j-1)*dtau);
end        
%calculate other u values
for j=1:M-1                
    for i=2:N-1
        u(j+1,i)=((1-2*a)*u(j,i)+a*u(j,i+1)+a*u(j,i-1)); 
    end
end
%期权价格
vprice=ones(M,N0);
for j=1:M
    for i=1:N0
        vprice(j,i)=E*u(j,i)*exp(-0.5*(k-1)*(Xminus+(i-1)*dx)-0.25*(k+1)^2*(j-1)*dtau);
    end
end
vprice
%股票价格的向量
Shares=ones(N0,1);
for n=0:N0-1
    Shares(n+1)=E*exp(Xminus+n*dx);
end
%exact value
[call,put]=blsprice(Shares,E,r,T,vol);

%compare
plot(vprice(M,:),'r-')
hold
plot(call,'b--')
legend('selfdef','matlab')