black scholes sur scilab

[invite6d4611d6]

Bonjour et merci pour ce forum, j'ai déjà pas mal avancé en lisant quelques posts.

Je fais appel à votre bon coeur pour m'aider à résoudre mon problème de programmation.
Je débute dans scilab et je fais face à des problèmes insondables pour moi. Mon exercice consiste à afficher la valeur d'un dérivé financier (un call).
En fait je veux remplir une matrice qui me permet d'obtenir la valeur d'une variable au temps suivant. Je l'ai rempli plutôt niaisement (peut être trop) et donc scilab se fache.

voila le bout de programme pour remplir cette matrice:

function spA=matriceA(nbpointsx,sigma,dx,dt,nbpointsx,r)
A=zeros(nbpointsx,nbpointsx);

alpha=1+dt/(dx*dx)*(sigma2/2-(r-sigma2/2)*dx-r*dx*dx);
 mu=(dt/(dx*dx))(-sigma2+(r-sigma2/2)*dx);
 omega=dt/(dx*dx)*sigma2/2;
 
for i=1:nbpointsx-2
     A(i,i)=alpha;
     A(i,i+1)=mu;
     A(i,i+2)=omega;  
end
A(nbpointsx,nbpointsx)=(1+sigma2*dt/2/dx/dx+(r-sigma2/2)*dt/dx-r*dt);
A(nbpointsx-1,nbpointsx-1)=A(nbpointsx,nbpointsx);

A(nbpointsx,nbpointsx-1)=-sigma2*dt/dx/dx-dt/dx*(r-sigma2/2);
A(nbpointsx-1,nbpointsx-2)=A(nbpointsx,nbpointsx-1);

A(nbpointsx,nbpointsx-2)=dt/dx/dx*sigma2/2;
A(nbpointsx-1,nbpointsx-3)=A(nbpointsx,nbpointsx-2);
spA=sparse(A);
endfunction

L'erreur serait à la ligne 10 ( A(i,i+1)=mu). Moi pas comprendre
Merci pour votre attention.
-----

[invite6d4611d6]

Bon en fait je ne sais pas si le problème viens vraiment de là ou plutôt du main...
Je suis vraiment bloqué. J'ai fait le travail de recherche depuis le début.
mvt brownien, lemme d'ito, arrivée à l'equation de black scholes. Bloquer ici c'est vraiment cruel.

// Résolution de l'équation de Black-Scholes
// Calcul de u

//constantes intervenant dans l'équation
sigma=0.3;sigma2=sigma*sigma;
r=0.05;
E=4;

//paramètres numériques
xmin=log(E)/2;xmax=log(3*E);
dt=1/365;dx=sqrt(dt/5);

//tmax est le temps final 
tmax=1;
//t0 est le temps entre chaque courbe;
//nbpointst est le nombre de points utilisés pour la figure;
nbpointsx=floor((xmax-xmin)/dx)+1;
nbpointst=floor(tmax/dt)+1;
x=logspace(xmin,xmax,nbpointsx);


getf('A.sci');

// calcul et représentation de la condition initiale 
  u=zeros(nbpointsx);
for i=1:nbpointsx
  u(i)=max(exp(x(i))-E,0);
  end;
 
xset("window",0);
plot2d(x,u,2,rect=[0,0,10,5]);
xtitle(['T=1'],'x','u' );


 
//Initialisation du temps et calcul de la matrice 
t=0;i=0;
spA=matriceA(nbpointsx,sigma,sigma2,dx,dt,nbpointsx,r);
A=I*spA;

while (i<=nbpointst) do
       i=i+1;
       u=A*u;
       
// calcul de la solution à l'instant suivant      
  //     b=u;
    //   [hand,rk]=lufact(spA);
      // u=lusolve(hand,b);ludel(hand);
       //end;
//tracé des courbes lorsque on a atteint l'écheance
  
          xset("window",0); xset('line style',3);
          plot2d(x,u,2,rect=[0,0,10,5]);
end;

[inviteaa261cd4]

Tu obtiens un message d'erreur décrivant la manière dont "scilab se fâche" ?

Dans ta fonction, tu ne listes pas ton sigma2 dans la liste des arguments, mais tu le rajoutes dans l'appel à la fonction que tu fais plus bas, donc tu as dû corriger ça entre temps.

En revanche, je ne vois pas où tu définis le I auquel tu multiplies ta matrice spA après l'avoir créée.

[inviteaa261cd4]

J'ai installé rapidement scilab pour tester.

