**1) Introduction and Motivation**.

. . . . . . . . . .Download Excel file
QMC_Black_Scholes.xls for European call option with three simulations

**2) The Basic Low Discrepancy Sequences**
(with animation).

**3) The Key Role of the Uniform Distribution [0, 1]
Numbers and the Moro's Inversion **

**4) Rate of Convergence and the
High-Dimensional Case**

**5) Halton, Faure, and Sobol Sequences**(with
animation).

. . . . . . . . . .Download Excel file
random_x_quasi-random.xls

. . . . . . . . . .Download Excel file
quasi-random-mult_dim.xls

**6) A Simple Hybrid Quasi-Monte Carlo
Approach** (with animation).

. . . . . . . . . .Download Excel
file hyb_qr_generator.xls

**7) Intuition and Definition of
Discrepancy and Low Discrepancy Sequences**

**9) Monte Carlo and Quasi-Monte Carlo
Internet Links**

Quasi-Monte Carlo simulation is the traditional Monte
Carlo simulation but using ** quasi-random sequences** instead
(pseudo) random numbers. These sequences are used to generate
representative samples from the probability distributions that we are
simulating in our practical problem.

The quasi-random sequences, also called

In reality, the low discrepancy sequences are ** totally
deterministic**, so the popular name

The essential characteristic of the Monte Carlo method
is the use of random sampling techniques to reach a solution of the
physical problem.

The use of these sequences of numbers in Monte Carlo simulations can be
view in two ways. First, the generation of (quasi or pseudo) random
numbers is a way to generate representative samples (a lot of scenarios),
to describe the uncertainties of our physical problem through
probabilities distributions. We will see that the Uniform [0, 1]
distribution permits to generate all the distributions that we need to
perform the simulations.

