- #1
tomatoyeungHK
- 4
- 0
Dear all:
I am studying the finite difference method for black scholes equation(heat equation form), and my Prof. asked me to use the Crank Nicolson and Compact Scheme to do it. But I cannot reach second order for the space, i.e., x . The error only reach 1st order and sometimes not stable, it is totally contradict with the unconditionally stable property of Crank Nicolson Method. Is there anyone can give e some help in solving this problem? Deadline is coming and I cannot find out any problems. Thanks a lot. My Matlab code is here:
I am studying the finite difference method for black scholes equation(heat equation form), and my Prof. asked me to use the Crank Nicolson and Compact Scheme to do it. But I cannot reach second order for the space, i.e., x . The error only reach 1st order and sometimes not stable, it is totally contradict with the unconditionally stable property of Crank Nicolson Method. Is there anyone can give e some help in solving this problem? Deadline is coming and I cannot find out any problems. Thanks a lot. My Matlab code is here:
Code:
clear
%the length of the interval , set: x_min = -1.4, x_max = 5;
%the length of the interval is 6.4
%Parameters
K = 100; sigma = 0.25; r = 0.05;
beta = sigma*sigma/2;
%Discretization of x and t
step_x = 640;
step_t = 640;
dx = 6.4 /(step_x);
%the T =0.25
dt = 0.25/(step_t);
%Creating Vectors
x = linspace(-1.4,5,(step_x-1));
t = linspace(0,0.25,step_t );
U = zeros(step_x -1, step_t);
%Initial value and Boundary Value
u0_C= K*(max(exp(x) -1,0));
%Setting the left hand matrix, i.e. coeff. of k+1.th layer vectors
%upper diagonal
upp = (-(r-beta)/(4*dx) - beta/(2*dx*dx));
%main diagonal
main = (beta/(dx*dx) + r/2 + 1/dt);
%lower diagonal
low = ((r-beta)/(4*dx) - beta/(2*dx*dx));
%main diagonal for right hand side,i.e. coeff. of k.th layer vectors
main_R = (beta/(dx*dx) + r/2 - 1/dt);
%Set up left hand matrix
l1 = [main, low, zeros(1,step_x -3)];
l2 = [main, upp, zeros(1,step_x -3)];
L_Mat = toeplitz(l1,l2);
%Set up right hand matrix
r1 = [-main_R, -low, zeros(1,step_x -3)];
r2 = [-main_R, -upp, zeros(1,step_x -3)];
R_Mat = toeplitz(r1,r2);
v = u0_C;
v = v';
for i = 1:step_t
%boundary conditions for call option
%Left points
uL0 = K*(max(exp(-1.4)-exp(-r*(i-1)*dt),0));
uL1 = K*(max(exp(-1.4)-exp(-r*(i*dt)),0));
%Right points
uR0 = K*(max(exp(5)-exp(-r*(i-1)*dt),0));
uR1 = K*(max(exp(5)-exp(-r*(i*dt)),0));
rhs = R_Mat*v;
%Updating the boundary conditions
rhs(1) = rhs(1) - low*uL0 - upp*uL1;
rhs(end) = rhs(end)-upp*uR0 - low*uR1;
v = L_Mat\rhs;
U(:,i)=v;
end
%U(:,i)
S=K*exp(x);
%U
[Call,Put] = blsprice(S,K,0.05,0.25,0.25);
plot(S,Call);
hold on
sol = U(:,end);
plot(S,sol,'r');
%finding the error at x=0,i.e. the error at exercise price K,
[x0,I0]=min(abs(x));
error_0=abs(sol(I0)-Call(I0));
error_0
%Norm_2 for the average error
error_2 = sqrt(6.4/step_x)*norm(sol'-Call,2);
error_2