mode(-1);  
clear;
//! ef_Q.sce A COMPLETER                                
//! Black et Scholes par Diff. Finies ou Éléments finis 

//-----------------------
//- Parametres financiers
//-----------------------
K=100; sigma=0.05; r=0.1; T=1;  	//- Diffusion + Transport

//- Fonction Payoff : 'put'
function y=P0(s); y=max(K-s,0); endfunction


//-----------------------
//- Paramètres numériques
//-----------------------
Smax=2*K; 
N=100; I=100-1;

METHODE ='DF'; CENTRAGE='CENTRE'; //- CENTRAGE= 'DROIT' ou 'CENTRE' , pour  le cas METHODE='DF'.
//METHODE='EF'; CENTRAGE='';
SCHEMA  ='CN';
MAILLAGE='UNIFORME'	//- 'UNIFORME' ou 'LOG' (cas METHODE='EF', sinon inutile.)
//MAILLAGE='LOG'		//- 'UNIFORME' ou 'LOG' (cas METHODE='EF', sinon inutile.)

//- Affichage des donnees:
printf('METHODE:%s  ',METHODE)
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=1000;   	//- Echelle d'erreur
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,Pex,1,rect=RECT);		//- Graphique BS
  plot2d(s,P,2,rect=RECT);		//- Graphique Schema
  xtitle('t='+string(t));		//- Titre et legendes
  legends(['Exacte',..
    'Approche: '+METHODE+'-'+SCHEMA+'-'+CENTRAGE+' (I+1='+string(I+1)+'; N='+string(N)+')',..
    'Erreur (scale = '+ string(err_scale)+')'],[1,2,5],1);

  //- Textes:
  errLI=norm(Pex-P,'inf');   
  select METHODE 
  case 'DF'; errL2=sqrt(h)*norm(Pex-P); 
  case 'EF'; errL2=0; // COMPLETER
  end
  printf('t=%5.2f; Erreur Linf=%8.5f, L2=%8.5f; min(P):%8.5f',t,errLI, errL2,min(P));  
endfunction 

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

case 'EF'  //-Maillage uniforme ou log.
  dt=T/N;
  select MAILLAGE           //- Maillage total  stotal = [s_0,...s_{I+1}]:
  case 'UNIFORME'
    stotal=(0:I+1)'/(I+1)*Smax;
  case 'LOG' 
    I1=(I+1)/2; 
    s1=K*log(1:I1)'/log(I1+1)
    s2=Smax- (Smax-K)*log(I1+1:-1:1)'/log(I1+1);
    stotal=[s1;s2];
    xbasc(); plot2d(stotal,zeros(stotal),-1,rect=RECT);
    scanf('')
  else; printf('MAILLAGE mal defini\n'); abort;
  end
  //- Definition des vecteurs s et h:
  h=stotal(2:I+2)-stotal(1:I+1); //- h=[h_0,..,h_{I}],  h_i=s_{i+1}-s_{i} 
  s=stotal(1:I+1)                //- s=[s_0,..,s_I] 
  
  //- Verification K dans maillage:
  i=find(abs(K-s)<1.e-10); 
  if i==[];
    printf('Warning: K pas sur maillage !\n'); //scanf('%c') 
  end


end

//-----------------------------
//- Initialisation des matrices
//-----------------------------
select METHODE
case 'DF' 
  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'  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);

case 'EF' 
  // COMPLETER M
  M=spzeros(I+1,I+1);



  A=spzeros(I+1,I+1);
  A(1,1)=r/2*h(1);
  for i=1:I;  j=i+1; A(j,j)=(s(j)*sigma)^2/2*(1/h(j-1)+1/h(j)) + r/2*(h(j-1)+h(j)); end
  for i=1:I;  j=i+1; A(j,j-1)=-(s(j)*sigma)^2/(2*h(j-1)) + r/2*s(j); end;
  for i=0:I-1;j=i+1; A(j,j+1)=-(s(j)*sigma)^2/(2*h(j)) - r/2*s(j); 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;

  //- Schema
  select METHODE
  case 'DF'		//-  P' + A P = 0
    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 );
    else; printf('SCHEMA non programmee'); abort;
    end
  case 'EF'   		//-  M P' + A P = 0
    select SCHEMA
    case 'EE'; P = P; // COMPLETER EE
    case 'EI'; P = P; // COMPLETER EI
    case 'CN'; P = P; // COMPLETER CN
    else; printf('SCHEMA non programmee'); abort;
    end
  else printf('METHODE not programmed');abort
  end


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

    ploot(t1); 

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

  end
end
printf('program finished normaly');