équation dynamique mathématique

Version imprimable

Bonjour à tous, j'aimerais programmer une équation dynamique en C++ qui est:
f(0,x)=[K-exp(x)]+ pour x<u;
f(i,x)=0 pour x>=u;
f(i+1,x)=pn*f(i,x+sigma*sqrt(delta))+(1-pn)*f(i,x-sigma*sqrt(delta));
voila mon code:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
  double partiepositive (double k,double e) // renvoie la partie positive de k-e
{
  if (k>e) return k-e;
  else 0;
}
 
double phi(double u, double K,double  T,double  sigma,double N,double r,double i,double x)//equiation dynamique
{
    double delta=T/N;
    double res;
    double mu=r-sigma*sigma/2;
    double pn=0.5+mu*sqrt(delta)/(2*sigma);
    res=partiepositive(K,exp(x));
     double w1=x+sigma*sqrt(delta);
   double w2=x-sigma*sqrt(delta);
 
    if((i==0)&&(x<u))
      return res;
    else
      {
    if((i==0)&&(x>=u))
      {return 0;}
    else
      {
        if(x>=u)
          return 0;
        else
          {   
        return pn*phi(u, K, T, sigma,N,r, i-1,w1)+(1-pn)*phi(u, K, T, sigma,N,r, i-1,w2);
          }
      }
      }
}

mon problème est que l'algorithme est très lent!!! à partir de i=27..
avez-vous des améliorations pour y remédier??
merci d'avance.

21/03/2007, 16h54
Ulmo

Tiens, un financier...

Pense à la balise code...

Et stocke plutôt les valeurs de f(.,.) dans un tableau. Actuellement tu calcules 2 fois f(i-2, x) au lieu d'une, par exemple (un fois pour x-sigma*sqrt(delta)+sigma*sqrt(delta) et une fois pour x+sigma*sqrt(delta)-sigma*sqrt(delta)).
21/03/2007, 17h09
souviron34

hou que c'est pas beau....

Bon c'est OK avec le langage, mais c'est pas esthéthique, pour du code de maths...

Déclare tes variables, et fait tes calculs.

Bon, ça c'est une chose (mais qui n'engage que moi)

Par contre, déjà tu calcules je sais combien de fois sqrt(delta)... Stockes le dans une variable.

Même chose pour sigma*sqrt(delta)

Au lieu de diviser par 2, mutiplies par 0.5

Ton K = exp, de même que toutes tes initialisations (et donc les calculs), tu les calcules à tous les coups, sans tester par rapport à u, alors que ton algo le dis

voilà pour l'instant
21/03/2007, 19h00
hengstepha

Citation:

Envoyé par Ulmo

Tiens, un financier...

Pense à la balise code...

Et stocke plutôt les valeurs de f(.,.) dans un tableau. Actuellement tu calcules 2 fois f(i-2, x) au lieu d'une, par exemple (un fois pour x-sigma*sqrt(delta)+sigma*sqrt(delta) et une fois pour x+sigma*sqrt(delta)-sigma*sqrt(delta)).

tu veux que je stocke tous les valeurs de f(,) ou bien ce que les f(,) que je calcule plusieurs fois.
21/03/2007, 19h05
hengstepha

Citation:

Envoyé par Ulmo

Tiens, un financier...

Pense à la balise code...

Et stocke plutôt les valeurs de f(.,.) dans un tableau. Actuellement tu calcules 2 fois f(i-2, x) au lieu d'une, par exemple (un fois pour x-sigma*sqrt(delta)+sigma*sqrt(delta) et une fois pour x+sigma*sqrt(delta)-sigma*sqrt(delta)).

je reformule...
je stocke tous les f(,) (je pense que c'est faisable mais faudrait un LONG tableau,ou des pointeurs de pointeurs.............. de pointeurs) ou seulement les f(,) que je calcule plusieurs fois( là c'est trop chaud pour un pauvre petit financier...)

ps: pas mal l'antithèse?
22/03/2007, 11h14
Ulmo

