1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ What is a "simple" implementation? What DID you find and why is it not OK? Implementation in pseudo-code or what language? $\endgroup$ – g g Apr 20 at 18:32
  • $\begingroup$ @gg: I edited the post to include the workings $\endgroup$ – SPaul Apr 21 at 5:07
  • 1
    $\begingroup$ 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"? $\endgroup$ – g g Apr 22 at 21:29
  • $\begingroup$ 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. $\endgroup$ – will Sep 7 at 9:53
  • $\begingroup$ 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). $\endgroup$ – will Sep 7 at 10:00
0
$\begingroup$

Is this simple enough in Python?

from scipy.stats import norm

uniform_rvs = generated_array_as_you_like(method='SOBOL_SEQUENCE')
normal_rvs = norm.ppf(uniform_rvs)
$\endgroup$
  • $\begingroup$ I edited the post to include the workings. Also, I didn't get what you meant by this : generated_array_as_you_like(method='SOBOL_SEQUENCE'). $\endgroup$ – SPaul Apr 21 at 5:09
  • 1
    $\begingroup$ Your workings are exactly what I have specified above. You have replaced generate_array_as_you_like(method='SOBOL_SEQUENCE') with isobol_seq.i4_sobol_generate(*args), i.e. you have generated a random uniform array with the method you like, and then generated normal random variables. If there is a problem it is not with your random number generator. The Normal distribution can contain any number from -inf to +inf $\endgroup$ – Attack68 Apr 21 at 5:48
  • $\begingroup$ 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 its very unlikely that the price would come down from 100 (initial price) to 8.11596295e-02. $\endgroup$ – SPaul Apr 21 at 7:19
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – will Sep 7 at 16:58
  • $\begingroup$ @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. $\endgroup$ – oliversm Sep 11 at 12:38
0
$\begingroup$

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:

enter image description here

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:

enter image description here

or with more paths and some transparency added: enter image description here

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.
$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – oliversm Sep 11 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.