Sobol (Quasi Random) sequence simplified

A number of major investment banks, hedge funds and Investment research firms use Sobol as their primary Random number generator for pricing OTC products especially for low dimension problems.

There are several technical papers already available online on Sobol sequence generation algorithm. Few such links have been provided towards the end of the article. This article however, does not intend to be too technical but rather aims at providing a simplified 'Dummies' guide to implementing Sobol sequence generation.

Before getting into the heart of the algorithm, let us briefly describe the following concepts.

Bitwise XOR operator

Represented by a sign. BitXor is used to perform a logical exclusion on two expressions. For example,

Primitive polynomials over GS(2)

[GS(2) simply means that the coeffecients a_k(s) in the polynomial below can either be 0 or 1]

To generate sobol quasi random sequence, we'll have to source in primitive polynomials over GS(2) as input to the algorithm. Various researchers have already generated and published these polynomials and can be used directly (and we recommend) while implementing the algorithm. One such primitive polynomial table is published at A011260.

So what if we have these primitive polynomials already available with us, it is nevertheless good to understand what these are. At this stage you could read on for further information or skip this section if you are in a hurry.

So what is a primitive polynomial?

Consider an irreducible polynomial over GS(2) of degree s_j

x^sj + a₁x^sj-1 + a2x^sj-2+.....+ a_sj-1 x + 1

This polynomial is said to be primitive if it has the order 2^sj - 1.

[Order of a polynomial??? The order of a polynomial P(x) with P(0) <> 0 is the smallest integer e for which P(x) completely divides x^e + 1 . For example, consider 1 + x + x² completely divides x³ - 1 and hence its order is 3 (3 = 2² - 1,S_j = 2). ]

Given below is a list of a few primitive polynomials over GS(2).

s_j  Primitive polynomials

1  1 + x

2  1 + x + x²

3  1 + x + x³ , 1 + x² + x³

4  1 + x + x⁴ , 1 + x³ + x⁴

5  1 + x² + x⁵ , 1 + x + x² + x³ + x⁵ , 1 + x³ + x⁵ , 1 + x + x³ + x⁴ + x⁵ ,

   1 + x² + x³ + x⁴ + x⁵ , 1 + x + x² + x⁴ + x⁵

Direction numbers vi(s) and initialisation number mi(s).

So now we understand what primitive polynomials are, what do we need next?

Direction numbers!!!!

To generate sobol low discrepancy quasi random sequence, we will now have to assume an initial set of initialisation numbers mⁱ(s) - m₁, m₂, m₃, . . . m_sj. (Remember s_j is also the degree of primitive polynomial we discussed previously). Each assumed m_i (i^th initialisation number), however, must satisfy two criterias - it should be an odd integer and it should be less than 2ⁱ. That is, m₁ < 2, m₂ < 2², m₃ < 2³ and so on.

However, this does not mean that any integer that satisfies these two criterias would generate a good quality low discrepancy sequence. In fact, the quality of the random number depends quite a bit on the assumed set of initialisation numbers. There are a few authors who have researched and arrived at a set of initial numbers which produce quality random numbers. See [3] for example.

The values for msj+1 and beyond are determined using the following recursive equation:

m_k = 2a₁m_k-1 2²a₂m_k-2 ... 2s^j-1a_sj-1m_k-sj+1 2^sj m_k-sj m_k-sj

where a_k are the coeffecients in the primitive polynomial as represented in x_sj + a₁x^{sj-1 +} a₂x^sj-2+.....+ a^sk-1 x + 1.

Example

Suppose we have a primitive polynomial x³ + x + 1. (a₁ = 0 and a₂ = 1) so that the recursive equation for mi simplifies to the following:

m_i = 4m_i-2 8m_i-3 m_i-3,

At this stage we assume the initialisation values m₁ = 1, m₂ = 3, and m₃ = 7.

We calculate m₄, m5 .. etc using the recursive equation above.

m₄ = 12 8 1 = 5,

m₅ = 28 24 3 = 7, and so on..

Now that we know how to calculate initialisation numbers m_i(s), we shall move on to understand how direction numbers are calculated.

