求教利用指数效用函数做期权套期保值的matlab程序-经管之家官网！

源于Michael Monoyio 2004年在Journal of Economic Dynamics & Control杂志上发表的一篇论文Option pricing with transaction costs using a Markov chain approximation，本人利用文中的算法编了一个matlab程序，利用最大化到期财富的效用来计算套期保值策略，可惜屡战屡败，编出的程序与结论相差太远，不知哪位高手编过类似相关问题的程序，烦请指点一二！附上本人程序求教。
%编写求t=M*deltat时刻价值函数的文件（就不采取买卖交易的情形），之后就可调用
%不带期权的
function f=cfv(n,y)
global p gama S M;
f=p*(-exp(-gama*y*S(M+1,n+1)))+(1-p)*(-exp(-gama*y*S(M+1,n)));
return;
%带期权的
function f=cfvw(n,y)
global p gama S K M;
if K<S(M+1,n)
f=p*(-exp(-gama*(y*S(M+1,n+1)+S(M+1,n+1)-K)))+(1-p)*(-exp(-gama*(y*S(M+1,n)+S(M+1,n)-K)));
elseif K>=S(M+1,n+1)
f=p*(-exp(-gama*y*S(M+1,n+1)))+(1-p)*(-exp(-gama*y*S(M+1,n)));
else
f=p*(-exp(-gama*(y*S(M+1,n+1)+S(M+1,n+1)-K)))+(1-p)*(-exp(-gama*y*S(M+1,n)));
end
return;
%编写一个求两个变量差的正部的函数
function f=Identifyf(a,b)
if a>b
f=a-b;
else
f=0;
end
return;
以上都是M文件。
global p gama S K M;
T=1;M=50;deltat=T/M;sigma=0.25;S0=15;K=15;r=0.15;mu=0.15;
lamda=0.005;gama=0.1;R=exp(r*deltat);
II=ones(M+1,M+1);S=II;
deltay=1*deltat;
p=0.5*(1+mu/sigma*deltat^0.5);
u=exp(sigma*deltat^0.5);d=exp(-sigma*deltat^0.5);
%将S二叉树划
for m=1:(M+1)
for n=1:m
S(m,n)=S0*u^(2*n-m-1);
end
end
%求t=M*deltat时的买卖NT边界上的价值函数
yu0=-3*II;yl0=-3*II;yuw=-3*II;ylw=-3*II;yu=-3*II;yl=-3*II;
for n=1:M
%不带期权
yu0(M,n)=(gama*S(M,n)*(u-d))^(-1)*log((p*(u-R*(1-lamda)))/((1-p)*(R*(1-lamda)-d)));
yl0(M,n)=(gama*S(M,n)*(u-d))^(-1)*log((p*(u-R*(1+lamda)))/((1-p)*(R*(1+lamda)-d)));
Ql0(M,n)=cfv(n,yl0(M,n));
Qu0(M,n)=cfv(n,yu0(M,n));
%带期权
yuw(M,n)=(gama*S(M,n)*(u-d))^(-1)*(log((p*(u-R*(1-lamda)))/((1-p)*(R*(1-lamda)-d)))+(Identifyf(S(M+1,n+1),K)-Identifyf(S(M+1,n),K))*gama);
ylw(M,n)=(gama*S(M,n)*(u-d))^(-1)*(log((p*(u-R*(1+lamda)))/((1-p)*(R*(1+lamda)-d)))+(Identifyf(S(M+1,n+1),K)-Identifyf(S(M+1,n),K))*gama);
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
[QQ,l]=max(Q);
Ql0(m,n)=QQ;yl0(m,n)=miny0+(l-1)*deltay;
clear QQ,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
[QQ,u]=max(Q);
Qu0(m,n)=QQ;yu0(m,n)=miny0+(u-1)*deltay;
clear QQ,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
[QQ,l]=max(Q);
Qlw(m,n)=QQ;ylw(m,n)=minyw+(l-1)*deltay;
clear QQ,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
[QQ,u]=max(Q);
Quw(m,n)=QQ;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;