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

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

Any help is appreciated.

• What is a "simple" implementation? What DID you find and why is it not OK? Implementation in pseudo-code or what language?
– g g
Apr 20, 2019 at 18:32
• @gg: I edited the post to include the workings Apr 21, 2019 at 5:07
• So maybe "sobol_seq.i4_sobol_generate" is incorrect? What did you do to test correctness? And how can we see from your code that "8.11596295e-02 should not be generated"?
– g g
Apr 22, 2019 at 21:29
• It looks like you're generating 1d sobol numbers for each time step. The sobol sequence is deterministic, which means that the numbers you're using for each time step are the same, so i would expect your paths to diverged in both directions and spread out (exponentially). You need to consider each time step as a different dimension when generating the sequence, and then use a different slice for each time step.
– will
Sep 7, 2019 at 9:53
• Please provide in your answer the final row of your paths array, I think you'll find that it has both very small and very large numbers. Also, if you were to plot each path (so each column) vs time, you'll see all of the lines are smooth, do not cross each other, and diverge away from 100 (drifting at r).
– will
Sep 7, 2019 at 10:00

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

def GBM(Ttm, TradingDaysInAYear, NoOfPaths, UnderlyingPrice, RiskFreeRate, Volatility):
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()

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.
• For the Sobol sequence generation in high dimensions, a colleague has made a C and Python wrapper to the generation with makes use of the Intel MKL, which he maintains: bitbucket.org/croci/mkl_sobol/src/master Sep 11, 2019 at 13:35
• @will Thanks for your answer and details. Can you please expand on brownian bridge to fill numbers in between ones from Sobol? if all i need is uniform random low discrepancy numbers with very high dimension (say 500), how can i generate numbers from sobol sequence (lets say 10 dimensions) and then fill rest of the dimensions. Can you please list out steps and also reference any link? Is there a brownian bridge for unform random variables or only for normal random variables? Jan 30, 2021 at 18:34
• @toing Is the variance contributed from each of your random numbers the same? or are you able to asign a different amount of variance to each required variable? IF you want more, then you can look at this code by S. Joe and F.Y.Kuo, it has direction numbers for sobol numbers up to just over 21k.
– will
Jan 31, 2021 at 11:24
• @will i found code in pytorch (i have listed as answer) that can help generate sobol sequences in higher dimensionality of upto ~1k. variance contributed that i need is same. Reason i am thinking of using brownian bridge or equivalent to fill in between numbers generated from sobol is cause as per some references, sobol sequences dont do well in higher dimensions. I was therefore thinking of generating sobol of lower dimensionality and then filling them in using some alternate techniques if possible. Jan 31, 2021 at 19:06

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 = (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.

• please can you give me a reference explaining why it's necessary to randomise the points before using them? I see no need for this.
– will
Sep 7, 2019 at 16:58
• @will This is known as randomised quasi-Monte Carlo. As the sequence is deterministic and the estimate from the first half of the sequence is not independent from the second half, we can't use the CLT. Instead all we can use to bound the error is the Koksma-Hlawka inequality, which is hard to estimate. To circumvent this we use randomised estimators. Some good material on this can be found in "Randomized Quasi-Monte Carlo: An Introduction for Practitioners" (L’Ecuyer) and "Quasi-Monte Carlo Sampling" (Owen). If you want to recover a confidence interval, then randomising will achieve this. Sep 11, 2019 at 12:38

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