// TP 4
// Méthode des differences finies pour BS.
// On passe en variable x=ln(S).
// La grille est uniforme en espace et en temps.

stacksize(1.e7);
clear
xbasc()

///////////////////////////
// Parametres financiers
///////////////////////////

// payoff : 'put' ou 'call'
CI='call';

T=1.;
r=0.1;
S0=100;
K=100;
sigma=0.05;

///////////////////////////
// Parametres numeriques
///////////////////////////

// Schéma en temps : 'EE', 'EI' ou 'CN'
ST='EE';
I=500; // nombres de pas en espace
N=100; // nombres de pas de temps
Smax=5*K; // valeur maximum du sous-jacent
xmax=log(5*K);
xmin=-log(5*K);
dx=(xmax-xmin)/I;
dt=T/N;

/////////////////////////
// Matrices
/////////////////////////

disp('Definition des matrices...');

// Matrices diagonales, diagonale inférieure, diagonale centrée,
// diagonale supérieure.
D1=sparse([1:(I-1);1:(I-1)]',ones(I-1,1),[I-1,I+1]);
D2=sparse([1:(I-1);2:I]',ones(I-1,1),[I-1,I+1]);
D3=sparse([1:(I-1);3:(I+1)]',ones(I-1,1),[I-1,I+1]);

// Matrice de diffusion Diff
// a compléter...

// Matrice d'advection Adv
// a compléter...

// Construction des matrices D et G
select ST
 case 'EE'
// a compléter...
 case 'EI'
// a compléter...
 case 'CN'
// a compléter...
end

// Extraction de la matrice carrée
GG=G(1:I-1,2:I);

// Condition initiale
Unew=zeros(I+1,1);

select CI
 case 'call'
  for i=1:I+1
	// a compléter...
  end
 case 'put'
  for i=1:I+1
	// a compléter...
  end
end
Uold=Unew;

// vecteur et fonctions conditions aux limites
CL=zeros(I+1,1);

// affichage
select CI
 case 'call'
  Rect=[0,-1,Smax,Smax-K*exp(-r*T)];
 case 'put'
  Rect=[0,-1,Smax,K];
end

S=exp(xmin:dx:xmax);

///////////////////////////////
// Boucle en temps
///////////////////////////////

disp('Boucle en temps...');

for theta=dt:dt:T
  
  // conditions aux limites au temps theta
  select CI
   case 'call'
	// a compléter...
   case 'put'
	// a compléter...
  end
  
  // calcul de U
  // a compléter...
  
  Unew(1)=CL(1);
  Unew(I+1)=CL(I+1);
  
  // affichage
  if (modulo(floor(theta/dt),N/50)==0) then
	xbasc();
	plot2d (S,Unew,rect=Rect);
  end
  
  Uold=Unew;
end // for en temps

xbasc();
plot2d (S,Unew,rect=Rect);
 
///////////////////////////////
// Affichage du prix
///////////////////////////////

// on interpole pour avoir le prix en S0
// a compléter...

disp(P,'prix calculé')

// calcul du vrai prix
d1=(1/(sigma*sqrt(T)))*(log(S0/K)+(r+sigma**2/2)*T);
d2=d1-sigma*sqrt(T);
select CI
 case 'call'
  prix=S0*cdfnor("PQ",d1,0,1)-K*exp(-r*T)*cdfnor("PQ",d2,0,1);
 case 'put'
  prix=-(S0-K*exp(-r*T))+(S0*cdfnor("PQ",d1,0,1)-K*exp(-r*T)*cdfnor("PQ",d2,0,1));
end

disp(prix,'prix exact');