Dérivées en code

À quelle vitesse cela change-t-il ?

La dérivée d'une fonction f en x est son taux de variation instantané : la pente de la courbe en ce point. En finance, presque toute sensibilité est une dérivée. Le delta d'une option est de combien son prix bouge par unité de mouvement du sous-jacent, son gamma est la vitesse de variation du delta lui-même, et la duration d'une obligation est une dérivée du prix par rapport au rendement. Le calcul différentiel donne des formules exactes pour les fonctions simples (la dérivée de x**2 est 2*x), mais les vraies fonctions de valorisation sont souvent trop compliquées ou n'existent que sous forme de code. On estime alors la dérivée numériquement.

# A function is just a Python function
def f(x):
    return x * x

# Slope between two nearby points approaches the true slope (2x = 6 at x = 3)
x = 3.0
for h in (1.0, 0.1, 0.001):
    slope = (f(x + h) - f(x)) / h
    print(h, round(slope, 4))
# 1.0 7.0
# 0.1 6.1
# 0.001 6.001

Différences finies : avant et centrée

La différence avant regarde un petit pas h devant : (f(x + h) - f(x)) / h. Son erreur diminue proportionnellement à h. La différence centrée regarde un pas de chaque côté : (f(x + h) - f(x - h)) / (2 * h). Les erreurs des deux côtés s'annulent, donc son erreur diminue comme h au carré et elle est exacte pour les polynômes du second degré. La dérivée seconde (gamma) combine trois points : (f(x + h) - 2 * f(x) + f(x - h)) / h**2. Passez la fonction elle-même en argument (f), une compétence des sujets sur les fonctions.

import math

def forward_diff(f, x, h=1e-6):
    return (f(x + h) - f(x)) / h

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def second_derivative(f, x, h=1e-4):
    return (f(x + h) - 2 * f(x) + f(x - h)) / h**2

print(round(central_diff(math.sin, 0.0), 6))            # 1.0  (cos 0)
print(round(central_diff(lambda x: x**3, 2.0), 4))      # 12.0 (3 x^2)
print(round(second_derivative(lambda x: x**3, 2.0), 3)) # 12.0 (6 x)

# A quadratic: central is exact even for a big step, forward is not
print(central_diff(lambda x: x * x, 3.0, h=0.5))  # 6.0
print(forward_diff(lambda x: x * x, 3.0, h=0.5))  # 6.5

Choisir la taille du pas

Un h plus petit réduit l'erreur de troncature de la formule, mais soustraire deux nombres presque égaux fait perdre des chiffres : cette erreur d'arrondi augmente quand h diminue, et pour une dérivée seconde elle augmente plus vite car on divise aussi par h**2. Un h trop petit (comme 1e-12) peut donner n'importe quoi. Un bon compromis est environ 1e-5 pour les dérivées premières en différences centrées et environ 1e-4 pour les dérivées secondes. Dans le code de valorisation réel, « bump and reprice » c'est exactement cela : relancer le valorisateur avec le spot augmenté et diminué d'un petit montant.

import math

def central_diff(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

# True derivative of exp at 0 is 1; watch the error as h shrinks
for h in (1e-1, 1e-3, 1e-5, 1e-8, 1e-12):
    error = abs(central_diff(math.exp, 0.0, h) - 1.0)
    print(f"h={h:g} error={error:.1e}")
# h=0.1 error=1.7e-03
# h=0.001 error=1.7e-07
# h=1e-05 error=1.2e-11
# h=1e-08 error=5.3e-10     (rounding error takes over)
# h=1e-12 error=3.3e-05     (much worse: too small)

Utiliser les dérivées : la méthode de Newton

Une dérivée indique dans quel sens aller pour atteindre une racine, un point où f(x) = 0. La méthode de Newton part d'une estimation x et la remplace sans cesse par x - f(x) / f'(x), en suivant la tangente jusqu'à zéro. Elle converge très vite près d'une racine, c'est pourquoi on l'utilise pour déduire une volatilité implicite d'un prix d'option. Arrêtez-vous quand abs(f(x)) passe sous une tolérance, abandonnez après un nombre maximal d'itérations, et levez une erreur si la dérivée est nulle ou si la boucle ne converge pas, plutôt que de renvoyer un mauvais nombre.

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def newton_root(f, x0, tol=1e-10, max_iter=50):
    x = x0
    for _ in range(max_iter):
        if abs(f(x)) < tol:
            return x
        slope = central_diff(f, x)
        if slope == 0:
            raise ValueError("flat slope: Newton cannot continue")
        x = x - f(x) / slope
    raise ValueError("did not converge")

# The square root of 2 is the positive root of x^2 - 2
root = newton_root(lambda x: x * x - 2, 1.0)
print(round(root, 8))           # 1.41421356

# x^2 + 1 never reaches zero
try:
    newton_root(lambda x: x * x + 1, 1.0)
except ValueError:
    print("no root found")      # no root found

• Avant (f(x+h) - f(x)) / h ; centrée (f(x+h) - f(x-h)) / (2h) ; seconde (f(x+h) - 2f(x) + f(x-h)) / h**2.
• Ni trop grand ni trop petit : environ 1e-5 pour les dérivées premières et 1e-4 pour les secondes.
• Newton : x = x - f(x) / f'(x) jusqu'à ce que abs(f(x)) soit minuscule ; levez une erreur si cela ne converge pas.

À suivre : les intégrales en code, additionner de nombreuses tranches fines pour trouver aires, totaux et espérances.

Lab : ce qui est noté

• Estimer une pente en regardant devant avec forward_diff25 pts
• Estimer une pente des deux côtés avec central_diff25 pts
• Estimer la courbure (gamma) avec second_derivative25 pts
• Trouver où une fonction s'annule avec newton_root25 pts