//The  CRR random walk /////////////
clear;clf(3);
n=100 ;T=1 ;delta_t=T/n ;
S0=140 ;sigma=0.2;
up=exp(sigma*sqrt(delta_t)) ;down=1/up;
function s=S(i,j);
  s=S0*up^j.*down^(i-j);
endfunction;
//computation of Call prices
K=S0*0.9;
xset("window",3);
plot2d(0:n,K*ones(1,n+1));
// Call option
function c=phi(S);
  c=max(S-K,0);
endfunction;
p=(1-down)/(up-down); //risk-neutral probability
CC(n+1,1:n+1)=phi(S(n,0:n)); //price wellknown at the end.
for i=n-1:-1 :0 //downward induction
  for j=0:i;      CC(i+1,j+1)=p*CC(i+2,j+1+1)+(1-p)*CC(i+2,j+1); // 
  end ;
end ;
function c=Call(i,j);//allows Call(0,0)
  c=CC(i+1,j+1)
endfunction;
//hedge factor A(i,j)=(Call(i+1,j+1)-Call(i+1,j))/(S(i+1,j+1)-S(i+1,j))
//NB: i<n
function a=A(i,j);
  a=(Call(i+1,j+1)-Call(i+1,j))/(S(i+1,j+1)-S(i+1,j));
endfunction;
/////////////////////////////////////////
colo=0;//clf(3);
/////////////////////////////////////////
// Simulations of M CRR trajectories /////////////////
M=1 ;colo=modulo(colo+1,8);
plot2d(0:n,K*ones(1,n+1));
//choice of M ramdom sequences of 0 and 1
//Probability q=0.5 for 1
q=0.5;
MdeltaJ=zeros(M,n);
MdeltaJ=int(q+rand(M,n));
MJ=zeros(M,n);
MJ=cumsum(MdeltaJ,"c");
// computation of resulting trajectories
Mttraj=zeros(M,n+1);
for m=1 :M
Mttraj(m,1)=S0;
for i=1 :n
Mttraj(m,i+1)=S(i,MJ(m,i));
end ;
end ;
//computation of Call-price trajectories
MCalltraj=zeros(M,n+1);
for m=1 :M
  MCalltraj(m,1)=Call(0,0);
  for i=1 :n
    MCalltraj(m,i+1)=Call(i,MJ(m,i));
  end ;
end ;
//computation of M Hedgefactor trajectories
MAtraj=zeros(M,n);
for m=1 :M
  MAtraj(m,0+1)=A(0,0);
  for i=1 :n-1
    MAtraj(m,i+1)=A(i,MJ(m,i));
  end ;
end ;
// plot of the M resulting stock trajectories
xset("window",3);
for m=1 :M
plot2d(0 :n,Mttraj(m,1 :n+1),style=colo) ;
end ;
//plot of the M K+Call trajectories in different colors
for m=1 :M
  plot2d(0:n,K+MCalltraj(m,1:n+1),style=colo) ;
end ;
//plot of M K*hedge-factor trajectories in different colors
for m=1 :M
  plot2d(0:n-1,K*MAtraj(m,1:n),style=colo) ;
end ;  
xs2eps(3,'dessin.');
//////////