QuantLib 金融计算——数学工具之随机数发生器
如果未做特别说明，文中的程序都是 Python3 代码。
QuantLib 金融计算——数学工具之随机数发生器
载入模块
import QuantLib as ql
import scipy
print(ql.__version__)
1.12
概述
随机模拟通常从产生均匀分布的随机数开始。假设 \(X \sim U [0, 1]\) 是均匀分布的随机变量。任意分布的随机数通常需要对 \(X\) 施加某种变换得到，一般情况下是用累积分布函数的逆函数 \(F^{−1}\)，\(F^{−1}(X)\) 的分布就是 \(F\)。其他的变换算法可能不需要 \(F^{−1}\)，比如用于生成正态分布的 Box Muller 变换算法。
均匀分布的随机数发生器主要分两种：
伪随机数
quantlib-python 提供了以下三种均匀分布的（伪）随机数发生器：
KnuthUniformRng，高德纳（Knuth）算法LecuyerUniformRng，L'Ecuyer 算法MersenneTwisterUniformRng，著名的梅森旋转（Mersenne-Twister）算法
随机数发生器的构造函数，
Rng(seed)
其中
seed，整数，默认值是 0，作为种子用于初始化相应的确定性序列；
随机数发生器的成员函数：
next()：返回一个SampleNumber对象，作为模拟的结果。
r = rng.next()
v = r.value(r)
用户通过反复调用成员函数 next() 获得一连串的随机数，需要注意的是 r 的类型是 SampleNumber，需要调用 value() 得到对应的浮点数。
例子 1，
def testingRandomNumbers1():
seed = 1
unifMt = ql.MersenneTwisterUniformRng(seed)
unifLec = ql.LecuyerUniformRng(seed)
unifKnuth = ql.KnuthUniformRng(seed)
print('{0:<25}{1:<25}{2:<25}'.format(
'Mersenne Twister', 'Lecuyer', 'Knut'))
for i in range(10):
print('{0:<25}{1:<25}{2:<25}'.format(
unifMt.next().value(),
unifLec.next().value(),
unifKnuth.next().value()))
testingRandomNumbers1()
Mersenne Twister Lecuyer Knut
0.41702199855353683 0.2853808990946861 0.4788952510312594
0.9971848082495853 0.2533581892659171 0.7694635535665499
0.7203244894044474 0.09346853100919404 0.47721285286866455
0.9325573613168672 0.6084968907396475 0.15752737762851
0.00011438119690865278 0.90342026007861 0.6065713927733087
正态分布（伪）随机数
随机模拟中最常见的分布是正态分布，quantlib-python 提供的正态分布随机数发生器有 4 类：
CentralLimitABCGaussianRngBoxMullerABCGaussianRngMoroInvCumulativeABCGaussianRngInvCumulativeABCGaussianRng
其中 ABC 特指一种均匀随机数发生器。
具体来讲 4 类发生器分为 12 种：
CentralLimitLecuyerGaussianRngCentralLimitKnuthGaussianRngCentralLimitMersenneTwisterGaussianRngBoxMullerLecuyerGaussianRngBoxMullerKnuthGaussianRngBoxMullerMersenneTwisterGaussianRngMoroInvCumulativeLecuyerGaussianRngMoroInvCumulativeKnuthGaussianRngMoroInvCumulativeMersenneTwisterGaussianRngInvCumulativeLecuyerGaussianRngInvCumulativeKnuthGaussianRngInvCumulativeMersenneTwisterGaussianRng
随机数发生器的构造函数：
rng = Rng(seed)
grng = Gaussianrng(rng)
正态分布随机数发生器接受一个对应的均匀分布随机数发生器作为源，以 BoxMullerMersenneTwisterGaussianRng 为例，需要配置一个 MersenneTwisterUniformRng 对象作为随机数的源，使用经典的 Box-Muller 算法得到正态分布随机数。
例子 2，
def testingRandomNumbers2():
seed = 12324
unifMt = ql.MersenneTwisterUniformRng(seed)
bmGauss = ql.BoxMullerMersenneTwisterGaussianRng(unifMt)
for i in range(5):
print(bmGauss.next().value())
testingRandomNumbers2()
-1.1756781173398896
0.14110041851886157
1.569582906805544
-0.026736779238941934
-0.8220676600472409
拟随机数
相较于之前描述的“伪”随机数，随机模拟中另一类重要的随机数成为“拟”随机数，也称为低偏差序列。因为收敛性更好，拟随机数通常用于高维随机变量的模拟。quantlib-python 提供的拟随机数有两类，
HaltonRsg: Halton 序列SobolRsg: Sobol 序列
HaltonRsg
HaltonRsg 的构造函数，
HaltonRsg(dimensionality,
seed,
randomStart,
randomShift)
其中，
dimensionality：整数，设置维度；seed，整数，默认值是 0，作为种子用于初始化相应的确定性序列；randomStart：布尔值，默认是True，是否随机开始；randomShift：布尔值，默认是False，是否随机平移。
HaltonRsg 的成员函数，
nextSequence()：返回一个SampleRealVector对象，作为模拟的结果；lastSequence()：返回一个SampleRealVector对象，作为上一个模拟的结果；dimension()：返回维度。
SobolRsg
SobolRsg 的构造函数，
SobolRsg(dimensionality,
seed,
directionIntegers=Jaeckel)
其中，
dimensionality：整数，设置维度；seed，整数，默认值是 0，作为种子用于初始化相应的确定性序列；directionIntegers，quantlib-python 的内置变量，默认值是SobolRsg.Jaeckel，用于 Sobol 序列的初始化。
SobolRsg 的成员函数，
nextSequence()：返回一个SampleRealVector对象，作为模拟的结果；lastSequence()：返回一个SampleRealVector对象，作为上一个模拟的结果；dimension()：返回维度。skipTo(n)：n是整数，跳转到抽样结果的第 n 个维度；nextInt32Sequence()：返回一个IntVector对象。
例子 3，
def testingRandomNumbers4():
dim = 5
haltonGen = ql.HaltonRsg(dim)
sobolGen = ql.SobolRsg(dim)
sampleHalton = haltonGen.nextSequence().value()
sampleSobol = sobolGen.nextSequence().value()
print('{0:<25}{1:<25}'.format(
'Halton', 'Sobol'))
for i in range(dim):
print('{0:<25}{1:<25}'.format(
sampleHalton[i],
sampleSobol[i]))
testingRandomNumbers4()
Halton Sobol
0.04081786540336907 0.5
0.8535710143553551 0.5
0.69400573329408 0.5
0.818105927979147 0.5
0.878826694887864 0.5
两类随机数的收敛性比较
最后用一个例子比较两类随机数的收敛性，分别产生正态分布的伪随机数和拟随机数，计算分布的四个统计指标：
- 均值（理论值等于 0.0）；
- 方差（理论值等于 1.0）；
- 偏度（理论值等于 0.0）；
- 超额峰度（理论值等于 0.0）。
def testingRandomNumbers5():
sobolGen = ql.SobolRsg(1)
seed = 12324
unifMt = ql.MersenneTwisterUniformRng(seed)
bmGauss = ql.BoxMullerMersenneTwisterGaussianRng(unifMt)
boxMullerStat = ql.IncrementalStatistics()
sobolStat = ql.IncrementalStatistics()
invGauss = ql.MoroInverseCumulativeNormal()
numSim = 10000
for j in range(numSim):
boxMullerStat.add(bmGauss.next().value())
currSobolNum = sobolGen.nextSequence().value()[0]
sobolStat.add(invGauss(currSobolNum))
stats = {
"BoxMuller Mean:": boxMullerStat.mean(),
"Sobol Mean:": sobolStat.mean(),
"BoxMuller Var:": boxMullerStat.variance(),
"Sobol Var:": sobolStat.variance(),
"BoxMuller Skew:": boxMullerStat.skewness(),
"Sobol Skew:": sobolStat.skewness(),
"BoxMuller Kurtosis:": boxMullerStat.kurtosis(),
"Sobol Kurtosis:": sobolStat.kurtosis()}
for k, v in stats.items():
print('{0:>25}{1:>30}'.format(k, v))
testingRandomNumbers5()
BoxMuller Mean: 0.005966482725988245
Sobol Mean: -0.0002364019095203635
BoxMuller Var: 1.0166044844467006
Sobol Var: 0.9986010126883317
BoxMuller Skew: 0.02100635339070779
Sobol Skew: -7.740573185322994e-05
BoxMuller Kurtosis: -0.0340476839897507
Sobol Kurtosis: -0.020768126049145776
直观上看 Sobol 序列的结果更加接近理论值，这证明使用拟随机数做模拟的收敛速度更好。
扩展阅读
★ 持续学习 ★ 坚持创作 ★
