Note 1: Normal (Gaussian) random numbers
A sequence of random numbers that has the mean value of 0 and standard deviation of 1 when its statistics is collected.
Note 2: Longer period
Ultimately, a computer can express only a
finite state. So, random numbers would return
to the initial value in a certain cycle,
which is called the period of random numbers.
One of the characteristics of Mersenne Twister is its very long cycle of 219937-1. A simulation usually needs to secure
the random number period targeting the cube
of the number to be used. Therefore, the
longer period is an important property for
the credit VaR calculation.
Note 3: Higher order of equidistribution
A property that guarantees the randomness of each sub-set of sequences in random combination when such sub-set is made out of a large set of pseudo-random numbers. This is a very important property in higher order Monte Carlo used for the risk management calculation.
Note 4: Pseudo-random numbers
Computer-generated random numbers are not actually the true random numbers. Thus, they are called "pseudo" random numbers.
Note 5: Physical random numbers
True random numbers (likely) obtained through physical phenomena. Instead of the conventional radioactive disintegration which produces a small amount of random numbers at a given time, presently the electric circuit are employed for the physical generation of random numbers.
Note 6: Antithetic variant method
A simple improvement method of Monte Carlo's convergency. When a normal random number is generated, change its sign to produce a new random number. This makes all the odd moments (e.g., mean, skewness) turned zero with a significant increase in convergency. A normal random numbers sequence generated by antithetic variant method is self-explanatory.
Note 7: Quadratic resampling
Generate random numbers for the desired number,
calculate the statistics of entire random
number sequences, then offset the random
numbers with the obtained statistics. This
method enables to adjust the 2nd moment which
the antithetic variant cannot eliminate.
Normally, this method is used with the antithetic
variant method. In NtRand, the quadratic resampling is implemented
based on the know-how in 1990. Regarding
the recent quadratic resampling method, some
participant employs the higher order of even-order
moment matching.
Note 8: Uniform/normal random numbers conversion
Uniform random numbers distribute constantly in [0,1] interval. As random numbers are usually obtained in the form of uniform random numbers, you need to convert them to normal random numbers. This process is called uniform/normal random numbers conversion.
Note 9: Box-Muller method
The simplest method among uniform/normal random numbers conversions. It is applicable to the conversion as far as it is objected for one sequence of pseudo-random numbers.
Note 10: Inverse function method (Moro's
algorithm)
One method of the uniform/normal random numbers
conversions. Performs conversion by applying
the approximate expression obtained by Taylor's
expansion while dividing it in intervals.
In Quasi Monte Carlo, the combined use of
Box-Muller and quasi-random number does not
bring an expected improvement due to the
problem of Box-Muller method in uniform/normal
random numbers conversion. Thus, the inverse
function method, rather complicated than
the Box-Muller, must be applied.
Note 11: Uniform random numbers
Random numbers uniformly distributed in the [0,1] interval.
Note 12: Multivariate correlation random
numbers
Multiple number of normal random number sequences which are correlational to each other. These make the key in multivariate simulation such as Monte Carlo VaR.
Note 13: Two seeds of random numbers
The original random number algorithm is extended
to prepare the seed of random numbers for
64 bits. This expansion is necessary to secure
a large initial-value space for the simulation
with repeatability (i.e., back testing, WHAT-IF
analysis) in the business practices. Generally,
such simulation is conducted by having the
key fields of the database associated with
the type of random numbers.
Note 14: Array formula in Excel
Excel performs multiple-calculation using array formula and returns multiple results at a time. To input the array formula, select the desired cell for output, type the expression, and then press CTRL+SHIFT+ENTER. For further information, see the topic about "Array Formula" in your Excel Online Help.
Note 15: Singular value decomposition (SVD)
One of the triangular decomposition methods
of matrix. To generate a multivariate correlation
random numbers, the covariance matrix entered
in the generation process shall be split
into two triangular matrix. Square rooting
and Cholesky decomposition are convenient
methods and covered in elementary textbooks
for their application in real symmetric positive
definite matrix. For the larger amount of
non-correlated data, the split Cholesky decomposition
assuming a "band matrix" for narrower
condition doubles the computing efficiency
compared to the standard Cholekey method.
However, in the computing using the higher
order economic time-series data, securing
the inter-sequence independency is difficult
and the positive definite matrix assumption
is too strict. Under the higher order environment,
solutions may not be stable because of the
increased possibility to encounter the underflow
error in the process of repeat operation
in the computer. In this case, the improved
SVD is an effective and practicable method.
Note 16: Microsoft Excel-specific Limitations
With regard to the Excel function, the maximum
number of data passable at one time is limited
to 32,767 records (both in row and column
directions) for an array formula, and to
32,767 records (in row) and 256 records (in
column) for a general expression. In Windows
95/98, the 16-bit limitation in addition
to the above restrictions restrains the maximum
passable data size at one time to 64 KB.
So, the Excel's standard "Numeric"
type (double-precision floating-point number)
data cannot be transferred more than 8,192
records. Thus, the 32-bit Windows NT/2000/XP operating
system is recommended for a higher order
multivariate Monte Carlo with NtRand. These limitations are caused from the Excel
and Windows, and not by NtRand.
Note 17: Different computing result by each
CPU, operating system, and compiler
The computing result depends on the computer's
CPU, operating system, and compiler. Thus
the same program would produce various results
on the different machines. For instance,
in triangular decomposition of a larger matrix,
the results are not identical on a Windows
machine and a UNIX workstation. The key reasons
for such difference are (1) the operation
mode of the floating-point processing unit
(FPU) on the CPU chip and (2) the microprogram
physically coded on the ROM on CPU chip.
This is the different question from the binary
underflow or round-off error, or the binary
coded decimal (BCD) problems often described
in the reference books. Assume that you run
a double data type program (most accurate numeric type among C/C++) on a Windows
machine. Its binary expression must be 64-bit in most compilers (e.g.,
Sun Pro CC, egcs, HP CC, Visual C++, Intel C/C++). However, the binary expression accuracy
of an Intel x86 FPU is settable to 24/53/64-bit "using the program,"
of which setting is usually specified in the operating system or compiler.
In addition, Pentium4 / Xeon has SSE II instruction set, enables to process
two double precision value at once by vectorized operation. This feature also
have sub-effect to accuracy, of course.
Different CPU ought to employ different microprogram algorithm, therefore
it is natural for the computers to generate various results in calculation
of no analytic solution functions (e.g., trigonometric function, exponential
function). That's why you always get different results on each operating
system such as Windows, Linux, and various version of UNIX.
Watch out, if you are using a personal computer or workstation for the
position evaluation of derivatives. A computer has a daemon inside. To
be frightened, to the knowledge of experts who are engaged in numerical
calculation, there are number of problems restrained from disclosure or
active opening due to the business operation reasons of computer manufacturers.
Numerical calculation is an abyss even to the system engineers. Please
study this paper (David Goldberg, "What Every Computer Scientist Should Know about Floating-Point Arithmetic") for more information.
Note 18: How to change random numbers automatically like Excel RAND, when
recalculate worksheet
NtRand addin functions gives stable number. This is by design. But you may want
to update random number every time when you recalculate the worksheet.
There is very simple way to achieve this. See following worksheet formula.
=NtRand(10000,0,RAND()*2147483647,RAND()*2147483647)
This is what you want. NtRand() function is called just once even if it returns multiple value, because it uses array formula. In this reason, RAND() function is never called 10000 or 20000 times
by this solution.
Note 19: How to use NtRand from VBA (Visual Basic for Application)
You may want to call NtRand addin functions from VBA (Visual Basic for Application). See following example. This code fragment implements same NtRandMultiNorm worksheet example written above.
Public Sub Test1()
Dim Result As Variant
Result = Application.Run("NtRandMultiNorm",
8, 0, 12345, 67890, True, True, True, Range(Cells(14, 3), Cells(16, 5)),
Range(Cells(12, 3), Cells(12, 5)))
... now 'Result' stores 3 x 8 array of multi-dimensional
normal random number sequence ...
... do what you want ...
End Sub
Remind that VBA is different from Visual Basic. The "Application.Run..."
technique works with VBA only. VBA can convert Excel internal data type
XLOPER inside. However, Visual Basic doesn't have this feature. It is not
simple matter to call internal Excel addin functions from Visual Basic.
Note 20: Enhancement of Monte Carlo by moment matching
The most serious headache of Monte Carlo
users may be its complexity. The theoretical
performance (computing error) of the crude Monte Carlo is represented by the expression below assuming
the count of simulation is n:
Here are typical two of the approaches for
improving the performance of crude Monte
Carlo:
- Quasi-random number Monte Carlo - Use a regular numerical sequence instead
of pseudo-random numbers, to approximate
more preferable probability distribution.
- Moment matching - Operate pseudo-random numbers to be generated
with statistical method, to approximate more
preferable probability distribution.
There is little to choose whichever approach
you should take. It only depends on your
purpose to use Monte Carlo. Pros and cons
of each approach is explained below:
- Quasi-random number Monte Carlo
Works best for specific purposes, such as
to evaluate the expected value or distribution
in a single variant or few dimension. For
instance, in the option pricing with a single
variant, this approach greatly improves the
accuracy with less times of repeat operations.
The quasi-random number Monte Carlo made
a boom among major players in financial communities
in the late 1980s and early 1990s. Today
the boom is almost gone, because the moment
matching with more times of operations will
do in almost all problems thanks to the current
improvement in computer performance. The
disadvantage on the other hand, is that this
approach is not the solution for computing
percentile numbers (e.g., VaR calculation)
because of the significant performance degradation
in higher order and the behavioral problem
of quasi-random numbers on the tail of the
distribution.
Assuming that the count of simulation is
n. The theoretical performance (computing
error) of the quasi-random number Monte Carlo
relies on the dimensional number d (following expression), which means the
more the variant increases, the inferior
the performance deteriorates than the crude
Monte Carlo.
- Moment matching
Versatile for almost all purposes. There
is no other powerful alternatives than this
approach to maintain a stable performance
in the higher dimension or the tail of the
distribution (e.g., percentile number calculation).
Regardless of the great improvement in performance,
the moment matching approach also has the
same theoretical performance limit as the
crude Monte Carlo's. Therefore its performance
is poor in the single variant environment
compared to the quasi-random number Monte
Carlo. In the financial community, the moment
matching approach has become widely used
as a standard in the higher order Monte Carlo
simulation since the middle of the 1990s
when the quantification of risk [i.e. value-at-risk
(VaR)] made presence as a significant subject.
What made the moment matching so popular?
Not to mention its user-friendliness and
good predictability (perspective), the progress
in computer performance as the infrastructure
has greatly contributed to its popularity.
Most people may think that "Quasi-random
number Monte Carlo ended its role because
now we have powerful computers capable of
finishing 10,000-time or more simulation
in less than one second." This idea
seems persuasive in a sense. Then, what if
you have to calculate a special option price
for which large number of simulations would
require very heavy computing load. This case
should use the quasi-random number Monte
Carlo than the moment matching.
NtRand is capable of using two moment matching
methods of antithetic variant method (Note 6) and quadratic resampling (Note 7) either individually or simultaneously. Especially
the simultaneous operation brings a remarkably
effective result, and thus this operation
was a sort of company secret in the financial
community (e.g., among the MBS/ABS derivative
players) in the 1980s. Today, almost in the
end of 1990s, this technology is out of the
closet while many theses are released one
after another regarding this subject.
Leave precise arguments to the papers by
experts. Let us run four types of 100-count
3-variate Monte Carlo simulation (n=100)
to assess the effects of the individual operation
and simultaneous operations of antithetic variant method (Note 6) and quadratic resampling (Note 7). The comparison tables indicate the first
moment (mean), second moment (standard deviation
or variance), third moment (skewness), fourth
moment (kurtosis), and correlation coefficients
of obtained random numbers.
- Input data
Suppose to generate drifting multivariate
normal random numbers having following statistics.
It is assessed that the more proximate the
result comes to these values, the more effective
the simulation case is.
|
|
Series 1 |
Series 2 |
Series 3 |
| 1st moment (mean) |
0.01000 |
0.02000 |
0.03000 |
| 2nd moment (standard deviation) |
0.29833 |
0.24083 |
0.21213 |
| 3rd moment (skewness) |
0.00000 |
0.00000 |
0.00000 |
| 4th moment (kurtosis) |
0.00000 |
0.00000 |
0.00000 |
|
|
| Correlation matrix |
Series 1 |
Series 2 |
Series 3 |
| Series 1 |
1.00000 |
0.54843 |
0.50011 |
| Series 2 |
0.54843 |
1.00000 |
0.68817 |
| Series 3 |
0.50011 |
0.68817 |
1.00000 |
|
- Case 1: Crude Monte Carlo - with no use of
moment matching
In only a hundred times of trial computing,
the obtained statistics from the evaluated
multivariate random numbers are substantially
far from those of the input data.
|
|
Series 1 |
Series 2 |
Series 3 |
| 1st moment (mean) |
0.05210 |
0.02710 |
0.03757 |
| 2nd moment (standard deviation) |
0.29159 |
0.22906 |
0.19863 |
| 3rd moment (skewness) |
-0.06096 |
0.39699 |
0.41137 |
| 4th moment (kurtosis) |
-0.30666 |
0.17643 |
0.55905 |
|
|
| Correlation matrix |
Series 1 |
Series 2 |
Series 3 |
| Series 1 |
1.00000 |
0.47414 |
0.45856 |
| Series 2 |
0.47414 |
1.00000 |
0.71127 |
| Series 3 |
0.45856 |
0.71127 |
1.00000 |
|
- Case 2: Only with antithetic variant method (Note 6) (odd-order moment matching)
Mean and skewness which make the odd-order
moment are adjusted so that they absolutely
match the input data. Higher odd-order moments
will be zeroed (not listed below).
|
|
Series 1 |
Series 2 |
Series 3 |
| 1st moment (mean) |
0.01000 |
0.02000 |
0.03000 |
| 2nd moment (standard deviation) |
0.28964 |
0.23389 |
0.21311 |
| 3rd moment (skewness) |
0.00000 |
0.00000 |
0.00000 |
| 4th moment (kurtosis) |
-0.13148 |
-0.21078 |
-0.04597 |
|
|
| Correlation matrix |
Series 1 |
Series 2 |
Series 3 |
| Series 1 |
1.00000 |
0.58045 |
0.52405 |
| Series 2 |
0.58045 |
1.00000 |
0.74021 |
| Series 3 |
0.52405 |
0.74021 |
1.00000 |
|
- Case 3: Only with quadratic resampling (Note 7) (1st- and 2nd-moment matching)
Mean, standard deviation, and correlation
coefficient are adjusted so that they absolutely
match the input data. Note that a slight
deviation may occur in very high order restricted
by the error limit (generated in the process
of inverse matrix computing and triangular
decomposition) caused by representation accuracy
of the computer.
|
|
Series 1 |
Series 2 |
Series 3 |
| 1st moment (mean) |
0.01000 |
0.02000 |
0.03000 |
| 2nd moment (standard deviation) |
0.29833 |
0.24083 |
0.21213 |
| 3rd moment (skewness) |
0.00470 |
0.40025 |
0.57564 |
| 4th moment (kurtosis) |
0.43710 |
1.71197 |
0.72973 |
|
|
| Correlation matrix |
Series 1 |
Series 2 |
Series 3 |
| Series 1 |
1.00000 |
0.54843 |
0.50011 |
| Series 2 |
0.54843 |
1.00000 |
0.68817 |
| Series 3 |
0.50011 |
0.68817 |
1.00000 |
|
- Case 4: Simultaneous use of antithetic variant method (Note 6) and quadratic resampling (Note 7) (moment matching with odd-order moment plus
2nd-moment)
All values except for kurtosis match the
input data. Note that a slight deviation
may occur in very high order restricted by
the error limit (generated in the process
of inverse matrix computing and triangular
decomposition) caused by representation accuracy
of the computer.
|
|
Series 1 |
Series 2 |
Series 3 |
| 1st moment (mean) |
0.01000 |
0.02000 |
0.03000 |
| 2nd moment (standard deviation) |
0.29833 |
0.24083 |
0.21213 |
| 3rd moment (skewness) |
0.00000 |
0.00000 |
0.00000 |
| 4th moment (kurtosis) |
-0.73001 |
1.19494 |
0.15166 |
|
|
| Correlation matrix |
Series 1 |
Series 2 |
Series 3 |
| Series 1 |
1.00000 |
0.54843 |
0.50011 |
| Series 2 |
0.54843 |
1.00000 |
0.68817 |
| Series 3 |
0.50011 |
0.68817 |
1.00000 |
|
Following table is another result of the
case 4, where the simulation was run by 10000
times (n=10000). The kurtosis shows an improvement.
|
|
Series 1 |
Series 2 |
Series 3 |
| 1st moment (mean) |
0.01000 |
0.02000 |
0.03000 |
| 2nd moment (standard deviation) |
0.29833 |
0.24083 |
0.21213 |
| 3rd moment (skewness) |
0.00000 |
0.00000 |
0.00000 |
| 4th moment (kurtosis) |
-0.03774 |
0.11285 |
0.00159 |
|
|
| Correlation matrix |
Series 1 |
Series 2 |
Series 3 |
| Series 1 |
1.00000 |
0.54843 |
0.50011 |
| Series 2 |
0.54843 |
1.00000 |
0.68817 |
| Series 3 |
0.50011 |
0.68817 |
1.00000 |
|
As discussed through above cases, the simultaneous
use of antithetic variant method (Note 6) and quadratic resampling (Note 7) is likely to serve for drastic improvement
in convergency just with a hundred-time Monte
Carlo simulation as if it would give a complete
solution. This must be right to a certain
extent, but a wrong conclusion. People are
often lured by this sort of logical development
named "Type I error" in statistics. To interpret the this
Type I error in the case 4, a quick improvement
in convergency by moment matching is quite
different a question from the one that this
method surely converges asymptotically to
the true solution. Regardless of the purposes
for calculating option or VaR, certain times
of trial computing should be necessary.
Note 21: Choice of the random number generator algorithm
NtRand has three pseudo random number generator algorithm. The selection is made by 2nd argument (Algorithm) of each random number functions.
- Algorithm=0: Mersenne Twister 2002
- Produce 53bit precision (0, 1) real number uniform random by revised Mersenne
Twister algorithm released 1/26/2002. At this moment, this is the 'standard'
Mersenne Twister. In this implementation, its random seed is 64 bit width.
I recommend this algorithm in most of the case.
- Algorithm=1: Mersenne Twister 1998
- Mersenne Twister used within Version 1.x NtRand. Its random seed is 64
bit width. Previous Mersenne Twister algorithm had a small problem. The
most significant bit of the seed is not well reflected to the state vector.
The problem is reported in Mersenne Twister home page, and TT800 problem report. The reports says Jeff Szuhay reveals the problem against MT's small cousin
TT800. However, this is not so 'small' problem and I also detected in 1998.
So I modified MT slightly when I released NtRand in 1998. In my resolution,
project 64bit random seed into 32bit space, and avoid this reflection problem.
(I didn't report the problem like Jeff does. Sorry, Prof. Matsumoto.) From
the first release in 1998, NtRand uses this algorithm. Now, the 'standard'
Mersenne Twister algorithm has released, I won't recommend my resolution.
Use this option just in case of you want to keep backward compatibility.
- Algorithm=2: Numerical Recipies ran2
- This is the recommended random number generator algorithm in the book, W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipies in C 2nd ed., Cambridge Univ. Press, 1992. This is very famous algorithm and it has long period (> 2 x 1018). However, unfortunately, the algorithm doesn't pass many random number
quality test (e.g. 2-dimensional random walk, n-block test, etc.). It has
clear correlation problem and I don't recommend this especially for multi-dimensional
Monte Carlo simulation. There are several other algorithms in Numerical
Recipies book, however, all of them are no good.
Version 2.01 (Released 1/6/2003)
Abstract
Install NtRand
