T=1
r=0.05
sigma=0.3
S0=140
K=145

function z=cdf(t)
    z=cdfnor("PQ",t,0,1)
endfunction

function y=CallBS(S,K,r,T,sigma)
    d1=(log(S/K)+T*(r+sigma^2/2))/(sigma*sqrt(T))
    d2=d1-sigma*sqrt(T)
    y=S*cdf(d1)-K*exp(-r*T)*cdf(d2)
endfunction


function y=PutBS(S,K,r,T,sigma)
    d1=(log(S/K)+T*(r+sigma^2/2))/(sigma*sqrt(T))
    d2=d1-sigma*sqrt(T)
    y=K*exp(-r*T)*cdf(-d2)-S*cdf(-d1)
endfunction

function delta=deltat(n)
    delta=T/n
endfunction
function re=R(n)
    re=exp(r*deltat(n))
endfunction
function u=up(n)
    u=exp(sigma*sqrt(deltat(n)))
endfunction
function dd=down(n)
    dd=exp(-sigma*sqrt(deltat(n)))
endfunction
function pp=p(n)
    pp=(R(n)-down(n))/(up(n)-down(n))
endfunction

function y=S(i,j,n)
    y=S0.*(up(n))^j.*(down(n))^(i-j)
endfunction
function phi=phi(S)
    phi=max(K-S,0)
endfunction
function z=C(i,j,n)
    z=(phi(S(n,j:(j+n-i),n))*binomial(p(n),n-i)')/R(n)^(n-i);
endfunction
Nmax=250
for n=1:Nmax
    Call(n)=C(0,0,n);
end
plot2d(10:Nmax,Call(10:Nmax))
plot2d(10:Nmax,ones(1,Nmax-9)*PutBS(S0,K,r,T,sigma))