Second, the simulation problem can be viewed as the solution of one (or
several) integral, the ** Monte Carlo Integral**. One of first
applications of Monte Carlo methods (40's) was the evaluation of
complicated integrals (not directly related to the probabilities
distributions).

Let us see the Monte Carlo integral in order to illustrate the main characteristic of the quasi-random numbers.

The evaluation of an integral using Monte Carlo simulations is the strongest application of quasi-random sequences, and hence was the initial motivation for research in this area. Quasi-random sequences are more evenly scattered throughout the region that the Monte Carlo integral is calculated, improving the accuracy of the integral evaluation. In order to see this, let the integral of one (any) function f(x) being evaluated using simulation. The idea is to use random points for the numerical evaluation of an integral, that is, using random points to determine the area under the function, see the picture below (based in Press et al., 1992).

The integral of the function f(x) is approximately the total area times
the fraction of points that fall under the curve of f(x). Naturally this
method for evaluation of an integral (using "random" points) is
competitive only for the multi-dimensional case and/or complicated
functions.

Note that the integral evaluation is better if the points are uniformly
scattered in the entire area or, for the multi-dimensional case, in the
*hypercube volume*.

The Monte Carlo simulation can be view as a problem of integral
evaluation. Recall that to calculate an expected value we have to evaluate
an integral (or a summation for discrete probability density case). Ripley
(1987, p.119) wrote: "*The object of any simulation study is the
estimation of one or more expectations of the form Ef(X)*".
In general, the Monte Carlo method solves multidimensional integrals, and
the expression for the Monte Carlo approximation for the multidimensional
integral over the *unit* hypercube is given by:

Referring the above expression, Sobol (1994, p.97) wrote: "*...our
estimate is an approximation of an integral over the n-dimensional unit
cube, derived by averaging the values of the integrand at independent
random points* G_{1}, ...,
G_{N}* uniformly distributed in
the cube*". In another topic is discussed the key
role of the uniform distribution in the [0, 1] interval generating
all the practical distributions that we need.

A finance example: the European call option solved by Monte Carlo simulation is given by the integral below:

Where **E _{t}[.]** is the payoff expectation under
risk-neutral (or martingale) measure conditional to the information at the
(current) time t;

In most cases are possible to write the integral in the interval [0, 1]
by using a function V = g(u) so that we can use f(u) instead f(V). In
other words, writing the integral with the primary uniform numbers that
the computer is generating, which are Uniform [0, 1].

The generation of the U[0, 1] numbers is the focus of this page and is
very important for the quality of the simulation.

Beyond the integration problem, there are other aspects in the
simulation that are important in several cases. For example, the *independence*
along the ** sample path** of a stochastic variable. For
problems that we need to simulate the entire path we have a

For while, let us see the different ways to generate these random points.

The ** pseudo-random** sequence of numbers
looks like random numbers because it looks unpredictable. However, pseudo
random numbers are generated with

Some risk-analysis software permit the user to choose the seed, allowing to reproduce the exact pseudo random sequence by using the same seed, which can be very useful for the modelers (see Vose, 2000, p.64).

In order to introduce the advantage of the quasi-random numbers over the
pseudo random numbers for integration problems, let us to consider the
evaluation of the European Call Option by Monte Carlo simulation in order
to compare with the very known closed-form formula of
Black-Scholes-Merton.

In the following Excel spreadsheet I use a VBA (Visual Basic for
Applications) for Excel program in order to evaluate the European call
using three types of simulation:

(a) the first simulation using the *quasi-random* numbers (the van
der Corput sequence, see topic 2) and *Moro's inversion* to get the
standard Normal distribution (see topic 3);

(b) the second simulation using the *pseudo random* numbers from
Excel with the *Excel's build in function* for the Normal inversion;
and

(c) the third simulation using the *pseudo random* numbers from
Excel but with the superior *Moro's inversion*.

**Download**
the *compressed* Excel file: **qmc_black_scholes.zip,
with 32 KB** (VBA code is not hidden).

** Alternative
download** the

By running the simulations, the reader can verify that in most cases the
simulation error with traditional Monte Carlo (pseudo random numbers) is
higher or much higher than the error with quasi-Monte Carlo simulation. In
many simulations this error is higher by a factor of 10 or more! In some
cases the error with pseudo-random sequence is lower, but this occurs by
chance in less frequent simulated cases.

Note also that for the quasi-random simulation, a higher number of
simulations always improves the accuracy of the simulation.

Even pseudo-random numbers from a reliable random
generator, which have "*the same relevant statistical properties
as a sequence of random numbers*" (Ripley, 1987, p.15), in many
situations exhibit a too slow rate of convergence for the main problem of
a Monte Carlo simulation. The main problem in finance and real options
(and in many other fields) is the evaluation of a (multidimensional)
integral of a function, e.g., the options payoff function with variable(s)
assuming values that obeys some probability density function(s). For this
type of application is possible even to sacrifice some statistical
properties (like independence), if not required in the practical problem,
in name of a faster convergence to the correct value. In this class of
problems enters the idea of low discrepancy sequences. For applications
requiring the independence property, we have many alternatives such as the
use of *hybrid *quasi-random sequences (see
this topic later).

The quasi-random numbers have the __low discrepancy property__ that
is a * measure of uniformity* for the distribution of the
points. This is very useful to represent an

The following table illustrates the nice statistical
properties of quasi-random sequence (used a *van der Corput sequence in
base 2*, see later the definition) compared with the (theoretical)
Uniform [0, 1] distribution, and with two typical pseudo-random sequences
(generate with Excel). In the table, the used QR sequence started with 0.5
instead zero, in order to see the minimum number different of zero.

Hence, quasi-random sequence presents a better
performance than typical pseudo-random sequences for all four
probabilistic moments, indicating that quasi-random sequence is more ** representative
of U[0, 1]** than pseudo-random numbers. However, the quasi-random
numbers were not developed by statisticians, but by number-theoreticians.

Instead statistics, Quasi-Monte Carlo simulation originates from a different mathematical area: the discipline of

Fortunately some technical details are not necessary to perform a good simulation, if the user owns knowledge of the main concepts and knows when and how to apply these quasi-random sequences. Some formal details are included in a more technical section below, for reference.

There are some problems when using QMC, as pointed out by Morokoff (2000):

First, quasi-Monte Carlo methods are valid for integration problems, but may not be directly applicable to simulations, due to the correlations between the points of a quasi-random sequence. This problem can be overcome in many cases....a second limitation: the improved accuracy of quasi-Monte Carlo methods is generally lost for problems of high dimension or problems in which the integrand is not smooth.

We will discuss later some ways to reduce or by-pass these problems, including the use of hybrid quasi-random sequences.

"Deterministic" Monte Carlo methods started in finance
applications with the paper of Barret & Moore & Wilmott (1992) in
*Risk*. Quasi Monte Carlo (QMC) methods began to be popular in
finance practice mainly after the work of Paskov & Traub from Columbia
University. They used an improved Sobol algorithm for evaluation of
derivatives, resulting in the software named Finder, which rose interest
of Wall Street firms. This applied research was reported at *Business
Week* in June 20, 1994, in article named "Suddenly, Number Theory
Makes Sense to Industry", on the superior results with QMC for
high-dimensional problems of derivative evaluation. In the beginning of
90's, the conventional wisdom was that, for the high-dimensional case, QMC
was inferior to Monte Carlo.

The basic nice characteristics of low discrepancy sequences are illustrated in the next topic, and later will be more formally defined.

**References:** There are many recommendable sources in the growing
literature on this subject that I used in this page. In general the
mathematically inclined reader should consider at least the theory on low
discrepancy sequences with the rigor of Niederreiter's book : *Random
Number Generation and Quasi-Monte Carlo Methods*, SIAM, CBMS 63, 1992.

For the practitioner, a good first reading on quasi-random approach is
the very didactic and comprehensive article of Galanti & Jung: "Low-Discrepancy
Sequences: Monte Carlo Simulation of Option Prices", *Journal of
Derivatives*, Fall 1997, pp.63-83.

A good complement is the paper: Boyle & Broadie & Glasserman: "Monte
Carlo Methods for Security Pricing" *Journal of Economic Dynamics
and Control*, June 1997, vol.21, n^{o} 8-9,
pp.1267-1321.

The Dupire's (eds.) book, published by Risk Books, has a collection of
many important articles in this subject.

The complete citations/references for quasi-Monte Carlo are found in
Monte Carlo bibliography page: the
Quasi-Monte Carlo topic.

The concept of low-discrepancy is associated with the property that the
successive numbers are added in a position as away as possible from the
others numbers, that is, *avoiding clustering* (groups of numbers
close to other). The sequence is constructed based in the schema that each
point is repelled from the others. So, if the idea for the points is *maximally
avoiding of each other*, the job for the numbers generated sequentially
is to fill in the larger "gaps" between the previous numbers of
the sequence.

The following picture illustrate the "filling in the gaps"
dynamic from the ** van der Corput (1935) sequence in base 2**,
which is the building block of several important low-discrepancy
sequences. This sequence starts from zero, and is confined in the interval
[0, 1) with the first 16 numbers (n from 0 to 15) given by:

We use the interval [0, 1) for the sequence definition, closed at zero
because this point is included in the sequence, but it is opened at 1
because the sequence never reach the number 1 (although for very large n
we get points as close 1 as we like). However, in several applications is
common to start from n = 1 instead n = 0 or even to discard the first n'
numbers.

Note: I think the sequence is better balanced with the (0, 1) interval
instead [0, 1) and, in addition, there are practical problems to work with
zero (for example the inversion to obtain the Normal distribution), so
that I don't include zero in my set of QR numbers for *practical
applications*.

Clearly the van der Corput sequence exhibits a cyclicality that can be
one problem in many applications: it moves upward and falls back in cycles
of two (however, a random permutation of the sequence can practically
eliminate this dependence, as we see later).

Let us see the static picture for the first 16 van der Corput (base2)
numbers. See in the picture below that there are no clustering of numbers
in this interval [0, 1).

There are van der Corput sequences for other base of prime numbers.
Different bases have different * cycle length*, the quantity
of numbers to cover the interval [0, 1) in each cycle. In base two, the
pairs (0, 1/2) and the pair (1/4, 3/4) are the two first cycles.

For other bases this cycle is larger. For example, in

**0, 1/3, 2/3, 1/9, 4/9, 7/9, 2/9, 5/9, 8/9, 1/27,
10/27, 19/27, .....**

The van der Corput sequence, for the number n and base b, is generated through by a three step procedure:

- The decimal-base number n is expanded in the base b. For example, n =
4 in base 2 is
**100**(4 =**1**x 2^{2}+**0**x 2^{1}+**0**x 2^{0}); - The number in base b is reflected. In the example,
**001**; and - Get the reflected base b number. In the example, the base 2 number
0.001 corresponds to the decimal number 1/8, that is, van der Corput =
1/8 = 0.125 (= 0 x 2
^{-1}+ 0 x 2^{-2}+ 1 x 2^{-3}).

With two equations we execute these steps. In the first equation we get
the inverted base b digits sequences **a _{j}** (e.g., for
base 2, a sequence of 0's and 1's; for base 3 a sequence of 0's, 1's and
2's, etc.), whereas in the second equation we get the corresponding van
der Corput sequence in base b for the required decimal number n:

Where m is the lowest integer that makes a_{j}(n) = 0 for all j
> m.

The reader would be interested in check the base 2 example presented
before (the numbers from 0 to 15 in base 2), in order to see that the
sequence 0, 1/2, 1/4, 3/4, etc. is reached with these two simple
equations. Even easier is to work with the algorithm code that I present
below.

The following VBA (Visual Basic for Applications) for Excel code, permits
the generation of van der Corput sequence in any base by a *user
defined function* **CorputBaseb(b, n)**, where b is the base (prime
number) and n is the number to be transformed (in base 2, for n = 0, 1,
2,... we get 0, 1/2, 1/4, ...). This code is a VBA version of the
pseudo-code presented in Wilmott (2001, pp.474-475):

There is a similar VBA code in the book of Jackson & Staunton (2001), but it is for base 2 (not any base) and I prefer to declare the integer variables "As Long" instead "As Integer" used in the book (due the "Integer" format limitation that does not permit to work with a reasonable number of simulations like 40,000 or more).

The next topic presents the explanation about the importance of the Uniform [0, 1] distribution for all simulations, and a relatively new technique to obtain the Normal distribution, which is largely used in finance and needs a special attention. The readers whom know this matter can skip this topic, jumping to item 4, but I recommend to see at least the last table of the section 3 (the accuracy problem with Excel).

The Uniform distribution in the interval [0, 1] is, for practical purposes, the only distribution that we need generate for our simulations. The reason is that the samples from the other distributions are derived using the Uniform Distribution. The Uniform distribution can be generate with either pseudo random number or quasi-random numbers, and algorithms are available to transform a Uniform distribution in any other distribution. See, for example, the excellent book of Law & Kelton, 2000, for algorithms using the Uniform distribution to generate the main distributions.

The main and direct way to do this transformation, is by the * cumulative
distribution function inversion*. In order to understand this
method let us see the picture below, a cumulative distribution of Standard
Normal Distribution, that has mean equal zero and variance equal 1, that
is, N(0, 1).

The red line is the cumulative distribution from a Standard Normal
Distribution, N(0, 1) , but could be any distribution with known
cumulative distribution that the reasoning is the same.

Note that the Y-axis has the same range of the Uniform [0, 1], because [0,
1] is the range where the probability distribution function is defined.
More important, the Y-axis space is equiprobable, so we have the __Uniform
Distribution for the range [0, 1] in the Y-axis__. This means that we
can associate each y-value drawn from U[0, 1] with any distribution with
known cumulative function, including our Standard Normal.

The blue line indicate that with a sample number drawn from the U[0, 1]
distribution (in the picture is for y ~ 0.3), we can associate a sample
number x for the distribution that we want simulate (in the picture, x ~
**-** 0.5).

The most important distribution for real options and finance applications in general is the Standard Normal Distribution. The cumulative distribution is the integral of its density function. The density distribution function for N(0, 1) looks like the picture below:

For Standard Normal distribution, the cumulative distribution (cumulative area from the picture above) is:

This equation describes the function showed in our first picture of this topic, the cumulative Normal distribution function.

In our simulation problem, we have the **Y(x)** samples (numbers
sampled from ~ U[0, 1]) but we want the corresponding x samples for any
distribution, e.g., the Normal distribution. In other words, we want **x(Y)**,
the inverse of the cumulative distribution. Hence, given one uniformly
distributed sample **u** we want to obtain the value **x(u)**,
that is, the correspondent sample for the desirable probability
distribution.

In terms of integral evaluation, let us work the simple example of *expected
value* calculus of a function g(x), which has Normal distribution with
density **f(x)** and cumulative distribution function **F(x)**.

The expected value of g(x) is given by:

However, substituting u = F(x) and letting U ~ Uniform[0, 1], we have the other expected value below.

For example, the integral for the evaluation of an European call option (C) over a stock that values V today (t = 0) with exercise price K, is given by:

Where, being F(x) the standard Normal cumulative distribution function,
was performed the substituition u = F(x), so that x = F^{-1}(u),
and u is drawn from the Uniform[0, 1] distribution.

The evaluation by simulation of the above integral with N samples u from
the Uniform distribution is performed using the equation:

Let us examine the case of Normal distribution inversion F^{-1}(u),
used in the equation above, due to its importance.

The known general way to obtain the inverse of a Normal distribution is
by using the ** Box-Muller algorithm** (see Press et al.,
1992). However, for low discrepancy sequences has been reporting that this
is not the better way because it damage the low discrepancy sequence
properties (alters the order of the sequence or scramble the sequence
uniformity), see Moro (1995), Galanti & Jung (1997), and Jackson &
Staunton (2001). In addition, has been reported that the Box-Muller
algorithm is slower than the

The traditional Normal *inversion* algorithm is given by Beasley &
Springer (1977), see Moro (1995). However, this algorithm is not very good
for the tails of the Normal distribution. The best way to obtain the
inversion from U[0, 1] to Normal distribution is by using an algorithm
presented in a famous short paper of Moro (1995). Moro presented a hybrid
algorithm: he uses the Beasley & Springer algorithm for the central
part of the Normal distribution and another algorithm for the * tails
of the distribution*.

He modeled the distribution tails using truncated *Chebyschev series*.
There is a C algorithm for this kind of series in Press et al. (1992),
eqs. 5.8.7 and 5.8.11, mainly. However, the Moro's C code presented in the
paper, is directly applied to the our problem of Normal tails.

The Moro's algorithm divides the domain for Y (where Y = uniform sample)
into two regions:

- The central region of the distribution, |Y| 0.42, is modeled as in Beasley & Springer;
- The tails of the distribution, |Y| > 0.42, are modeled with Chebyschev series.

The Moro's algorithm is very simple, the only job is to write some
constants and a few lines of code. The paper of Joy & Boyle & Tan
(1996) presents the Moro's algorithm but with some differences, mainly in
the constants (apparently, Joy et al. look for more precision with the
Chebyschev series), used also in Galanti & Jung (1997). However, the
tests that I performed with the original Moro's constants were very
satisfactory.

There is a VBA version of this algorithm in the book of Jackson &
Staunton (2001): "Advanced Modelling in Finance Using Excel and VBA",
which uses the constants from the Moro's original work. With simple short
algorithms codes in C and VBA available in the literature, I'll not repeat
it here. The interested reader can look directly in these references.

In the next two tables, I present some comparison between the Moro's
inversion and the popular Excel inversion. The Excel weakness is noted
only in the second table.

In the first following table, I present a result of a limited test (only
one quasi-random sample with 1,000 points) about the accuracy of Moro's
inversion (using the Moro's original constants) and the inversion using
the Excel internal function named **NORMSINV()**. I used the program
*Best Fit* to check the Normal distribution hypothesis after the
inversion (for both inversions, Normal distribution was the best fitting
for both chi-squared and A-D tests).

For this test, using a van der Corput quasi-random sequence in base 2, the
table shows that both inversions (Moro and Excel) presented practically
the same (acceptable) accuracy. *In this case*, the Excel 97
inversion looks not too bad with QRN as I supposed before, see the results
in the table below. Of course, are necessary new tests (including other
bases, other sequences, more simulations, etc.) to make a more conclusive
statement.

However, the Excel inversion *weak point* is in the * tails
accuracy*. Let us see the next table, which we take uniform numbers
(

On the contrary, note that the Moro's inversion get the same result of
the column "exact" for the precision of 5 decimal digits. This
nice performance of the Moro's inversion for the tails can make a
significant difference for real options valuation.

The "tail performance" is important for a large number of
simulations (so that non-zero low discrepancy numbers get closer of the
limits of the interval). For typical (real) options problems, mainly for
out-of-money options, the performance of the tails is important because
the options exercise is related to occur at the upside (more extreme
value), highlighting the weight of the tail in the (real) option value.

Moro's inversion became popular in computational finance community. More recent finance Monte Carlo's books (for advanced finance researchers) also mentioned Moro's inversion: Jackel (2002, p.10) and mainly Glasserman (2004, pp.67-69).

So, I strongly recommend the Moro's inversion instead the Excel
build-in inversion for problems that requires some accuracy. McCullough &
Wilson (2001) argue that the statistical problems with Excel 1997 were
not fixed in the new version Excel 2000.

In addition, perhaps the Moro's inversion is faster, but this test was not
performed here. **Last Time News**: a small part of the problems (like the big error in the tails) were fixed with Excel **XP**. However, most problems remains with the XP version: see the paper of McCullough (in the version I read, table 5 showed the tails problem in the Normal inversion).

The practitioner of simulation is always concerned with the
computational time and the numerical error of the simulation. Both are
related with the simulation *rate of convergence* to the true value.
The standard error in the simulation can be viewed as the standard
deviation of the (random) sample divided by an increasing function of **N**,
the number of simulations or paths.

One way to reduce the error is by looking the numerator, by using
techniques of reduction of standard deviation of the sample ("variance
reduction techniques"), which are not scope of this section.

The other way to reduce the error is by increasing the number of
simulations N. However, there is a trade-off between the computational
time and the numerical error. The natural and simplest way to reduce the
simulation numerical error is by increasing N, but the computational time
increases directly with the number of simulations N.

The rates of convergence to the true value, which are function of N, are
given below for both Monte Carlo and Quasi-Monte Carlo simulations:

- The traditional Monte Carlo (MC) using (pseudo) random numbers, has a
convergence rate of only
**O(N**(think with Central Limit Theorem arguments or see Niederreiter, 1992, Theorem 1.1, p.4). Note that this rate is independent of the dimension d, it depends only of the number of simulations N;^{-½}) - Quasi-Monte Carlo (QMC) rate of convergence can be much faster with
errors approaching size of
**O(N**in optimal cases (see e.g. Morokoff, 2000, pp.765-766, and Birge, 1995, p.2).^{-1})

One example: for MC is necessary to increase 100 times the number of simulations N to reduce the error by a factor of 10, whereas the QMC requires less (in general much less) that 100 times, and only 10 times in optimal cases; - The theoretic upper bound rate of convergence (or maximum error) for
the multi-dimensional low discrepancy sequences (Halton and others, see
next topic) is of
**O((ln N)**, see Niederreiter (p.32), where^{d}/N)**d**is the number of dimensions from the practical problem. So, unlike the pseudo-random bound convergence, the quasi-random sequence performance decreases with the dimension d.

The different criteria from items 2 and 3, which can lead to very different results for the rates of convergence of Quasi-Monte Carlo simulations, see the table below (d = number of dimensions; N = number of simulations):

Note that for dimensions higher than 5 in the table, the worst-case for
Quasi-Monte Carlo (QMC) is *much* inferior to the crude Monte Carlo
with pseudo random numbers. However the best QMC case is always better
than MC. The worst-case bound is not very reliable for practical purposes
and, in high-dimensional problems, the reader needs to be aware of the too
large required value for N that makes the QMC's worst-case better than the
crude Monte Carlo.

Let us to explain better the problem of high-dimensionality.

In finance applications is more common to get results closer of the best
rate of convergence than the theoretic worst-case. This occurs sometimes
because the use of some technique for the *reduction of effective
dimension*, sometimes because the favorable *smooth payoff function*
and even the discount effect. For instance, in the "performance table" from the Excel file QMC_Black_Scholes.xls, the absolute error function = 7 N^{ - 0.85} (fits with R2 = 0.9995), shows a convergence closer to best case O(N^{ - 1}) than the traditional MC error of O(N^{ - ½}).

Here I refer "smooth function" as, e.g., a function with second
derivative *Lipschitz-continuous*. In this case is possible to show
that the worse limit for rate of convergence decreases to **O((ln N) ^{(d-1)/2}/N^{3/2})**,
see Caflisch.& Morokoff & Owen (1997, p.306) using a hybrid of
quasi-Monte Carlo with scrambled (or randomized) nets. The effective
dimension of a function is related to its

Boyle at al (1997) wrote: "*The reason for the difference.....
may be that the integrands typically found in finance applications behave
better than those used by numerical analysts to compare different
algorithms. Another important consideration is that finance applications
typically involve discounting, and this may effectively reduce
dimensionality*".

Press et al. (1992, p.314) wrote that "*Quasi-random sampling
does better when the integrand is smooth .... than when it has step
discontinuities.*" and their Figure 7.7.2 suggests than this *smoothness
effect* is stronger for quasi-random than pseudo-random sequences.

Paskov & Traub (1996) found better results for QMC than MC for a
problem with **360 dimensions**. They evaluated a collateralized
mortgage obligation (CMO) using an ** improved Sobol sequence**.
Their Sobol algorithm is different of the one reported in Press et
al.(1992). They report that the improvements include the development of a
table of

Another good result for QMC in high-dimensional problems was found by
Caflisch.& Morokoff & Owen (1997) also for 360 dimensions, but
working with a mortgage-backed securities. They reached good results for
QMC methods using techniques like * Brownian bridge* to
reduce the effective dimensionality of the problem.

An introductory and didactic explanation of Brownian bridge with quasi-random sequences is founded at Jung (1998). However, the Brownian bridge can demand additional computation complexity/time (reordering the path) for path-dependent problems, so that it is not a panacea. Willard (1997, p.58) wrote about the Brownian bridge: "

In most real options problems, we are interested in American options and path-dependence caused by the technical uncertainty, so that we tend to prefer other techniques to avoid multi-dimensional problems.

Caflisch.& Morokoff & Owen (1997) concluded the paper suggesting
the use of ** hybrid Quasi-Monte Carlo**, for example using
quasi-random sequences for a small number of dominant dimensions and
pseudo random (with a good variance reduction technique such as Latin
Hypercube sampling) for the remaining dimensions. This point can be
exploited by the practitioner. The idea is to get the accuracy of
quasi-random approach for the dimensions with higher impact on the results
without the disadvantages of higher dimensionality behavior for these
sequences.

In a later topic, I will present a different but very simple hybrid
quasi-random algorithm, which permits to joint the nice quasi-random
vector properties, with the independence property of pseudo random numbers
to model the *sequence of QR vectors*. This approach uses
independent (pseudo) random permutations of quasi-random vectors, and
apparently avoiding all the problems of high-dimensional degradation of
low discrepancy sequences that we discuss above and that we will see
better in the next topic with some graphs.

The van der Corput sequence is the basic one dimensional low discrepancy
sequence. However, How to generate multi-dimensional low discrepancy
sequences? This is the object of this topic.

The van der Corput base2 sequence is also the **first dimension of the
Halton sequence**, the most basic low discrepancy sequence of our
interest, which can be viewed as the *building block* of other low
discrepancy sequences used in finance. **Halton** (1960), **Faure** (1982), **Sobol** (1967), and
**Niederreiter** (1987), are the best known low-discrepancy sequences,
but this is a growing research area and new sequences are being proposed
(see Boyle et al., 1997, p.1296).

The construction process of new LD sequences involves sub-dividing the
unit hypercube into sub-volumes (boxes) of constant volume, which have
faces parallel to the hypercube's faces. The idea is to put a number in
each of these sub-volumes before going to a finer grid.

Faure and Sobol sequence are like the Halton sequence, but with only one base and with a reordering of elements (we see this later). The main difference is related with the multi-dimensional construction of the sequences. The Faure and Sobol sequences change of dimension with a permutation of the quasi-random vectors, that is, reorder the basic QR vector within each dimension.

In our more basic real option problem, with one source of uncertainty,
we can think the **number of dimensions d as the number of discrete time
intervals of one sample path**, so d = T/D
t, where T is typically the expiration of our real options or the horizon
of interest. The **number of iterations N is the number of sample paths**.
So, across the paths at one specific time instant, we want a good *uniform
sample numbers* in order to generate for example a good Normal
distribution of a Brownian motion, that probably will result in a
Log-normal distribution for prices.

The main challenge for the low discrepancy sequences is to avoid the ** multi-dimensional
clustering** caused by the correlations between the dimensions. We
wish to generate LD sequences with no correlation for every pair of
dimensions (different time instants t and s).

Below are discussed the Halton, Faure, and Sobol approaches, showing the each one has a different way to build a multi-dimensional sequence.

The Halton sequence uses one different prime base for each dimension. For the first dimension uses base 2, for the second dimension Halton uses base 3, for the third dimension uses base 5, and so on. Higher base means higher cycle and higher computational time.

For the case with two dimensions, Halton sequence works with base 2 in the first dimension and base 3 in the second dimension. The two dimensional Halton sequence with 1,000 x 1,000 is displayed below.

Compare this image with the same plotting but using (pseudo) random numbers of Excel 97:

The following animation shows the "filling in the gaps" process for Halton two-dimensional plotting. See the number of iterations in blue, at the upper right area of this picture (from 50 iterations until 1,000 iterations). Note the nice property of low discrepancy, which increasing finer grid points are generated overlaying the dimensions 1 x 2.

In the previous charts including the animation above, the Halton number sequence started from n = 16 (instead the more common start of n = 1). The sequence preserves its basic properties starting from a different number, so that is not necessary to start from n = 0 or n =1. In addition, there are some advantages to cut the first n* numbers of a sequence to improve the uniformity in higher dimensions, as reported by many authors (e.g., Galanti & Jung, 1997, p.69).

**Download:** The following Excel spreadsheet permits to compare the
quasi-random with random sequences in a two-dimensional plotting. The VBA
for Excel code is available without password.

Download the
*compressed* file ** random_x_quasi-random.zip**
with 200 KB (contains the Excel-VBA file).

The ** alternative**
is to download the

The major problem for the pure quasi-random sequences is their
degradation when the dimension is large. The generation process of
uniformly distributed points in [0, 1)^{d} becomes increasingly
harder as d increases because the space to fill becomes too large.

The high-dimensional Halton sequences exhibit long cycle lengths, due the
large prime number base. For example, in the dimension 50, is used the
base 229, the 50th prime number. The long cycle length means that the
high-dimensional sequence needs several numbers for an entire walk in the
interval [0, 1).

As pointed out by Galanti & Jung (1997): "*This implies that
the speed at which increasing finer grid points are generated decreases
with increasing dimension*"

The Halton sequence for 20 x 21 dimensions is showed below and we can see
a visible degradation.

The figure below indicates that Halton sequence while still looks uniform in the 2D area, exhibits correlation between dimensions. So, Halton sequence becomes unsatisfactory after ~ dimension 14.

In practice, due the correlation, many people can prefer to avoid the use of Halton sequence for more than 6 or 8 dimensions.

**Download:** The following Excel with VBA functions and charts,
permits the reader to reproduce the results above in multi-dimension by
seeing the 2D projections of these dimensions. The VBA for Excel code is
available without password.

Download the
*compressed* file ** quasi-random-mult_dim.zip**
with 750 KB (contains the Excel-VBA file).

The alternative
is to download the *non-compressed* file version
** quasi-random-mult_dim.xls**
with 1,850 KB.

Recall, there are some problems when using QMC, as pointed out by Morokoff (2000):

First, quasi-Monte Carlo methods are valid for integration problems, but may not be directly applicable to simulations, due to the correlations between the points of a quasi-random sequence. This problem can be overcome in many cases....a second limitation: the improved accuracy of quasi-Monte Carlo methods is generally lost for problems of high dimension or problems in which the integrand is not smooth.

These problems have been solved by different means, and later will be presented a simple practical solution with hybrid-low discrepancy sequences. However, let us see some improvements over the Halton sequence with Faure and Sobol sequences.

The Faure sequence is like the Halton sequence, but with two differences. First, it uses only one base for all dimensions. Second, it uses a permutation of the vector elements for each dimension.

The base of a Faure sequence is the smallest primer number that is
larger than or equal to the number of dimensions in the problem, or equal
2 for one dimensional problem. So, the Faure sequence works with van der
Corput sequences of long cycle for high-dimensional problem. Long cycles
have the problem of higher computational time (compared with shorter cycle
sequences).

As occurred with high-dimensional Halton sequence, there is the problem of
low speed at which the Faure sequence generates increasing finer grid
points to cover the unit hypercube. However, this problem is not too
severe as the Halton sequence. For example, if the dimension of the
problem is 50, the last Halton sequence (in dimension 50) uses the 50th
prime number that is 229, whereas the Faure sequence uses the first prime
number after 50, that is a base 53, which is much smaller than 229. So,
the "filling in the gaps" in high-dimensions is faster with
Faure sequence when compared with the Halton one.

By reordering the sequence within each dimension, Faure sequences
prevents some problems of correlation for sequential high-dimensions that
occurred with the Halton sequence. Faure approach makes a link between the
low discrepancy sequences theory and the combinatorial theory for the
vector reordering.

The algorithm for the Faure sequence uses the two equations given in the
topic of van der Corput sequence (see the equations presented before the
VBA code), but before using that second equation, there are a
combinatorial rearrange of the **a _{j}**. This is performed
using a

In order to use the recursive formulae above, start the first dimension
using the van der Corput sequence with the specific Faure's base b,
reorder the numbers with the equation above for the dimension d = 2, and
so on.

Remember from the basic number theory that mod means modulo and "x
mod b" is the remainder from dividing x by b.

Galanti & Jung (1997) use a slight different equation because they
use the first dimension (the basic van der Corput sequence) to generate
all the dimensions, and the equation above uses the previous dimension (d
- 1) to generate the next dimension (d).

There are some variations for this class of sequences. Tesuka (1993)
formulated the * Generalized Faure Sequence*, based on Halton
sequence but using polynomials for reordering.

As reported by Galanti & Jung (1997, p.71), among others, the Faure
sequence suffers from the *problem of start-up* (n_{0})
and, in particular for high-dimension and low values for n_{0},
the Faure numbers exhibits clustering about zero. In order to reduce this
problem, Faure suggests discarding the first n = (b^{4} - 1)
points, where b is the base. The start-up problem is reported also for
other sequences with the same practical suggestion of discarding some
initial numbers from the sequence.

The Faure sequence exhibits high-dimensional degradation at
approximately the 25th dimension.

In relation to computational time, as reported by Galanti & Jung
(1997, p.80), Faure sequence is ** significantly slower** than
Halton and Sobol sequences.

The Sobol sequence, like the Faure sequence, has the same base for all dimensions and proceeds a reordering of the vector elements within each dimension. The Sobol sequence is simpler (and faster) than the Faure sequence in the aspect that Sobol sequence uses base 2 for all dimensions. So, there is some computational time advantage due the shorter cycle length.

However, the simplicity of Sobol sequence compared with the Faure one,
ends at this point because the reordering task is more complex. Sobol
reordering is based on a set of "direction numbers", {v_{i}}.
The v_{i} numbers are given by the equation **v _{i} =
m_{i}/2^{i}**, "

In other words, Sobol sequence uses the coefficients of irreducible
primitive polynomials of modulo 2 for its complex reordering algorithm.

There is a C code for an improved Sobol algorithm (due Antonov &
Saleev, 1979) in the "algorithm bible" Press et al (1992).

The paper of Galanti & Jung (1997) presents a patient description of Sobol algorithm, including simple numerical examples, and the interested reader should read this paper for an introduction of the algorithm, before look for more advanced sources.

While Sobol sequence "fills in the gaps" at higher rate in high-dimensional problems (due its shorter cycle), there is the problem that at high-dimensions the Sobol points tend to repeat themselves across dimensions, resulting in high-dimensional clustering. The traditional solution is again to discard the first n points. Boyle & Broadie & Glasserman (1995) discard the first 64 points, but it seems arbitrary.

Paskov & Traub from Columbia University used an improved Sobol
algorithm for evaluation of derivatives, resulting in the software named
** Finder**. See the article at

Apparently they obtained an even better Sobol algorithm version (improvement over Antonov & Saleev, 1979) with specially designed "direction numbers", reaching a better reordering for the Sobol's vectors in high-dimensional problems.

The Sobol sequence (even Antonov & Saleev version) appears to resist
more to the high-dimensional degradation. As reported by Galanti &
Jung (1997, p. 75): "*Surprisingly, the Sobol sequence does not
show any degradation through the 260th dimension*". So, the Sobol
sequence outperforms both Faure and Halton sequences.