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

• 绝对精度：如果当前计算结果和前一个计算结果的差小于精度，则停止计算。
• 最大计算次数：如果达到最大计算次数，则停止计算。

## 常见积分方法

• TrapezoidIntegralMidPoint
• SimpsonIntegral
• GaussLobattoIntegral
• GaussKronrodAdaptive
• GaussKronrodNonAdaptive

myIntegrator = ql.XXXintegrator(absoluteAccuracy,
maxEvaluations)


myIntegrator(f, a, b)


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)

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


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

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

$s = \sigma \sqrt{\tau}$

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

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


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


Call Value:  2.611902550625855


## 高斯积分

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

• 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}$$

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


$\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


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


## 扩展阅读

《QuantLib 金融计算》系列合集

posted @ 2018-09-24 16:31  xuruilong100  阅读(1254)  评论(0编辑  收藏  举报