请教循环向量化

喜新不厌旧

偶有个程序，循环比较大，运行不出来，有人建议循环改成MEX方式或向量化，可惜偶不会。
希望这里的高手能给予指教，小女子多谢啦
QQ 5795687
zhe9264@hotmail.com

function y=MRJ(dP,dt,lambda,sigma,eta,Pm,dPhi)
% Real Options, Mean-Reversion with Jumps
% Solution of the Finite Difference Problem
format compact
filename = input('Enter test number: ', 's');
%%%%%%%%%%%%
% Real Option Parameters :
rho = 0.10; % Exogeneous discount rate
lambda = 0.15; % Annual frequency of jumps
sigma = 0.22; % Volatility
eta = 0.03; % Reversion speed
Pbar = 20; % Average oil price (US$/bbl)
T = 10; % Lease period (Expiry) (in years)
q = 1/3; % Economic quality of the reserve
% Finite Difference Grid :
Pm = 45; % Truncation of the space grid
dP = 0.1; % Space Grid interval
dt = 1e-4; % Time Grid interval
D = 5+dP/2; % Development cost (US$)
% Jump characteristics
% Reference to function pdf_phi
m1 = 1/2;
m2 = 2;
dPhi=5e-2; % Step size for discretisation of PDF of Phi
% Number of steps : L=3/dPhi
% Grid Parameters :
m = Pm/dP % Number of space steps, better be an integer
n = T/dt % Number of time steps, better be an integer
L=3/dPhi % Has to be an integer
%%%%%%%%%%%%
% Initialisation
t_stor=100;
k = (m1+m2)/2 -1;
A = dt/(rho*dt+1);
for i = 0:m
P(i+1) = dP*i; % Prices scale
payoff(i+1)= max( q*P(i+1) - D , 0); % Payoff constraint
pplus(i+1) = A*1/2*(sigma^2*i^2+i*eta*Pbar-i^2*eta*dP-i*lambda*k);
pzero(i+1) = A*(1/dt-sigma^2*i^2 -lambda);
pminus(i+1)= A*1/2*(sigma^2*i^2-i*eta*Pbar+i^2*eta*dP+i*lambda*k);
pjump(i+1) = A*lambda;
if pplus(i+1)<0
  disp('pplus is negative at i=')
  disp(i)
end
if pzero(i+1)<0
  disp('pzero is negative at i=')
  disp(i)
end
if pminus(i+1)<0
  disp('pzero is negative at i=')
  disp(i)
end
end

vold = payoff;
%%%%%%%%%%%%
% Time loop
h = waitbar(0,'Please wait...');
time(n+1)=n*dt;
Stor(n/t_stor+1,:) = vold;
thres(n+1)=D/q;
for t = n-1:-1:0
time(t+1)=t*dt;
F(1)=0;
% Space loop
for i = 1:m
    % Expectation term
    % Probability density function of Phi
    % Parameters of the PDF
    sigma1=1/8; m1=1/2;
    sigma2=2/7; m2=2;
    temp = 0;
   
    for l = 1:L
     Phil = (l-1/2)*dPhi;
     lb = floor(i*Phil);
     ub = ceil(i*Phil);
     prob = dPhi*1/(2*sigma1*sqrt(2*pi))*exp(-1/2*((Phil-m1)/sigma1)^2)...
     +1/(2*sigma2*sqrt(2*pi))*exp(-1/2*((Phil-m2)/sigma2)^2);

     if ub<=m & lb~=ub
        vold2 = interp1([lb ub],[vold(lb+1) vold(ub+1)],i*Phil);
     elseif ub <= m & lb == ub
        vold2 = vold(lb);
     else
        vold2 = vold(m+1);
     end
     
     expt(i+1)=cumsum(prob*vold2);
    end %end for l = 1:L
   
   clear Phil lb ub prob % Elimination of unuseful variables
   % Calculation of the Option value at time t
   vold(m+2)=max( q*(m+1)*dP - D , 0);
   F(i+1) = pplus(i+1)*vold(i+2) + pzero(i+1)*vold(i+1) + pminus(i+1)*vold(i) + pjump(i+1)*expt(i+1);
   
   if F(i+1)<payoff(i+1)
      F(i+1)=payoff(i+1); % Boundary condition
   end
   
   %%%%%%%%%%
   % Threshold
   if (F(i+1)==payoff(i+1)) & (F(i)~=payoff(i))
      thres(t+1)=i*dP;
   end
   
