//! DF pour Black et Scholes 
mode(-1); 
clear

//------------------------
//- DONNEES
//------------------------
K=100; sigma=0.2; r=0.1; T=1;  Smin=20; Smax=200;  

//- CDI: Condition initiale
//COMPLETER definition de P0, ug,ud
function y=P0(s); y=max(K-s,0); endfunction
//function y=ug(t); y=K*exp(-r*t)-Smin; endfunction
function y=ug(t); y=K-Smin; endfunction
function y=ud(t); y=0; endfunction


//------------------------
//- DONNEES NUMERIQUES
//------------------------
I=20-1;
N=10;

//- CHOIX SCHEMA:
//SCHEMA='EE'; // 'EE' ou 'EI' ou 'CN' 
//SCHEMA  ='EE-AMER'
//SCHEMA  ='EI-AMER-PSOR'
SCHEMA  ='EI-AMER-PUL'
//SCHEMA  ='EI-AMER-NEWTON'
//SCHEMA  ='EI-AMER-SPLIT'


//- CHOIX CENTRAGE
CENTRAGE='CENTRE'; // 'CENTRE', 'DROIT', 'GAUCHE' 


//- Paramètres pour la fenetre graphique:
Xmin=Smin; Xmax=Smax; Ymin=-20; Ymax=K; 
err_scale=0; //- Echelle pour le graphe d'erreur.
deltan=N/10; // Eventuellement, Affichage uniquement tous les deltan pas.


//- Affichage des données:
printf('sigma=%5.2f, r=%5.2f, Smax=%5.2f\n',sigma,r,Smax);
printf('Maillage I= %5i, N=%5i\n',I,N);
printf('CENTRAGE : %s\n',CENTRAGE)
printf('SCHEMA: %s\n',SCHEMA)

//-----------------------------------------------
//- FORMULE de Black et Scholes.
//-----------------------------------------------
function y=Normal(x)
  y=cdfnor("PQ",x,zeros(x),ones(x)); //- N'accepte pas x=-%inf
endfunction

function P=BS(t,s,K,r,sigma)
  //- Formule de Black et scholes pour le put europeen vanille.
  if t==0 
    P=max(K-s,0);
  else
    P=ones(s)*K*exp(-r*t); 
    i=find(s>0); 
    tau=sigma^2*t;
    dm=(log(s(i) /K) + r*t - 0.5*tau) / sqrt(tau);
    dp=(log(s(i) /K) + r*t + 0.5*tau) / sqrt(tau);
    P(i)=K*exp(-r*t)*(Normal(-dm)) - s(i).*(Normal(-dp));
  end
endfunction

//--------------------
//- Fonction graphique
//--------------------
function ploot(t,P)
  Pex=BS(t,s);
  clf;
  a=gca();
  a.data_bounds=[Xmin,Xmax,Ymin,Ymax]
  xtitle('t='+string(t));
  sgraph=[Smin;s;Smax];
  //Pgraph=[ug(t);P;ud(t)];
  Pgraph=[ug(t);P;ud(t)];
  Pexgraph=BS(t,sgraph);

  //- Pour un simple graphique:
  //plot(sgraph,Pgraph,'blue')

  //- Pour un graphique plus complet, comparaison avec solution BS :
  xset('thickness',2);
  //plot2d(sgraph,Pexgraph,style=1);
  plot2d(sgraph,[P0(Smin);P0(s);P0(Smax)],style=5);
  plot2d(sgraph,Pgraph,style=2);
  xset('thickness',1);
  plot2d(sgraph,zeros(sgraph),style=0);


  Title_schema=SCHEMA+'-'+CENTRAGE+' (I='+string(I)+'; N='+string(N)+')';
  //Title_ex='Exacte';
  Title_pa='Payoff'
  //legend(Title_ex,Title_pa,Title_schema);
  legend(Title_pa,Title_schema);

  //- Pour rajouter graphique d'erreur:
  //plot2d(sgraph,err_scale*(Pgraph-Pexgraph),style=5);  //- Graphique d'Erreur
  //Title_err='Erreur (scale='+string(err_scale)+')';
  //legend(Title_ex,Title_schema,Title_err);
endfunction

//--------------------
//- Maillage
//--------------------
//COMPLETER: dt, h (pas d'espace), s (maillage: vecteur colonne), 
dt=T/N; 
h=(Smax-Smin)/(I+1);  
s=Smin + (1:I)'*h; 

//- COEFFICIENT CFL
//COMPLETER
cfl=dt/h^2 * (sigma*Smax)^2; 
printf('CFL : %5.3f\n',cfl); 


//--------------------------
//- Initialiser la matrice A
//--------------------------
//COMPLETER:  matrice A , fonction q(t)

select CENTRAGE

case 'CENTRE';

  A=zeros(I,I);
  alpha=sigma^2/2 * s.^2 /h^2
  bet=r*s/(2*h)
  for i=1:I;   A(i,i) = 2*alpha(i) + r; end;
  for i=2:I;   A(i,i-1) = -alpha(i) + bet(i); end;
  for i=1:I-1; A(i,i+1) = -alpha(i) - bet(i); end;

  function y=q(t)
     y=zeros(s);
     y(1) = (-alpha(1) + bet(1))* ug(t); 
     y($) = (-alpha($) - bet($))* ud(t); 
  endfunction

