mode(-1);  
clear;
//! tp.sce                                 
//! Opt Cout Transaction;       Fév 2010   


//-----------------------------------
//- Parametres financiers; Payoff
//-----------------------------------
test=1

select test
case 1
  printf('Test 1: solution exacte connue \n');
  K=100; sigma=0.2; r=0.1; T=1; Smax=200;
  eps=sigma^2/2;
  eps=0.0;
  alist=[sigma^2 - eps, sigma^2 + eps];
  function y=P0(s); y=max(K-s,0); endfunction
  function y=Pexact(t,s)
    sigma_ex=sqrt(sigma^2-eps);
    y=BS(t,s,K,r,sigma_ex);
  endfunction

case 2
  printf('test 2: comparaison avec/sans cout de transaction\n');
  K=100; sigma=0.4; r=0.1; T=0.2; Smax=200;
  //- Sans cout: retirer les commentaires. 
  alist=[sigma^2, sigma^2]; file_name='P.dat' 
  //- Avec cout: retirer les commentaires.
  //eps=sigma^2/5; alist=[sigma^2 - eps, sigma^2 + eps]; file_name='Pc.dat' 

  function y=P0(s);
     y=zeros(s);
     i=find(40<s & s<=50); y(i)=(s(i)-40);
     i=find(50<s & s<=60); y(i)=10;
     i=find(60<s & s<=70); y(i)=70-s(i);
     y=8*y;
  endfunction
  function y=Pexact(t,s)
    y=zeros(s);
  endfunction

end

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

I=100-1;
N=100; 

//CENTRAGE='DROIT'	//- 'DROIT' ou 'CENTRE' 
CENTRAGE='CENTRE'
MAILLAGE='UNIFORME'     //- Pas d'autres maillages programmés ici.
SCHEMA  ='EI'           //- Choix: 'EE' , 'EI' ,  'CN',  Cout transaction:  'EE-CT',  'EI-SPLIT-CT', 'EI-CT'
//SCHEMA  ='EE-CT' 
//SCHEMA  ='EI-SPLIT-CT' 
//SCHEMA  ='EI-CT' 

//- 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);

//-----------------------------------
//- Paramètres 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,P)

  Pbs=BS(t,s,K,r,sigma);     //-Solution BS (pour le coefficient sigma)
  Pex=Pexact(t,s);           //-Solution exacte si disponible

  xbasc();
  plot2d(s,Pex,1,rect=RECT);        //- Graphique sol exacte 
  plot2d(s,P,2,rect=RECT);          //- Graphique Schema
  //plot2d(s,Pbs,0,rect=RECT);        //- Graphique sol BS
  plot(s,Pbs,'k--');        //- Graphique sol BS
  //plot2d(s,err_scale*(P-Pex),5,rect=RECT); //- Graphique erreur

  xtitle('t='+string(t));               // Titre
  xset('font size',1);
  titre_num=SCHEMA+' ('+CENTRAGE+'; I+1='+string(I+1)+'; N='+string(N)+')';
  legend(['Exacte', titre_num,'Black et Scholes']);
  //legends(['Exacte',titre_num,'Black et Scholes'],[1,2,0],1);

             //'Erreur (scale = '+ string(err_scale)+')'],[1,2,5],1);
  plot2d(s,zeros(s),0,rect=RECT);   //- Graphique Maillage
  xset('font size',1);

  //- Textes:
  errLI=norm(Pex-P,'inf');
  errL2=sqrt(h)*norm(Pex-P);
  printf('t=%5.2f; Erreur Linf=%8.5f, L2=%8.5f; ',t,errLI, errL2);
endfunction

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

//-Maillage uniforme  
dt=T/N; 
h=Smax/(I+1);  
s=(0:I)'*h;
printf('DF: Maillage uniforme.\n'); 

//-----------------------------------
//- Initialisation des matrices A1,A2
//-----------------------------------
// COMPLETER A1,A2
// correspondant a sigma^2-eps et sigma^2+eps

for k=1:2
  alpha=alist(k);

  J=diag(sparse(ones(I,1)),1);
  a=sparse(alpha/2*s.^2);
  Lap=diag(a)*(2*speye(J)-J-J')/h^2;
  select CENTRAGE
  case 'DROIT'  then; b=sparse(-r*s);  D=diag(b) * (J-speye(J))/h;
  case 'CENTRE' then; 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);

  select k
  case 1; A1=A;
  case 2; A2=A;
  end
end


//- 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;

  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 'EE-CT'; 
      // COMPLETER
      P = min((speye(A1)-dt*A1)*P, (speye(A2)-dt*A1)*P);

  case 'EI-SPLIT-CT'; 
      // COMPLETER
      P = min((speye(A1)+dt*A1)\P, (speye(A2)+dt*A1)\P);

  case 'EI-CT'; 
      // COMPLETER en suivant la structure suivante:
      if n==0; getf Newton.sci; B1=speye(A1)+dt*A1; B2=speye(A2)+dt*A2; end;
      b=P; x0=P; eta=1.e-5; kmax=100;
      [P,k,err]=Newton(B1,B2,b,b,x0,eta,kmax)
      //- Verification eventuelle:
      //v=max(B1*P-b,B2*P-b); eps=norm(v,'inf'); printf('verif: k=%i, eps=%f\n',k,eps); 

  else; printf('SCHEMA non programmé'); abort;
  end


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

    ploot(t1); 

    printf('\n'); 

  end
end
//------------------- SAUVEGARDE POUR COMPARAISON 
if exists('file_name');
  printf('Sauvarge solution dans fichier %s\n',file_name);
  unix('rm -f '+file_name'); write(file_name,P);
end
//------------------- FIN
printf('programme terminé normalement');