Finite Differences Method – Matlab Code

In this Matlab code we implement the following finite differences methods: implicit, explicit and Crank-Nicolson to estimate option values. We also plot our concluded values vs. the ones generated by the Black-Scholes closed-form solution.


function FiniteDifferences()

S0 = 10; sigma = 0.20;r = 0.04; dT = 0.002;X0 = log(S0);T =0.5
delX = [sigma*sqrt(dT),sigma*sqrt(3*dT),sigma*sqrt(4*dT)];
figure;
for i =1:3
dX = delX(i);
Path = ceil(log(18/2) / dX) ;

Time = (T/dT)+1;

X = zeros(Path,1);
P = zeros(Path,Time);

X(1) = log(2);

P(1,Time) = 10 - exp(X(1));
for j =2:Path
X(j) = X(j-1) + dX;
P(j,Time) = max(10 - exp(X(j)),0);
end
% Explicit different method
[p_d, p_m, p_u] = compute_explicit_prob(r,sigma,dX,dT);
A = GetExplicitA(p_u, p_m, p_d,Path);

PutValue = zeros(Path,3);

% Main Block<span></span>
for j = Time-1:-1:1
P(2:Path-1,j) = A*P(:,j+1);
P(1,j) = P(2,j) -(exp(X(1)) - exp(X(2)));
P(Path,j) = P(Path-1,j) ;
end
Y = zeros(size(X,1),1);
for m =1:size(X,1)
[~, Y(m)] = blsprice(exp(X(m)),S0,0.04,0.5,0.20,0);
end
subplot(3,3,(i-1)*3+1);
plot(exp(X),P(:,1),exp(X),Y);
legend('Explicit Difference Price','BS Price');

P(:,1:Time-1) = zeros(Path,Time-1);
[p_d, p_m, p_u] = compute_implicit_prob(r,sigma,dX,dT);
A_I = GetImplicitA(p_u, p_m, p_d,Path);
for j = Time-1:-1:1
b = [exp(X(2)) - exp(X(1)); P(2:Path-1,j+1);0 ];
P(:,j) = A_I\b;
end

subplot(3,3,(i-1)*3+2);
plot(exp(X),P(:,1),exp(X),Y);
legend('Implicit Difference Price','BS Price');

P(:,1:Time-1) = zeros(Path,Time-1);
[p_d, p_m, p_u] = compute_crank_nicol_prob(r,sigma,dX,dT);
Z = zeros(Path,1);
A_I = GetImplicitA(p_u, p_m, p_d,Path);
for j = Time-1:-1:1

Z(2:Path-1) = GetCrankZ(p_u, p_m, p_d,Path,P(:,j+1));
Z(1) = -(exp(X(1)) - exp(X(2)));
Z(Path) =Z(Path-1) ;
P(:,j) = A_I\Z;
end

subplot(3,3,(i-1)*3+3);
plot(exp(X),P(:,1),exp(X),Y);
legend('Crank Nicolson Price','BS Price');

end

end

function [p_u, p_m, p_d] = compute_explicit_prob(r,sigma,dX,dT)
t1 = (sigma^2)/(2*(dX^2));
t2 = (r - ((sigma^2)/2))/(2*dX);
t3 = dT * ((sigma^2)/(dX^2));
t4 = r*dT;
p_u = dT*(t1 +t2);
p_m = 1 - t3 -t4;
p_d = dT*(t1 -t2);
end

function [p_u, p_m, p_d] = compute_implicit_prob(r,sigma,dX,dT)
v = r - ((sigma^2)/2);
t1 = (sigma^2)/((dX^2));
t2 = (v)/(dX);
t3 = dT * ((sigma^2)/(dX^2));
t4 = r*dT;
p_u = -0.5* dT*(t1 +t2);
p_m = 1 + t3 + t4;
p_d = -0.5* dT*(t1 -t2);
end

function [p_u, p_m, p_d] = compute_crank_nicol_prob(r,sigma,dX,dT)
v = r - ((sigma^2)/2);
t1 = (sigma^2)/((dX^2));
t2 = (v)/(dX);
t3 = dT * ((sigma^2)/(2*dX^2));
t4 = r*dT*0.5;
p_u = -0.25* dT*(t1 +t2);
p_m = 1 + t3 + t4;
p_d = -0.25* dT*(t1 -t2);
end

function [A] = GetExplicitA(p_u, p_m ,p_d, N)
A = zeros(N-2,N);
A(1,1:3)= [p_u p_m p_d];
for i = 2:N-2
A(i,:) = circshift(A(i-1,:),[0,1]);
end

end
function [A] = GetImplicitA(p_u, p_m ,p_d, N)
A = zeros(N,N);
A(1,1:2) = [ 1 -1];
A(2,1:3)= [p_u p_m p_d];
for i = 3:N-1
A(i,:) = circshift(A(i-1,:),[0,1]);
end
A(N,N-1:N) = [ 1 -1];
end

function [A] = GetCrankZ(p_u, p_m ,p_d, N,P)
A = zeros(N-2,N);
A(1,1:3)= [-p_u -(p_m-2) -p_d];
for i = 2:N-2
A(i,:) = circshift(A(i-1,:),[0,1]);
end
A = A*P;
end

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s