Octave：如何对这些 FOR 循环进行矢量化？

Question

我正在编写一个 Octave 脚本来计算欧式期权的价格。

第一部分使用蒙特卡罗模拟 n 个时间段内的标的资产价格。这会重复 nIter 次。

Octave 使设置初始矩阵变得非常容易。但是我还没有找到以向量化的方式完成任务，避免FOR循环的方法：

%% 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

请注意，n 中的价格等于 n-1 中的价格乘以一个因子，因此以下方法不起作用...

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

...因为它需要 S0 并将其乘以所有因子，产生与预期随机游走不同的结果。

Answer 1

由于迭代之间的依赖关系以获取每个新列相对于前一列的结果，看来您至少需要一个循环，但以矢量化方式在列中执行所有操作，这可能会加快它为你准备。两个嵌套循环的 vectorized 替换看起来像这样 -

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

我突然想到，可以使用cumprod 来处理依赖关系，这让我们得到cumulative product，这基本上是在这里完成的，因此会导致一个无循环的解决方案！这是实现 -

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

---

在 MATLAB 上进行基准测试

基准代码 -

%// 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

运行时结果 -

-------------------------------------- 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.

现在，运行时确实表明第一个提议的方法可能更适合这里！