10 YEARS OF MATHEMATICA

Quasi-Random Numbers

Roman E. Mäder

MathConsult Dr. R. Mäder
Samstagernstrasse 58a
8832 Wollerau
Switzerland
mailto:maeder@mathconsult.ch
http://www.mathconsult.ch/

Daniel Egloff

Computer Services Kaisha (Schweiz) AG
Industriestrasse 50a
8304 Wallisellen
Switzerland
http://www.csk.ch/

A first look at a Mathematica package for the generation of Niederreiter sequences and other quasi-random numbers. Such sequences can be used in Monte-Carlo simulation; Their lower discepancies compared to good pseudo-random numbers lead to much improved convergence.

The package is a joint development of CSK (Schweiz) AG and MathConsult Dr. R. Mäder.

© 1998 MathConsult Dr. R. Mäder. All rights reserved.

Monte-Carlo Simulation

What is Monte-Carlo Simulation?

An Obvious Example: Expected Values

The expected value of f(X), where X is a random variable with probability density function p(x) and support [a, b], where \(TraditionalForm
      \`\(-∞\) ≤ a < b ≤ ∞\), is

\(TraditionalForm
    \`E(f) = ∫\_a\%b f\ \((x)\)\ p\ \((x)\) dx\)

The expectation is approximated by the sample mean

M = (1/n) (∑_(i = 1)^n) f((Subscript[x, i]))

where the point set (Subscript[x, i]) is generated in a suitable way to sample the distribution of X.

Pricing a European Call Option

The assumption underlying the Black-Scholes option pricing formula is that stock prices are lognormally distributed, i.e., the returns (Subscript[r, i]) = log ((Subscript[s, i])/(Subscript[s, i - 1])) are normally distributed.

The payout function of a European call option with strike price 1.

call[x_] := Max[x - 1, 0]

Interest rate and volatility

r = 0.04; σ = 0.2;

Generate 10000 lognormally distributed random numbers.

\(SeedRandom[2001]; n = 10000; 

    xs = RandomArray[LogNormalDistribution[r - σ\^2\/2, σ], n]; 
    \)

Evaluate f at these points.

fx = call/@xs;

Take the average

Plus@@(fx/Length)[fx]

Discount to present value

fr = %  (E^(-r))

Compare with Black-Scholes option pricing formula

bs = BlackScholes[1, 1, σ, r, 1]

-Log[10, Abs[fr - bs]]

How fast (slow) does it converge?

confr = 
      Table[With[{fx1 = Take[fx, i]}, 
          (E^(-r)) Plus@@(fx1/Length)[fx1]], {i, 500,  n, 
          500}];

ListPlot[confr - bs, PlotJoined -> True, AxesOrigin -> {1, 0}, 
      Ticks -> {Table[{i, 500  i}, {i, 2, Length[confr], 2}], Automatic}, 
      AxesLabel -> {"n", "error"}];

The expected error is of order (1/(Sqrt[n])).

Why use good random numbers?

What's wrong with taking equidistant values?
Numerical integration works this way.

Convert uniform values to lognormally distributed values (the inverse of the lognormal CDF).

UniformToLogNormal[μ_, σ_] := 
      Function[{y}, (E^(μ + (Sqrt[2]) σ InverseErf[-1 + 2 y]))]

Generate 1000 equidistant values between 0 and 1 (exclusive) and convert to lognormal

\(m = 1000; 

    ys = UniformToLogNormal[r - σ\^2\/2, σ]/@
        Range[1. /m, 1 - 1. /m, 1. /m]; \)

Proceed as before

fy = call/@ys;

Plus@@(fy/Length)[fy]

fe = %  (E^(-r))

-Log[10, Abs[fe - bs]]

We need something in between

Uniform distribution is good, but what about correlation in several dimensions?

Low-Discrepancy Sequences

Discrepancy

The discrepancy (Subscript[D, n]) is a measure for the deviation of a sequence of points from the uniform distribution.

A sequence of points ((Subscript[x, i])) is uniformly distributed if in the limit the fraction of points lying in any measurable set is equal to the area of that set.

The discrepancy of pseudo-random numbers is of the order O( (((((log log n)/n)))^(1/2 ))).
A lower bound for the discrepancy of any sequence in the s-dimensional hypercube is ((((Subscript[C, s])(log n))^((s - 1)/2))/n).(Subscript[C, s]) is a constant (depending on the dimension).

Low-discrepancy sequences are sequences whose discrepancy is less than (Subscript[C, s]) (((log(N))^s)/N).

No correlation between dimensions

This condition is implied by low discrepancy.

Note that the (sequential) autocorrelation of successive values is quite different from that of Brownian motion.

Effect on integration error

For integrands of bounded variation a low-discrepancy sequence gives an error estimate of O((((log(n))^s)/n)), which is significantly better than the Monte-Carlo error estimate.

Niederreiter Sequences

Niederreiter H.: Low-Discrepancy and Low-Dispersion Sequences, Journal of Number Theory 30,  (1988) p. 51-70.

Basic idea

For a prime base p and an exponent d construct rational numbers of the form

((Subscript[x, i])/(p^d)), 0 ≤ (Subscript[x, i]) < (p^d)

The sequence (Subscript[(((Subscript[x, i]))), i = 1, …, (p^d)]) is constructed as a linear recurring sequence using some power of an irreducible polynomial in (Subscript[𝔽, p]) as a generator. For random vectors, use different polynomials for each component.

We identify the p-adic expansion \(TraditionalForm
      \`\((n)\)\_p = \ a\_\(d - 1\)\ \(a\_\(d - 2\)\) …a\_0\) of an integer n in base p with at most d digits as an elements of \(TraditionalForm
      \`\((a\_0, …a\_\(d - 1\))\) ∈ 
        𝔽\^d\) such that the first component corresponds to the least significant digit.

Given the square matrix C ∈ (𝔽^((d^2))) of generating elements, the i-th element of a one-dimensional Niederreiter sequence in base p is calculated as follows: Represent i - 1 as  a(i) ∈ (𝔽^d) and define y(i)  =  C . a(i). Then reflect y(i) at the decimal point to get the i-th element \(TraditionalForm
      \`0. \ \(X\_1\) …\ X\_d ∈ \(\([0, 1\))\)\) of the Niederreiter sequence as a fraction in base p.

The Example of Bratley-Fox-Niederreiter

Choose a base , exponent, and an irreducible polyonial

\(p\  = \ 3; 
d\  = \ 4; 
irredpoly\  = \ {1, 0, 1}; \)

Construct the matrix of generating elements

\(c = GeneratingElementsImpulseResponse[irredpoly, p, d]; 

    c\  // \ TableForm\)

Convert an integer to its list of digits to the given base (minus 1).

ilist = IndexToList[8, p]

Multiply the first columns of the matrix with the digit list.

dlist = Mod[TakeColumns[c, Length[ilist]]  .  ilist, p]

Convert the result into a fraction.

NestList[#*1/p&, 1/p, Length[dlist] - 1]

dlist . %

Implementation

Note: this is work in progress.

Features

An Iterator-Style Interface

Iterators are functions the deliver a new element form a set or sequence every time they are called.

We use an object-oriented approach where independent streams of quasi-random numbers can be created (and destroyed) and when used as a function, produce as many sequence elements as desired.

Create a stream of quasi-random points with base 5 and two dimensions

s2  =  MakeQRStream[5,  2];

Called as a function with no arguments it delivers the next element, repeatedly.

s2[]

With an argument n it produces a list of n values.

s2[6]

Statistical Properties

Here we generate 2000 values and plot them.

qvals = s2[2000];

ListPlot[qvals, AspectRatio -> 1];

Compare with random numbers.

\(SeedRandom[2002]; 
rvals = Table[{Random[], Random[]}, {2000}]; \)

ListPlot[rvals, AspectRatio -> 1];

Correlations

CorrelationMatrix[qvals] // MatrixForm

CorrelationMatrix[rvals] // MatrixForm

The first four moments of both coordinates and their theoretical values

TableForm[
      Through[{Mean, Variance, Skewness, Kurtosis}[#]]&/@
        Transpose[qvals], 
      TableHeadings -> {None, {Mean, Variance, Skewness, Kurtosis}}]

TableForm[
      Through[{Mean, Variance, Skewness, Kurtosis}[#]]&/@
        Transpose[rvals], 
      TableHeadings -> {None, {Mean, Variance, Skewness, Kurtosis}}]

What it should be:

TableForm[
      Through[{Identity, N}[
          Through[{Mean, Variance, Skewness, Kurtosis}[
              UniformDistribution[0, 1]]]]], 
      TableHeadings -> {None, {Mean, Variance, Skewness, Kurtosis}}]

Implementation Details

Each instance of a QRStream has its own private instance variables, created as local variables in a Module, then held unevaluated in the object.

As much as possible is precomputed; constants are inserted literally into the instance.

SetAttributes[QRStream, HoldAll]

Options[MakeQRStream] = {
    FirstIndex->Automatic,
    ElementType->Real,
    Converter->Identity
}

MakeQRStream[b_Integer?PrimeQ, dim_Integer?Positive,
             opts___?OptionQ] :=
    Module[{genPolys, genMats, nextIndex, d, invBase, maxi,
             type, conv},
        {nextIndex,type,conv} = {FirstIndex,ElementType,Converter} /.
            Flatten[{opts}] /. Options[MakeQRStream];
        If[ nextIndex === Automatic, nextIndex = base^2 ];
        Switch[ type,
            Rational, invBase = {1/b},
            Real,    invBase = {1.0/b},
            _,       Return[$Failed]
        ];
        genPolys = GeneratingPolys[base, dim];
        d = 0; maxi=1; (* not yet initialised *)
        With[{conv=conv, genPolys=genPolys},
          QRStream[{genPolys, genMats, nextIndex, b, d, invBase, maxi},
                   conv]
        ]
    ]

The iterator is implemented as a subvalue.

QRStream[{genPolys_, genMats_, nextIndex_, b_, d_,
          invBase_, maxi_}, conv_][n_]:=
    With[{last = nextIndex+n-1},
        If[ last >= maxi, (* extend precision *)
            d = Ceiling[Log[b, last+1]];
            maxi = b^d;
            genMats = GeneratingElementsImpulseResponse[genPolys, b, d];
            Do[AppendTo[invBase, invBase[[-1]]/b], {d-Length[invBase]}];
        ];
        Table[conv[quasiRandom[genMats, b, invBase, IndexToList[nextIndex++, b]]], {n}]
    ]

Here is the internal representation of a QR stream.

s = MakeQRStream[2, 2]

At first use its variables get initialized.

s[]

To see the values of the symbols, we replace the head QRStream with List, which is not HoldAll.

List@@s

After generation of many elements, the precision has been extended.

s[100];

List@@s

Gaussian and Lognormal QR numbers

The value of the Converter option is a function that is applied to the generated quasi-random numbers.

Gauss[μ,σ] gives the inverse of the CDF of the Gaussian (normal) distribution. It converts uniformly distributed numbers into Gaussian random numbers.

Gauss[0, 1]

In this easy way, we can create a stream of quasi-random Gaussian numbers.

g1  =  MakeQRStream[3, Converter -> Gauss[0, 1]];

gvals = g1[1000];

ListPlot[gvals];

Statistical Properties

Through[{Min, Max}[gvals]]

Bin counts and bell curve

{maxbin, nbins} = {Max[Through[{Min, Max}[Abs[gvals]]]], 50};

dx = 2.0  maxbin/nbins;

bins  =  BinCounts[ gvals,  {-maxbin, maxbin, dx} ];

mid  =  
      Drop[ Range[ -maxbin, maxbin, dx ]  +  maxbin/nbins,  
        -1 ];

buckets = Transpose[ {bins,  mid} ];

bars=BarChart[ buckets,PlotRange->All,BarLabels->None,BarEdges->False];

Superimpose the bell curve of the normal distribution. The challenge is to scale it properly to fit the bar chart.

nsample = Length[gvals];

nd = NormalDistribution[Mean[gvals], (Sqrt[Variance])[gvals]];

nplot = 
      Plot[With[{
            x = (i - ((nbins + 1)/2))/((nbins - 1)/(2  maxbin))}, 
          nsample (CDF[nd, x + dx/2] - CDF[nd, x - dx/2])], {i, 1, 
          nbins}, DisplayFunction -> Identity];

Show[{bars, nplot}];

An Application

Multinormal Distribution

Construct a vector stream with a given covariance matrix.

Setup

cor0=(({{1, 0.2, 0.7}, {0, 1, -0.5}, {0, 0, 1}}));

Given a correlation matrix:

(corr = 
          cor0 + Transpose[cor0] - DiagonalMatrix[Transpose[cor0, {1, 1}]])
         // MatrixForm

And a list of standard deviations

v = DiagonalMatrix[{1, 0.5, 1.25}];

Transform

Construct the covariance matrix;

(Σ = v . corr . v) // MatrixForm

Cholesky decomposition:

(l = Chop[Transpose[CholeskyDecomposition[Σ]]]) // 
      MatrixForm

QR Stream

A Gaussian quasi-random stream:

dim = Length[Σ];

nStream = MakeQRStream[7, dim, Converter -> Gauss[0, 1]];

Convert to desired covariances:

cStream[n_] := l . #&/@nStream[n]

Run a simulation

t = cStream[1000];

Statistical Properties

Means and standard deviations

TableForm[Through[{Mean, StandardDeviation}[#]]&/@Transpose[t], 
      TableHeadings -> {Automatic, {Mean, StandardDeviation}}]

Covariance

(Σt = CovarianceMatrix[t]) // MatrixForm

RMS error:

(Sqrt[Plus@@Flatten[(((Σt - Σ))^2)]/(dim^2)
            ])

Plots

ScatterPlot3D[t, AxesLabel -> {1, 2, 3}];

Cross-sections:

Function[{xyz}, 
        ListPlot[#[[xyz]]&/@t, AspectRatio -> 1, 
          AxesLabel -> xyz]]/@{{1, 2}, {1, 3}, {2, 3}};

Outlook

Look for this and other packages to be available in purely electronic form this Fall at

http://www.mathconsult.ch/

For β testing and volume purchases/site licenses, please contact MathConsult (maeder@mathconsult.ch).

CSK (Schweiz) AG will offer this functionality in its products for risk-management.

Scratch

"FileName" /. NotebookInformation[EvaluationNotebook[]]

SetDirectory[
      Replace[%, {
          FrontEnd`FileName[dir_List, __] :> ToFileName[dir], 
          name_String :> "."}]]

Needs["Statistics`MultiDescriptiveStatistics`"]

Needs["Statistics`ContinuousDistributions`"]

Needs["Statistics`DataManipulation`"]

Needs["Graphics`Graphics3D`"];

Needs["Graphics`Graphics`"]

Needs["Graphics`MultipleListPlot`"]

Needs["LinearAlgebra`MatrixManipulation`"];

Needs["LinearAlgebra`Cholesky`"]

Needs["Niederreiter`"];

Black-Scholes

Dabei bedeuten die Parameter:

s	Preis der Aktie
x	Ausübungspreis der Option
σ	Volatilität
r	risikofreier Zinssatz
t	verbleibende Laufzeit

Module[{y}, 
      With[{nd = NormalDistribution[0,  1]}, 
        
CND[x_?NumericQ | x_DirectedInfinity] = 
          CDF[nd, x]; 
        
Derivative[n_][CND][x_?NumericQ] := 
          D[CDF[nd, y], {y, n}] /. y -> x; ]]