Nullstellensuche: Bisektion und Newton

f(x) = 0 lösen

Viele quantitative Fragen laufen darauf hinaus, das x zu finden, an dem eine Funktion null wird: die Rendite, die eine Anleihe bewertet, die Volatilität, die zu einer Option passt, der Zins, der den Wert eines Projekts null macht. Eine Nullstelle von f ist ein x mit f(x) = 0. Meist gibt es keine Formel, also nähert sich ein Programm an. Zwei Verfahren leisten fast die ganze Arbeit. Die Bisektion ist langsam, aber sicher; das Newton-Verfahren ist schnell, braucht aber Sorgfalt. Beide erhalten die Funktion selbst als Argument, was Python erlaubt, weil Funktionen Werte sind.

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

Bisektion: langsam und sicher

Ist f stetig und haben f(a) und f(b) entgegengesetzte Vorzeichen, liegt dazwischen eine Nullstelle. Man wertet die Mitte aus: Eine der beiden Hälften hat an ihren Enden weiter entgegengesetzte Vorzeichen, dort ist die Nullstelle also eingeschlossen. Jeder Schritt halbiert das Intervall, egal wie f aussieht; der Fehler nach n Schritten ist (b - a) / 2**n, und man weiß vorher, wie viele Schritte nötig sind: ceil(log2((b - a) / tol)). Diese Garantie ohne Ableitung macht die Bisektion zur sicheren Standardwahl.

# 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: der Tangente folgen

Das Newton-Verfahren ersetzt f durch seine Tangente an der aktuellen Näherung und springt dorthin, wo diese Gerade null schneidet: x_new = x - f(x) / df(x). In der Nähe einer einfachen Nullstelle wird der Fehler je Schritt ungefähr *quadriert*, die Zahl der richtigen Stellen verdoppelt sich also: von 1 auf 2, 4, 8, 16 in vier Schritten. Es braucht die Ableitung df (bei der impliziten Volatilität ist das Vega, das Sie schon haben) und einen guten Startwert.

# 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))

Wenn Newton scheitert

Die Tangente ist ein lokales Bild. Ist die Ableitung null oder winzig, ist der Sprung riesig oder undefiniert; liegt der Start weit von der Nullstelle, können die Iterierten wegspringen oder kreisen (für atan(x) ab x = 2 wachsen sie unbegrenzt). Ein robuster Löser schützt also die Division, begrenzt die Iterationen und löst einen Fehler aus, statt eine falsche Zahl zu liefern. Viele Bibliotheken starten mit einer Einschließung, nutzen Newton darin und fallen auf Bisektion zurück, wenn ein Schritt sie verlässt.

• Eine Nullstelle von f ist ein x mit f(x) = 0; Löser erhalten f (und df) als Funktionsargumente.
• Bisektion: die Hälfte mit Vorzeichenwechsel behalten; ceil(log2((b - a) / tol)) Schritte, garantiert, ohne Ableitung.
• Newton: x - f(x) / df(x); die Stellen verdoppeln sich je Schritt nahe einer einfachen Nullstelle, aber es braucht df, einen guten Start und Schutz (Ableitung null, keine Konvergenz).

Als Nächstes: Binomialbäume, eine Option bewerten, indem die Aktie auf- und abspringt und man rückwärts rechnet.

Lab: Was bewertet wird

• Mit bisection eine Nullstelle durch Halbieren einfangen25 Pkt.
• Mit bisection_steps die Halbierungen zählen25 Pkt.
• Mit newton_step einmal der Tangente folgen25 Pkt.
• Mit newton gegen eine Nullstelle konvergieren25 Pkt.