用于蒙特卡洛模拟的 RNG 中某处出错？

Question

因此，对于 Monte-Carlo 类，我创建了一个统一的随机数生成器来模拟正态分布并将其用于 MC 期权定价，但我在某处出现了严重错误。它使用一个简单的线性同余生成器 (lcg)，它生成一个随机向量，该向量被输入到一个逆正态分布的数值近似中（beasley-springer-morrow 算法）生成标准正态分布值（具体过程详解见here）。

这是我目前的代码。

Rng：

def lcrm_generator(num, a, seed, c):

    mod = 4294967296 # 2^32
    val = int(seed)  #int() is to ensure seeds with a leading zero still get accepted by the program
    rng = []
    for i in range(num):
        val = (a * val + c) % mod
        rng.append(val/(mod-1))

    return rng

逆正态逼近器：

def bsm_algorithm(u):

    # These are my necessary initial constants
    a0 = 2.50662823884; a1 = -18.61500062529; a2 = 41.39119773534; a3 = -25.44106049637;

    b0 = -8.47351093090; b1 = 23.08336743743; b2 = -21.06224101826; b3 = 3.13082909833;

    c0 = 0.3374754822726147; c1 = 0.9761690190917186; c2 = 0.1607979714918209; c3 = 0.0276438810333863;
    c4 = 0.0038405729373609; c5 = 0.0003951896511919; c6 = 0.0000321767881768; c7 = 0.0000002888167364;
    c8 = 0.0000003960315187;

    x = [0]*len(u)
    for i in range(len(u)):
        y = u[i] - 0.5
        if abs(y) < 0.42:
            r = y**2
            x[i] = y*(((a3*r+a2)*r+a1)*r+a0)/((((b3*r+b2)*r+b1)*r+b0)*r+1)
        else:
            r = u[i]
            if y > 0:
                r = 1 - u[i]
            r = log(-log(r))
            x[i] = c0+r*(c1+r*(c2+r*(c3+r*(c4+r*(c5+r*(c6+r*(c7+r*c8)))))))
            if y < 0:
                x[i] = -x[i]

    return x

将这两者与以下结合并绘制直方图显示数据看起来正确正常，

a=lcrm_generator(100000,301,"0",21)
b = bsm_algorithm(a)
plt.hist(b, bins=100)
plt.show()

及期权定价函数：

def LookbackOP(S,K,r,sigma,intervals,sims,Call_Put=1):

    ## My objects that will determine the option prices.
    path = [0]*intervals
    values = [0]*sims

    ## Objects to hold the random nums used for simulation.
    randuni = [0]*sims
    randnorm = [0]*sims
    for i in range(sims):
        randuni[i] = lcrm_generator(intervals,301,i,21)
        randnorm[i] = bsm_algorithm(randuni[i]) 

    # Generating the simulation one by one.
    for i in range(sims):
        path[0] = 1

        ## Below is to generate one whole simulated path.

        ################## MY INCORRECT WAY ##################
        for j in range(1,intervals):
            path[j] = path[j-1]*exp((r - .5*sigma**2)*(1/intervals) + sqrt(1/intervals)*randnorm[i][j])

        ################## CORRECT BUILT-IN WAY ##################
            # path[j] = path[j-1]*exp((r - .5*sigma**2)*(1/intervals) + sqrt(1/intervals)*np.random.normal(0,1))

        ## For each separate simulation, price the option either as a call or a put.    
        if Call_Put == 1:
            values[i] = max(S*max(path)-K,0)
        elif Call_Put == 0:
            values[i] = max(K-S*min(path),0)
        else:
            print("Error: You inputted something other than '1 = Call', or '0 = Put'")
    plt.hist(values,bins=30)
    plt.show()

    ## To get expected return of option by takeing their avg.
    option_value = np.mean(values)
    print(option_value)
    return option_value

在最后一段代码中，指出了我的错误，这似乎可以通过简单地使用 numpy's 普通 rng 来修复。使用一个与另一个产生截然不同的答案，我很想相信 numpy 而不是自己，但我的代码看起来很正常，所以我哪里出错了。

Answer 1

首先，

a=lcrm_generator(100000,301,"0",21)

这看起来很奇怪 - 为什么你需要一个角色作为种子？无论如何，好的参数都在这里的表格中：https://en.wikipedia.org/wiki/Linear_congruential_generator。但我相信问题不在于 LCG，但您的高斯可能存在系统性差异。

我跑了代码

from scipy.stats import norm
q = lcrm_generator(100000, 301, "0", 21)
g = bsm(q)

param = norm.fit(g)
print(param)

对于 100 000、1 000 000 和 10 000 000 个样本，我的输出是

(0.0009370998731855792, 0.9982155665317061)
(-0.0006429485100838258, 0.9996875045073682)
(-0.0007464875819397183, 1.0002711142625116)

在 1 000 000 和 10 000 000 个样本之间没有任何改进。基本上，你对高斯反函数使用了一些近似，这些只是近似的产物，没有什么可做的。

我相信 Numpy 正在使用一种精确的正态采样方法（我认为是 Box-Muller）