QuantLib 金融計算——數學工具之數值積分

目錄

• QuantLib 金融計算——數學工具之數值積分
  • 概述
  • 常見積分方法
  • 高斯積分

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

QuantLib 金融計算——數學工具之數值積分

載入模塊

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

print(ql.__version__)

1.12

概述

quantlib-python 提供了許多方法計算標量函數 \(f : R \to R\) 在閉區間上的積分：

\[ \int_a^b f(x) dx \]

對於主要的積分方法，必須提供兩個參數：

• 絕對精度：若是當前計算結果和前一個計算結果的差小於精度，則中止計算。
• 最大計算次數：若是達到最大計算次數，則中止計算。

對於某些特殊的數值積分，例如高斯積分，還須要提供其餘額外參數。

常見積分方法

首先討論最普通最多見的一類數值積分，quantlib-python 提供了下列方法：

• TrapezoidIntegralMidPoint
• SimpsonIntegral
• GaussLobattoIntegral
• GaussKronrodAdaptive
• GaussKronrodNonAdaptive

這些方法在通常的數值分析教科書中都有詳細的討論。在 quantlib-python 中，上述數值積分器對象的構造方式是相同的，以下：

myIntegrator = ql.XXXintegrator(absoluteAccuracy,
                                maxEvaluations)

計算閉區間 \([a, b]\) 上的積分值：

myIntegrator(f, a, b)

其中 f 是一個「單參數」函數，返回一個浮點數。

例子 1，標準正態密度函數上的積分

def testIntegration1():
    absAcc = 0.00001
    maxEval = 1000

    a = 0.0
    b = scipy.pi

    numInt1 = ql.TrapezoidIntegralMidPoint(absAcc, maxEval)
    numInt2 = ql.SimpsonIntegral(absAcc, maxEval)
    numInt3 = ql.GaussLobattoIntegral(maxEval, absAcc)
    numInt4 = ql.GaussKronrodAdaptive(absAcc, maxEval)
    numInt5 = ql.GaussKronrodNonAdaptive(absAcc, maxEval, absAcc)

    analytical = norm.cdf(b) - norm.cdf(a)

    print('{0:<30}{1}'.format('Analytical:', analytical))
    print('{0:<30}{1}'.format('Midpoint Trapezoidal:', numInt1(norm.pdf, a, b)))
    print('{0:<30}{1}'.format('Simpson:', numInt2(norm.pdf, a, b)))
    print('{0:<30}{1}'.format('Gauss Lobatto:', numInt3(norm.pdf, a, b)))
    print('{0:<30}{1}'.format('Gauss Kronrod Adpt:', numInt4(norm.pdf, a, b)))
    print('{0:<30}{1}'.format('Gauss Kronrod Non Adpt:', numInt5(norm.pdf, a, b)))

testIntegration1()

Analytical:                   0.4991598418317367
Midpoint Trapezoidal:         0.4991643496589137
Simpson:                      0.4991598398355923
Gauss Lobatto:                0.49916005276697556
Gauss Kronrod Adpt:           0.49915984183173506
Gauss Kronrod Non Adpt:       0.4991598418317367

全部結果幾乎是一致的。

下面是一個更復雜的例子，直接從歐式看漲期權的積分形式近似計算期權價格。

敲訂價格爲 \(K\) 的看漲期權的積分形式爲：

\[ e^{-r \tau} E(S - K)^+ = e^{-r \tau} \int_{K}^{\infty} (x-K)f(x)dx \]

其中 \(f(x)\) 是對數正態分佈的密度函數，均值爲：

\[ \log(S_0) + (r + \frac{1}{2} \sigma^2)\tau \]

方差爲：

\[ s = \sigma \sqrt{\tau} \]

一般 quantlib-python 提供的數值積分方法不接受額外參數，若是計算涉及額外參數，須要作特殊的轉換，將額外參數和積分函數「綁定」成爲一個單參數函數。

Python 的語言機制很是靈活，能夠經過構造實現「函數體」來綁定積分區間和積分函數，積分區間做爲類的參數。或者，能夠更簡單地編寫一個返回函數的函數，

例子 2，積分上限採用 \(10 \times K\)

def callFunc(spot,
             strike,
             r,
             vol,
             tau):
    mean = scipy.log(spot) + (r - 0.5 * vol * vol) * tau
    stdDev = vol * scipy.sqrt(tau)

    def inner_func(x):
        return (x - strike) * \
               lognorm.pdf(
                   x, stdDev, loc=0, scale=scipy.exp(mean)) * \
               scipy.exp(-r * tau)

    return inner_func

其中，內部函數 inner_func 做爲對象被返回，inner_func 是一個單參數函數。

