Apprendre le Modèle de Heston de a à z sur Python — Partie I

Nous allons essayer de développer l’intuition derrière le modèle de Heston afin de pouvoir l'implémenter. Les démonstrations mathématiques sont moins importantes pour nous dans ce document.

Nous fournissons également la formule de pricing de base et les différents problèmes qui pourraient survenir lors du pricing avec ce premier.

Il vous suffit d‘avoir les notions de probabilité et les intégrales. Le reste n'est que des ajustements à ce que vous savez déjà.

Ce document est également disponible sur GitHub. Donc, si vous souhaitez obtenir les fichiers .ipynb ou .py, vous pouvez les trouver ici!

Table des matières

1. Modèle de Black & Scholes
2. De Black & Scholes à Heston
3. Modèle de Heston: composant par composant
4. Version simple de Heston
5. Code: étape par étape
6. Code complet

Modèle de Black & Scholes

Nous allons commencer par présenter le modèle Black & Scholes et puis nous allons faire les parallèles avec le Heston.

Le modèle de pricing de Black and Scholes donne la formule suivante pour un call:

𝑁 est la fonction de distribution cumulative pour une variable aléatoire normale. 𝑁 () est simplement la fonction NORM.DIST dans Excel.

Cette fonction prend les entrées suivantes: NORM.DIST (x, moyenne, écart-type) et calcule, pour toute entrée x, l'aire sous la courbe en cloche en partant par la gauche.

Supposons que nous avons 0 comme moyenne et 1 comme écart type pour notre distribution.

Cela signifie qu'à 0, la zone du graphique en-dessous (à partir du point le plus éloigné à gauche) sera exactement la moitié.

Cette zone représente la probabilité qu'une valeur aléatoire choisie parmi une distribution de moyenne 0 et d'écart type 1 soit inférieure ou égale à 𝑥.

Par conséquent, plus le nombre 𝑥 est grand, plus la probabilité qu'un nombre aléatoire sélectionné soit inférieur à 𝑥 est élevée.

De Black & Scholes à Heston

Il faut d'abord avoir une vue d’ensemble avant de s'aventurer dans les détails.

La formule Heston semble intimidante mais en quelques étapes, nous allons vous montrer que ce n’est pas vraiment si loin du modèle Black & Scholes.

Notre exemple n’est qu’une simple option d’achat. Disons que vous avez une option européenne expirant dans 2 ans avec un prix d’exercice à 100 et le prix de l’action (disons Apple) est à 110 aujourd’hui.

Nous supposons également que le taux d’intérêt sans risque est de 1%. On s’attend à ce que la valeur de l’option soit d’environ 10 si le taux d’intérêt est suffisaments bas et qu’il y ait pas de volatilité dans le marché sous le Black & Scholes.

Nous pouvons utiliser une notation générale et définir:

• 𝐶 la valeur du call
• 𝑆𝑡 la valeur de l’action
• 𝐾 le prix d’exercice
• 𝑇 la maturité (ainsi 𝑇−𝑡 = τ est l’échéance.)
• 𝑟 le taux sans risque

La formule de pricing du modèle Heston est:

Dans le modèle de Black Scholes, les fonctions P1 et P2 sont en fait les fonctions N(d1) et N(d2).

Dans le Heston, on remplace N(d1) par:

et N(d2) avec:

Ne soyez pas intimidé par la formule. N'oubliez pas que l'utilisation est la même que N(d1) et N(d2).

Notre formule complète est donc:

Modèle de Heston: composant par composant

Nous allons décomposer cela morceau par morceau afin de mieux comprendre le modèle. Après cela, nous allons commencer par écrire les éléments les plus simples et ensuite écrire l'intégrale complète.

Nous avons deux intégrales. Il y a un Re encapsulant chacun des intégrales. Cela signifie simplement que nous avons des valeurs complexes et que lorsque nous calculons le prix, nous ne considérerons que la partie réel.

Alors ne vous inquiétez pas si vous n’avez pas traité de nombres complexes car vous n’avez pas à le faire! On note que iest defini comme i = √-1.

Nous revenons donc au P1 ou P2. Nous pouvons les raccourcir en Pj:

Chaque intégrale a une fonction: fj. Ces fonctions sont définies comme:

et les fonctions C et D sont définies comme:

Nous définissons également les éléments finaux x, d, g et u qui sont utilisés dans fj:

Nous avons également un élément BRS qui est défini comme:

Version simple de Heston

Au lieu de faire deux intégrales, vous pouvez tout mettre sous une intégrale unique. (Juste pour être sûr, j'ai montré comment faire cela dans l'annexe de ce document jupyter notebook)

et nous définissons donc une fonction f:

Code: étape par étape

1. Définition de f (fHeston)