En lui demandant de m'afficher les valeurs de alpha, mu et omega avant de tenter les affectations, j'ai remarqué que mu n'est pas un nombre. Ça a du coup attiré mon attention sur le fait que tu as oublié un signe multiplicatif entre tes deux parenthèses dans le calcul de mu.

À partir de là, ton code ne bloque plus que sur le I que je t'ai signalé plus tôt.

[A voir en vidéo sur Futura]

[invite6d4611d6]

Merci pour ton aide

je modifie et je vois ce qui se passe.
en fait I=eye(nbpointsx,nbpointsx);
j'ai lu qu'il faut faire ca a partir d'une matrice creuse pour retrouver la matrice "normale".

Bon il me dit que les coordonnées de u ne sont pas des nombres, je regarde pourquoi.

[invite6d4611d6]

J'ai un peu change mon equation de base pour la rendre plus simple (equation de la chaleur) puis je fais le changement de variable inverse.
Ca marche mais je trouve que ce n'est pas precis.
Je souhaiterais savoir si c'est normal.

voici les fichiers

clear;
stacksize(10^8);
// Tracé de la solution à l'équation de black-scholes
// Calcul de C

// ne converge pas pour sigma>0.5
// Euler explicite
//constantes intervenant dans l'équation
sigma=0.3;sigma2=sigma*sigma;
r=0.05;
E=4;
k=r*2/sigma2;
T=1;

//paramètres numériques

Tau=sigma2*T/2;
rho=-((1-k)*(1-k)/2)*((1-k)*(1-k)/2);
xmin=-1; xmax=2*E;dx=(xmax-xmin)/600;

//dt dépend de dx
dt=0.1*dx*dx;//On peut changer la valeur du coefficient pour avoir une meilleure approximation


//nbpointst est le nombre de points utilisés pour la figure;
nbpointsx=floor(xmax-xmin)/dx;
nbpointst=floor(Tau/dt);


alpha=dt/dx/dx;

x=linspace(xmin,xmax,nbpointsx)';
i=0;
for i=1:nbpointsx

end;
S=E*exp(x);
exec('A.sci');

// calcul et représentation de la condition initiale 
   
   u=max(exp((k+1)*x/2)-exp((k-1)*x/2),0);
   C=zeros(nbpointsx);
   i=0;
   for i=1:nbpointsx
     C(i)=E*exp(-x(i)*(1-k)/2)*u(i);
     end;
// affichage des valeurs du call C suivant le prix de l'action S 
xset("window",0);
plot2d(S(1:nbpointsx),C(1:nbpointsx),2,rect=[0,0,6,0.5]);
xtitle(['T=1'],'S','C' );


 

A=matriceA(nbpointsx,alpha);

//calcul d'un taux

i=0;
while (i*dt<=Tau) do
  i=i+1;
  u=A*u;
end;

//changement de variable x->C
for i=1:nbpointsx
  C(i)=E*exp(-(1-k)/2*x(i)+rho*sigma2*T/2)*u(i);
  end;

//tracé de la courbe explicite lorsque on a atteint l'écheance 

          xset("window",0); xset('line style',3);
          plot2d(S(1:nbpointsx),C(1:nbpointsx),2,rect=[0,0,6,0.5]);

ainsi que ceci
a appeler A.sce

function A=matriceA(nbpointsx,alpha)
  A=zeros(nbpointsx,nbpointsx);
  
  A(1,1)=1-2*alpha; A(1,2)=alpha;
  A(nbpointsx,nbpointsx-1)=alpha; A(nbpointsx,nbpointsx)=1-2*alpha;
for i=2:nbpointsx-1
     A(i,i)=1-2*alpha;
     A(i,i+1)=alpha;
     A(i,i-1)=alpha;  
end

endfunction

Merci à vous

[inviteaa261cd4]

Tu as tenté d'augmenter le nombre de points en réduisant dx (donc en augmentant ton dénominateur ligne 19 ?

Je sais pas ce que tu es censé obtenir au final, mais en passant de 600 à 1200, le "démarrage" de S -> C (première courbe affichée) est déjà plus net. (Au prix d'une durée de calcul multipliée par 16, mais ça reste encore raisonnable sur mon système)

[invite6d4611d6]

Oui c'est mieux en effet. Je n'ai pas un pc puissant alors j'avais arreté le programme avant qu'il affiche la courbe.
J'aurais peut être du utiliser une autre méthode approximation: Crank.
Je dois le rendre malheureusement.
merci Kiv