Recherche de racines : bissection et Newton

Résoudre f(x) = 0

Beaucoup de questions quantitatives reviennent à trouver le x où une fonction s'annule : le rendement qui valorise une obligation, la volatilité qui correspond à une option, le taux qui rend nulle la valeur d'un projet. Une racine de f est un x tel que f(x) = 0. Il n'existe en général pas de formule, donc un programme s'en approche. Deux méthodes font presque tout le travail. La dichotomie est lente mais sûre ; la méthode de Newton est rapide mais demande de l'attention. Toutes deux reçoivent la fonction elle-même en argument, ce que Python permet puisque les fonctions sont des valeurs.

def bisection(f, a, b, tol=1e-10):
    if f(a) * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    while b - a > tol:
        mid = (a + b) / 2
        if f(a) * f(mid) <= 0:       # the root is in the left half
            b = mid
        else:
            a = mid
    return (a + b) / 2

root = bisection(lambda x: x * x - 2, 0, 2)
print(root)                          # 1.414213562355144
print(root * root)                   # 1.9999999999492266

Dichotomie : lente et sûre

Si f est continue et que f(a) et f(b) sont de signes opposés, une racine se trouve entre les deux. On évalue le milieu : l'une des deux moitiés a encore des signes opposés à ses extrémités, la racine y est donc piégée. Chaque étape divise l'intervalle par deux, quelle que soit f, donc l'erreur après n étapes est (b - a) / 2**n et on peut dire à l'avance combien d'étapes il faut : ceil(log2((b - a) / tol)). Cette garantie, sans dérivée, fait de la dichotomie le choix sûr par défaut.

# How many halvings does a given tolerance cost? About 3.3 per decimal digit
import math

def bisection_steps(a, b, tol):
    if b - a <= tol:
        return 0
    return math.ceil(math.log2((b - a) / tol))

for tol in (1e-2, 1e-4, 1e-8, 1e-12):
    print(tol, bisection_steps(0, 2, tol))       # 1e-2: 8, 1e-4: 15, 1e-8: 28, 1e-12: 41

Newton : suivre la tangente

La méthode de Newton remplace f par sa tangente au point courant et saute là où cette droite coupe zéro : x_new = x - f(x) / df(x). Près d'une racine simple, l'erreur est à peu près *élevée au carré* à chaque étape, donc le nombre de chiffres exacts double : de 1 à 2, 4, 8, 16 en quatre étapes. Elle exige la dérivée df (pour la volatilité implicite, c'est le véga, que vous avez déjà) et un bon point de départ.

# Newton for sqrt(2): print every iterate and the error that remains
import math

def newton_step(f, df, x):
    return x - f(x) / df(x)

f = lambda x: x * x - 2
df = lambda x: 2 * x
x = 1.0
for step in range(6):
    print(step, x, abs(x - math.sqrt(2)))
    x = newton_step(f, df, x)

# Newton on the implied volatility of a call, using vega as the derivative
import math

def norm_cdf(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def norm_pdf(x):
    return math.exp(-0.5 * x * x) / math.sqrt(2 * math.pi)

def bs_call(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))

def vega(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_pdf(d1) * math.sqrt(t)

market = 8.0
sigma = 0.5                                              # a rough start
for step in range(1, 6):
    error = bs_call(100, 100, 0.05, sigma, 0.5) - market
    sigma -= error / vega(100, 100, 0.05, sigma, 0.5)
    print(step, round(sigma, 10))

Quand Newton échoue

La tangente est une image locale. Si la dérivée est nulle ou minuscule, le saut est énorme ou indéfini ; si le départ est loin de la racine, les itérés peuvent s'éloigner ou osciller (pour atan(x) depuis x = 2, ils croissent sans limite). Un solveur robuste protège donc la division, limite le nombre d'itérations et lève une erreur plutôt que de renvoyer un nombre faux. Beaucoup de bibliothèques partent d'un encadrement, utilisent Newton à l'intérieur et reviennent à la dichotomie quand un pas en sort.

• Une racine de f est un x avec f(x) = 0 ; les solveurs prennent f (et df) comme arguments fonctions.
• Dichotomie : garder la moitié où le signe change ; ceil(log2((b - a) / tol)) étapes, garanties, sans dérivée.
• Newton : x - f(x) / df(x) ; les chiffres doublent à chaque étape près d'une racine simple, mais il faut df, un bon départ et des garde-fous (dérivée nulle, non-convergence).

À suivre : les arbres binomiaux, valoriser une option en faisant monter et descendre l'action puis en remontant à rebours.

Lab : ce qui est noté

• Piéger une racine en coupant en deux avec bisection25 pts
• Compter les divisions par deux avec bisection_steps25 pts
• Suivre la tangente une fois avec newton_step25 pts
• Converger vers une racine avec newton25 pts