g05xe:: Random Number Generators (NAG Toolbox)

The Brownian bridge algorithm (see Glasserman (2004)) is a popular method for constructing a Wiener process at a set of discrete times,

t_{0} < t_{1} < t_{2} <, \dots, < t_{N} < T

, for

N \geq 1

. To ease notation we assume that

T

has the index

N + 1

so that

T = t_{N + 1}

. Inherent in the algorithm is the notion of a bridge construction order which specifies the order in which the

N + 2

points of the Wiener process,

X_{t_{0}}, X_{T}

and

X_{t_{i}}

, for

i = 1, 2, \dots, N

, are generated. The value of

X_{t_{0}}

is always assumed known, and the first point to be generated is always the final time

X_{T}

. Thereafter, successive points are generated iteratively by an interpolation formula, using points which were computed at previous iterations. In many cases the bridge construction order is not important, since any construction order will yield a correct process. However, in certain cases, for example when using quasi-random variates to construct the sample paths, the bridge construction order can be important.

Supported Bridge Construction Orders

The time indices interval is processed in levels

L^{i}

, for

i = 1, 2, \dots

. Each level

L^{i}

contains

n_{i}

points

L_{1}^{i}, \dots, L_{n_{i}}^{i}

where

n_{i} \leq 2^{i - 1}

. The number of points at each level depends on the value of

N

. The points

L_{j}^{i}

for

i \geq 1

and

j = 1, 2, \dots n_{i}

are computed as follows: define

L_{0}^{0} = N + 1

and set

\begin{matrix} L_{j}^{i} = J + (K - J) / 2 & where \\ J = \max \{L_{k}^{p} : 1 \leq k \leq n_{p}, ​ 0 \leq p < i ​ and ​ L_{k}^{p} < L_{j}^{i}\} & ​ and ​ \\ K = \min \{L_{k}^{p} : 1 \leq k \leq n_{p}, ​ 0 \leq p < i ​ and ​ L_{k}^{p} > L_{j}^{i}\} \end{matrix}

By convention the maximum of the empty set is taken to be to be zero. Figure 1 illustrates the algorithm when

N + 1

is a power of two. When

N + 1

is not a power of two, one must decide how to round the divisions by

2

. For example, if one rounds down to the nearest integer, then one could get the following:

Figure 2

From the series of bisections outlined above, two ways of ordering the time indices

L_{j}^{i}

are supported. In both cases, levels are always processed from coarsest to finest (i.e., increasing

i

). Within a level, the time indices can either be processed left to right (i.e., increasing

j

) or right to left (i.e., decreasing

j

). For example, when processing left to right, the sequence of time indices could be generated as:

\begin{matrix} N + 1 & L_{1}^{1} & L_{1}^{2} & L_{2}^{2} & L_{1}^{3} & L_{2}^{3} & L_{3}^{3} & L_{4}^{3} & \dots \end{matrix}

while when processing right to left, the same sequence would be generated as:

\begin{matrix} N + 1 & L_{1}^{1} & L_{2}^{2} & L_{1}^{2} & L_{4}^{3} & L_{3}^{3} & L_{2}^{3} & L_{1}^{3} & \dots \end{matrix}

nag_rand_bb_make_bridge_order (g05xe) therefore offers four bridge construction methods; processing either left to right or right to left, with rounding either up or down. Which method is used is controlled by the bgord argument. For example, on the set of times

\begin{matrix} t_{1} & t_{2} & t_{3} & t_{4} & t_{5} & t_{6} & t_{7} & t_{8} & t_{9} & t_{10} & t_{11} & t_{12} & T \end{matrix}

the Brownian bridge would be constructed in the following orders:

bgord = 1

(processing left to right, rounding down)

\begin{matrix} T & t_{6} & t_{3} & t_{9} & t_{1} & t_{4} & t_{7} & t_{11} & t_{2} & t_{5} & t_{8} & t_{10} & t_{12} \end{matrix}

bgord = 2

(processing left to right, rounding up)

\begin{matrix} T & t_{7} & t_{4} & t_{10} & t_{2} & t_{6} & t_{9} & t_{12} & t_{1} & t_{3} & t_{5} & t_{8} & t_{11} \end{matrix}

bgord = 3

(processing right to left, rounding down)

\begin{matrix} T & t_{6} & t_{9} & t_{3} & t_{11} & t_{7} & t_{4} & t_{1} & t_{12} & t_{10} & t_{8} & t_{5} & t_{2} \end{matrix}

bgord = 4

(processing right to left, rounding up)

\begin{matrix} T & t_{7} & t_{10} & t_{4} & t_{12} & t_{9} & t_{6} & t_{2} & t_{11} & t_{8} & t_{5} & t_{3} & t_{1} \end{matrix} .

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example

function g05xe_example


fprintf('g05xe example results\n\n');

