请教循环向量化
偶有个程序，循环比较大，运行不出来，有人建议循环改成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(,,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 fori = 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); 给你一个例子吧，向量化上面一个循环：
将
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 另外，请将运行条件也贴一下，便于测试用。
下面一个循环，你可以按照上面的方法。
一般的，循环用矩阵运算代替时，原来的循环变量变成了矩阵。所以，运算* / ^ 等时需要用
点运算 如： .* ./ .^ 来代替。
加减不用。
还有什么问题再讨论吧。 谢谢你费心啦！！
就是说如果将变量的下标去掉，并将乘除与平方前面加点，就变成矩阵运算了
不知道理解的对不对$考虑$
我没有什么运行条件呀，用这个循环就是为画出两个图，一个是T为横轴，P为纵轴; 一个是P为横轴，F为纵轴
我就
T=linspace(1,10,100001)
plot(T,thres);
和
plot(P,F)
也不知道对不对 $frage$ $害羞$ 这样理解没有错。但有一点要主意，就是如果运算没有前后的依赖关系，这样的替换就没有问题了。假如前后的运算有依赖关系，比如你的第二个for 循环，里面的内循环依赖外循环，而且在条件判断上；所以，简化起来不是那么简单。
还有，你的第二个循环 i 从 1开始，可是你的 F 从1 + 1开始索引，是不是有什么问题？
评论
围观........
页:
[1]