[call,put]=americanoption(s,k,r,T,vol,d)

function [call,put]=americanoption(s,k,r,T,vol,d)

%author:QuantHR
%Date:2005-11-01;


h=1-exp(-r*T);
a=2*r/(vol)^2;
b=2*(r-d)/(vol)^2;
r1=0.5*(-(b-1)+sqrt((b-1)^2+4*a/h));
r2=0.5*(-(b-1)-sqrt((b-1)^2+4*a/h));

maxiter = 5000;
tol = 1e-5;
o = optimset('MaxIter', maxiter, 'TolFun', tol);

[ss, fval, exitFlag] = fminbnd(@corp1,s/20,s*20,o,s,k,T,r,vol,d,r1);
[sss, fval, exitFlag] = fminbnd(@corp2,s/20,s*20,o,s,k,T,r,vol,d,r2);

a1=ss/r1*(1-exp(-d*T)*norm1(ss,k,r,T,vol,d));
a2=-sss/r2*(1-exp(-d*T)*norm2(sss,k,r,T,vol,d));

if s<ss
	call=blsprice(s,k,r,T,vol,d)+a1*(s/ss)^(r1);
else
	call=s-k;
end

if s>sss
	[cc,pp]=blsprice(s,k,r,T,vol,d);
	put=pp+a2*(s/sss)^(r2);
else
	put=k-s;
end



function cd=corp1(ss,s,k,T,r,vol,d,r1)
cd=abs(ss-k-blsprice(ss,k,r,T,vol,d)-(1-exp(-d*T)*norm1(ss,k,r,T,vol,d))*ss/r1);

function cd=corp2(sss,s,k,T,r,vol,d,r2)
[c,p]=blsprice(sss,k,r,T,vol,d);
cd=abs(k-sss-p+(1-exp(-d*T)*norm2(sss,k,r,T,vol,d))*sss/r2);

function nd=norm1(s,k,r,T,vol,d)
d1=(log(s./k)+(r-d+vol.^2/2).*T)./(vol.*sqrt(T));
nd=normcdf(d1);

function nd=norm2(s,k,r,T,vol,d)
d1=(log(s./k)+(r-d+vol.^2/2).*T)./(vol.*sqrt(T));
nd=normcdf(-d1);