Nous commençons donc par décomposer cette fonction monstreuse. Nous allons commencer par définir 𝑓 comme fHeston dans Python avec les inputs s, St, K, r, T, sigma, kappa, theta, volvol, rho, les paramètres du Heston étant{𝜎, 𝜅, 𝜃, 𝜈, ρ} .

Parce que notre fonction utilisera des nombres complexes, nous pouvons simplement définir i comme une variable globale. Chaque nombre complexe ne sera essentiellement qu'une opération sur i.

Vous noterez également que dans la mise en œuvre, vous pouvez simplement traiter i comme n'importe quelle autre variable et tout devrait bien se passer. Alors ne vous inquiétez pas si vous ne savez pas comment travailler avec des nombres complexes!

L'expression ρ.𝜎.𝑖.𝑠 est beaucoup utilisée. Nous pouvons définir cela à l'avance et l'appeler simplement prod.

import numpy as npi = complex(0,1)
# Sera utilisée dans le pricer Heston
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
    # Beaucoup utilise
    prod = rho * sigma *i *s

Nous étendons notre fonction. Nous divisons et calculons d et g tels que nous les avons définis dans la section ci-dessus:

# Sera utilisée dans le pricer Heston
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
    # Beaucoup utilise
    prod = rho * sigma *i *s 
    
    # Calcule d
    d1 = (prod - kappa)**2
    d2 = (sigma**2) * (i*s + s**2)
    d = np.sqrt(d1 + d2)
    
    # Calcule g
    g1 = kappa - prod - d
    g2 = kappa - prod + d
    g = g1/g2

Ensuite, décomposez la grande exponentielle. Nous traitons d’abord la première partie:

# Sera utilisée dans le pricer Heston
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
    # Beaucoup utilise
    prod = rho * sigma *i *s 
    
    # Calcule d
    d1 = (prod - kappa)**2
    d2 = (sigma**2) * (i*s + s**2)
    d = np.sqrt(d1 + d2)
    
    # Calcule g
    g1 = kappa - prod - d
    g2 = kappa - prod + d
    g = g1/g2
    
    # Calcule le premier exponentiel
    exp1 = np.exp(np.log(St) * i *s) * np.exp(i * s* r* T)
    exp2 = 1 - g * np.exp(-d *T)
    exp3 = 1- g
    mainExp1 = exp1*np.power(exp2/exp3, -2*theta*kappa/(sigma **2))

On traite ensuite la deuxième intégrale pour obtenir la fonction fHeston complète.

import numpy as npi = complex(0,1)
u = 1
# Sera utilisée dans le pricer Heston
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
    # To be used a lot
    prod = rho * sigma *i *s 
    
    # Calcule d
    d1 = (prod - kappa)**2
    d2 = (sigma**2) * (i*s + s**2)
    d = np.sqrt(d1 + d2)
    
    # Calcule g
    g1 = kappa - prod - d
    g2 = kappa - prod + d
    g = g1/g2
    
    # Calcule le premier exponentiel
    exp1 = np.exp(np.log(St) * i *s) * np.exp(i * s* r* T)
    exp2 = 1 - g * np.exp(-d *T)
    exp3 = 1- g
    mainExp1 = exp1*np.power(exp2/exp3, -2*theta*kappa/(sigma**2))
    
    # Calcule le deuxieme exponentiel
    exp4 = theta * kappa * T/(sigma **2)
    exp5 = volvol/(sigma **2)
    exp6 = (1 - np.exp(-d * T))/(1 - g * np.exp(-d * T))
    mainExp2 = np.exp((exp4 * g1) + (exp5 *g1 * exp6))
    
    return (mainExp1 * mainExp2)

2. Calculer l’intègral

Fixer la limite de l’intègral

Vu que notre intégrale monte à l'infini, nous ne pouvons pas l’evaluer en pratique. Ainsi, nous allons faire sélectionner une limite suffisamment grande pour permettre à notre intégrale de converger.

Par convergence, nous entendons simplement qu'à un moment donné, les zones supplémentaires sous la courbe que nous calculons sont si petites qu'elles sont insignifiantes.

Dans notre cas, nous avons sélectionné 100. Vous pouvez essayer ceci avec différentes valeurs pour tester.

Sélection d'un schéma d'intégration

Pour calculer notre intégrale, nous utiliserons un schéma d'intégration très simple. C'est le schéma que nous avons tous utilisé au lycée.

Nous divisons simplement notre zone en rectangles et trouvons le point médian. Nous utilisons ce point médian pour calculer l'aire de chaque rectangle et les additionner.

Dans notre cas, nous initialisons P, le prix final, à 0, nous divisons notre zone en 1000 rectangles et fixons la limite à 100.

Donc, en substance, la largeur de chaque rectangle sera du = 100/1000.

Nous calculons également la première partie avant l'intégrale:

