Python中使用准随机标准正态数的montecarlo模拟使用sobol序列给出错误值

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中有一个bug。Anaconda，python 3.7，64位，Windows 10 x64，通过pip安装sobol_seq

pip install sobol_seq

# Name                    Version                   Build  Channel
sobol-seq                 0.1.2                    pypi_0    pypi

简单代码

^{pr2}$

产出

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

来自http://people.sc.fsu.edu/~jburkardt/py_src/sobol/sobol.html的代码，sobol_库.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()

这是图片