end % end for  i = 1:m

vold = F; % New initial guess
% Reduction of vectors sizes

if t/t_stor==floor(t/t_stor)
  tt=t/t_stor;
  Stor(tt+1,:) = vold; % Storage
end

waitbar((n+1-t)/(n+1),h)
end %end for t = n-1:-1:0

close(h)
F_init=F; % Initial value of the option
save(filename)

T=linspace(1,10,100001)
plot(T,thres);

recbio

给你一个例子吧，向量化上面一个循环：
将
for i = 0:m
P(i+1) = dP*i; % Prices scale
payoff(i+1)= max( q*P(i+1) - D , 0); % Payoff constraint
pplus(i+1) = A*1/2*(sigma^2*i^2+i*eta*Pbar-i^2*eta*dP-i*lambda*k);
pzero(i+1) = A*(1/dt-sigma^2*i^2 -lambda);
pminus(i+1)= A*1/2*(sigma^2*i^2-i*eta*Pbar+i^2*eta*dP+i*lambda*k);
pjump(i+1) = A*lambda;
if pplus(i+1)<0
  disp('pplus is negative at i=')
  disp(i)
end
if pzero(i+1)<0
  disp('pzero is negative at i=')
  disp(i)
end
if pminus(i+1)<0
  disp('pzero is negative at i=')
  disp(i)
end
end

替换为：

i = 0:m;
P = dP*i; % Prices scale
payoff = max( q.*P - D , 0); % Payoff constraint
tmpA = sigma*sigma*i.^2;
tmpB = i*eta*Pbar -i.^2*eta*dP - i*lambda*k;
pplus = A/2*(tmpA + tmpB);
pzero = A*(1/dt - tmpA -lambda);
pminus = A/2*(tmpA - tmpB);
pjump = i - i + A*lambda ;

for i = 0:m;       %% 原来的循环显示可以照旧
if pplus(i+1)<0
  disp('pplus is negative at i=')
  disp(i)
end
if pzero(i+1)<0
  disp('pzero is negative at i=')
  disp(i)
end
if pminus(i+1)<0
  disp('pzero is negative at i=')
  disp(i)
end
end

recbio

另外，请将运行条件也贴一下，便于测试用。

下面一个循环，你可以按照上面的方法。

一般的，循环用矩阵运算代替时，原来的循环变量变成了矩阵。所以，运算* / ^ 等时需要用
点运算 如： .* ./ .^ 来代替。

加减不用。

还有什么问题再讨论吧。

喜新不厌旧

谢谢你费心啦！！

就是说如果将变量的下标去掉，并将乘除与平方前面加点，就变成矩阵运算了

不知道理解的对不对$考虑$

我没有什么运行条件呀，用这个循环就是为画出两个图，一个是T为横轴，P为纵轴; 一个是P为横轴，F为纵轴

我就


T=linspace(1,10,100001)
plot(T,thres);

和


plot(P,F)

也不知道对不对 $frage$ $害羞$

recbio

这样理解没有错。但有一点要主意，就是如果运算没有前后的依赖关系，这样的替换就没有问题了。假如前后的运算有依赖关系，比如你的第二个for 循环，里面的内循环依赖外循环，而且在条件判断上；所以，简化起来不是那么简单。

还有，你的第二个循环 i 从 1开始，可是你的 F 从1 + 1开始索引，是不是有什么问题？

ssjypm568

围观........[img][/img]