QuantLib 金融計算——數學工具之求解器

目錄

• QuantLib 金融計算——數學工具之求解器
  • 概述
    • 調用方式
  • 非 Newton 算法（不須要導數）
  • Newton 算法（須要導數）

若是未作特別說明，文中的程序都是 Python3 代碼。

QuantLib 金融計算——數學工具之求解器

載入模塊

import QuantLib as ql
import scipy
from scipy.stats import norm

print(ql.__version__)

1.12

概述

QuantLib 提供了多種類型的一維求解器，用以求解單參數函數的根，

\[ f(x)=0 \]

其中 \(f : R \to R\) 是實數域上的函數。

QuantLib 提供的求解器類型有：

• Brent
• Bisection
• Secant
• Ridder
• Newton（要求提供成員函數 derivative，計算導數）
• FalsePosition

這些求解器的構造函數均爲默認構造函數，不接受參數。例如，Brent 求解器實例的構造語句爲 mySolv = Brent()。

調用方式

求解器的成員函數 solve 有兩種調用方式：

solve(f,
      accuracy,
      guess,
      step)

solve(f,
      accuracy,
      guess,
      xMin,
      xMax)

• f：單參數函數或函數對象，返回值爲一個浮點數。
• accuracy：浮點數，表示求解精度 \(\epsilon\)，用於中止計算。假設 \(x_i\) 是根的準確解，
  • 當 \(|f(x)| < \epsilon\)；
  • 或 \(|x - x_i| < \epsilon\) 時中止計算。
• guess：浮點數，對根的初始猜想值。
• step：浮點數，在第一種調用方式中，沒有限定根的區間範圍，算法須要本身搜索，肯定一個範圍。step 規定了搜索算法的步長。
• xMin、xMax：浮點數，左右區間範圍

根求解器在量化金融中最經典的應用是求解隱含波動率。給按期權價格 \(p\) 以及其餘參數 \(S_0\)、\(K\)、\(r_d\)、\(r_f\)、\(\tau\)，咱們要計算波動率 \(\sigma\)，知足

\[ f(\sigma) = \mathrm{blackScholesPrice}(S_0 , K, r_d , r_f , \sigma , \tau, \phi) - p = 0 \]

其中 Black-Scholes 函數中 \(\phi = 1\) 表明看漲期權；\(\phi = −1\) 表明看跌期權。

非 Newton 算法（不須要導數）

下面的例子顯示瞭如何加一個多參數函數包裝爲一個單參數函數，並使用 QuantLib 求解器計算隱含波動率。

例子 1

# Black-Scholes 函數
def blackScholesPrice(spot,
                      strike,
                      rd,
                      rf,
                      vol,
                      tau,
                      phi):
    domDf = scipy.exp(-rd * tau)
    forDf = scipy.exp(-rf * tau)
    fwd = spot * forDf / domDf
    stdDev = vol * scipy.sqrt(tau)

    dp = (scipy.log(fwd / strike) + 0.5 * stdDev * stdDev) / stdDev
    dm = (scipy.log(fwd / strike) - 0.5 * stdDev * stdDev) / stdDev

    res = phi * domDf * (fwd * norm.cdf(phi * dp) - strike * norm.cdf(phi * dm))

    return res


# 包裝函數
def impliedVolProblem(spot,
                      strike,
                      rd,
                      rf,
                      tau,
                      phi,
                      price):
    def inner_func(v):
        return blackScholesPrice(spot, strike, rd, rf, v, tau, phi) - price

    return inner_func


def testSolver1():
    # setup of market parameters
    spot = 100.0
    strike = 110.0
    rd = 0.002
    rf = 0.01
    tau = 0.5
    phi = 1
    vol = 0.1423

    # calculate corresponding Black Scholes price

    price = blackScholesPrice(spot, strike, rd, rf, vol, tau, phi)
    # setup a solver
    mySolv1 = ql.Bisection()
    mySolv2 = ql.Brent()
    mySolv3 = ql.Ridder()

    accuracy = 0.00001
    guess = 0.25

    min = 0.0
    max = 1.0

    myVolFunc = impliedVolProblem(spot, strike, rd, rf, tau, phi, price)

    res1 = mySolv1.solve(myVolFunc, accuracy, guess, min, max)
    res2 = mySolv2.solve(myVolFunc, accuracy, guess, min, max)
    res3 = mySolv3.solve(myVolFunc, accuracy, guess, min, max)

    print('{0:<35}{1}'.format('Input Volatility:', vol))
    print('{0:<35}{1}'.format('Implied Volatility Bisection:', res1))
    print('{0:<35}{1}'.format('Implied Volatility Brent:', res2))
    print('{0:<35}{1}'.format('Implied Volatility Ridder:', res3))