def testIntegration4():
    spot = 100.0
    r = 0.03
    tau = 0.5
    vol = 0.20
    strike = 110.0

    a = strike
    b = strike * 10.0

    ptrF = callFunc(spot, strike, r, vol, tau)

    absAcc = 0.00001
    maxEval = 1000
    numInt = ql.SimpsonIntegral(absAcc, maxEval)

    print("Call Value: ", numInt(ptrF, a, b))


testIntegration4()

與標準 Black-Scholes 公式得出的結果幾乎一致。

Call Value:  2.611902550625855

高斯積分

一般，一個 n 點高斯求積經過選取合適的 \(x_i\) 和 \(w_i\)（\(i = 1, ..., n\)）產生 2n − 1 階（或較低階）多項式的準確積分值構造出來，

\[ \int_{-1}^1 f(x)dx \approx \sum_{i=1}^n w_if(x_i) \]

存在不一樣類型的權重函數和區間形式，quantlib-python 提供了以下幾種：

• GaussLaguerreIntegration：計算 \(\int_0^{\infty} f(x)dx\) 的廣義 Gauss Laguerre 積分；權重函數爲 \(w(x,s) := x^s e^{-x} , s>-1\)
• GaussHermiteIntegration：計算 \(\int_{-\infty}^{\infty} f(x)dx\) 的 Gauss Hermite 積分；權重函數爲 \(w(x,\mu) = |x|^{2\mu} e^{-x^2} , \mu > -0.5\)
• GaussJacobiIntegration：計算 \(\int_{-1}^1 f(x)dx\) 的Gauss Jacobi 積分；權重函數爲 \(w(x,\alpha, \beta) = (1-x)^\alpha(1+x)^\beta , \alpha,\beta > 1\)
• GaussHyperbolicIntegration：計算 \(\int_{-\infty}^{\infty} f(x)dx\) 的高斯雙曲積分；權重函數爲 \(w(x) = \frac{1}{\cosh(x)}\)
• GaussLegendreIntegration：計算 \(\int_{-1}^1 f(x)dx\) 的 Gauss Legendre 積分；權重函數爲 \(w(x)=1\)
• GaussChebyshevIntegration：計算 \(\int_{-1}^1 f(x)dx\) 的第一類 Gauss Chebyshev 積分；權重函數爲\(w(x) = \sqrt{(1-x^2)}\)
• GaussChebyshev2ndIntegration：計算 \(\int_{-1}^1 f(x)dx\) 的第二類 Gauss Legendre 積分；權重函數爲 \(w(x, \lambda) = (1+x^2)^{\lambda - 1/2}\)

例子 3

def testIntegration2():
    gLagInt = ql.GaussLaguerreIntegration(16)  # [0,\infty]
    gHerInt = ql.GaussHermiteIntegration(16)  # (-\infty, \infty)
    gChebInt = ql.GaussChebyshevIntegration(64)  # (-1, 1)
    gChebInt2 = ql.GaussChebyshev2ndIntegration(64)  # (-1, 1)

    analytical = norm.cdf(1) - norm.cdf(-1)

    print('{0:<15}{1}'.format("Laguerre:", gLagInt(norm.pdf)))
    print('{0:<15}{1}'.format("Hermite:", gHerInt(norm.pdf)))
    print('{0:<15}{1}'.format("Analytical:", analytical))
    print('{0:<15}{1}'.format("Cheb:", gChebInt(norm.pdf)))
    print('{0:<15}{1}'.format("Cheb 2 kind:", gChebInt2(norm.pdf)))

Laguerre:      0.49999230923944715
Hermite:       0.9999999834745512
Analytical:    0.6826894921370859
Cheb:          0.6827380724493052
Cheb 2 kind:   0.682595292164792

一般 quantlib-python 提供的高斯積分方法只針對固定的區間，例如 \([-1,1]\)，若是須要計算其餘區間上的積分，須要作特殊的轉換，將積分區間和積分函數「綁定」成爲一個單參數函數。區間 \([−1, 1]\) 向 \([a, b]\) 的轉換至關簡單

\[ \int_a^b f(x)dx = \frac{b-a}{2} \int_{-1}^1f \left(\frac{b-a}{2}x + \frac{b+a}{2}\right) dx \]

相似以前的作法，

def Func(f, a, b):
    t1 = 0.5 * (b - a)
    t2 = 0.5 * (b + a)

    def inner_func(x):
        return t1 * f(t1 * x + t2)

    return inner_func

例子 4

def testIntegration3():
    a = -1.96
    b = 1.96

    gChebInt = ql.GaussChebyshevIntegration(64)

    analytical = norm.cdf(b) - norm.cdf(a)
    f = Func(norm.pdf, a, b)

    print('{0:<15}{1}'.format("Analytical:", analytical))
    print('{0:<15}{1}'.format("Chebyshev:", gChebInt(f)))


testIntegration3()

Analytical:    0.950004209703559
Chebyshev:     0.9500271929144378