Just started to code in matlab because I am studying the book An Introduction to Financial Option Valuation by Higham. One of the example codeblocks he gives (this is the source) is:
V = zeros(M,1);
Vanti = zeros(M,1);
for i = 1:M
samples = randn(N,1);
% standard Monte Carlo
Svals = S*cumprod(exp((r-0.5*sigma^2)*Dt+sigma*sqrt(Dt)*samples));
Smax = max(Svals);
if Smax < B
V(i) = exp(-r*T)*max(Svals(end)-E,0);
end
% antithetic path
Svals2 = S*cumprod(exp((r-0.5*sigma^2)*Dt-sigma*sqrt(Dt)*samples));
Smax2 = max(Svals2);
V2 = 0;
if Smax2 < B
V2 = exp(-r*T)*max(Svals2(end)-E,0);
end
Vanti(i) = 0.5*(V(i) + V2);
end
I am trying to get this loop more efficient, so I am trying to remove the for loop. This is what I wrote so far:
V = zeros(M,1);
Vanti = zeros(M,1);
samples = randn(N,M);
Svals = S*cumprod(exp((r-0.5*sigma^2)*Dt+sigma*sqrt(Dt)*samples));
Svals2 = S*cumprod(exp((r-0.5*sigma^2)*Dt-sigma*sqrt(Dt)*samples));
Send = Svals(end,:);
Send2 = Svals2(end,:);
Smax = max(Svals);
Smax2 = max(Svals2);
V2 = zeros(M,1);
for i = 1:M
if Smax(i) < B
V(i) = exp(-r*T)*max(Send(i)-E,0);
end
if Smax2(i) < B
V2(i) = exp(-r*T)*max(Send2(i)-E,0);
end
end
Vanti = 0.5*(V + V2);
aM = mean(V); bM = std(V);
conf = [aM - 1.96*bM/sqrt(M), aM + 1.96*bM/sqrt(M)]
aManti = mean(Vanti); bManti = std(Vanti);
confanti = [aManti - 1.96*bManti/sqrt(M), aManti + 1.96*bManti/sqrt(M)]
toc
This already made the code significantly quicker, because there aren't any randn variables generated inside the loop. I don't know however, how I am able to remove the other part of the loop. Is it even possible?
samples
the other way around: accessing a contiguous column is faster than accessing a contiguous row (since array storage in matlab is column-major). Also: you can try usingprofile
to see which part of your code takes the most amount of time (though it's not always reliable, especially for fast programs).Svals
is scalar, save forsamples
, you can calculate the terms pre-loop. Like:alpha1 = cumprod(exp((r-0.5*sigma^2)*Dt+sigma*sqrt(Dt)*samples),2)
and thenSvals=S*alpha1(i,:);
.