mode(-1);  
clear;
//! amer.sce                                   
//! Put Americain par Diff. Finies; Jan 2008   
//!                                            
//-----------------------
//- Parametres financiers:
//-----------------------
K=100; sigma=0.2; r=0.1; T=1;  
S=100;  //- valeur cherchée : P(T,S)

//- Fonction Payoff : 'put'
function y=P0(s); y=max(K-s,0); endfunction
//function y=P0(s); y=zeros(s); i=find(50<=s & s<= 100); y(i)=K; endfunction


//-----------------------
//- Paramètres numériques
//-----------------------
Smax=2*K; 

N=10; I=20-1;

CENTRAGE='DROIT'	//- 'DROIT' ou 'CENTRE' 
MAILLAGE='UNIFORME'     //- Pas d'autres maillages programmés ici.
SCHEMA  ='EE' 		//- EE ou EI ou CN : pour l'option europeenne.


//- Variantes: 
//SCHEMA  ='AMER-EE'
//SCHEMA  ='AMER-EI-PSOR';
//SCHEMA  ='AMER-EI-PUL'
//SCHEMA  ='AMER-EI-N'
//SCHEMA  ='AMER-EI-SPLIT'
//SCHEMA  ='AMER-CN-N'; CENTRAGE='CENTRE';
//SCHEMA  ='AMER-CN-SPLIT'; CENTRAGE='CENTRE';

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

//-----------------------------------
//- Parametres d'affichage graphiques
//-----------------------------------
err_scale=0;   	//- Echelle d'erreur: NON UTILISÉE 
RECT=[0,-10,Smax,K]; 
deltan=N/10; 		//- Affichage tous les deltan pas.

//-----------------------------
//- Formule de Black et Scholes
//-----------------------------
function y=Normal(x)    
//- y = int_{-infty}^{x} exp(-x^2/2)/sqrt(2 pi).
y=cdfnor("PQ",x,zeros(x),ones(x));
endfunction

function P=BS(t,s,K,r,sigma)
//- Une formule de Black et scholes pour le put europeen vanille.
if t==0 then
  P=max(K-s,0);
else
  P=K*exp(-r*t)*ones(s); 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)*(1-Normal(dm)) - s(i).*(1-Normal(dp));
end
endfunction


//------------------------------------
//- Fonction affichage des graphiques:
//------------------------------------
function ploot(t)
  xbasc(); 
  
  //Pex=BS(t,s);				//- Solution Black et scholes
  //plot2d(s,err_scale*abs(P-Pex),5,rect=RECT); //- Graphique erreur
  plot2d(s,P0(s),1,rect=RECT);		//- Graphique Payoff
  plot2d(s,P,2,rect=RECT);		//- Graphique Schema
  plot2d(s,zeros(s),0,rect=RECT);	//- Graphique maillage
  xtitle('t='+string(t));		//- Titre et legendes
  legends(['Payoff',..
    'Schema: '+SCHEMA+'-'+CENTRAGE+' (I+1='+string(I+1)+'; N='+string(N)+')'],..
    [1,2],1);
  //  'Erreur (scale = '+ string(err_scale)+')'],

  //- Textes:
  //errLI=norm(Pex-P,'inf');   
  //select METHODE 
  //case 'DF'; errL2=sqrt(h)*norm(Pex-P); 
  //case 'EF'; errL2=sqrt(sum(h.*((Pex-P)^2))); 
  //end
  //- Estimation de P(t,S) (Cas MAILLAGE UNIFORME ici)
  hval=h; 
  Sval=S; // Valeur estimée en P(t,Sval); cas d'un maillage uniforme uniquement.
  j0=floor(Sval/hval); j1=j0+1;
  S0=j0*hval; S1=S0+hval;
  Pval=P(1+j1)*(Sval-S0)/hval + P(1+j0)*(S1-Sval)/hval; // interpolation P1.

  //printf('t=%5.2f; Erreur Linf=%8.5f, L2=%8.5f; min(P):%8.5f',t,errLI, errL2,min(P));  
  //printf('t=%5.2f; S=%8.5f, P(t,S)=%8.5f, min(P):%8.5f',t,S,Pval,min(P));  
  printf('t=%5.2f; S=%10.5f, P(t,S)=%10.5f',t,S,Pval);  

endfunction 

//--------------------------
//- Maillage espace et temps 
//--------------------------

//-Maillage uniforme  
dt=T/N; 
h=Smax/(I+1);  
s=(0:I)'*h;
printf('DF: Maillage uniforme.\n'); 
//- K dans maillage ?
if (K/h-floor(K/h))>1.e-10; printf('Warning: K pas sur maillage;\n'); end

//-----------------------------
//- Initialisation des matrices
//-----------------------------
  J=diag(sparse(ones(I,1)),1);
  a=sparse(sigma^2/2*s.^2);
  Lap=diag(a)*(2*speye(J)-J-J')/h^2;
  select CENTRAGE
  case 'DROIT' ; b=sparse(-r*s);  D=diag(b) * (J-speye(J))/h;
  case 'CENTRE'; b=sparse(-r*s);  D=diag(b) * (J-J')/(2*h);
  else; printf('this CENTRAGE not programmed'); abort
  end
  A=Lap + D + r*speye(D);

//- Initialiser P et graphique
P=P0(s); 
ploot(0)
printf('\n'); 
//printf('.. Appuyer sur la touche Return'); scanf('');

//- Boucle sur n
for n=0:N-1
 
  t1=(n+1)*dt;

  //- Schema
  select SCHEMA
  case 'EE'; P = P - dt*A*P;
  case 'EI'; P = (speye(A) + dt*A)\(P);
  case 'CN'; P = (speye(A) + dt/2*A)\ ( P - dt/2*A*P );

  case 'AMER-EE'; 
    // COMPLETER

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

  case 'AMER-EI-PUL'
    B=speye(A)+dt*A; b=P; g=P0(s);
    if n==0
      getf(PWD+'/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)); //scanf('')
    end;
    c=U\b; 
    P=p_descente(L,c,g);
    //- Verification:  
    //v=min(B*P-b,P-g); eps=norm(v,'inf'); printf('verif: eps=%f\n',eps);

  case 'AMER-EI-N' 
    if n==0; getf(PWD+'/Newton.sci'); end;
    // COMPLETER
    

    //- Verification: 
    //v=min(B*P-b,P-g); eps=norm(v,'inf'); printf('verif: eps=%f\n',eps);

  case 'AMER-EI-SPLIT'
    // COMPLETER
    
   

  case 'AMER-CN-N' 
    if n==0; getf(PWD+'/Newton.sci'); end;
    // COMPLETER
    
    

  case 'AMER-CN-SPLIT'
    // COMPLETER
    
   
  else; printf('SCHEMA non programmé'); abort;
  end


  if modulo(n+1,deltan)==0 

    ploot(t1); 

    printf('\n'); 
    //scanf('');

  end
end
printf('program finished normaly');