إيجاد الجذور: التنصيف ونيوتن

حل المعادلة f(x) = 0

كثير من الأسئلة الكمية تعود إلى إيجاد x الذي تنعدم عنده دالة: العائد الذي يسعّر سندًا، والتقلب الذي يطابق خيارًا، والمعدل الذي يجعل قيمة مشروع صفرًا. الجذر لدالة f هو x يحقق f(x) = 0. وغالبًا لا توجد صيغة، فيقترب منه البرنامج تدريجيًا. وتقوم طريقتان بمعظم العمل. التنصيف بطيء لكنه مضمون؛ وطريقة نيوتن سريعة لكنها تحتاج حذرًا. وكلتاهما تستقبل الدالة نفسها وسيطًا، وهذا ممكن في بايثون لأن الدوال قيم.

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

التنصيف: بطيء ومضمون

إذا كانت f متصلة وكان لـ f(a) وf(b) إشارتان متعاكستان فهناك جذر بينهما. نقيّم المنتصف: أحد النصفين ما زال طرفاه متعاكسَي الإشارة، فيُحصر الجذر فيه. كل خطوة تنصّف المجال أيًّا كانت f، فيكون الخطأ بعد n خطوة هو (b - a) / 2**n ويمكن معرفة عدد الخطوات مسبقًا: ceil(log2((b - a) / tol)). وهذا الضمان دون حاجة إلى مشتقة هو ما يجعل التنصيف الخيار الآمن الافتراضي.

# 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

نيوتن: اتبع المماس

تستبدل طريقة نيوتن بـ f مماسها عند التخمين الحالي وتقفز إلى حيث يقطع هذا المستقيم الصفر: x_new = x - f(x) / df(x). وقرب جذر بسيط يتربع الخطأ تقريبًا في كل خطوة فيتضاعف عدد الأرقام الصحيحة: من 1 إلى 2 فـ4 فـ8 فـ16 في أربع خطوات. وتحتاج إلى المشتقة df (وهي في التقلب الضمني فيغا التي لديك أصلًا) وإلى بداية جيدة.

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

متى تفشل طريقة نيوتن

المماس صورة محلية. فإن كانت المشتقة صفرًا أو صغيرة جدًا كانت القفزة هائلة أو غير معرَّفة؛ وإن كانت البداية بعيدة عن الجذر فقد تبتعد التقديرات أو تدور (في atan(x) من x = 2 تنمو بلا حد). لذلك يحمي المحلّ المتين القسمة ويحدّ عدد التكرارات ويطلق خطأً بدل أن يعيد رقمًا خاطئًا. وتبدأ مكتبات كثيرة من مجال محصور وتستخدم نيوتن داخله وتعود إلى التنصيف عندما تخرج خطوة منه.

• جذر f هو x يحقق f(x) = 0؛ وتستقبل المحلّات f (وdf) كوسائط دالة.
• التنصيف: أبقِ النصف الذي تتغير فيه الإشارة؛ ceil(log2((b - a) / tol)) خطوة مضمونة ودون مشتقة.
• نيوتن: x - f(x) / df(x)؛ تتضاعف الأرقام كل خطوة قرب جذر بسيط، لكنه يحتاج df وبداية جيدة وحمايات (مشتقة صفرية وعدم تقارب).

التالي: الأشجار ثنائية الحدّ، تسعير خيار بجعل السهم يصعد ويهبط ثم العمل رجوعًا.

المختبر: ما الذي يُقيَّم

• احصر جذرًا بالتنصيف باستخدام bisection25 نقطة
• احسب عدد التنصيفات باستخدام bisection_steps25 نقطة
• اتبع المماس مرة واحدة باستخدام newton_step25 نقطة
• تقارب نحو جذر باستخدام newton25 نقطة