使用sobol序列使用准随机标准正态数的Python中的蒙特卡罗模拟给出了错误的值

Question

我正在尝试使用准随机标准正态数进行蒙特卡罗模拟。我知道我们可以使用Sobol序列生成统一数字，然后使用概率积分变换将它们转换为标准正态数。我的代码给出了模拟资产路径的不切实际的值：

import sobol_seq
import numpy as np
from scipy.stats import norm

def i4_sobol_generate_std_normal(dim_num, n, skip=1):
    """
    Generates multivariate standard normal quasi-random variables.
    Parameters:
      Input, integer dim_num, the spatial dimension.
      Input, integer n, the number of points to generate.
      Input, integer SKIP, the number of initial points to skip.
      Output, real np array of shape (n, dim_num).
    """

    sobols = sobol_seq.i4_sobol_generate(dim_num, n, skip)

    normals = norm.ppf(sobols)

    return normals

def GBM(Ttm, TradingDaysInAYear, NoOfPaths, UnderlyingPrice, RiskFreeRate, Volatility):
    dt = float(Ttm) / TradingDaysInAYear
    paths = np.zeros((TradingDaysInAYear + 1, NoOfPaths), np.float64)
    paths[0] = UnderlyingPrice
    for t in range(1, TradingDaysInAYear + 1):
        rand = i4_sobol_generate_std_normal(1, NoOfPaths)
        lRand = []
        for i in range(len(rand)):
            a = rand[i][0]
            lRand.append(a)
        rand = np.array(lRand)

        paths[t] = paths[t - 1] * np.exp((RiskFreeRate - 0.5 * Volatility ** 2) * dt + Volatility * np.sqrt(dt) * rand)
    return paths

GBM(1, 252, 8, 100., 0.05, 0.5)

array([[1.00000000e+02, 1.00000000e+02, 1.00000000e+02, ...,
        1.00000000e+02, 1.00000000e+02, 1.00000000e+02],
       [9.99702425e+01, 1.02116774e+02, 9.78688323e+01, ...,
        1.00978615e+02, 9.64128959e+01, 9.72154915e+01],
       [9.99404939e+01, 1.04278354e+02, 9.57830834e+01, ...,
        1.01966807e+02, 9.29544649e+01, 9.45085180e+01],
       ...,
       [9.28295879e+01, 1.88049044e+04, 4.58249200e-01, ...,
        1.14117599e+03, 1.08089096e-02, 8.58754653e-02],
       [9.28019642e+01, 1.92029616e+04, 4.48483141e-01, ...,
        1.15234371e+03, 1.04211828e-02, 8.34842557e-02],
       [9.27743486e+01, 1.96094448e+04, 4.38925214e-01, ...,
        1.16362072e+03, 1.00473641e-02, 8.11596295e-02]])

不应生成像8.11596295e-02这样的值，因此我认为代码中存在错误。如果我使用numpy库rand = np.random.standard_normal(NoOfPaths)的标准正常抽奖，那么价格与Black Scholes价格相匹配。因此我认为问题在于随机数生成器。价值8.11596295e-02是指路径中的价格，价格不太可能从100（初始价格）下降到8.11596295e-02。

参考文献：1，2，3。

Answer 1

看来sobol_seq有一个错误。 Anaconda，python 3.7,64bit，Windows 10 x64，通过pip安装sobol_seq

pip install sobol_seq

# Name                    Version                   Build  Channel
sobol-seq                 0.1.2                    pypi_0    pypi

简单的代码

print(sobol_seq.i4_sobol_generate(1, 1, 0))
print(sobol_seq.i4_sobol_generate(1, 1, 1))
print(sobol_seq.i4_sobol_generate(1, 1, 2))
print(sobol_seq.i4_sobol_generate(1, 1, 3))

产出产量

[[0.5]]
[[0.5]]
[[0.5]]
[[0.5]]

来自http://people.sc.fsu.edu/~jburkardt/py_src/sobol/sobol.html的代码，sobol_lib.py表现得很合理（好吧，除了第一点）。

好吧，封闭的代码看起来可能有效，将种子与采样数组保持在一起。慢，但......

import sobol_seq
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def i4_sobol_generate_std_normal(dim_num, seed, size=None):
    """
    Generates multivariate standard normal quasi-random variables.
    Parameters:
      Input, integer dim_num, the spatial dimension.
      Input, integer n, the number of points to generate.
      Input, integer seed, initial seed
      Output, real np array of shape (n, dim_num).
    """

    if size is None:
        q, seed = sobol_seq.i4_sobol(dim_num, seed)
        normals = norm.ppf(q)
        return (normals, seed)

    if isinstance(size, int) or isinstance(size, np.int32) or isinstance(size, np.int64) or isinstance(size, np.int16):
        rc = np.empty((dim_num, size))

        for k in range(size):
            q, seed = sobol_seq.i4_sobol(dim_num, seed)
            rc[:,k] = norm.ppf(q)

        return (rc, seed)
    else:
        raise ValueError("Size type is not recognized")

    return None


seed = 1
x, seed = i4_sobol_generate_std_normal(1, seed)
print(x)

x, seed = i4_sobol_generate_std_normal(1, seed)
print(x)

seed = 1
x, seed = i4_sobol_generate_std_normal(1, seed, size=10)
print(x)

x, seed = i4_sobol_generate_std_normal(1, seed, size=1000)
print(x)

hist, bins = np.histogram(x, bins=20, range=(-2.5, 2.5), density=True)
plt.bar(bins[:-1], hist, width = 0.22, align='edge')

plt.show()

这是图片