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

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