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

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=(1/dx^2)*D1-(2/dx^2)*D2+(1/dx^2)*D3;
Diff=Diff*sigma^2/2;

// Matrice d'advection
// centre
Adv=-1/(2*dx)*D1+1/(2*dx)*D3;
// upwind
// $$$ Adv=-(1/dx)*D2+(1/dx)*D3;
Adv=Adv*(r-sigma^2/2);

select ST
 case 'EE'
  D=Diff+Adv+(1/dt-r)*D2;
  G=(1/dt)*D2;
 case 'EI'
  D=(1/dt)*D2;
  G=-Diff-Adv+(1/dt+r)*D2;
 case 'CN'
  D=Diff/2+Adv/2+(1/dt-r/2)*D2;
  G=-Diff/2-Adv/2+(1/dt+r/2)*D2;
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
	Unew(i)=max(0,exp(xmin+(i-1)*dx)-K);
  end
 case 'put'
  for i=1:I+1
	Unew(i)=max(0,K-exp(xmin+(i-1)*dx));
  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'
	CL(I+1)=Smax-K*exp(-r*theta);
   case 'put'
	CL(1)=K*exp(-r*theta);
  end
  
  // calcul de U
  Unew(2:I)=GG\(D*Uold-G*CL);
  
  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
i=floor((log(S0)-xmin)/dx);
P=Unew(i+1)*(xmin+(i+1)*dx-log(S0))/dx+Unew(i+2)*(log(S0)-(xmin+i*dx))/dx;

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

///////////////////////////////
// Calcul de l'erreur
///////////////////////////////

erreur=0;

for i=1:I+1
  s=S(i);
  d1=(1/(sigma*sqrt(T)))*(log(s/K)+(r+sigma**2/2)*T);
  d2=d1-sigma*sqrt(T);
  select CI
   case 'call'
	prix=s*cdfnor("PQ",d1,0,1)-K*exp(-r*T)*cdfnor("PQ",d2,0,1);
   case 'put'
	prix=-(s-K*exp(-r*T))+(s*cdfnor("PQ",d1,0,1)-K*exp(-r*T)*cdfnor("PQ",d2,0,1));
  end
  if (abs(prix-Unew(i))> erreur)
	erreur=abs(prix-Unew(i));
  end
end

disp(erreur,'erreur');