Pour calculer f(n, x), tu as besoin d'une valeur de f(n, .) (en x).
Donc de 2 valeurs de f(n-1, .).
Donc de 3 valeurs de f(n-2, .) puisque les incréments sont les mêmes vers le haut et vers le bas.
...
Donc de N+1 valeurs de f(0, .) que tu connais explicitement.

Il suffit donc d'un tableau de taille n+1, mais il faut faire attention à l'ordre de remplissage pour n'écraser que les valeurs de l'étape précédente dont on n'a plus besoin.

Au passage, tu as une formule explicite avec la formule du binôme, mais ce n'est peut-être pas le but...
23/03/2007, 15h51
hengstepha

Citation:

Envoyé par Ulmo

Pour calculer f(n, x), tu as besoin d'une valeur de f(n, .) (en x).
Donc de 2 valeurs de f(n-1, .).
Donc de 3 valeurs de f(n-2, .) puisque les incréments sont les mêmes vers le haut et vers le bas.
...
Donc de N+1 valeurs de f(0, .) que tu connais explicitement.

Il suffit donc d'un tableau de taille n+1, mais il faut faire attention à l'ordre de remplissage pour n'écraser que les valeurs de l'étape précédente dont on n'a plus besoin.

Au passage, tu as une formule explicite avec la formule du binôme, mais ce n'est peut-être pas le but...

