cevfds.m

用matlab实现障碍期权的复制

function price=CEVFDS(Nx,Nt,Y,S,K,r,q,sigma,T,Smax,Smin,type)
h=(Smax-Smin)/Nx; k=T/Nt; mu=k/h^2;
D1=(diag(1:Nx-1).*h+diag(ones(1,Nx-1)).*Smin);
D2=(diag(1:Nx-1).*h+diag(ones(1,Nx-1)).*Smin).^(2*Y);
T1=diag(ones(Nx-2,1),1)+diag(-ones(Nx-2,1),-1);
T2=diag(ones(Nx-1,1)*(-2))+diag(ones(Nx-2,1),-1)+...
diag(ones(Nx-2,1),1);
B=(1+r*k)*eye(Nx-1)-0.5*mu*sigma^2*D2*T2-...
(r-q)*k/(2*h)*D1*T1;
S0=Smin+h:h:Smax-h;
Tvec=0:k:T;
% The boundry value changes depending on the type of
% option for which a price is required. Hence, the use
% of the "switch" function.
switch type
case 'call'
V0n=max((S0-K),0)';
VmSmin=(0*Tvec)';
VmSmax=(Smax-K*exp(-r*Tvec))';
case 'put'
V0n=max((K-S0),0)';
VmSmin=K*exp(-r*Tvec)';
VmSmax=zeros(1,Nt+1)';
case 'doc'
V0n=max((S0-K),0)';
VmSmin=(0*Tvec)';
VmSmax=(Smax-K*exp(-r*Tvec))';
case 'dop'
V0n=max((K-S0),0)';
VmSmin=(0*Tvec)';
VmSmax=zeros(1,Nt+1)';
case 'uoc'
V0n=max((S0-K),0)';
VmSmin=(0*Tvec)';
VmSmax=zeros(1,Nt+1)';
case 'uop'
V0n=max((K-S0),0)';
VmSmin=K*exp(-r*Tvec)';
VmSmax=zeros(1,Nt+1)';
end
bm=zeros(1,Nx-1)';
sol=zeros(Nt,Nx-1);
Vm=V0n;
for i=1:Nt;
bm(1,1)=VmSmin(i+1,1)*((r-q)*k/(2*h)*(Smin+h)-...
sigma^2*mu/2*(Smin+h)^(2*Y));
bm(Nx-1,1)=VmSmax(i+1,1)*((r-q)*k/(2*h)*(Smin+(Nx-1)*h)-...
sigma^2*mu/2*(Smin+(Nx-1)*h)^(2*Y));
Vmplus1=B\(Vm-bm);
sol(i,:)=Vmplus1';
Vm=Vmplus1;
end;
sol=flipud(sol);
sol=[sol;V0n'];
sol=[flipud(VmSmin) sol flipud(VmSmax)];
price=interp1([Smin S0 Smax],sol(1,:),S);