//The  CRR random walk as a function of n /////////////
clear;
T=1 ;Nmax=500;
function delta_t=delta_t(n);
  delta_t=T./n ;//mind the "./"
endfunction;
S0=140 ;sigma=0.2;
function up=up(n);
  up=exp(sigma*sqrt(delta_t(n)));
endfunction;
function down=down(n)
  down=exp(-sigma*sqrt(delta_t(n)));
endfunction;
function s=S(n,i,j);
  s=S0*up(n)^j.*down(n)^(i-j);
endfunction;
//computation of Call prices
K=S0*0.9;
//plot2d(0:Nmax,K*ones(1,Nmax+1));
// Call option
function c=phi(S);
  c=max(S-K,0);
endfunction;
function p=p(n); //risk-neutral probability
  p=(1-down(n))./(up(n)-down(n));// mind the "./"
endfunction; 
// backward induction computation.
//CC=zeros(Nmax,Nmax+1,Nmax+1);
//for n=1:Nmax
//  CC(n,n+1,1:n+1)=phi(S(n,n,0:n)); //price wellknown at the end.
//  for i=n-1:-1 :0 //downward induction
//    for j=0:i; 
//      CC(n,i+1,j+1)=p(n)*CC(n,i+2,j+1+1)+(1-p(n))*CC(n,i+2,j+1); 
//    end ;
//  end ;
//end;
//function c=CallFromMatrix(n,i,j);//allows CallFromMatrix(n,0,0)
//  c=CC(n,i+1,j+1)
//endfunction;
//disp(CallFromMatrix(Nmax,0,0));

// the binomial expectation computation
function c=Call(n)
  c=binomial(p(n),n)*phi(S(n,n,0:n))'
endfunction;
//price at time=0
C=zeros(1,Nmax);
for n=1:Nmax;
  C(1,n)=Call(n);
end;
//disp(max(abs(CallFromMatrix(1:Nmax,0,0)-C(1:Nmax)')),'error=')
plot(20:Nmax,C(1,20:Nmax));
//// Normal distribution and erf (error function///
function y=N(x);
  y=(erf(x/sqrt(2))+1)/2;
endfunction;
function BS=BlackScholes(S,K,r,T,sigma);
  d1=(log(S/K)+T*(r+(sigma^2)/2))/(sigma*sqrt(T));
  d2=d1-sigma*sqrt(T);
  BS=S*N(d1)-K*exp(-r*T)*N(d2);
endfunction;
CallBS=BlackScholes(S0,K,0,T,sigma);
plot(1:Nmax,CallBS,'-r');