Octave: how can these FOR loops be vectorized?

Question

I am writing an Octave script to calculate the price of an European option.

The first part uses Monte Carlo to simulate the underlying asset price over n number of time periods. This is repeated nIter number of times.

Octave makes it very easy to setup initial matrices. But I haven't found the way to complete the task in a vectorized way, avoiding FOR loops:

%% Octave simplifies creation of 'e', 'dlns', and 'Prices'
e = norminv(rand(nIter,n));
dlns = cat(2, ones(nIter,1), exp((adj_r+0.5*sigma^2)*dt+sigma*e.*sqrt(dt)));
Prices = zeros(nIter, n+1);

for i = 1:nIter              % IS THERE A WAY TO VECTORIZE THESE FOR LOOPS?
  for j = 1:n+1
    if j == 1
      Prices(i,j)=S0;
    else
      Prices(i,j)=Prices(i,j-1)*dlns(i,j);
    end
  endfor
endfor

Note that the price in n is equal to price in n-1 times a factor, hence the following does not work...

Prices(i,:) = S0 * dlns(i,:)

...since it takes S0 and multiplies it by all the factors, yielding different results than the expected random walk.

Answer 1

Because of the dependency between iterations to obtain results for each new column with respect to the previous column, it seems you would need at least one loop there, but do all operations within a column in a vectorized fashion and that might speed it up for you. The vectorized replacement for the two nested loops would look something like this -

Prices(:,1)=S0;
for j = 2:n+1
    Prices(:,j) = Prices(:,j-1).*dlns(:,j);
endfor

It just occurred to me that the dependency can be taken care of with cumprod that gets us cumulative product which is essentially being done here and thus would lead to a no-loop solution! Here's the implementation -

Prices = [repmat(S0,nIter,1) cumprod(dlns(:,2:end),2)*S0]

---

Benchmarking on MATLAB

Benchmarking Code -

%// Parameters as told by OP and then create the inputs
nIter= 100000;
n = 100;
adj_r = 0.03;
sigma = 0.2;
dt = 1/n;
S0 = 60;
e = norminv(rand(nIter,n));
dlns = cat(2, ones(nIter,1), exp((adj_r+0.5*sigma^2)*dt+sigma*e.*sqrt(dt)));

disp('-------------------------------------- With Original Approach')
tic
Prices = zeros(nIter, n+1);
for i = 1:nIter
    for j = 1:n+1
        if j == 1
            Prices(i,j)=S0;
        else
            Prices(i,j)=Prices(i,j-1)*dlns(i,j);
        end
    end
end
toc, clear Prices

disp('-------------------------------------- With Proposed Approach - I')
tic
Prices2(nIter, n+1)=0; %// faster pre-allocation scheme
Prices2(:,1)=S0;
for j = 2:n+1
    Prices2(:,j)=Prices2(:,j-1).*dlns(:,j);
end
toc, clear Prices2

disp('-------------------------------------- With Proposed Approach - II')
tic
Prices3 = [repmat(S0,nIter,1) cumprod(dlns(:,2:end),2)*S0];
toc, clear Prices3

Runtimes results -

-------------------------------------- With Original Approach
Elapsed time is 0.259054 seconds.
-------------------------------------- With Proposed Approach - I
Elapsed time is 0.020566 seconds.
-------------------------------------- With Proposed Approach - II
Elapsed time is 0.067292 seconds.

Now, the runtimes do suggest that the first proposed approach might be a better fit here!