case 'DROIT' ;

  A=zeros(I,I);
  alpha=sigma^2/2 * s.^2 /h^2
  bet=r*s/(h)
  for i=1:I;   A(i,i) = 2*alpha(i) + bet(i) + r; end;
  for i=2:I;   A(i,i-1) = -alpha(i) ; end;
  for i=1:I-1; A(i,i+1) = -alpha(i) - bet(i); end;

  function y=q(t)
     y=zeros(s);
     y(1) = (-alpha(1))* ug(t); 
     y($) = (-alpha($) - bet($))* ud(t); 
  endfunction

case 'GAUCHE';

  A=zeros(I,I);
  alpha=sigma^2/2 * s.^2 /h^2
  bet=r*s/(h)
  for i=1:I;   A(i,i) = 2*alpha(i) - bet(i) + r; end;
  for i=2:I;   A(i,i-1) = -alpha(i) + bet(i) ; end;
  for i=1:I-1; A(i,i+1) = -alpha(i) ; end;

  function y=q(t)
     y=zeros(s);
     y(1) = (-alpha(1) + bet(1))* ug(t); 
     y($) = (-alpha($))* ud(t); 
  endfunction

else 

  printf('CENTRAGE not programmed !'); abort

end
 
//--------------------
//- Initialiser P et graphique
//--------------------
P=P0(s);
ploot(0,P);
printf('appuyer sur la touche Enter'); 


//--------------------
//- BOUCLE PRINCIPALE
//--------------------
//- demarrage compteur temps
timer(); 

for n=0:N-1

  t=n*dt;

  //- Schema
  select SCHEMA 
  case 'EE'; 
    P = P - dt*(A*P+q(t));

  case 'EI'; 
    t1=(n+1)*dt;
    P = (eye(A) + dt*A)\(P-dt*q(t1));

  case 'CN';
    t1=(n+1)*dt;
    P = (eye(A) + dt/2*A) \ ( (eye(A) - dt/2*A) * (P - dt*(q(t)+q(t1))/2) );

  case 'EE-AMER'; 
    // COMPLETER
    P = max( P - dt*(A*P+q(t)), P0(s) );
    
  case 'EI-AMER-PSOR'; 
    // A TESTER
    if n==0; exec('PSOR.sci'); end;
    B=speye(A)+dt*A; b=P-dt*q(t+dt); x0=P; g=P0(s); eta=1.e-5; kmax=100;
    [P,k,err]=PSOR(B,b,g,x0,eta,kmax)

  case 'EI-AMER-PUL'
    if n==0
      B=speye(A)+dt*A; 
      exec('PUL.sci');
      printf('dec U.L de B:..');  [U,L]=uldecomp(B);
      printf('(verification): norme de B-U*L=%f\n',norm(B-U*L)); 
      input('appuyer sur Enter');
    end;
    b=P-dt*q(t+dt); g=P0(s);
    c=U\b; 
    P=p_descente(L,c,g);
    //- Verification:  
    v=min(B*P-b,P-g); eps=norm(v,'inf'); printf('  Verification(eps=0?): eps=%f\n',eps);

  case 'EI-AMER-NEWTON' 
    if n==0; exec('Newton.sci'); end;
    // COMPLETER
    B=speye(A)+dt*A; b=P-dt*q(t+dt); x0=P; g=P0(s); eta=1.e-5; kmax=100;
    [P,k,err]=Newton(B,b,g,x0,eta,kmax)
    //- Verification: 
    v=min(B*P-b,P-g); eps=norm(v,'inf'); printf('  Verification(eps=0?): eps=%f\n',eps);

  case 'EI-AMER-SPLIT'
    // COMPLETER
    P = (speye(A) + dt*A)\(P-dt*q(t+dt));
    P = max(P,P0(s));


  else
    printf('SCHEMA not programmed'); abort;

  end

  if modulo(n+1,deltan)==0; // Affichage tous les deltan pas.

   //- Graphiques:
   t1=(n+1)*dt; 
   ploot(t1,P);
 
   //- Calculs d'erreurs:
   //COMPLETER
   Pex=BS(t1,s);             //- Appel de la solution Black et Scholes
   errLI=norm(Pex-P,'inf');  //- Calcul erreur Linfty
   //printf('t=%5.2f; Err.Linf=%8.5f (difference avec Formule B&S)',t1,errLI);  
   printf('\n');
   //input('');
  end

end

//- ESTIMATION VALEUR DU PRIX EN Sval 

Seval=95.0;
j=find(s<=Seval); j=j($);
Peval=((s(j+1)-Seval)*P(j) +  (Seval-s(j))*P(j+1))/(s(j+1)-s(j));
printf('I=%5i, N=%5i, Seval=%5.2f: Peval=%15.12f\n',I,N,Seval,Peval);




t_total=timer()
printf('temps total = %5.2f\n',t_total);
printf('programme terminé normalement');