a ok en faite on oublie la forme itérative ( on utilise pas l'équation dynamique). et on essai de trouver la forme générale?? si j'ai bien compris c'est bien cela?
23/03/2007, 15h59
Ulmo

Quand on peut c'est mieux, on gagne du temps. Ici il suffit de compter combien de fois il y aura f(0, x+k*h) dans f(n, x).
23/03/2007, 17h32
salseropom
Citation:
Envoyé par hengstepha

Code:

1 2 3 4 5 double partiepositive (double k,double e) // renvoie la partie positive de k-e { if (k>e) return k-e; else 0; }
Salut, ta fonction partiepositive marche, mais elle renvoie des warning. Je te propose la chose suivante
Code:

1 2 3 4 5 double partiepositive (double k,double e) // renvoie la partie positive de k-e { return (k>e) ? k-e : 0.; }
et peut-être que le mieux serait de faire une macro
Code:

1 2 #define MAX(a,b) ((a>b) ? (a) : (b))
et ensuite de faire partout, à la place de partiepositive(a,b) de mettre
Code:

1 2 MAX(k-e,0.);
je ne pense pas que ça ira bcp plus vite mais au moins tu auras enlevé qq warning lors de la compilation et je pense que le code sera plus lisible...:roll:
23/03/2007, 18h49
hengstepha

Citation:

Envoyé par Ulmo

Quand on peut c'est mieux, on gagne du temps. Ici il suffit de compter combien de fois il y aura f(0, x+k*h) dans f(n, x).

ok
voila mais changement : j'ai sépare le cas pair et impair mais je ne trouve pas de bon résulats...

double fact(double n)
{
if (n<=1)
return 1;
else
return n*fact(n-1);
}

double binom(double n, double p)
{
return fact(n)/(fact(p)*(fact(n-p)));
}

double phiimpaire(double u, double K,double T,double sigma,int N,double r,double x)
{
int k=-1;
double * px=new double[N+1];
double delta=T/N;
double mu=r-sigma*sigma/2;
double pn=0.5+mu*sqrt(delta)/(2*sigma);
double j=sigma*sqrt(delta);
double res;

for (int i=0; i<floor(N/2)+1;i++)
{
px[i]=binom(N,i)*pow(pn,N-i)*pow(1-pn,i)*phiint(u,K,x+(N-2*i)*j);
}
for (int i=(int)(floor(N/2)+1);i<N+1;i++)
{
k++;
px[i]=binom(N,i)*pow(pn,N-i)*pow(1-pn,i)*phiint(u,K,x-(1+2*k)*j);
}
for(int i=0; i<N+1;i++)
{
res+=px[i];
}
delete [] px;
return res;

}

double phipaire(double u, double K,double T,double sigma,int N,double r,double x)
{
int k=0;
double *px=new double[N+1];
double delta=T/N;
double mu=r-sigma*sigma/2;
double pn=0.5+mu*sqrt(delta)/(2*sigma);
double j=sigma*sqrt(delta);
double res;

for (int i=0; i<N/2;i++)
{
px[i]=binom(N,i)*pow(pn,N-i)*pow(1-pn,i)*phiint(u,K,x+(N-2*i)*j);
}
px[N/2]=binom(N,N/2)*pow(pn,N-N/2)*pow(1-pn,N/2)*phiint(u,K,x);
for (int i=N/2+1;i<N+1;i++)
{
k++;
px[i]=binom(N,i)*pow(pn,N-i)*pow(1-pn,i)*phiint(u,K,x-(2*k)*j);
}
for(int i=0; i<N+1;i++)
{
res+=px[i];
}
delete [] px;
return res;

}
si tu peux reperer l'erreur sa serait cool :D
23/03/2007, 19h08
hengstepha

Citation:

Envoyé par Ulmo

Quand on peut c'est mieux, on gagne du temps. Ici il suffit de compter combien de fois il y aura f(0, x+k*h) dans f(n, x).

j'ai remarqué que ma fonction binomiale marcche pour des petites valeurs maisne marchait pas pour les grandes par exemple n=400, c'est normal?

La première chose à faire est d'éliminer les calculs redondants (sqrt) ou à priori pas forcéments utiles (exp). Par exemple :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
double calc(double u, double K,double  SigmaDel,double pn, double i,double x)//equiation dynamique
{
    if(x>=u) return 0;
    if(i==0)
    {
      double e = exp(x);
      return K>e? K-e:0;
    }
    return pn*calc(u, K,SigmaDel,pn,i-1,x+SigmaDel)+(1-pn)*calc(u, K,SigmaDel,pn, i-1,x-SigmaDel);
}
double phi(double u, double K,double  T,double  sigma,double N,double r,double i,double x)//equiation dynamique
{
    double sq = sqrt(T/N);
    double pn=0.5+(r-sigma*sigma/2)*sq/(2*sigma);
    return calc(u,K,sigma*sq,pn,i,x);
}

Ensuite, il faut abandonner la structure recursive pour une structure itérative

Citation:
Envoyé par diogene
La première chose à faire est d'éliminer les calculs redondants (sqrt) ou à priori pas forcéments utiles (exp). Par exemple :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 double calc(double u, double K,double SigmaDel,double pn, double i,double x)//equiation dynamique { if(x>=u) return 0; if(i==0) { double e = exp(x); return K>e? K-e:0; } return pn*calc(u, K,SigmaDel,pn,i-1,x+SigmaDel)+(1-pn)*calc(u, K,SigmaDel,pn, i-1,x-SigmaDel); } double phi(double u, double K,double T,double sigma,double N,double r,double i,double x)//equiation dynamique { double sq = sqrt(T/N); double pn=0.5+(r-sigma*sigma/2)*sq/(2*sigma); return calc(u,K,sigma*sq,pn,i,x); }
Ensuite, il faut abandonner la structure recursive pour une structure itérative
salut a toi
j'ai essayé ce que tu as donné mais pour un N >27 c'est très long....

Citation:
Envoyé par diogene
La première chose à faire est d'éliminer les calculs redondants (sqrt) ou à priori pas forcéments utiles (exp). Par exemple :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 double calc(double u, double K,double SigmaDel,double pn, double i,double x)//equiation dynamique { if(x>=u) return 0; if(i==0) { double e = exp(x); return K>e? K-e:0; } return pn*calc(u, K,SigmaDel,pn,i-1,x+SigmaDel)+(1-pn)*calc(u, K,SigmaDel,pn, i-1,x-SigmaDel); } double phi(double u, double K,double T,double sigma,double N,double r,double i,double x)//equiation dynamique { double sq = sqrt(T/N); double pn=0.5+(r-sigma*sigma/2)*sq/(2*sigma); return calc(u,K,sigma*sq,pn,i,x); }
Ensuite, il faut abandonner la structure recursive pour une structure itérative
merci du conseil mais le probleme pour passer à la forme itérative c'est que ici elle dépend de calc(u, K,SigmaDel,pn,i-1,x+SigmaDel) et de
calc(u, K,SigmaDel,pn,i-1,x-SigmaDel) pas très simple..

26/03/2007, 15h38
Ulmo

Citation:

Envoyé par hengstepha

j'ai remarqué que ma fonction binomiale marcche pour des petites valeurs maisne marchait pas pour les grandes par exemple n=400, c'est normal?

Tu as des des erreurs de débordement : 400! ne tient pas dans un double. Tu ne peux aller que jusqu'à 170.
De toute façon, à partir de N = 20 tu peux remplacer efficacement l'arbre par une variable gaussienne. Tu as alors un calcul presque explicite (la formule de Black & Scholes).
26/03/2007, 16h29
souviron34

Citation:

Envoyé par Ulmo

Tu as des des erreurs de débordement : 400! ne tient pas dans un double. Tu ne peux aller que jusqu'à 170.
De toute façon, à partir de N = 20 tu peux remplacer efficacement l'arbre par une variable gaussienne. Tu as alors un calcul presque explicite (la formule de Black & Scholes).

Complètement faux.. Il y a tout un thread sur le sujet :

http://www.developpez.net/forums/sho...d.php?t=278583

Ne confondez pas nombre de chiffres et valeur....

Sur des sytèmes 32 bits (double = 64 bits), on peut aller jusqu'à :

factorielle(1.97969 e+308)

(voir post bas de page 2..)
26/03/2007, 17h12
Ulmo

Citation:

Envoyé par souviron34

Complètement faux.. Il y a tout un thread sur le sujet :

http://www.developpez.net/forums/sho...d.php?t=278583

Ne confondez pas nombre de chiffres et valeur....

Sur des sytèmes 32 bits (double = 64 bits), on peut aller jusqu'à :

factorielle(1.97969 e+308)

(voir post bas de page 2..)

Plus précisemment jusqu'à :
factorielle(171) = 1.97969 e+308 ;)