Direction numbers v_k are simply defined by the equation v_k = m_k/2^k

So that, if m₁ = 1, m₂ = 3 , m₃ = 7, m₄ = 5, m₅ = 7 then the direction numbers v₁ = 0.5, v₂ = 0.75, v₃ = 0.875, v₄ = 0.3125 and v₅ = 0.21875.

The Gray code

They gray code for a number n is defined by the equation:

G(n) = n[n/2]   [.] represents the integral part function

Using the above equation we get, G(1) = 1, G(2) = 3, G(3) = 2, G(4) = 6, G(5) = 7 and so on..

The Gray code algorithm

The Sobol low discrepancy sequence can now be generated using the equation:

xⁿ = g₁v₁ g₂v₂ g₃v₃ ........

where ....g₃g₂g₁ is the binary representation of the G(n)

Using this we get, x¹ = 0.5

x² = 0.5 XORoperator0.75 = 0.25 [G(2) = 3]

x³ = 0.75 [G(3) = 2]

x⁴ = 0.75 0.875 = 0.125 [G(4) = 6] .. etc and so on.

A still faster method

So now we know how to generate the sequence. But there is still a more efficient way of implementing the alorithm. This was shown by Antonov and Saleev [2]. Before that let us define one more sequece {c_i}

c_i = index of the first 0 digit from the right in the binary representation of i = (. . . i₃i₂i₁)2. We have c₀ = 1, c₁ = 2, c₂ = 1, c₃ = 3, c₄ = 1, c₅ = 2, etc.

The sequence can now be generated using the following equation :

xⁿ⁺¹ = xⁿ v_cn

As we will show, the above equation also generates the same sequence as that generated using the equation mentioned in the previous section.

The starting point of the sequence is assumed to be x⁰ = 0.

x₁ = v_c0 = v₁ = 0.5

x₂ = 0.5 v_c1 = 0.5 v₂ = 0.5 0.75 = 0.25

x₃ = 0.25 v_c2 = 0.25 v₁ = 0.25 0.5 = 0.75

x₄ = 0.75 v_c3 = 0.75 v₃ = 0.75 0.875 = 0.125

Higher dimension problems

We showed you how to generate sobol sequence for one dimension. So how do we generate this sequence for more dimensions?

The answer is simple. For each dimension use one primitive polynomial, assume a set of initialisation numbers based on the criteria above, generate rest of the initialisation numbers using the recursive equation described before, calculate direction numbers vi(s) and finally use the equation stated in the previous section to generate the sequence.

Using Sobol low discrepancy sequence for Monte Carlo simulations

A lot of Monte Carlo based pricers use sobol for spot path generations. It would thus be justified to make a final comment, before closing this article, on how to use sobol within the monte carlo closed pricer.

Lets suppose we need to simulate spot paths for 'm' underlyings with 'n' time steps. Simply generate the sobol sequences for mxn dimensions and use random numbers from a unique dimension for each time step and underlying.

Random Number Generator (xll) Excel addin

DeltaQuants' Random Number Generator library can be downloaded from here. The library has been implemented in C++. The library documentation as well as other excel libraries can be downloaded from here.

References

• Bratley, P. and Fox, B. L. (1988), "Algorithm 659: Implementing Sobol’s quasirandom sequence generator".
• Antonov, I.A. and Saleev, V.M. (1979) "An economic method of computing LPτ-sequences". Zh. Vych. Mat. Mat. Fiz. 19: 243–245 (in Russian); U.S.S.R Comput. Maths. Math. Phys. 19: 252–256 (in English).
• http://web.maths.unsw.edu.au/~fkuo/sobol/
• SOBOL’, I. M. Points which uniformly fill a multidimensional cube. Math. Cybern. 2 (1985) Znanie, Moscow (in Russian).
• S. Joe and F. Y. Kuo, Remark on Algorithm 659: Implementing Sobol's quasirandom sequence generator, ACM Trans. Math. Softw. 29, 49-57 (2003).

Related Articles

• Risk analysis of Cliquets
• Payoff catalog of some of the most commonly traded equity structured products
• Risk analysis of Shark notes
• Risk analysis of Autocallable notes
• Greeks plot for option strategies