## 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
```