testSolver1()

# Input Volatility:                  0.1423
# Implied Volatility Bisection:      0.14229583740234375
# Implied Volatility Brent:          0.14230199334812577
# Implied Volatility Ridder:         0.1422999996313447

Newton 算法（須要導數）

Newton 算法要求爲根求解器提供 \(f(\sigma)\) 的導數 \(\frac{\partial f}{\partial \sigma}\)（即 vega）。下面的例子顯示瞭如何將導數添加進求解隱含波動率的過程。爲此咱們須要一個類，一方面提供做爲一個函數對象，另外一方面要提供成員函數 derivative。

例子 2

class BlackScholesClass:
    def __init__(self,
                 spot,
                 strike,
                 rd,
                 rf,
                 tau,
                 phi,
                 price):
        self.spot_ = spot
        self.strike_ = strike
        self.rd_ = rd
        self.rf_ = rf
        self.phi_ = phi
        self.tau_ = tau
        self.price_ = price
        self.sqrtTau_ = scipy.sqrt(tau)
        self.d_ = norm

        self.domDf_ = scipy.exp(-self.rd_ * self.tau_)
        self.forDf_ = scipy.exp(-self.rf_ * self.tau_)
        self.fwd_ = self.spot_ * self.forDf_ / self.domDf_
        self.logFwd_ = scipy.log(self.fwd_ / self.strike_)

    def blackScholesPrice(self,
                          spot,
                          strike,
                          rd,
                          rf,
                          vol,
                          tau,
                          phi):
        domDf = scipy.exp(-rd * tau)
        forDf = scipy.exp(-rf * tau)
        fwd = spot * forDf / domDf
        stdDev = vol * scipy.sqrt(tau)

        dp = (scipy.log(fwd / strike) + 0.5 * stdDev * stdDev) / stdDev
        dm = (scipy.log(fwd / strike) - 0.5 * stdDev * stdDev) / stdDev

        res = phi * domDf * (fwd * norm.cdf(phi * dp) - strike * norm.cdf(phi * dm))

        return res

    def impliedVolProblem(self,
                          spot,
                          strike,
                          rd,
                          rf,
                          vol,
                          tau,
                          phi,
                          price):
        return self.blackScholesPrice(
            spot, strike, rd, rf, vol, tau, phi) - price

    def __call__(self,
                 x):
        return self.impliedVolProblem(
            self.spot_, self.strike_, self.rd_, self.rf_,
            x,
            self.tau_, self.phi_, self.price_)

    def derivative(self,
                   x):
        # vega
        stdDev = x * self.sqrtTau_
        dp = (self.logFwd_ + 0.5 * stdDev * stdDev) / stdDev
        return self.spot_ * self.forDf_ * self.d_.pdf(dp) * self.sqrtTau_


def testSolver2():

    # setup of market parameters
    spot = 100.0
    strike = 110.0
    rd = 0.002
    rf = 0.01
    tau = 0.5
    phi = 1
    vol = 0.1423

    # calculate corresponding Black Scholes price
    price = blackScholesPrice(
        spot, strike, rd, rf, vol, tau, phi)
    solvProblem = BlackScholesClass(
        spot, strike, rd, rf, tau, phi, price)

    mySolv = ql.Newton()

    accuracy = 0.00001
    guess = 0.10
    step = 0.001

    res = mySolv.solve(
        solvProblem, accuracy, guess, step)

    print('{0:<20}{1}'.format('Input Volatility:', vol))
    print('{0:<20}{1}'.format('Implied Volatility:', res))


testSolver2()

# Input Volatility:   0.1423
# Implied Volatility: 0.14230000000000048

導數的使用明顯提升了精度。