acpricecn.m

finite difference method with Crank-Nikolson approach(SOR) for American call option

function Price=acpricecn(CurrentPrice,StrikePrice,IntRate,TimetoMaturity,Std,DeltaT,DeltaS,OverRel,Tolerance)
%finite difference method with Crank-Nikolson approach(SOR) for American call option
%CurrentPrice---price of underlying asset
%StrikePrice---strike price
%IntRate---interest rates
%TimetoMaturity---time to maturity
%Std----sigma of underlying asset
%DeltaT---the length of time
%DeltaS---the length of Undelying asset
%OverRel---Overrelax Paramater
%Tolerance---the tolerance parameter.
%q=0---installment rates
%DivRate=0---divident rate
%还要修改,它不同看跌期权,因为f_{i,M}<>0

q=1;
DivRate=0.04;
MaxPrice=2*max(CurrentPrice,StrikePrice)*exp(IntRate*TimetoMaturity);
NumOftime=round(TimetoMaturity/DeltaT);%取整
DeltaT=TimetoMaturity/NumOftime;
NumOfPrice=round(MaxPrice/DeltaS);
DeltaS=MaxPrice/NumOfPrice;
OldValue=zeros(NumOfPrice-1,1);
NewValue=zeros(NumOfPrice-1,1);
%y = linspace(a,b,n) generates a row vector y of n points linearly spaced ;
%between and including a and b;
AssetPrice=linspace(0,MaxPrice,NumOfPrice+1)';%生成0--MaxPrice的向量值
VectorI=0:NumOftime;
VectorJ=0:NumOfPrice;

Payoff=max(AssetPrice(2:NumOfPrice)-StrikePrice,0);%Call option
PastValue=Payoff; %Boundary condition_1 for call
BoundValue=0;%NumOfPrice*DeltaS-StrikePrice*exp(-IntRate*DeltaT*(NumOftime-VectorI));%Boundary condition_2 for Call

A=0.25*DeltaT*(Std^2*VectorJ-(IntRate-DivRate)).*VectorJ;
B=-DeltaT*0.5*(Std^2*VectorJ.^2+IntRate);
C=0.25*DeltaT*(Std^2*VectorJ+(IntRate-DivRate)).*VectorJ;
M2=diag(A(3:NumOfPrice),-1)+diag(1+B(2:NumOfPrice))+diag(C(2:NumOfPrice-1),1);
%standard American option
%Aux=zeros(NumOfPrice-1,1);
%installment American option
Aux=q*ones(NumOfPrice-1,1);%where q is installment rates

for i=NumOftime:-1:1
      Aux(1)=q;
     D=M2*PastValue+Aux;
    OldValue=PastValue;
    Error=realmax;
    while Tolerance<Error
         NewValue(1)=max(Payoff(1),OldValue(1)+OverRel/(1-B(2))*(D(1)-(1-B(2))*OldValue(1)+C(2)*OldValue(2)));
         for k=2:NumOfPrice-2
             NewValue(k)=max(Payoff(k),OldValue(k)+OverRel/(1-B(k+1))*(D(k)+A(k+1)*NewValue(k-1)-(1-B(k+1))*OldValue(k)+C(k+1)*OldValue(k+1)));
         end
         NewValue(NumOfPrice-1)=max(Payoff(NumOfPrice-1),OldValue(NumOfPrice-1)+OverRel/(1-B(NumOfPrice))*(D(NumOfPrice-1)+A(NumOfPrice)*NewValue(NumOfPrice-2)-(1-B(NumOfPrice))*OldValue(NumOfPrice-1)));
         Error=norm(NewValue-OldValue);
         OldValue=NewValue;
     end
     PastValue=NewValue;
 end
 NewValue=[BoundValue(1);NewValue;0];
 LocUp=ceil(CurrentPrice/DeltaS);%取整>=
 LocDown=floor(CurrentPrice/DeltaS);%取整<=
 if LocUp==LocDown
     Price=NewValue(LocDown+1,1);
 else
     Price=NewValue(LocDown+1,1)+(CurrentPrice-LocDown*DeltaS)*(NewValue(LocUp+1,1)-NewValue(LocDown+1,1))/DeltaS;
 end
 %appricecn(96,100,0.05,3/12,0.2,1/400,0.8,1.5,0.001)