Info

selects the Knuth-TAOCP RNG (see, Knuth, 1998), as the default uniform random number generator, and the Box-Muller algorithm as the default normal RNG. The selected uniform random number generator is then used to obtain all other non-uniform variates. Table 16.1 shows the names and the parameters of some the uniform R's various non-uniform RNG functions useful for parametric and nonparametric Monte Carlo methods. For instance, to generate a sample of size 20 from a binomial distribution with 10 independent Bernoulli trials, with probability of success fixed to 0.4 across trial, we can use the following command.

 52555335345465633124

Modern languages offer many other functions and specialized libraries, that can facilitate Monte Carlo studies. Consider the case of bootstrapping. For instance, to obtain samples of size 12, with replacement, from a vector (in this case the sequence of integers between 1 and 12), in R we can type

It is important to always carefully look for the available options in a software instruction reference before implementing a Monte Carlo study.

7. Quality of Random Number Generators

Three types of considerations should guide the practitioner in the choice of the appropriate RNG: theoretical, empirical, and practical. These three aspects should always be jointly considered when selecting a suitable random number generator for a Monte Carlo study.

7.1 Theoretical Arguments: Period Length

From a theoretical point of view one of the most critical issues concerning a random number generator is its period length. The finite period of the generators implies that there is a limit on how many random numbers can be actually generated with a particular RNG. Arguments based on the "drawing without replacement" analogy and on more refined spatial considerations, have been used to justify the need to choose a generator with a period much larger than the number of random numbers actually needed for a Monte Carlo study. The recommended minimum length of the period, p, is taken to be a function of the number of random draws, n, used in the Monte Carlo experiment. Knuth, for instance, (1998, p. 185) recommended p > 1000 • n. Other authors have made more conservative suggestions. MacLaren (1992) recommended p > n2, while Hellekalak (1998) and Ripley (1987, p. 26) recommended a much more conservative p ^ 200 • n2. For instance, if 8 million numbers are required, p should be at least 8 • 109 w 233 following Knuth's recommendation, and 12, 800 • 1012 w 254, following Ripley's suggestion. Empirical evidence supports these recommendations (see, e.g., L'Ecuyer, 1998, and references therein). It is worth highlighting that the above mentioned arguments relegate the use of a single LCG only to exploratory or small Monte Carlo projects. Any serious or intensive application would require generators with much larger periods.

Another approach, suggested by Knuth (1998), is to set the parameters of a LCG (see, Equation 6.1) to different values every time a few million numbers have been generated. For instance, in GAUSS, the functions rndcon and rnd-mult can be used to reset the constant shift parameter, c, and the multiplier, a, respectively of the uniform RNG. All parameters must be integers in the range [0, 231 - 1]. The default values are c = 0, a = 397204094, and m = 232. We can immediately see, that some choices will produce numbers that do not appear random at all. For instance, setting seed = 312479559, a = 1, and c = 0 produces a sequence with constant value seed/m = 3 124 7 9 5 59/232. In R,6

seed <- 312479559 a <- 1 c <-0 m <- 2~32

x <- ( a * seed + c ) %% m u <- x / m u returns 0.07275482 (by default R displays only 7 digits7). Other choices are much harder to assess. For instance, in GAUSS, if we change the default multiplier, from a = 1 to a = 65539, using the following set of commands:

seed = 312479559; a = 2~16 + 3; rndmult a; rndseed seed;

we obtain a well-known "bad" RNG. As Figure 16.3 clearly illustrates, though the generated uniform random sequence, when plotted in pairs in the unit square does not arouse e any suspicion, it appears extremely correlated when plotted in triplets in the unit cube. In fact, all the points lie on 15 planes in the three-dimensional space. This should highlight the fact that setting the parameter of a LCG requires considerable care.

For a theorem on how to choose a and c in order to guarantee that the maximum period is achieved and to obtain parameter values that are known to produce reasonable generators, consult Knuth (1998).

7.2 Empirical Support

Theoretical arguments alone are not enough to help in the selection of the appropriate RNG method. In fact, even if supported by sound theoretical arguments, a particular implementation of a RNG in a software application still requires testing, as they might be implemented incorrectly.

Many random number generators implemented in various software application are known to fail even simple test of randomness (see, e.g., Sawitzki, 1985; Park and Miller, 1988; and McCullough, 1999, 2000). It is advisable to perform many tests before concluding that it is safe to use a particular generator. The quality of non-uniform random variates depends closely on the quality of the uniform random numbers on which they are based. For instance, the Box-Muller method is not recommended by the literature on random number generation because of its slowness and sensitivity to the underlying uniform random number generator (see, e.g., Ripley, 1987).

6Note that R uses double-precision arithmetic and not integer arithmetic as required to implement the LCG. Even using the function as.integer the RNG cannot be correctly reproduced since 232 overflows in R whereas it is ignored in computer system implementing two's complement arithmetic.

7The version of GAUSS used for this chapter actually produces 0.072754817 with those settings.

Figure 16.3. Plots of 100,000 draws from GAUSS's LCG with a = 215 +3, c = 0, and m = 232. (a) Two-dimensional plot of pseudo-random pairs (Ui, Ui+i). (b) Three-dimensional plot of pseudo-random triplets (Ui, Ui+i, Ui+2).

(a) Two-dimensional plot of pseudo-random pairs (Ui, Ui+i).

(b) Three-dimensional plot of pseudo-random triplets (Ui, Ui+i, Ui+2).

There are several collection of tests for uniform RNG in wide use today (see, e.g., Gentle, 2003). A popular collection is the DIEHARD battery of randomness tests introduced by Marsaglia (1996). This battery of tests provides a wide range of statistical tests for evaluating the stream of output of RNGs. The DIEHARD program, provided as an MSDOS executable or as C source code, is freely available from http:fstat.fsu.edu/~geo. The DIEHARD battery of tests8 requires as input a specially formatted binary file of 10 to 11 megabytes size. The RNG should produce 32-bit positive integers that need to be saved in a text file in hexadecimal form, 8 hex 'digits' per integer, 10 integers per line, and with no intervening spaces. Consider testing a random number generator in GAUSS for Windows. The rndKMi function in GAUSS returns a matrix of random integers, between 0 and 232, and the state (seed) of the random number generator. A file containing 3 millions uniform random numbers can be created by typing the following commands in the GAUSS for Windows command prompt:

output file = rndKMi.dat reset;

screen off; print x; output off;

8The version available at the moment of writing was: DOS, Jan 7, 1997.

We can create an hex file satisfying the stated conditions using the following Perl9 script:

open( IN, "rndKMi.dat" ); open( OUT, ">rndKMi.hex" ); while ( \$line = <IN> ) { chomp \$line;

printf OUT "%08x", \$line; # prints random 32-bit integer if ( \$. % 10 == 0 ) { printf OUT "\n" } ; # new line every 10 no.

which saves the output in a file named rndKMi.hex. The MSDOS auxiliary program, asc2bin.exe, provided in the test suit, can then be used to convert the hexadecimal file into a binary file, that can be directly fed to the main DIEHARD program. The following script runs the auxiliary asc2bin utility and the main diehard programs providing them with the necessary input arguments. The results from the 18 DIEHARD tests is saved in the file rnd-KMi.out.

# converts hex file to binary file using utility program asc2bin open( INPUT, "| asc2bin" );

select(INPUT);

print "rndKMi.hex\n"; # hex input file name print "rndKMi.bin\n"; # binary output file name close(INPUT);

# performs DIEHARD test open( INPUT, "| diehard" ); select(INPUT); print "rndKMi.bin\n"; print "rndKMi.out\n"; print "111111111111111\n"; close(INPUT);

# binary input file name

# test output file name

# asks for all 18 tests

These tests are general and are not able to detect all possible problems a generator might encounter in practice. With a specific application in mind it is possible to design ad hoc tests that are more suitable to assess the usefulness of a generator for a particular application.

9Perl, an acronym for Practical Extraction and Report Language is a cross-platform, high-level, open-source programming language with network and object-oriented programming support. Perl is derived mainly from the C programming language and to a lesser extent from sed, awk, grep, and the Unix shell. Perl's process, file, and text manipulation facilities make it particularly well-suited for a wide range of tasks such as data preparation and processing. Perl can also automate what would otherwise be tedious, repetitive, and error prone activities. Perl can be downloaded from the WWW's Comprehensive Perl Archive Network (CPAN) located at http://www.cpan.org/, where the source package, modules, documentation, and links to binaries distributions of Perl are available for various platforms, including Win32, Mac, and Unix/Linux. For more on Perl and its uses in econometrics and statistics, see Baiocchi (2003) and Baiocchi (2004).

For non-uniform variates a chi-squared or a Kolmogorov-Smirnov goodness-of-fit test cab be used to test the quality of the generator. For more information of generating and testing non-uniform random variates, consult Devroye (1986), Ripley (1987), and Gentle (2003).

8. Practical Issues

Several important practical issue can play a role in the choice of a generator. We will briefly mention only two here for the sake of time and space: execution time and available software options.

Execution time considerations are not usually thought to be relevant for economic applications (see, e.g., Geweke, 1996). The assumption is that computations using the pseudo-random numbers are generally much more time consuming than generating them. However, for some applications considerable savings can be afforded by using the appropriate generation method. For instance, Train (1999) and Bhat (2000, 2001) by introducing an alternative type of random number generator in the estimation of mixed logit models, found that the accuracy of the results is greatly improved. For four and five dimension integrals, the quasi-random approach based on Halton sequences required 125 draws to achieve the same level accuracy as 2,000 draws obtained using standard pseudo-random number sequences. They found that, in terms of execution time, using quasi-random numbers can take only 10 per cent of the time required with pseudo-random numbers to deliver the same degree of accuracy. Figure 16.4 illustrates the difference between the two random number generation approaches.10 The figure clearly illustrates that quasi-random pairs provide a much better coverage in [0,1]2. For implementation in practice tables of parameters for good random number generators (see, e.g., Knuth, 1998, page 106) are needed. Portability is an issue of practical relevance. For repro-ducibility purposes a portable implementations a generator should be preferable. Portability issues are often related to hardware structure, differences in the supporting software standard, etc.

Ideally, the software and algorithms used in a Monte Carlo should be both subject to peer review and based on openly published and freely available algorithms and source code. Freely available and open-source software would also guarantee the highest degree of quality and reproducibility. Also, though most university have a portfolio of software applications useful to economists, there will always be researchers with no access to all these applications. Finally, vendors of proprietary software rarely describe the algorithms used to implement econometric and statistical procedures, nor do they provide sufficiently detailed information about their reliability. This is a serious omission

10The Halton sequences were obtained using R's rnorm.halton function provided by the fOptions library.

Figure 16.4.

APPLICATIONS OF SIMULATION METHODS 5,000 draws in two dimensions from uniform distribution.

Figure 16.4. fa) Plot of pseudo-random pairs {Ui, t^+i). (b) Plot of quasi-random pairs (Ui, Ui^i).

that makes the use of "black box" packages less attractive for academic research.

9. Testing the Random Number Generators

The results of Marsaglia's DIEHARD batteries of tests of the uniform random number generators in GAUSS for Windows version 6.0, LIMDEP version 7.0,11 and R version 2.0.0, are shown in Table 16.2. Care is required when analyzing the output of the tests. The "p-values" of the DIEHARD test are actually values of the CDF of the test statistics. For instance, consider the output line taken from a DIEHARD test output file.

chisquare for 99 degrees of freedom=100.946; p-value= .573288

The "p-value" is the value at 100.946 of the cumulative chi-squared distribution function with 3 degrees of freedom, as the following calculation in R shows.

When testing the random numbers we want to reject not only when the observed frequencies are too different from the expected ones, but also when they agree too closely.12 The standard practice of doubling the one-sided p-value to

11 The results for LIMDEP are taken from McCullough (1999, p. 199).

12Fisher (1936) suggested that the null in a goodness-of-fit test could also be rejected if the p-value is too large. Sometimes the agreement between observed and expected frequencies is "too good" to have arisen by chance only. According to Fisher (1936):

obtain a two-sided p-value is inappropriate for an asymmetric distribution. The CDF scale can be conveniently used to assess the results of the tests. A large CDF-value corresponds to a small p-value and indicates a significant result. A small CDF-value indicates a close match between observed and expected frequencies. In theory, if the CDF-value is less that 0.025 or greater than 0.975 we reject the null at the 5 per cent significance level, however, in practice, because of the large number of tests, unless stated differently, a test is considered failed if all CDF-values are either zero or one. Table 16.2 summarizes all the results of the DIEHARD tests.

LIMDEP uses L'Ecuyers (1999) multiple part random number generator. This generator has a period of roughly 2191.

GAUSS implements two generators. The functions rndi (and rndus) is based on the multiplicative (c = 0) congruential generator with period 232. GAUSS' rndKMi (and rndKMu) function is reportedly (see, Ford and Ford, 2001) based on the "KISS" and on the "Monster" uniform generators introduced by Marsaglia (2000) and has period of greater then 108888 and dimension greater than 920. The KISS generator is used to initialize the 920 seeds that the Monster generator number generator requires. The initialization of this generator requires 920 random numbers. Marsaglia's KISS generator has a period of about 2124 and itself requires 6 seeds to be initialized. In order to simplify the use of rndkmi all but one of these seeds are chosen according to Marsaglia's suggestions (Marsaglia, 2000). In practice, for the the KISS generator, only one seed requires setting.

Results in table 16.2 show that R passes all the tests13 and, given its long period, appears to be suitable for large projects. According to McCullough (1998) results, LIMDEP passes almost all tests, and because of its shorter period is useful for small projects. The results for the GAUSS rndi (and rndus) function agrees with Vinod's (2000) findings. The rndKMi (and rndKMu), which was not tested before, fails 5 tests. We agree with Vinod's conclusions that GAUSS fails too many tests to be trusted for serious Monte Carlo simulations.

These test results should be interpreted with extreme caution. The DIEHARD tests were designed with generators with a 232 period in mind. The observed inverse relationship between period length and number of failed tests suggests that further testing is required, before a more definite conclusion on the safety of a generator can be reached.

Fictitious data can seldom survive a careful scrutiny, and, since most men underestimate the frequency of large deviations arising by chance, ...

13For details on the generator used by R see footnote 2.

Table 16.2. Marsaglia's DIEHARD tests.

Test

GAUSS rndi

GAUSS rndKMi

LimDep