Python exemple Loi normale centrée réduite Algorithme

[Sonny Sachs]

Bonjour,

Tout d'abord, je précise que je suis sous python 3.x

Je dois appliquer un nombre aléatoire de la loi normale centrée réduite à une variable.

Cette variable, dans mon script, est telle que :

x=1+2
N(x)

N=loi normale
x=variable

N est une loi normale centrée réduite tel que N(0,1)

Soit :

mu=0
sigma=1

J'ai cherché une fonction toute faite, j'ai trouvé dans un premier temps :

random.gauss(0,1)

Mais ça me dit que le module random ne trouve pas gauss... Pourtant elle existe : http://docs.python.org/3.1/library/random.html

Donc je me suis penché sur deux lib apparemment fonctionnelles : NumPy et SciPy.

http://www.scipy.org/install.html

Sauf que cela ne m'intéresse pas vraiment car je dois les installer sur mon ordinateur, or j'aime bien python pour sa portabilité donc préciser sur mon script de devoir installer une librairie me semble paradoxal.

Donc bon, en gros, je vais devoir me retaper l'algo de la loi normale standard. C'est pas bien grave !

Par souci d'efficacité, je cherche un code déjà fait, et je suis tombé sur un code de Tyrtamos que voici :

from __future__ import division
 
from math import *
 
def pgaussred(x):
    """pgaussred(x): probabilité qu'une variable aléatoire distribuée selon une loi normale réduite soit inférieure à x"""
    if x==0:
        return 0.5
    if x>=7.56:
        return 1.0
    if x<=-7.56:
        return 0.0
    u=abs(x)
    n=int(u*2000)
    du=u/n
    k=1/sqrt(2*pi)
    u1=0
    f1=k
    p=0.5
    for i in xrange(0,n):
        u2=u1+du
        f2=k*exp(-0.5*u2*u2)
        p=p+(f1+f2)*du*0.5
        u1=u2
        f1=f2
    if x>0:
        return p
    else:
        return 1.0-p

# exemple d'utilisation

print pgaussred(-8)

Que vous trouverez sur cette page :

http://python.jpvweb.com/mesrecettes...id=loi_normale

Sauf que :

Il a bien mis mu=0 (x==0) mais j'ai l'impression que son sigma=7.56 et non 1 (enfin je crois, si j'ai bien capté son code, sinon d'où sort ce 7.56 ?)

De plus, quand je tente de l'associer à mon code, j'ai une erreur. Je crois que c'est parce que sa fonction pgaussred() requiert un nombre et que moi c'est une variable (mais cette variable est normalement un nombre à ce stade du script ?). En effet, comme j'ai précisé plus haut, ça fait :

x=1+2
pgaussred(x)
#erreur

Par ailleurs j'ai pas vraiment bien saisi tout son code... Il aurait dû (pu) mieux le commenter...

Merci !

[fred]

Pour la fonction gauss, normalement si tu importes bien ton module random,

>>> from random import gauss
>>> gauss(0,1)
0.08263973494981278

ça doit fonctionner, non?

[tyrtamos]

Bonjour,

Si tu hésites avec random.gauss (qui fonctionne très bien), le code dont tu parles n'est pas celui que tu cherches. En effet, sont différents:
- générer une variable aléatoire distribuée selon la loi normale réduite,
- calculer la probabilité qu'une variable aléatoire distribuée selon la loi normale centrée réduite soit inférieure à une valeur donnée.

Si c'est le 1er cas, j'ai un code qui fait ça en bas de page, mais utilise plutôt random.gauss: il est probable qu'il soit plus rapide.


Je réponds tout de même à tes questions concernant le code dont tu parles dans ton message.

On a la loi normale réduite (formule mathématique) qui est une "loi de densité de probabilité", et on cherche la probabilité qu'une variable aléatoire soit inférieure à une valeur donnée x. Comme la probabilité est une surface et qu'il n'y a aucune formule exacte pour la calculer, on fait un calcul approché: on calcule 2000 tranches fines verticales par écarts-type et on additionne les surfaces correspondantes.

Comme indiqué dans le texte de présentation, les 7.56 correspondent au fait que, au delà de cette valeur de x, la probabilité est tellement proche de 1 qu'il n'est pas utile de continuer à calculer. D'ailleurs, avec mon code, on trouve 0.9999999999999981 pour x=7.55.

