[TOC]
如果未做特别说明，文中的程序都是 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
