g05xc:: Random Number Generators (NAG Toolbox)

Fix two times

t_{0} < T

and let

{(t_{i})}_{1 \leq i \leq N}

be any set of time points satisfying

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

. Let

{(X_{t_{i}})}_{1 \leq i \leq N}

denote a

d

-dimensional Wiener sample path at these time points, and let

C

be any

d

d

matrix such that

C C^{T}

is the desired covariance structure for the Wiener process. Each point

X_{t_{i}}

of the sample path is constructed according to the Brownian bridge interpolation algorithm (see Glasserman (2004) or Brownian Bridge in the G05 Chapter Introduction). We always start at some fixed point

X_{t_{0}} = x \in ℝ^{d}

. If we set

X_{T} = x + C \sqrt{T - t_{0}} Z

where

Z

is any

d

-dimensional standard Normal random variable, then

X

will behave like a normal (free) Wiener process. However if we fix the terminal value

X_{T} = w \in ℝ^{d}

, then

X

will behave like a non-free Wiener process.

Implementation

Suppose we require the increments of

P

Wiener sample paths each of dimension

d

. The main input to the Brownian bridge increments generator is then an array of quasi-random points

Z^{1}, Z^{2},\dots, Z^{P}

where each point

Z^{p} = (Z_{1}^{p}, Z_{2}^{p}, \dots, Z_{D}^{p})

has dimension

D = d (N + 1)

D = d N

depending on whether a free or non-free Wiener process is required. When nag_rand_bb_inc (g05xd) is called, the

p

th sample path for

1 \leq p \leq P

is constructed as follows: if a non-free Wiener process is required set

X_{T}

equal to the terminal value

w

, otherwise construct

X_{T}

X_{T} = X_{t_{0}} + C \sqrt{T - t_{0}} [\begin{matrix} Z_{1}^{p} \\ ⋮ \\ Z_{d}^{p} \end{matrix}]

where

C

is the matrix described in Brownian Bridge Algorithm. The array times holds the remaining time points

t_{1}, t_{2},\dots t_{N}

in the order in which the bridge is to be constructed. For each

j = 1,\dots, N

set

r = times (j)

, find

q = max \{t_{0}, times (i) : 1 \leq i < j, times (i) < r\}

and

s = min \{T, times (i) : 1 \leq i < j, times (i) > r\}

and construct the point

X_{r}

X_{r} = \frac{X_{q} (s - r) + X_{s} (r - q)}{s - q} + C \sqrt{\frac{(s - r) (r - q)}{(s - q)}} [\begin{matrix} Z_{j d - a d + 1}^{p} \\ ⋮ \\ Z_{j d - a d + d}^{p} \end{matrix}]

where

a = 0

a = 1

depending on whether a free or non-free Wiener process is required. The function nag_rand_bb_make_bridge_order (g05xe) can be used to initialize the times array for several predefined bridge construction orders. Lastly, the scaled Wiener increments

\frac{X_{t_{1}} - X_{t_{0}}}{t_{1} - t_{0}}, \frac{X_{t_{2}} - X_{t_{1}}}{t_{2} - t_{1}},\dots, \frac{X_{t_{N}} - X_{t_{N - 1}}}{t_{N} - t_{N - 1}}, \frac{X_{T} - X_{t_{N}}}{T - t_{N}}

are computed.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

The efficient implementation of a Brownian bridge algorithm requires the use of a workspace array called the working stack. Since previously computed points will be used to interpolate new points, they should be kept close to the hardware processing units so that the data can be accessed quickly. Ideally the whole stack should be held in hardware cache. Different bridge construction orders may require different amounts of working stack. Indeed, a naive bridge algorithm may require a stack of size

\frac{N}{4}

or even

\frac{N}{2}

, which could be very inefficient when

N

is large. nag_rand_bb_inc_init (g05xc) performs a detailed analysis of the bridge construction order specified by times. Heuristics are used to find an execution strategy which requires a small working stack, while still constructing the bridge in the order required.

Example

The following example program calls nag_rand_bb_init (g05xa) and nag_rand_bb (g05xb) to generate two sample paths from a two-dimensional free Wiener process. It then calls nag_rand_bb_inc_init (g05xc) and nag_rand_bb_inc (g05xd) with the same input arguments to obtain the scaled increments of the Wiener sample paths. Lastly, the program prints the Wiener sample paths from nag_rand_bb (g05xb), the scaled increments from nag_rand_bb_inc (g05xd), and the cumulative sum of the unscaled increments side by side. Note that the cumulative sum of the unscaled increments is identical to the output of nag_rand_bb (g05xb).

function g05xc_example