% Get information required to set up the bridge
[bgord,t0,tend,ntimes,intime,nmove,move] = get_bridge_init_data();

% Make the bridge construction bgord
[times, ifail] = g05xe( ...
                        t0, tend, intime, move, 'bgord', bgord);

% Initialize the Brownian bridge generator
[rcomm, ifail] = g05xa( ...
                        t0, tend, times);

% Get additional information required by the bridge generator
[npaths,d,start,a,term,c] = get_bridge_gen_data();

% Generate the Z values
[z] = get_z(npaths, d, a, ntimes);

% Call the Brownian bridge generator routine
[z, b, ifail] = g05xb( ...
                       npaths, start, term, z, c, rcomm, 'a', a);

% Display the results
for i = 1:npaths
  fprintf('Weiner Path %d, %d time steps, %d dimensions\n', i, ntimes+1, d);
  w = transpose(reshape(b(:,i), d, ntimes+1));

  ifail = x04ca('G', ' ', w, '');

  fprintf('\n');
end



function [bgord,t0,tend,ntimes,intime,nmove,move] = get_bridge_init_data()
  % Set the basic parameters for a Wiener process
  t0 = 0;
  ntimes = int64(10);

  % We want to generate the Wiener process at these time points
  intime = 1.71*double(1:ntimes) + t0;
  tend = t0 + 1.71*double(ntimes + 1);

  % We suppose the following 3 times are very important and should be
  % constructed first. Note: these are indices into intime
  nmove= int64(3);
  move = [int64(3), 5, 4];
  bgord = int64(3);

function [npaths,d,start,a,term,c] = get_bridge_gen_data();
  % Set the basic parameters for a non-free Wiener process
  npaths = int64(2);
  d = 3;
  a = int64(0);
  start = zeros(d, 1);
  term  = zeros(d, 1);

  % As a = 0, term need not be initialized

  % We want the following covariance matrix
  c = [ 6,   1,  -0.2; 
        1,   5,   0.3;  
       -0.2, 0.3, 4  ];

  % Cholesky factorize of the covariance matrix c
  [c, info] = f07fd('l', c);

function [z] = get_z(npaths, d, a, ntimes)
  idim = d*(ntimes+1-a);

  % We now need to generate the input pseudorandom points

  % First initialize the base pseudorandom number generator
  state = initialize_prng(int64(6), int64(0), [int64(1023401)]);

  % Scrambled quasi-random sequences preserve the good discrepancy
  % properties of quasi-random sequences while counteracting the bias
  % some applications experience when using quasi-random sequences.
  % Initialize the scrambled quasi-random generator.
  [iref, state] = initialize_scrambled_qrng(int64(1), int64(2), ...
                                            idim, state);

  % Generate the quasi-random points from N(0,1)
  xmean = zeros(idim, 1);
  std   = ones(idim, 1);
  [z, iref, ifail] = g05yj( ...
                            xmean, std, npaths, iref);
  z = z';

function [state] = initialize_prng(genid, subid, seed)
  % Initialize the generator to a repeatable sequence
  [state, ifail] = g05kf( ...
                          genid, subid, seed);


function [iref, state] = initialize_scrambled_qrng(genid,stype,idim,state)
  iskip = int64(0);
  nsdigits = int64(32);
  [iref, state, ifail] = g05yn( ...
                                genid, stype, int64(idim), ...
                                iskip, nsdigits, state);

g05xe example results

Weiner Path 1, 11 time steps, 3 dimensions
              1          2          3
  1     -2.1275    -2.4995    -6.0191
  2     -6.1589    -1.3257    -3.7378
  3     -5.1917    -3.1653    -6.2291
  4    -11.5557    -5.9183    -5.9062
  5     -9.2492    -5.7497    -4.2989
  6     -6.7853   -13.9759    -0.8990
  7    -12.7642   -15.6386    -3.6481
  8    -12.5245   -11.8142     3.3504
  9    -15.1995   -15.5145     0.5355
 10    -16.0360   -14.4140     0.0104
 11    -22.6719   -14.3308    -0.2418

Weiner Path 2, 11 time steps, 3 dimensions
              1          2          3
  1     -0.0973     3.7229     0.8640
  2      0.8027     8.5041    -0.9103
  3     -3.8494     6.1062     0.1231
  4     -6.6643     4.9936    -0.1329
  5     -6.8095     9.3508     4.7022
  6     -7.7178    10.9577    -1.4262
  7     -8.0711    12.7207     4.4744
  8    -12.8353     8.8296     7.6458
  9     -7.9795    12.2399     7.3783
 10     -6.4313    10.0770     5.5234
 11     -6.6258    10.3026     6.5021

NAG Toolbox: nag_rand_bb_make_bridge_order (g05xe)

▸▿ Contents

Purpose

Syntax

Description

Supported Bridge Construction Orders

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example