Edit : D'ailleurs 171! ne tient pas dans un double, on s'arrête bien à 170!.
Et à 1754! avec des long double.
27/03/2007, 01h24
souviron34

milles excuses :oops: :hola:

j'avais affiché ceci en faisant une estimation rapide par dychotomie (pas en faisant le calcul complet)...

Si le sujet est toujours d'actualité, tu peux essayer ce code :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
double valeur(double K, double x)
{
     double e = exp(x);
     return K>e? K-e:0;
}
// Les coefficients du binôme
long cnp(int n,int p)
{
 long C = 1;
 int i;
 if(p<0) return 0;
 for(i=1; i<=p; i++)C= C*(n-i+1)/i;
 return C ;
}
double calc2(double u, double K,double  SigmaDel,double pn, double n,double x)
{
  double del = (x-u)/SigmaDel/2;
  double expo = pow(pn,n);
  double fact = (1-pn)/pn;
  int q = n/2;
  int i;
  long coef;
  int m0;
  int i0;
  double result=0;
  if(x>=u) return 0;
  if(n==0)return valeur(K,x);
  i0 = ceil(n/2.0+del);
  i0 = i0<0 ? 0 : i0;
  for(i = i0; i <= n; i++)
  {
   m0 = ceil(n-i+2*del);
   m0 = m0<0? 0 : m0;
   coef =0;
   if(i<n-q  && m0<=i) coef = cnp(n,i)-cnp(n,m0-1);
   if(i>=n-q && m0<=n-i) coef = cnp(n,n-i)-cnp(n,m0-1);
   if(coef!=0)result += expo*valeur(K,x+(n-2*i)*SigmaDel)* coef;
   expo = expo*fact;
  }
  return result;
 
}
double phi(double u, double K,double  T,double  sigma,double N,double r,long i,double x)
{
    double sq = sqrt(T/N);
    double pn=0.5+(r-sigma*sigma/2)*sq/(2*sigma);
    return calc2(u,K,sigma*sq,pn,i,x);
}

D'après qq tests, il semble coller et être beaucoup plus rapide. On peut encore optimiser un peu ce code .