Enfin, mon code est fait pour Python 2. Il y a 2 corrections à faire pour Python 3: xrange => range et print est une fonction avec parenthèses. Voilà la version corrigée pour Python 3:

#!/usr/bin/python
# -*- coding: utf-8 -*-
# Python 3
 
from math import *

def pgaussred(x):
    """pgaussred(x): probabilité qu'une variable aléatoire distribuée selon une loi normale réduite soit inférieure à x"""
    if x==0:
        return 0.5
    if x>=7.56:
        return 1.0
    if x<=-7.56:
        return 0.0
    u=abs(x)
    n=int(u*2000)
    du=u/n
    k=1/sqrt(2*pi)
    u1=0
    f1=k
    p=0.5
    for i in range(0,n):
        u2=u1+du
        f2=k*exp(-0.5*u2*u2)
        p=p+(f1+f2)*du*0.5
        u1=u2
        f1=f2
    if x>0:
        return p
    else:
        return 1.0-p

print(pgaussred(1.0))
# affiche: 0.841344741027

A ta disposition si tu as d'autres questions.

[Sonny Sachs]

Bonjour,

Voici mon code (qui ne fonctionne pas, mais qui devrait ?!!) :

#!/usr/local/bin/python
# -*- coding: Utf-8 -*-

from math import log, exp, sqrt
from random import gauss

S=float(input("S="))
K=float(input("K="))
r=float(input("r="))
v=float(input("v="))
T=float(input("T="))

d1=float((log(S/K)+(r+v*v/2)*T)/(v*sqrt(T)))
d2=float(d1-v*sqrt(T))

N=float(gauss(0,1))

c=float(S*N(d1)-K*exp(-r*T)*N(d2))
#TypeError: 'float' object is not callable

p=float(-S*N(-d1)+K*exp(-r*T)*N(-d2))
#TypeError: 'float' object is not callable

print(c,p)

J'ai cette erreur j'ignore pourquoi, il n'y a pas de faute dans les formules et toutes les variables se sont déjà vues attribuer une valeur.

Ensuite, je ne sais pas si gauss(0,1) correspond bien à ma demande, mais je crois que si. Voir plus bas la formule en question.

J'ai une question bête :

Dans la formule que je dois appliquer, on voit qu'il y a 2 fois N pour calculer c et 2 fois N pour calculer p. Je me demande si les valeurs doivent changer pour chaque N ou si ça doit rester la même pour les 4 ?

Voici la formule que je dois appliquer :

http://fr.wikipedia.org/wiki/Modèle_..._Black-Scholes

Voici un exemple d'application de la formule en plusieurs langages dont en Python :

http://www.espenhaug.com/black_scholes.html

Son code est donc :

from math import *
# Cumulative normal distribution

def CND(X):

    (a1,a2,a3,a4,a5) = (0.31938153, -0.356563782, 1.781477937, 

     -1.821255978, 1.330274429)
    L = abs(X)

    K = 1.0 / (1.0 + 0.2316419 * L)

    w = 1.0 - 1.0 / sqrt(2*pi)*exp(-L*L/2.) * (a1*K + a2*K*K + a3*pow(K,3) +

    a4*pow(K,4) + a5*pow(K,5))
    if X<0:

        w = 1.0-w

    return w
    
# Black Scholes Function

def BlackScholes(CallPutFlag,S,X,T,r,v):

    d1 = (log(S/X)+(r+v*v/2.)*T)/(v*sqrt(T))

    d2 = d1-v*sqrt(T)
    if CallPutFlag=='c':

        return S*CND(d1)-X*exp(-r*T)*CND(d2)

    else:

        return X*exp(-r*T)*CND(-d2)-S*CND(-d1)

