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

Least Squares Monte Carlo – Matlab Code

In this Matlab code, we value American options using the Least Squared Monte Carlo method developed by Longstaff and Schwartz. Three different polynomials are used in this example to solve for the future continuation value. Finally, antithetic variates are also used to provide a faster convergence.

% inputs
path = 100000;
timestep = 20;
S0 = 40;
r = 0.06;
sigma =0.20;
K = 40;

method = ['Hermite ';'Laguerre ';'Monomials'];
cellmethod = cellstr(method);
W = zeros(path,timestep);
S = zeros(path,timestep);
PutValue = zeros(27,4);
Display = zeros(81,6);
counter =1;
for methodi=1:3
name = cellmethod{methodi};
count = 1;

zero_start = zeros(path,1);
for terms =2:4
for K = 36:4:44
for T = [0.5 1 2]
dt = T /(timestep-1);

drift = r - ((sigma*sigma)/2);
temp = sqrt(dt)* randn(path/2,timestep-1);
dW = [temp;-temp];
W = cumsum(dW,2);
W = [zero_start W];
Time = (0:dt:T);
Time = repmat(Time,path,1);
drift_t = (drift*Time);
vol = sigma*W;
S = S0*exp(drift_t+vol);
Y = max(K - S(:,timestep),0);
for t=timestep-1:-1:2
Y = Y*exp(-r*dt);
X = S(:,t);
ECY = GetECV(name,terms,X,Y);
EV = max(K - S(:,t),0);
YIndex = (EV > ECY) & (ECY > 0);
Y(YIndex) = EV(YIndex);
end

PutValue(count,methodi) = mean(exp(-r*dt)*Y);
Display(counter,1) = methodi;
Display(counter,2) = terms;
Display(counter,3) = K;
Display(counter,4) = T;
Display(counter,5) = PutValue(count,methodi);
[~,val]= binprice(40, K, 0.06, T, T/20, 0.20, 0, 0, 0, 0);
Display(counter,6) = val(1,1);
counter = counter+1;
count= count+1;
end
end
end
end
count =1;
% price options for different strike prices and time frames
for K = 36:4:44
for T = [0.5 1 2]

[~,option_value] = binprice(40, K, 0.06, T, T/20, 0.20, 0, 0, 0, 0);
PutValue(count,4) = option_value(1,1);
count = count+1;
end
end

PutValue(:,4)=[PutValue(1:count-1,4);PutValue(1:count-1,4);PutValue(1:count-1,4)];

% functions

function [ ECY ] = GetECV( method,terms,X,Y )
switch method
case 'Laguerre'
L = getLaguerreL(terms,X);
case 'Hermite'
L = getHermite(terms,X);
case 'Monomials'
L = getMonomials(terms,X);
end
A = L' * L;
b = Y'*L;
b = b';
a = A\b;
ECY = L*a;
end

function [L] = getLaguerreL(terms,X)
m = size(X,1);
terms = terms+1;
L = zeros(m,terms);
L(:,1) = exp(-X/2);
L(:,2) = (1-X).*exp(-X/2);
if(terms > 2)
L(:,3) = (1-(2.*X) + ((X.*X) /2)).*exp(-X/2);
end
if(terms > 3)
L(:,4) = (1-(3.*X) + (3*(X.*X) /2) - ((X.*X.*X)/6)).*exp(-X/2);
L(:,5) = (1-(4.*X) + (3*(X.*X)) - ((2*X.^3)/3) + ((X.^4)/24)).*exp(-X/2);
end
end

function [L] = getHermite(terms,X)
m = size(X,1);
L = zeros(m,terms);
L(:,1) = 1;
L(:,2) = 2*X;
for i=3:terms
L(:,i) = (2*(L(:,i-1).*X)) - (2*(i-1).*L(:,i-2));
end
end

function [L] = getMonomials(terms,X)
m = size(X,1);
L = zeros(m,terms);
L(:,1) = 1;
for i=2:terms
L(:,i) = X.*L(:,i-1);
end
end

Building Systematic Trading Strategies – Part I

This series of posts is to give a broad overview about building quantitative strategies. This post is not about the intricate details on how to build an infrastructure to simulate such strategies or how to implement technically, but it should just serve as a primer for all those aspiring to build quantitative portfolios.

The Idea

Every strategy starts off with an idea. Unlike a qualitative assessment (say, Apple’s stock is going to increase in the next month) that involves an analyst’s research on primarily following a few stocks, a quantitative strategy is more involved and needs to subscribe to the following rules.

  1. It should be a general rule that can be applied to broad set of stocks (what we call, Universe)
  2. The strategy should not suffer from forward bias. In simple words, it means you cannot use future data to predict future data. The strategy can however use all data up to that time, to predict future directions.
  3. One should be able to express the rule in a mathematical form.
  4. The output of the rule should be a set of dollar amounts that can be assigned to each stock in the Universe

Say we have a simple idea that stocks that have gone up recently will come down and vice-versa. This behavior pattern is what is popularly called Mean reversion anomaly.  Let us how we “quantify” this idea.

Quantification of the idea

There are different ways on how we can “quantify” this idea. There is no such correct way in building an idea. Most of it involves trial and error,  and to see how different factors impact an idea.

Let us start with the mean reversion idea.  One way to state it mathematically is below

Amount Invested for Stock I and Day D  (weight[i,d] )=  – 5 day returns of Stock I on Day D ; for all Top 500 stocks by volume (returns5[i,d]).

In this expression we quantify “recent returns” as  the last 5 day returns. The universe is the top 500 stocks traded by volume in  US equity markets.

weight[i,d] = -returns5[i,d]

One immediate concern with this strategy is that we observe, the sum of dollar invested will not equal our notional value that we plan to invest. So to normalize to the notional we can do the following operation.

weight[i,d] = (weight[i,d]*notional value)/ (absolute sum of all weight[i,d])

That is it! This is the hello world version of a quantitative strategy. Lets see how it performs

Following up lets look at ways to improve the idea and make it more tradable!