fprintf('g05xc 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
[rcommb, ifail] = g05xa( ...
                         t0, tend, times);
[rcommd, ifail] = g05xc( ...
                         t0, tend, times);

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

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

% Call the Brownian bridge generator routine
[zb, bb, ifail] = g05xb( ...
                         npaths, start, term, z, c, rcommb, 'a', a);
[zd, bd, ifail] = g05xd( ...
                         npaths, dif, z, c, rcommd, 'a', a);

% Display the results
unscaled = zeros(1, d);
for n = 1:npaths
  fprintf('Weiner Path %d, %d time steps, %d dimensions\n', n, ntimes+1, d);
  fprintf('     Output of g05xb    Output of g05xd    Sum of g05xd\n');

  cum = start;
  unscaled = bd(1:d, n)*(intime(1)-t0);
  cum = cum + unscaled;
  fprintf('%2d %8.4f %8.4f  %8.4f %8.4f  %8.4f %8.4f\n', ...
          1, bb(1:d,n), bd(1:d,n), cum(1:d));

  jd1 = 1;
  for j = 2:ntimes
    jd1 = jd1 + d;
    unscaled = bd(jd1:j*d, n)*(intime(j)-intime(j-1));
    cum = cum + unscaled;
    fprintf('%2d %8.4f %8.4f  %8.4f %8.4f  %8.4f %8.4f\n', ...
            j, bb(jd1:j*d,n), bd(jd1:j*d,n), cum(1:d));
  end

  j = ntimes+1;
  jd1 = ntimes*d + 1;
  unscaled = bd(jd1:j*d, n)*(tend-intime(j-1));
  cum = cum + unscaled;
  fprintf('%2d %8.4f %8.4f  %8.4f %8.4f  %8.4f %8.4f\n\n', ...
          j, bb(jd1:j*d,n), bd(jd1:j*d,n), cum(1:d));
end



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

  % We want to generate the Wiener process at these time points
  intime = 1:n + t0;
  tend = t0 + n + 1;

  nmove = int64(0);
  move  = zeros(nmove, 1, 'int64');
  bgord = int64(3);

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

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

  % 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)]);

  % Generate the pseudorandom points from N(0,1)
  [state, z, ifail] = g05sk( ...
                             int64(idim*npaths), 0, 1, state);

  z = transpose(reshape(z, npaths, idim));

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

g05xc example results

Weiner Path 1, 11 time steps, 2 dimensions
     Output of g05xb    Output of g05xd    Sum of g05xd
 1  -2.2323   1.6656   -2.2323  -0.3344   -2.2323   1.6656
 2  -5.2301   1.2812   -2.9978  -0.3844   -5.2301   1.2812
 3  -0.9025  -1.2421    4.3276  -2.5234   -0.9025  -1.2421
 4  -3.6799  -0.3972   -2.7774   0.8449   -3.6799  -0.3972
 5  -6.5789  -2.0358   -2.8990  -1.6386   -6.5789  -2.0358
 6 -11.2879  -1.1972   -4.7090   0.8385  -11.2879  -1.1972
 7  -8.8959  -1.6751    2.3919  -0.4779   -8.8959  -1.6751
 8  -9.7103  -2.0523   -0.8144  -0.3772   -9.7103  -2.0523
 9  -8.5720  -3.3306    1.1383  -1.2783   -8.5720  -3.3306
10  -9.8245  -3.2035   -1.2524   0.1271   -9.8245  -3.2035
11  -4.9941  -8.3506    4.8304  -5.1471   -4.9941  -8.3506

Weiner Path 2, 11 time steps, 2 dimensions
     Output of g05xb    Output of g05xd    Sum of g05xd
 1  -1.4101   0.0576   -1.4101  -1.9424   -1.4101   0.0576
 2  -3.5738   0.2519   -2.1637   0.1943   -3.5738   0.2519
 3  -5.2528   1.7232   -1.6790   1.4713   -5.2528   1.7232
 4  -0.8540   1.0897    4.3988  -0.6335   -0.8540   1.0897
 5   0.4905  -0.9098    1.3445  -1.9995    0.4905  -0.9098
 6   2.3322   1.3415    1.8417   2.2514    2.3322   1.3415
 7   3.0105  -4.3312    0.6783  -5.6728    3.0105  -4.3312
 8   2.6776  -3.4437   -0.3329   0.8875    2.6776  -3.4437
 9   0.6546  -2.7291   -2.0230   0.7146    0.6546  -2.7291
10  -1.3175  -3.8166   -1.9721  -1.0875   -1.3175  -3.8166
11  -3.0214  -3.5439   -1.7039   0.2727   -3.0214  -3.5439

NAG Toolbox: nag_rand_bb_inc_init (g05xc)

▸▿ Contents

Purpose

Syntax

Description

Brownian Bridge Algorithm

Implementation

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example