Apparemment son X est égal à mon K (simple convention d'écriture, ça n'a pas d'importance).
Le "CallPutFlag" est donc c et p (avec son code on affiche l'un ou l'autre, avec le mien ça doit retourner les deux)

Merci !

[fred]

J'ai regardé que le 1er code en diagonale, mais il me semble que tu as oublié l'opérateur * entre N et ta parenthèse

c=float(S*N*(d1)-K*exp(-r*T)*N*(d2))

À tester.

[tyrtamos]

Bonjour,

+1 pour fred. Il manque aussi les 2 signes '*' dans le calcul de p, et ça marche! A noter que tous les float() ne sont pas utiles.

je suggère aussi, après l'importation de random, d'exécuter seed() pour qu'on ne retrouve pas toujours la même séquence aléatoire à chaque lancement:

from random import gauss, seed
seed()

[Sonny Sachs]

Ok d'accord j'avais oublié mes '*' ... Évidemment ! Merci à vous.

voici donc le code fonctionnel, j'y ai juste rajouté quelques print() qui vont me servir à vous exprimer le problème auquel je me trouve désormais confronté :

#!/usr/local/bin/python
# -*- coding: Utf-8 -*-

from math import log, exp, sqrt
from random import gauss, seed

seed()

S=float(input("S=")) #100
K=float(input("K=")) #105
r=float(input("r=")) #0.1
v=float(input("v=")) #0.2
T=float(input("T=")) #0.5

d1=float((log(S/K)+(r+v*v/2)*T)/(v*sqrt(T)))
d2=float(d1-v*sqrt(T))
print(d1,d2)

N=float(gauss(0,1))

print(N)

c=float(S*N*(d1)-K*exp(-r*T)*N*(d2))
p=float(-S*N*(-d1)+K*exp(-r*T)*N*(-d2))

print(c,p)

Output code Python :

d1 = 0.07926550931782529
d2 = -0.06215584691948424

N = 1.3601173885812377

c = 19.224742896876926
p = 19.224742896876926

Vérification du calcul "à la main" (via Google et l'IDLE) :

d1 = (log(100/105)+(0.1+0.2*0.2/2)*0.5)/(0.2*sqrt(0.5)) = 0.274433098
d2 = 0.274433098-0.2*sqrt(0.5) = 0.133011742

N = 1.3601173885812377 <-- je reprend le résultat du code Python

c = 100*1.3601173885812377*(0.274433098)-105*exp(-0.1*0.5)*1.3601173885812377*(0.133011742) = 19.2568386
p = -100*1.3601173885812377*(-0.274433098)-105*exp(-0.1*0.5)*1.3601173885812377*(-0.133011742) = 55.3954071

Donc, avec les inputs suivants :

S=100
K=105
r=0.1
v=0.2
T=0.5

J'obtiens des résultats différents pour les 3 procédés suivants :

1. Via ce code Python :

d1 = 0.07926550931782529
d2 = -0.06215584691948424

N = 1.3601173885812377

c = 19.224742896876926
p = 19.224742896876926

2. Via calcul "à la main" (Google et IDLE) :

d1 = 0.274433098
d2 = 0.133011742

N = 1.3601173885812377

c = 19.2568386
p = 55.3954071

3. Résultats corrects (vérifiés sur 2 tools différents "publics", donc a priori corrects, qui donnent les même résultats) :

d1 = ?
d2 = ?

N = ?

c = 5.6945
p = 5.5736

J'ignore les valeurs de d1 et d2 ainsi que N cependant les résultats sont cohérents.
Ceux issus du code Python semblent clairement incorrects : 19.225 et 19.225 (non seulement trop élevé mais surtout ils sont identiques !)
Ceux issus du calcul "à la main" sembles très clairement incorrects : 19.257 et 55.395 (c'est bien trop élevé !)

Par ailleurs, si je relance mon code Python, avec seed() la valeur de N changera. Donc les résultats de c et de p seront différents... Tandis que sur les logiciels en ligne "publics" les résultats ne changent jamais ! Comme s'ils ne prenaient pas en compte N...

EDIT : J'ai de nouveaux éléments, j'ai pratiquement trouvé la solution mais je cale sur un point !

Je vous présente la chose sous forme d'exemple.

ex:

S=42
K=40
r=0.1
v=0.2
T=0.5

d1=(log(S/K)+(r+v*v/2)*T)/(v*sqrt(T))
d1=(log(42/40)+(0.1+0.2*0.2/2)*0.5)/(0.2*sqrt(0.5))=0.7692626281060315

d2=d1-v*sqrt(T)
d2=0.7692626281060315-0.2*sqrt(0.5)=0.627841271868722

c=S*N*(d1)-K*exp(-r*T)*N*(d2)
p=K*exp(-r*T)*N*(-d2)-S*N*(-d1)

K*exp(-r*T)=40*exp(-0.1*0.5)=38.04917698002856

c=42*N*(0.7692626281060315)-38.04917698002856*N*(0.627841271868722)
p=38.04917698002856*N*(-0.627841271868722)-42*N*(-0.7692626281060315)

En utilisant l'approximation polynomiale, on obtient :

N(0.7692626281060315)=0.7791
N(-0.7692626281060315)=0.2209
N(0.627841271868722)=0.7349
N(-0.627841271868722)=0.2651

# comment arrive-t-on a ce résultat ???????!!!!!!!!!

Ce qui nous donne c et p :

c=42*0.7791-38.04917698002856*0.7349=4.759859837377011
p=38.04917698002856*0.2651-42*0.2209=0.8090368174055715

Soit :

c=4.76
p=0.81

Auriez-vous une idée de comment appliquer cette approximation polynomiale en Python ?

[tyrtamos]

Bonjour,

Intrigué par ton problème que je ne comprenais pas, je suis aller voir le lien que tu donnes plus haut: la formule de Black-Scholes.

J'ai vérifié tes calculs, et voilà ce que j'ai trouvé: N(x) est une fonction (ce qui explique l'absence initiale des '*' avant la parenthèse), et correspond bien à la fonction de répartition de la loi normale centrée réduite, et plus précisément à la probabilité qu'une variable distribuée selon la LNCR soit inférieure à x, et pas au random.gauss() qui donne une valeur aléatoire.

=> Il faut donc prendre ma formule (ou une table).

Voilà ma proposition de code modifié (Python 3). J'ai renommé ma formule N(x) pour coller à ton problème:

from math import log, exp, sqrt, pi

def N(x):
    """probabilité qu'une variable aléatoire distribuée selon une loi normale
       centrée réduite soit inférieure à x
    """
    if x==0:
        return 0.5
    if x>=7.56:
        return 1.0
    if x<=-7.56:
        return 0.0
    u=abs(x)
    n=int(u*2000)
    du=u/n
    k=1/sqrt(2*pi)
    u1=0
    f1=k
    p=0.5
    for i in range(0,n):
        u2=u1+du
        f2=k*exp(-0.5*u2*u2)
        p=p+(f1+f2)*du*0.5
        u1=u2
        f1=f2
    if x>0:
        return p
    else:
        return 1.0-p

S=float(input("S="))
K=float(input("K="))
r=float(input("r="))
v=float(input("v="))
T=float(input("T="))

d1 = (log(S/K)+(r+v*v/2)*T)/(v*sqrt(T))
d2 = d1-v*sqrt(T)
print("d1 =", d1)
print("d2 =", d2)

c = S*N(d1)-K*exp(-r*T)*N(d2)
p = -S*N(-d1)+K*exp(-r*T)*N(-d2)
print("c =", c)
print("p =", p)

Avec les valeurs:

S = 100
K = 105
r = 0.1
v = 0.2
T = 0.5

On trouve:

d1 = 0.07926550931782529
d2 = -0.06215584691948424
c = 5.694450705766563
p = 5.573540278341532

Et, bien sûr, on obtient le même résultat à chaque exécution puisqu'il n'y a plus de valeur aléatoire.

Ok?

[Sonny Sachs]

Ok ! Merci beaucoup ! Dans mes tests j'avais pourtant tenté d'utiliser les tables mais les résultats étaient pas vraiment optimums.

Bon, j'ai une dernière petite question, je dois calculer la dérivée seconde de N(x).

Soit N'(x)

J'ignore s'il est vraiment bon d'utiliser une apostrophe dans le code ex: print("N\'(x)") on peut utiliser Np(d2) sinon.

J'ai par exemple ceci à calculer :

θc = -(S*Np(d1)*v)/(2*sqrt(T))-r*K*exp(-r*T)*N(d2)
print("θc=", θc)

θp = -(S*Np(d1)*v)/(2*sqrt(T))+r*K*exp(-r*T)*N(-d2)
print("θp=", θp)

Où N'(x)=(1/(sqrt(2*pi))*exp(-x*x/2)

J'ai tenté de mettre une deuxième fonction, tel que :

def N'(x):

Et de seulement changer la valeur de k.

k par N(x) = k=1/sqrt(2*pi)
k par N'(x) = k=(1/(sqrt(2*pi))*exp(-x*x/2)

Mais ça n'a pas marché... J'obtiens des valeurs anormales de θ.

[tyrtamos]

Bonjour,

Tu parles de dérivée seconde de N(x), mais je crois que tu veux la dérivée première. En effet, si N(x) est la fonction de répartition de la loi normale centrée réduite, elle est une intégrale de la loi de densité de probabilité. La dérivée de N(x) est donc la loi de densité elle-même, et on la retrouve bien avec ta formule:

def Np(x):
    return 1/(sqrt(2*pi))*exp(-x*x/2)

Et non: on ne peut mettre ni apostrophe ni lettre grecque dans les noms de variable ou de fonction en python. Tu devrais apprendre un peu plus le langage Python: tu gagnerais du temps.

Bref, si je fais ton calcul comme je le comprends, voilà ce que ça donne:

from math import log, exp, sqrt, pi

def N(x):
    """fonction de répartition de la loi normale centrée réduite
       => probabilité qu'une variable aléatoire distribuée selon
       cette loi soit inférieure à x
    """
    if x==0:
        return 0.5
    if x>=7.56:
        return 1.0
    if x<=-7.56:
        return 0.0
    u=abs(x)
    n=int(u*2000)
    du=u/n
    k=1/sqrt(2*pi)
    u1=0
    f1=k
    p=0.5
    for i in range(0,n):
        u2=u1+du
        f2=k*exp(-0.5*u2*u2)
        p=p+(f1+f2)*du*0.5
        u1=u2
        f1=f2
    if x>0:
        return p
    else:
        return 1.0-p

def Np(x):
    """densité de probabilité de la loi normale centrée réduite"""
    return 1/(sqrt(2*pi))*exp(-x*x/2)

S=float(input("S="))
K=float(input("K="))
r=float(input("r="))
v=float(input("v="))
T=float(input("T="))

d1 = (log(S/K)+(r+v*v/2)*T)/(v*sqrt(T))
d2 = d1-v*sqrt(T)
print("d1 =", d1)
print("d2 =", d2)

Tc = -(S*Np(d1)*v)/(2*sqrt(T))-r*K*exp(-r*T)*N(d2)
print("Tc=", Tc)

Tp = -(S*Np(d1)*v)/(2*sqrt(T))+r*K*exp(-r*T)*N(-d2)
print("Tp=", Tp)

Avec ces valeurs:

S = 100
K = 105
r = 0.1
v = 0.2
T = 0.5

Cela donne:

d1 = 0.07926550931782529
d2 = -0.06215584691948424
Tc= -10.370647296772137
Tp= -0.3827383395146393

Mais je ne sais pas si ces valeurs sont correctes pour toi.

Et ce serait encore mieux si tu disais à quoi ça correspond: il ne faut jamais perdre une occasion d'apprendre quelque chose...

[Sonny Sachs]

Yes, ça fonctionne. C'était tout bête pour le coup j'aurai pu trouver tout seul...

Sinon, oui c'est la dérivée première bien entendu, petite erreur.

Par contre es-tu sûr pour les lettres grecques ? Moi ça fonctionne :

>>> Δ=1
>>> print("Δ=",Δ)
1

Et pour l'apostrophe je peux toujours l'utiliser dans l'output tel que :

>>> Np=1
>>> print("N\'=",Np)
N'=
>>>

Bref, ce ne sont que des détails.

Merci à toi Tyrtamos !

PS : concernant ce à quoi cela correspond, n'as-tu pas vu à quoi servait Black & Scholes ?

Je pourrai te refiler le petit pricer une fois fini si cela t'intéresses, je doute que ça intéresse grand monde ici ! C'était juste un petit souci de Python et je ne devrais plus rencontrer de problème avec désormais.

Merci encore !

[tyrtamos]

Par contre es-tu sûr pour les lettres grecques ? Moi ça fonctionne

Effectivement, ça fonctionne sous Python 3 (et moi je travaille surtout sous Python 2 qui ne le supporte pas).

PS : concernant ce à quoi cela correspond, n'as-tu pas vu à quoi servait Black & Scholes ?

Pour le 1er code, bien sûr! J'ai lu avec beaucoup d'intérêt la page de wikipedia sur le sujet. Je savais que les mathématiques étaient bien utilisées dans le monde financier, mais je suis toujours surpris de voir le haut niveau pratiqué. Avec un peu d'inquiétude aussi: le fait qu'une partie de plus en plus grande des cours de la bourse résulte des combats ente ordinateurs qui luttent au 1/1000 de seconde près.

Par contre, je n'ai pas retrouvé les calculs de θc et θp, et c'est pour eux que je posais la question.

Je pourrai te refiler le petit pricer une fois fini si cela t'intéresses

Avec plaisir! Merci d'avance