projet.sci

//PROJET EDP CHHOR DUBOIS

mode(-1);
clear

//initialisation
K=100; sigma=0.3; r=0.1; T=0.2; Smax=4*K;
compteur=0;

//Donnees du schema
M=20-1; // longueur de l'intervalle en espace
N=20;  // longueur de l'intervalle de temps
METHODE='BrSc-CN'; // 'BrSc-EE', 'BrSc-EI', 'BrSc-CN', 'H-EI'
CENTRAGE='CENTRE'; // 'DROIT' 'GAUCHE' ou 'CENTRE'
PAYOFF='payoff2';  // 'payoff1', 'payoff2'

//Parametres graphiques:
err_scale=10;
Xmin=0; Xmax=Smax; Ymin=-20; Ymax=K;
deltan=1; // Eventuellement, Affichage tous les deltan pas.

//Affichage des donnees:
printf('sigma=%5.2f, r=%5.2f, Smax=%5.2f\n',sigma,r,Smax);
printf('Maillage M+1= %5i, N=%5i\n',M+1,N);
printf('CENTRAGE : %s\n',CENTRAGE)
printf('METHODE  : %s\n',METHODE)
printf('PAYOFF  : %s\n',PAYOFF)

//Fonctions externes
printf('chargement des fonctions Normal(), Log(), BS();\n');
getf ('I:\Master2\EDP\func.sci');

//Fonction graphique
function ploot(t)
  xbasc();
  RECT=[Xmin,Ymin,Xmax,Ymax]
  Pex=BS(t,s);
  plot2d(s,PF,1,rect=RECT)
  plot2d(s,P,2,strf='000');   
  plot2d(s,zeros(s),0,strf='000');
  //plot2d(s,err_scale*(P-Pex),6,strf='000');
  xtitle('t='+string(t));
  xset('font size',3);
  legends(['Exacte',..
  METHODE+'-'+CENTRAGE+' (M+1='+string(M+1)+'; N='+string(N)+')',..
  //'Erreur (scale = '+ string(err_scale)+')'
  ],[1,2,6],1);
  xset('font size',1);
endfunction

//- Maillage
dt=T/N;
h=Smax/(M+1);
s=(0:M)'*h;
cfl=dt/h^2*(sigma*Smax)^2;
printf('nombre CFL : %5.3f\n',cfl);
M1=find(s==80); M2=find(s==90); M3=find(s==100);

//- Verification K dans maillage:
if (K/h-floor(K/h))<>0 then
  printf('attention K pas sur maillage;'); abort
end

//- CDI

//payoff classique
function y=P0(s);
  y=max(K-s,0);
endfunction

//payoff de l'application demandée
function y=P1(s);
  dim=size(s); m=dim(1)
  for k=1:m
    a=s(k);
    if ( a > (K/2) ) then
      y(k)= max(K-a,0);
    else
      y(k)=0;
    end;
  end; 
endfunction


//--------------------------
//- Initialiser la matrice A
//--------------------------
J=diag(ones(M,1),1);

