// Mod`ele de Ho et Lee TP0910MMDFA.sce
clear ;Nmax=8 ;Tmax=Nmax;delta_t=Tmax/Nmax;
pi=0.5 ;delta=1.01;r=0.02;

//Une fonction Z0 avec taux actuariel initial independant de T
//
function z0=Z0(k); z0=(1+r)^(-k*delta_t) ; endfunction;
//  
// une fonction Z0 plus réaliste: smirk de maturité
// 
rr=zeros(Nmax+1);ZZ0=zeros(Nmax+1);
for k=0:Nmax;
    rr(k+1)=r*(1+(1/(30*Nmax))*k*(k-Nmax*3/2));
    ZZ0(k+1)=(1+rr(k+1))^(-(k)*delta_t);
end;    
//function z0=Z0(k); z0=ZZ0(k+1) ; endfunction;
xset("window",1);
plot(0 :Nmax,Z0(0 :Nmax)) ;
//Introduction de la fonction èta du modèle
function ee=eta(l,x);
if x==0
ee=(1^l)./(pi+(1-pi)*delta^l);
else ee=delta^l./(pi+(1-pi)*delta^l);
end ;
endfunction;
//calcul des valeurs du zero coupons: la fonctions 
//Z(n,j,k)=ZZ(n+1,j+1,k+1) où t=n*delta_t, j est le nombre 
//de up et T=k*delta_t. 
ZZ=ones(Nmax+1,Nmax+1,Nmax+1);
for k=0 :Nmax
ZZ(0+1,1,k+1)=Z0(k);
end ;
for n=1 :Nmax
for k=n :Nmax
ZZ(n+1,0+1,k+1)=eta(k-n,0)*ZZ(n-1+1,0+1,k+1)/ZZ(n-1+1,0+1,n+1);
for j=1 :n
ZZ(n+1,j+1,k+1)=eta(k-n,1)*ZZ(n-1+1,j-1+1,k+1)/ZZ(n-1+1,j-1+1,n+1);
end ;
end ;
end ;
// pour s'affranchir de "+1" dans les 3 arguments de ZZ, 
//on introduit la fonction Z:
function z=Z(n,j,k) //t=n*delta_t et T=k*delta_t
z=ZZ(n+1,j+1,k+1);
endfunction;

// Calcul du taux actuariel A(t,T) associé au zéro coupon Z(t,T): 
function a=A(n,j,k)
a=exp(-log(Z((n),j,k))'./((k)-(n)))-1 // attention: n<k
endfunction;

//évolution de la courbe des taux
xset("window",5);
//                  t=0
AbscT=zeros(Nmax+1);
OrdoZ=zeros(Nmax+1);
for kk=1:Nmax+1 ; //kk commence toujours à 1
   AbscT(kk)=kk-1;
   OrdoZ(kk)=Z(0,0,kk-1);
end;
plot2d(AbscT,OrdoZ);
//                  t=delta_t
AbscT=zeros(Nmax);
OrdoZ=zeros(Nmax);
for kk=1:Nmax ; //kk commence toujours à 1
   AbscT(kk)=kk-1;  
   OrdoZ(kk)=Z(1,0,kk+1-1);//apres un "down"
end;
plot(AbscT,OrdoZ,"--r");
for kk=1:Nmax ; //kk commence toujours à 1
   AbscT(kk)=kk-1;
   OrdoZ(kk)=Z(1,1,kk+1-1);//après un "up"
end;
plot(AbscT,OrdoZ,"--b");
//                  t=2*delta_t
AbscT=zeros(Nmax-1);
OrdoZ=zeros(Nmax-1);
//j=0
for kk=1:Nmax-1 ; //kk commence toujours à 1
   AbscT(kk)=kk-1;  
   OrdoZ(kk)=Z(2,0,kk+2-1);
end;
plot(AbscT,OrdoZ,"-+r");
//j=1 , AbscT inchangé
OrdoZ=zeros(Nmax);
for kk=1:Nmax-1 ; //kk commence toujours à 1
   OrdoZ(kk)=Z(2,1,kk+2-1);//j=0
end;
plot(AbscT,OrdoZ,"-+g");
//j=2 , AbscT inchangé
OrdoZ=zeros(Nmax);
for kk=1:Nmax-1 ; //kk commence toujours à 1 
   OrdoZ(kk)=Z(2,2,kk+2-1);//j=0
end;
plot(AbscT,OrdoZ,"-+b");

xset("window",6);
// Par boucles
//                  t=n*delta_t
for n=0:2
    AbscT=zeros(Nmax+1-n);
    OrdoZ=zeros(Nmax+1-n);
    for j=0:n;
        for kk=1:Nmax+1-n ; //kk commence toujours à 1
          AbscT(kk)=kk-1;   //k=n+kk-1
          OrdoZ(kk)=Z(n,j,n+kk-1);
        end;
    plot(AbscT,OrdoZ);
    end;
end;    
// taux actuariel
xset("window",7);
// Par boucles
//                  t=n*delta_t
for n=0:2
AbscT=zeros(Nmax+1-n);
OrdoA=zeros(Nmax+1-n);
    for j=0:n;
        for kk=1:Nmax-n ; //kk commence toujours à 1
          AbscT(kk)=kk;   //k=n+kk
          OrdoA(kk)=A(n,j,n+kk);
        end;
     plot(AbscT,OrdoA);
    end;
end;    
////////////////////////////////////////////////////