// TP 4
// Méthode des differences finies pour BS.
// 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='CN';
I=1000; // nombres de pas en espace
N=100; // nombres de pas de temps
Smax=5*K; // valeur maximum du sous-jacent
dS=Smax/I;
dt=T/N;

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

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

vec_S=dS:dS:(Smax-dS);
vec_SS=vec_S^2;

DS1=sparse([1:(I-1);1:(I-1)]',vec_S,[I-1,I+1]);
DS2=sparse([1:(I-1);2:I]',vec_S,[I-1,I+1]);
DS3=sparse([1:(I-1);3:(I+1)]',vec_S,[I-1,I+1]);

DSS1=sparse([1:(I-1);1:(I-1)]',vec_SS,[I-1,I+1]);
DSS2=sparse([1:(I-1);2:I]',vec_SS,[I-1,I+1]);
DSS3=sparse([1:(I-1);3:(I+1)]',vec_SS,[I-1,I+1]);

D2=sparse([1:(I-1);2:I]',ones(I-1,1),[I-1,I+1]);

// Matrice de diffusion
Diff=(1/dS^2)*DSS1-(2/dS^2)*DSS2+(1/dS^2)*DSS3;
Diff=Diff*sigma^2/2;

// Matrice d'advection
// centre
Adv=-1/(2*dS)*DS1+1/(2*dS)*DS3;
// upwind
// $$$  Adv=-(1/dS)*DS2+(1/dS)*DS3;
Adv=r*Adv;

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
Pnew=zeros(I+1,1);

select CI
 case 'call'
  for i=1:I+1
	Pnew(i)=max(0,(i-1)*dS-K);
  end
 case 'put'
  for i=1:I+1
	Pnew(i)=max(0,K-(i-1)*dS);
  end
end
Pold=Pnew;

// 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=0:dS:Smax;

///////////////////////////////
// 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 P
  Pnew(2:I)=GG\(D*Pold-G*CL);
  
  Pnew(1)=CL(1);
  Pnew(I+1)=CL(I+1);
  
  // affichage
  if (modulo(floor(theta/dt),N/50)==0) then
	xbasc();
	plot2d (S,Pnew,rect=Rect);
  end
  
  Pold=Pnew;
end // for en temps

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

// on interpole pour avoir le prix en S0
i=floor(S0/dS);
P=Pnew(i+1)*((i+1)*dS-S0)/dS+Pnew(i+2)*(S0-i*dS)/dS;

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

disp(abs(P-prix),'erreur');