Monte Carlo simulations in Python using quasi random standard normal numbers using sobol sequences gives erroneous values

Question

I'm trying to perform Monte Carlo Simulations using quasi-random standard normal numbers. I understand that we can use Sobol sequences to generate uniform numbers, and then use probability integral transform to convert them to standard normal numbers. My code gives unrealistic values of the simulated asset path:

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

Values like 8.11596295e-02 should not be generated, hence I think there is something wrong in the code. If I use standard normal draws from the numpy library rand = np.random.standard_normal(NoOfPaths) then the price matches with the Black Scholes price. Hence I think the problem is with the random number generator. The value 8.11596295e-02 refers to a price in a path, and it's very unlikely that the price would come down from 100 (initial price) to 8.11596295e-02.

References: 1, 2, 3.

Answer 1

For future reference, we added Sobol' sequence in SciPy 1.7: scipy.stats.qmc.Sobol

Answer 2

It appears there is a bug in sobol_seq. Anaconda, python 3.7, 64bit, Windows 10 x64, installed sobol_seq via pip

pip install sobol_seq

# Name                    Version                   Build  Channel
sobol-seq                 0.1.2                    pypi_0    pypi

Simple code

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

produced output

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

Code from http://people.sc.fsu.edu/~jburkardt/py_src/sobol/sobol.html, sobol_lib.py behaves reasonable (well, except first point).

Well, enclosed code looks like it might work, keeping seed together with sampled array. Slow, though...

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

Here is the picture