# Pricer Heston
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
    P, iterations, maxNumber = 0,1000,100
    ds = maxNumber/iterations
    
    element1 = 0.5 * (St - K * np.exp(-r * T))

Nous définissons s1 et s2 tels que définis dans la formule. Vous observez que nous incrémentons s1 à chaque exécution de la boucle.

Par exemple, à j = 1, nous faisons s1 = du*(2 * 1+1)/2 qui i donne s1 = du * 1.5.Lorsque j = 2, alors s1 = du*(2 * 2+1)/2 qui est s1 = 2.5 * du. Cela vérifie notre règle du point médian pour l'intégration.

Vous noterez également que nous ne partons pas de 0. En effet, se rapprocher de 0 peut signifier diviser par 0. Bien sûr, cela posera un problème.

# Pricer Heston
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
    P, iterations, maxNumber = 0,1000,100
    ds = maxNumber/iterations
    
    element1 = 0.5 * (St - K * np.exp(-r * T))
    
    # Calcule l'integral complexe
    # Utilisant j au lieu de i dans la boucle
    for j in range(1, iterations):
        s1 = ds * (2*j + 1)/2
        s2 = s1 - i
        
        numerator1 = fHeston(s2,  St, K, r, T, 
                             sigma, kappa, theta, volvol, rho)
        numerator2 = K * fHeston(s1,  St, K, r, T, 
                             sigma, kappa, theta, volvol, rho)
        denominator = np.exp(np.log(K) * i * s1) *i *s1

Donc, à ce stade, nous avons la largeur de tous nos rectangles, et la hauteur est simplement la valeur du fHeston lorsque nous nous déplaçons le long de l’axe horizontale.

La somme des résultats de fHeston divisée par le dénominateur correspond à la longeur de nos réctangles. Nous pouvons maintenant calculer les superficies de nos réctangles et les additionner.

# Pricer Heston
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
    P, iterations, maxNumber = 0,1000,100
    ds = maxNumber/iterations
    
    element1 = 0.5 * (St - K * np.exp(-r * T))
    
    # Calcule l'integral complexe
    # Utilisant j au lieu de i dans la boucle
    for j in range(1, iterations):
        s1 = ds * (2*j + 1)/2
        s2 = s1 - i
        
        numerator1 = fHeston(s2,  St, K, r, T, 
                             sigma, kappa, theta, volvol, rho)
        numerator2 = K * fHeston(s1,  St, K, r, T, 
                             sigma, kappa, theta, volvol, rho)
        denominator = np.exp(np.log(K) * i * s1) *i *s1
      
        P += ds *(numerator1 - numerator2)/denominator
    
    element2 = P/np.pi
    
    return np.real((element1 + element2))

Code complet

Nous avons donc réussi à implémenter le modèle Heston en moins de 50 lignes de code!

import numpy as npi = complex(0,1)
u = 1
# Sera utilisée dans le pricer Heston
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
    # To be used a lot
    prod = rho * sigma *i *s 
    
    # Calcule d
    d1 = (prod - kappa)**2
    d2 = (sigma**2) * (i*s + s**2)
    d = np.sqrt(d1 + d2)
    
    # Calcule g
    g1 = kappa - prod - d
    g2 = kappa - prod + d
    g = g1/g2
    
    # Calcule le premier exponentiel
    exp1 = np.exp(np.log(St) * i *s) * np.exp(i * s* r* T)
    exp2 = 1 - g * np.exp(-d *T)
    exp3 = 1- g
    mainExp1 = exp1*np.power(exp2/exp3, -2*theta*kappa/(sigma**2))
    
    # Calcule le deuxieme exponentiel
    exp4 = theta * kappa * T/(sigma **2)
    exp5 = volvol/(sigma **2)
    exp6 = (1 - np.exp(-d * T))/(1 - g * np.exp(-d * T))
    mainExp2 = np.exp((exp4 * g1) + (exp5 *g1 * exp6))
    
    return (mainExp1 * mainExp2)
# Pricer Heston
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
    P, iterations, maxNumber = 0,1000,100
    ds = maxNumber/iterations
    
    element1 = 0.5 * (St - K * np.exp(-r * T))
    
    # Calcule l'integral complexe
    # Utilisant j au lieu de i dans la boucle
    for j in range(1, iterations):
        s1 = ds * (2*j + 1)/2
        s2 = s1 - i
        
        numerator1 = fHeston(s2,  St, K, r, T, 
                             sigma, kappa, theta, volvol, rho)
        numerator2 = K * fHeston(s1,  St, K, r, T, 
                             sigma, kappa, theta, volvol, rho)
        denominator = np.exp(np.log(K) * i * s1) *i *s1
      
        P += ds *(numerator1 - numerator2)/denominator
    
    element2 = P/np.pi
    
    return np.real((element1 + element2))