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

Question

I am 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.

References: https://stats.stackexchange.com/questions/27450/best-method-for-transforming-low-discrepancy-sequence-into-normal-distribution, https://stackoverflow.com/questions/9412339/recommendations-for-low-discrepancy-e-g-sobol-quasi-random-sequences-in-pytho, https://github.com/naught101/sobol_seq

Any help is appreciated.

Answer 1

It is happening because you're using the same (psuedo/quasi) random numbers for each time step.

in your code here:

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) # <-- this is the issue
        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

on this line rand = i4_sobol_generate_std_normal(1, NoOfPaths) you are settign the random numbers to be the same at each time point. The impact of this is that a path that starts with an unlikely probability will get that same unlikely probability in every time step, and so all of your paths diverge:

You need to generate all of the random numbers you need for all of your different dimensions - where each time step counts as a different dimension. I have rewritten your code, see here:

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

from matplotlib import pyplot

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(s0, r, vol, t, n_paths=8192, dt_days=3, sbol_skip=0):
    dt = dt_days / 365.25

    time_points = []
    t_ = 0
    while t_ + dt < t:
        t_ += dt
        time_points.append(t_)

    time_points.append(t)
    time_points = np.array(time_points)

    n_time_points = len(time_points)

    rand = i4_sobol_generate_std_normal(n_time_points, n_paths, skip=sbol_skip)

    paths = np.zeros((n_paths, n_time_points))
    paths[:,0] = s0
    for t_i in range(1, n_time_points):
        dt_ = time_points[t_i] - time_points[t_i-1]
        paths[:, t_i] = paths[:, t_i-1] * np.exp((r - 0.5*vol**2) * dt_ + np.sqrt(dt_) * vol * rand[:, t_i])

    return time_points, paths

n_paths = 1024
time_points, paths = GBM(100, 0, 0.2, 0.3, n_paths=n_paths)

fig = pyplot.figure()
ax = fig.add_subplot(1,1,1)

for i in range(n_paths):
    ax.plot(time_points, paths[i, :], c=(1,0,0,0.05))

ax.legend()

pyplot.show()

And i added some code to plot the paths, now they look like this:

or with more paths and some transparency added:

Now, there is an issue with this -> sobol_seq is limited to 40 dimensions. There are 2 ways to fix this:

1. Write your own sobol sequence code and include more direction numbers so that you can go to higher dimensions. This is not trivial, but by no means achievable. There is code here, credit to Stephen Joe, and Frances Kuo, complete with direction numbers to take you up to 21k dimensions.
2. You use a Brownian Bridge. Using this technique, you can re-attribute the low discrepancy nature of the sobol sequence to the points which control the most variance in the MC. For the purposes of simply generating paths, this means that you use the first dimension for the total time step, then the second point to generate the bridging point in the middle, then you recursively do the same and fill the gaps in. once you run out of sobol numbers, you have divided the space up into 40 sub time periods. You then simply fill the gaps in of these spaces with random numbers, as the gains you'll get now from low discrepancy are extremely small.

Answer 2

As a heads up, if you're using low discrepancy points, then you should be randomising these before transforming them to the normal distribution. There are many ways to do this, but for coding simplicity I have used a uniform translation (% 1 in Python). As an example this would look like

import numpy as np
from scipy.stats import norm
import sobol_seq as ss
randomisation = np.random.random()  # Ensures the sequence is randomised.
uniforms = np.concatenate(ss.i4_sobol_generate(1, 30))  # Some definitions start with 0!
uniforms = (randomisation + uniforms) % 1  # Digital shifting would also be an option
normals = norm.ppf(uniforms)  # The inverse transform method.

As for the small values, given that the log-process of a Geometric Brownian motion $S_t$ is a drifting Brownian motion with mean $\log(S_0) + (\mu - \tfrac{\sigma^2}{2})t$ and variance $\sigma^2 t$ you can compute how likely it is that the minimum of this process drops below a given level. In your case about 4-orders of magnitude in base-10 and about 9 in base-$\rm{e}$). The expression for this probability can be found here. You can evaluate this for the parameters you have and see how frequently you expect to see a value this small, and compare it to your simulations.

Answer 3

You can now use pytorch that can generate scrambled sobol numbers (low discrepancy random numbers) in a QMC.

https://pytorch.org/docs/stable/generated/torch.quasirandom.SobolEngine.html