select CENTRAGE
case 'DROIT' ; b=-r*s;  A1=diag(b) * (J-eye(J))/h;
case 'GAUCHE'; b=-r*s;  A1=diag(b) * (eye(J)-J')/h;
case 'CENTRE'; b=-r*s;  A1=diag(b) * (J-J')/(2*h);
else printf('CENTRAGE not programmed'); abort
end

a=sigma^2/2*s.^2;
A2=diag(a)*(2*eye(J)-J-J')/h^2;
A=A2+A1+r*eye(A2);
B=(speye(A)+dt*A);
B1=(2*speye(A)+dt*A);
B2=(2*speye(A)-dt*A);


//Initialiser P et graphique
select PAYOFF
  case 'payoff1' ; PF=P0(s); P=P0(s); 
  case 'payoff2' ; PF=P1(s); P=P1(s); 
end
ploot(0);
printf('appuyer sur la touche Return');  scanf('%c');
timer(); // demarrage compteur temps


getf ('I:\Master2\EDP\fonction.sci');


//Boucle sur P
for n=0:N-1

  t1=(n+1)*dt;

  //Schema
 
select METHODE
 

case 'BrSc-EE';
      P=max(P-dt*A*P,PF);
      
      
case 'BrSc-EI';
     [U,L]=decomp(B);
      C=monte(U,P);//Calcul de C tq UC=P;
      P2=P;
      P=descenteP(L,C,PF);//P tel que min( L*X-C,X-g)=0
       
 
case 'BrSc-CN';
      C=B2*P;
      [U,L]=decomp(B1);
      Y=monte(U,C);
    P2=P;
      P=descenteP(L,C,PF);// P tel que min(L*X-C,X-g)=0
      

case 'H-EI';
       compteur=0;
       eps=0.1;//seuil d'arret
       g=PF;//notre foncction payoff
       P=g;
       P=B\P;
       err=eps+1;
       dim=size(P);
       m=dim(1);
       while ( err > eps)     
               Pav=P;               
               a1=B*P-Pav;
               a2=P-g;
               C=B;           
            for i=1:m
             if ( a1(i) > a2(i) ) then
              C(i,:)=0;//on modifie la matrice B suivant la règle
              C(i,i)=1;
              P(i)=g(i)
               end;
                 end;
                 P=C\P;       
          compteur=compteur+1;
          err=norm(P-Pav);
       end;
      
      
case 'H-CN'; // A revoir
       compteur=0;
       eps=0.1;
       g=PF;
       P=PF;
       //P=B\P;
       err=eps+1;
       dim=size(P);
       m=dim(1);
       while ( err > eps)     
               Pav=P;
               b=B2*P;
               P=B1\P;               
               a1=B1*P-b;
               a2=P-g;
            for i=1:m     
                  if ( a1(i) > a2(i) ) then
              B1(i,:)=0;
              B1(i,i)=1;
              //P(i)=g(i)
              b(i)=g(i)
               end;
                 end;       
          compteur=compteur+1;
          err=norm(P-Pav);
       end;
      

case 'FT-EI';//A revoir
  I=size(B);
  Imax=I(1);
   g=PF;
   P=PF;
   P = B\P;
   k=1;
   Pav=P;
      P = B\P;
      P(1)=g(1);
      P(Imax)=0;
      i=1;
      //j=0;
  while (i<=Imax)
        if (i<k) then P(i)=g(i);  end;
    if (i>=k)
             P(i)= B(i,:)\P;          
              if (i==k) then
              if (P(i)-g(i)<0) then i=k+1;
              end;              
                  P(i)=B(i,:)*P   - Pav(i);
              else
                  if (P(i-1)<0) then
                      i=k-1;             
                      P(i)=B(i,:)*P   - Pav(i);
                  end;
              end;
              i=i+1;
         end;
          //k=k+1;               
end;       
           

else; printf('METHODE non programmee'); abort;
end


  if modulo(n+1,deltan)==0 then
    //- Graphiques:
    ploot(t1)
    //-Solution Black et scholes
    Pex=BS(t1,s);
    //- Textes:
    errLI=norm(Pex-P,'inf');
    errL2=sqrt(h)*norm(Pex-P);
    //printf('t=%5.2f; Err.Linf=%8.5f, Err.L2=%8.5f; ',t1,errLI, errL2);
    //scanf('%c');
  end
 
  P(M+1)=0;
 
end

//on aurait voulu mettre en évidence les points S=80, 90, 100
plot2d4(s(M1),P(M1),5); // affichage du point d'abscisse M1=80
plot2d4(s(M2),P(M2),5); // affichage du point d'abscisse M2=90
plot2d4(s(M3),P(M3),5); // affichage du point d'abscisse M3=100
t_total=timer();
printf('\n')
if (compteur<>0) then  printf('Convergence en %5.2f iteration\n ',compteur);end
printf('temps total = %5.2f\n',t_total);
printf('programme  termine'' normalement');