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 o 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 au} E(S - K)^+ = e^{-r au} int_{K}^{infty} (x-K)f(x)dx ]

其中 (f(x)) 是对数正态分布的密度函数，均值为：

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

方差为：

[s = sigma sqrt{ au} ]

通常 quantlib-python 提供的数值积分方法不接受额外参数，如果计算涉及额外参数，需要做特殊的转换，将额外参数和积分函数“绑定”成为一个单参数函数。

Python 的语言机制非常灵活，可以通过构造实现“函数体”来绑定积分区间和积分函数，积分区间作为类的参数。或者，可以更简单地编写一个返回函数的函数，

例子 2，积分上限采用 (10 imes 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|^{2mu} e^{-x^2} , mu > -0.5)
GaussJacobiIntegration：计算 (int_{-1}^1 f(x)dx) 的Gauss Jacobi 积分；权重函数为 (w(x,alpha, eta) = (1-x)^alpha(1+x)^eta , alpha,eta > 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} ight) 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