12. SciPy#
En plus de ce qui est inclus dans Anaconda, ce cours nécessitera les bibliothèques suivantes :
!pip install --upgrade quantecon
Nous utilisons les importations suivantes.
import numpy as np
import quantecon as qe
12.1. Vue d’ensemble#
SciPy s’appuie sur NumPy pour fournir des outils courants de programmation scientifique tels que
etc., etc.
Comme NumPy, SciPy est stable, mature et largement utilisé.
De nombreuses routines de SciPy sont de fines surcouches autour de bibliothèques Fortran standard de l’industrie telles que LAPACK, BLAS, etc.
Il n’est pas vraiment nécessaire d”« apprendre » SciPy dans son ensemble.
Une approche plus courante consiste à avoir une idée de ce que contient la bibliothèque, puis à consulter la documentation au besoin.
Dans ce cours, nous cherchons uniquement à mettre en évidence quelques parties utiles du package.
12.2. SciPy versus NumPy#
SciPy est un package qui contient divers outils construits sur NumPy, en utilisant son type de données tableau et les fonctionnalités connexes.
Note
Dans les anciennes versions de SciPy (scipy < 0.15.1), importer le package importait aussi les symboles NumPy dans l’espace de noms global, comme on peut le voir dans cet extrait du fichier d’initialisation de SciPy :
from numpy import *
from numpy.random import rand, randn
from numpy.fft import fft, ifft
from numpy.lib.scimath import *
Cependant, il est préférable d’utiliser explicitement les fonctionnalités de NumPy.
import numpy as np
a = np.identity(3)
Les versions plus récentes de SciPy (1.15+) n’importent plus automatiquement les symboles NumPy.
Ce qui est utile dans SciPy, c’est la fonctionnalité de ses sous-packages
scipy.optimize,scipy.integrate,scipy.stats, etc.
Explorons quelques-uns des principaux sous-packages.
12.3. Statistiques#
Le sous-package scipy.stats fournit
de nombreux objets de variables aléatoires (densités, fonctions de répartition, échantillonnage aléatoire, etc.)
certaines procédures d’estimation
certains tests statistiques
12.3.1. Variables aléatoires et lois de probabilité#
Rappelons que numpy.random fournit des outils pour générer des variables aléatoires
rng = np.random.default_rng()
rng.beta(5, 5, size=3)
array([0.52619744, 0.56542371, 0.50532167])
Ceci génère un tirage de la loi ayant la fonction de densité ci-dessous lorsque a, b = 5, 5
Parfois, nous avons besoin d’accéder à la densité elle-même, ou à la fonction de répartition, aux quantiles, etc.
Pour cela, nous pouvons utiliser scipy.stats, qui fournit toutes ces fonctionnalités ainsi que la génération de nombres aléatoires dans une seule interface cohérente.
Voici un exemple d’utilisation
from scipy.stats import beta
import matplotlib.pyplot as plt
q = beta(5, 5) # Beta(a, b), avec a = b = 5
obs = q.rvs(2000) # 2000 observations
grid = np.linspace(0.01, 0.99, 100)
fig, ax = plt.subplots()
ax.hist(obs, bins=40, density=True)
ax.plot(grid, q.pdf(grid), 'k-', linewidth=2)
plt.show()
L’objet q qui représente la loi possède des méthodes utiles supplémentaires, notamment
q.cdf(0.4) # Fonction de répartition
np.float64(0.26656768000000003)
q.ppf(0.8) # Fonction quantile (fonction de répartition inverse)
np.float64(0.6339134834642708)
q.mean()
np.float64(0.5)
La syntaxe générale pour créer ces objets qui représentent des lois de probabilité (de type rv_frozen) est
name = scipy.stats.distribution_name(shape_parameters, loc=c, scale=d)
Ici, distribution_name est l’un des noms de lois de probabilité de scipy.stats.
Les paramètres loc et scale transforment la variable aléatoire originale
\(X\) en \(Y = c + d X\).
12.3.2. Syntaxe alternative#
Il existe une manière alternative d’appeler les méthodes décrites ci-dessus.
Par exemple, le code qui génère la figure ci-dessus peut être remplacé par
obs = beta.rvs(5, 5, size=2000)
grid = np.linspace(0.01, 0.99, 100)
fig, ax = plt.subplots()
ax.hist(obs, bins=40, density=True)
ax.plot(grid, beta.pdf(grid, 5, 5), 'k-', linewidth=2)
plt.show()
12.3.3. Autres trouvailles dans scipy.stats#
Il existe une variété de fonctions statistiques dans scipy.stats.
Par exemple, scipy.stats.linregress implémente la régression linéaire simple
from scipy.stats import linregress
x = rng.standard_normal(200)
y = 2 * x + 0.1 * rng.standard_normal(200)
gradient, intercept, r_value, p_value, std_err = linregress(x, y)
gradient, intercept
(np.float64(2.0074526915071473), np.float64(0.004899837480739758))
Pour voir la liste complète, consultez la documentation.
12.4. Racines et points fixes#
Une racine ou un zéro d’une fonction réelle \(f\) sur \([a,b]\) est un \(x \in [a, b]\) tel que \(f(x)=0\).
Par exemple, si nous traçons la fonction
avec \(x \in [0,1]\), nous obtenons
f = lambda x: np.sin(4 * (x - 1/4)) + x + x**20 - 1
x = np.linspace(0, 1, 100)
fig, ax = plt.subplots()
ax.plot(x, f(x), label='$f(x)$')
ax.axhline(ls='--', c='k')
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$f(x)$', fontsize=12)
ax.legend(fontsize=12)
plt.show()
L’unique racine est d’environ 0,408.
Considérons quelques techniques numériques pour trouver les racines.
12.4.1. Bisection#
L’un des algorithmes les plus courants pour la recherche numérique de racines est la bissection.
Pour comprendre l’idée, rappelez-vous le jeu bien connu où
Le joueur A pense à un nombre secret entre 1 et 100
Le joueur B demande s’il est inférieur à 50
Si oui, B demande s’il est inférieur à 25
Si non, B demande s’il est inférieur à 75
Et ainsi de suite.
C’est la bissection.
Voici une implémentation simpliste de l’algorithme en Python.
Elle fonctionne pour toutes les fonctions continues croissantes suffisamment bien comportées avec \(f(a) < 0 < f(b)\)
def bisect(f, a, b, tol=10e-5):
"""
Implements the bisection root finding algorithm, assuming that f is a
real-valued function on [a, b] satisfying f(a) < 0 < f(b).
"""
lower, upper = a, b
while upper - lower > tol:
middle = 0.5 * (upper + lower)
if f(middle) > 0: # la racine est entre lower et middle
lower, upper = lower, middle
else: # la racine est entre middle et upper
lower, upper = middle, upper
return 0.5 * (upper + lower)
Testons-la en utilisant la fonction \(f\) définie dans (12.2)
bisect(f, 0, 1)
0.408294677734375
Sans surprise, SciPy fournit sa propre fonction de bissection.
Testons-la en utilisant la même fonction \(f\) définie dans (12.2)
from scipy.optimize import bisect
bisect(f, 0, 1)
0.4082935042806639
12.4.2. La méthode de Newton-Raphson#
Un autre algorithme de recherche de racines très courant est la méthode de Newton-Raphson.
Dans SciPy, cet algorithme est implémenté par scipy.optimize.newton.
Contrairement à la bissection, la méthode de Newton-Raphson utilise l’information sur la pente locale afin d’accélérer la convergence.
Étudions cela en utilisant la même fonction \(f\) définie ci-dessus.
Avec une condition initiale appropriée pour la recherche, nous obtenons la convergence :
from scipy.optimize import newton
newton(f, 0.2) # Démarre la recherche à la condition initiale x = 0.2
np.float64(0.40829350427935673)
Mais d’autres conditions initiales conduisent à un échec de la convergence :
newton(f, 0.7) # Démarre la recherche à x = 0.7 à la place
np.float64(0.7001700000000279)
12.4.3. Méthodes hybrides#
Un principe général des méthodes numériques est le suivant :
Si vous avez une connaissance spécifique d’un problème donné, vous pourriez être en mesure de l’exploiter pour gagner en efficacité.
Sinon, le choix de l’algorithme implique un compromis entre vitesse et robustesse.
En pratique, la plupart des algorithmes par défaut pour la recherche de racines, l’optimisation et les points fixes utilisent des méthodes hybrides.
Ces méthodes combinent généralement une méthode rapide avec une méthode robuste de la manière suivante :
Tenter d’utiliser une méthode rapide
Vérifier les diagnostics
Si les diagnostics sont mauvais, passer à un algorithme plus robuste
Dans scipy.optimize, la fonction brentq est une telle méthode hybride et constitue un bon choix par défaut
from scipy.optimize import brentq
brentq(f, 0, 1)
0.40829350427936706
Ici, la solution correcte est trouvée et la vitesse est meilleure que la bissection :
with qe.Timer(unit="milliseconds"):
brentq(f, 0, 1)
0.0470 ms elapsed
with qe.Timer(unit="milliseconds"):
bisect(f, 0, 1)
0.0973 ms elapsed
12.4.4. Recherche de racines multivariée#
Utilisez scipy.optimize.fsolve, une surcouche pour une méthode hybride dans MINPACK.
Voir la documentation pour plus de détails.
12.4.5. Points fixes#
Un point fixe d’une fonction réelle \(f\) sur \([a,b]\) est un \(x \in [a, b]\) tel que \(f(x)=x\).
SciPy dispose également d’une fonction pour trouver des points fixes (scalaires)
from scipy.optimize import fixed_point
fixed_point(lambda x: x**2, 10.0) # 10.0 est une supposition initiale
array(1.)
Si vous n’obtenez pas de bons résultats, vous pouvez toujours revenir au chercheur de racines brentq, car
le point fixe d’une fonction \(f\) est la racine de \(g(x) := x - f(x)\).
12.5. Optimisation#
La plupart des packages numériques ne fournissent que des fonctions de minimisation.
La maximisation peut être réalisée en rappelant que le maximiseur d’une fonction \(f\) sur le domaine \(D\) est le minimiseur de \(-f\) sur \(D\).
La minimisation est étroitement liée à la recherche de racines : pour les fonctions lisses, les optima intérieurs correspondent aux racines de la dérivée première.
Le compromis vitesse/robustesse décrit ci-dessus est également présent dans l’optimisation numérique.
À moins que vous n’ayez des informations préalables que vous pouvez exploiter, il est généralement préférable d’utiliser des méthodes hybrides.
Pour la minimisation univariée (c’est-à-dire scalaire) sous contrainte, une bonne option hybride est fminbound
from scipy.optimize import fminbound
fminbound(lambda x: x**2, -1, 2) # Recherche dans [-1, 2]
np.float64(0.0)
12.5.1. Optimisation multivariée#
Les optimiseurs locaux multivariés incluent minimize, fmin, fmin_powell, fmin_cg, fmin_bfgs et fmin_ncg.
Les optimiseurs locaux multivariés sous contrainte incluent fmin_l_bfgs_b, fmin_tnc, fmin_cobyla.
Voir la documentation pour plus de détails.
12.6. Intégration#
La plupart des méthodes d’intégration numérique fonctionnent en calculant l’intégrale d’un polynôme approximant.
L’erreur résultante dépend de la qualité de l’ajustement du polynôme à l’intégrande, qui dépend à son tour du degré de « régularité » de l’intégrande.
Dans SciPy, le module pertinent pour l’intégration numérique est scipy.integrate.
Un bon choix par défaut pour l’intégration univariée est quad
from scipy.integrate import quad
integral, error = quad(lambda x: x**2, 0, 1)
integral
0.33333333333333337
En fait, quad est une interface vers une routine d’intégration numérique très standard de la bibliothèque Fortran QUADPACK.
Elle utilise la quadrature de Clenshaw-Curtis, basée sur un développement en termes de polynômes de Tchebychev.
Il existe d’autres options pour l’intégration univariée — une option utile est fixed_quad, qui est rapide et fonctionne donc bien à l’intérieur des boucles for.
Il existe également des fonctions pour l’intégration multivariée.
Voir la documentation pour plus de détails.
12.7. Algèbre linéaire#
Nous avons vu que NumPy fournit un module d’algèbre linéaire appelé linalg.
SciPy fournit également un module d’algèbre linéaire du même nom.
Ce dernier n’est pas un sur-ensemble exact du premier, mais dans l’ensemble il possède davantage de fonctionnalités.
Nous vous laissons explorer l”ensemble des routines disponibles.
12.8. Exercices#
Les premiers exercices concernent l’évaluation d’une option d’achat européenne sous l’hypothèse de neutralité au risque. Le prix satisfait
où
\(\beta\) est un facteur d’actualisation,
\(n\) est la date d’échéance,
\(K\) est le prix d’exercice et
\(\{S_t\}\) est le prix de l’actif sous-jacent à chaque instant \(t\).
Par exemple, si l’option d’achat consiste à acheter des actions Amazon au prix d’exercice \(K\), le propriétaire a le droit (mais non l’obligation) d’acheter 1 action Amazon au prix \(K\) après \(n\) jours.
Le gain est donc \(\max\{S_n - K, 0\}\)
Le prix est l’espérance du gain, actualisée à la valeur actuelle.
Exercice 12.1
Supposons que \(S_n\) suit la loi log-normale de paramètres \(\mu\) et \(\sigma\). Soit \(f\) la densité de cette loi. Alors
Tracez la fonction
sur l’intervalle \([0, 400]\) lorsque μ, σ, β, n, K = 4, 0.25, 0.99, 10, 40.
Indication
Depuis scipy.stats, vous pouvez importer lognorm puis utiliser lognorm.pdf(x, σ, scale=np.exp(μ)) pour obtenir la densité \(f\).
Solution
Voici une solution possible
from scipy.integrate import quad
from scipy.stats import lognorm
μ, σ, β, n, K = 4, 0.25, 0.99, 10, 40
def g(x):
return β**n * np.maximum(x - K, 0) * lognorm.pdf(x, σ, scale=np.exp(μ))
x_grid = np.linspace(0, 400, 1000)
y_grid = g(x_grid)
fig, ax = plt.subplots()
ax.plot(x_grid, y_grid, label="$g$")
ax.legend()
plt.show()
Exercice 12.2
Afin d’obtenir le prix de l’option, calculez l’intégrale de cette fonction numériquement en utilisant quad de scipy.integrate.
Solution
P, error = quad(g, 0, 1_000)
print(f"The numerical integration based option price is {P:.3f}")
The numerical integration based option price is 15.188
Exercice 12.3
Essayez d’obtenir un résultat similaire en utilisant Monte-Carlo pour calculer le terme d’espérance dans le prix de l’option, plutôt que quad.
En particulier, utilisez le fait que si \(S_n^1, \ldots, S_n^M\) sont des tirages indépendants de la loi log-normale spécifiée ci-dessus, alors, par la loi des grands nombres,
Posez M = 10_000_000
Solution
Voici une solution :
rng = np.random.default_rng()
M = 10_000_000
S = np.exp(μ + σ * rng.standard_normal(M))
return_draws = np.maximum(S - K, 0)
P = β**n * np.mean(return_draws)
print(f"The Monte Carlo option price is {P:3f}")
The Monte Carlo option price is 15.186671
Exercice 12.4
Dans ce cours, nous avons abordé le concept d”appels récursifs de fonctions.
Essayez d’écrire une implémentation récursive de la fonction de bissection maison décrite ci-dessus.
Testez-la sur la fonction (12.2).
Solution
Voici une solution raisonnable :
def bisect(f, a, b, tol=10e-5):
"""
Implements the bisection root-finding algorithm, assuming that f is a
real-valued function on [a, b] satisfying f(a) < 0 < f(b).
"""
lower, upper = a, b
if upper - lower < tol:
return 0.5 * (upper + lower)
else:
middle = 0.5 * (upper + lower)
print(f'Current mid point = {middle}')
if f(middle) > 0: # Implique que la racine est entre lower et middle
return bisect(f, lower, middle)
else: # Implique que la racine est entre middle et upper
return bisect(f, middle, upper)
Nous pouvons la tester comme suit
f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1
bisect(f, 0, 1)
Current mid point = 0.5
Current mid point = 0.25
Current mid point = 0.375
Current mid point = 0.4375
Current mid point = 0.40625
Current mid point = 0.421875
Current mid point = 0.4140625
Current mid point = 0.41015625
Current mid point = 0.408203125
Current mid point = 0.4091796875
Current mid point = 0.40869140625
Current mid point = 0.408447265625
Current mid point = 0.4083251953125
Current mid point = 0.40826416015625
0.408294677734375