Qlw(M,n)=cfvw(n,ylw(M,n));
Quw(M,n)=cfvw(n,yuw(M,n));
end
miny0=min(yl0(M,1:M));maxy0=max(yu0(M,1:M));
minyw=min(ylw(M,1:M));maxyw=max(yuw(M,1:M));
t0=ceil((maxy0-miny0)/deltay)+1;tw=ceil((maxyw-minyw)/deltay)+1;
%求t=M*deltat时刻,最大NT区域的价值函数
II0=ones(M,t0);Q0=II0;Q1=II0;IIw=ones(M,tw);Qw=IIw;Q2=IIw;
%不带期权
for n=1:M
for t=1:t0
y=miny0+(t-1)*deltay;
if yl0(M,n)<=y<=yu0(M,n)
Q0(n,t)=cfv(n,y);
elseif yl0(M,n)>y
Q0(n,t)=exp(gama*R*(1+lamda)*S(M,n)*(yl0(M,n)-y))*Ql0(M,n);
else
Q0(n,t)=exp(gama*R*(lamda-1)*S(M,n)*(y-yu0(M,n)))*Qu0(M,n);
end
end
end
%带期权
for n=1:M
for t=1:tw
y=minyw+(t-1)*deltay;
if ylw(M,n)<=y<=yuw(M,n)
Qw(n,t)=cfvw(n,y);
elseif ylw(M,n)>y
Qw(n,t)=exp(gama*R*(1+lamda)*S(M,n)*(ylw(M,n)-y))*Qlw(M,n);
else
Qw(n,t)=exp(gama*R*(lamda-1)*S(M,n)*(y-yuw(M,n)))*Quw(M,n);
end
end
end
%依次用后向递推探测任一时刻的NT下界、NT上界以及最大NT区域内的价值函数
%不带期权
for m=(M-1):(-1):1
R1=exp(r*(M+1-m)*deltat);
for n=1:m
%NT下界
for t=1:t0
Q(t)=exp(gama*R1*(1+lamda)*t*deltay*S(m,n))*(p*Q0(n+1,t)+(1-p)*Q0(n,t));
end
[,l]=max(Q);
Ql0(m,n)=;yl0(m,n)=miny0+(l-1)*deltay;
clear ,y;
%NT上界
for t=1:t0
Q(t)=exp(gama*R1*(lamda-1)*(t0+1-t)*deltay*S(m,n))*(p*Q0(n+1,t)+(1-p)*Q0(n,t));
end
[,u]=max(Q);
Qu0(m,n)=;yu0(m,n)=miny0+(u-1)*deltay;
clear ,y;
%最大NT区域的价值函数
for t=1:t0
y=miny0+(t-1)*deltay;
if yl0(m,n)<=y<=yu0(m,n)
Q1(n,t)=p*Q0(n+1,t)+(1-p)*Q0(n,t);
elseif yl0(m,n)>y
Q1(n,t)=exp(gama*R1*(1+lamda)*S(m,n)*(yl0(m,n)-y))*Ql0(m,n);
else
Q1(n,t)=exp(gama*R1*(lamda-1)*S(m,n)*(y-yu0(m,n)))*Qu0(m,n);
end
end
end
Q0=Q1;
end
%带期权
for m=(M-1):(-1):1
R1=exp(r*(M+1-m)*deltat);
for n=1:m
%NT下界
for t=1:tw
Q(t)=exp(gama*R1*(1+lamda)*t*deltay*S(m,n))*(p*Qw(n+1,t)+(1-p)*Qw(n,t));
end
[,l]=max(Q);
Qlw(m,n)=;ylw(m,n)=minyw+(l-1)*deltay;
clear ,l;
%NT上界
for t=1:tw
Q(t)=exp(gama*R1*(lamda-1)*(tw+1-t)*deltay*S(m,n))*(p*Qw(n+1,t)+(1-p)*Qw(n,t));
end
[,u]=max(Q);
Quw(m,n)=;yuw(m,n)=minyw+(u-1)*deltay;
clear Q,u;
%最大NT区域的价值函数
for t=1:tw
y=minyw+(t-1)*deltay;
if ylw(m,n)<=y<=yuw(m,n)
Q2(n,t)=p*Qw(n+1,t)+(1-p)*Qw(n,t);
elseif ylw(m,n)>y
Q2(n,t)=exp(gama*R1*(1+lamda)*S(m,n)*(ylw(m,n)-y))*Qlw(m,n);
else
Q2(n,t)=exp(gama*R1*(lamda-1)*S(m,n)*(y-yuw(m,n)))*Quw(m,n);
end
end
end
Qw=Q2;
end
clear tw t0 y Q0 Q1 Q2 Qw;
yl=ylw-yl0;yu=yuw-yu0;