Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_rand_bb_make_bridge_order (g05xe)

## Purpose

nag_rand_bb_make_bridge_order (g05xe) takes a set of input times and permutes them to specify one of several predefined Brownian bridge construction orders. The permuted times can be passed to nag_rand_bb_init (g05xa) or nag_rand_bb_inc_init (g05xc) to initialize the Brownian bridge generators with the chosen bridge construction order.

## Syntax

[times, ifail] = g05xe(bgord, t0, tend, intime, move, 'ntimes', ntimes, 'nmove', nmove)
[times, ifail] = nag_rand_bb_make_bridge_order(bgord, t0, tend, intime, move, 'ntimes', ntimes, 'nmove', nmove)

## Description

The Brownian bridge algorithm (see Glasserman (2004)) is a popular method for constructing a Wiener process at a set of discrete times, t0 < t1 < t2 < ,, < tN < T ${t}_{0}<{t}_{1}<{t}_{2}<,\dots ,<{t}_{N}, for N1$N\ge 1$. To ease notation we assume that T$T$ has the index N + 1$N+1$ so that T = tN + 1$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$N+2$ points of the Weiner process, Xt0,XT${X}_{{t}_{0}},{X}_{T}$ and Xti${X}_{{t}_{i}}$, for i = 1,2,,N$\mathit{i}=1,2,\dots ,N$, are generated. The value of Xt0${X}_{{t}_{0}}$ is always assumed known, and the first point to be generated is always the final time XT${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

nag_rand_bb_make_bridge_order (g05xe) accepts as input an array of time points t1 , t2 ,, tN , T ${t}_{1},{t}_{2},\dots ,{t}_{N},T$ at which the Wiener process is to be sampled. These time points are then permuted to construct the bridge. In all of the supported construction orders the first construction point is T$T$ which has index N + 1$N+1$. The remaining points are constructed by iteratively bisecting (sub-intervals of) the time indices interval [0,N + 1] $\left[0,N+1\right]$, as Figure 1 illustrates:
Figure 1
The time indices interval is processed in levels Li${L}^{\mathit{i}}$, for i = 1,2,$\mathit{i}=1,2,\dots$. Each level Li${L}^{i}$ contains ni${n}_{i}$ points L1i,,Lnii${L}_{1}^{i},\dots ,{L}_{{n}_{i}}^{i}$ where ni2i1${n}_{i}\le {2}^{i-1}$. The number of points at each level depends on the value of N$N$. The points Lji${L}_{j}^{i}$ for i1$i\ge 1$ and j = 1,2,ni$j=1,2,\dots {n}_{i}$ are computed as follows: define L00 = N + 1${L}_{0}^{0}=N+1$ and set
 Lji = J + (K − J) / 2 where J = max {Lkp : 1 ≤ k ≤ np, ​0 ≤ p < i​ and ​Lkp < Lji} ​ and ​ K = min {Lkp : 1 ≤ k ≤ np, ​0 ≤ p < i​ and ​Lkp > Lji}
$Lji = J+ (K-J)/2 where J= max{ Lkp : 1≤k≤np , ​ 0≤p Lji }$
By convention the maximum of the empty set is taken to be to be zero. Figure 1 illustrates the algorithm when N + 1$N+1$ is a power of two. When N + 1$N+1$ is not a power of two, one must decide how to round the divisions by 2$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 Lji${L}_{j}^{i}$ are supported. In both cases, levels are always processed from coarsest to finest (i.e., increasing i$i$). Within a level, the time indices can either be processed left to right (i.e., increasing j$j$) or right to left (i.e., decreasing j$j$). For example, when processing left to right, the sequence of time indices could be generated as:
 N + 1 L11 L12 L22 L13 L23 L33 L43 ⋯
$N+1 L11 L12 L22 L13 L23 L33 L43 ⋯$
while when processing right to left, the same sequence would be generated as:
 N + 1 L11 L22 L12 L43 L33 L23 L13 ⋯
$N+1 L11 L22 L12 L43 L33 L23 L13 ⋯$
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 parameter. For example, on the set of times
 t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 T
$t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 T$
the Brownian bridge would be constructed in the following orders:
bgord = 1${\mathbf{bgord}}=1$ (processing left to right, rounding down)
 T t6 t3 t9 t1 t4 t7 t11 t2 t5 t8 t10 t12
$T t6 t3 t9 t1 t4 t7 t11 t2 t5 t8 t10 t12$
bgord = 2${\mathbf{bgord}}=2$ (processing left to right, rounding up)
 T t7 t4 t10 t2 t6 t9 t12 t1 t3 t5 t8 t11
$T t7 t4 t10 t2 t6 t9 t12 t1 t3 t5 t8 t11$
bgord = 3${\mathbf{bgord}}=3$ (processing right to left, rounding down)
 T t6 t9 t3 t11 t7 t4 t1 t12 t10 t8 t5 t2
$T t6 t9 t3 t11 t7 t4 t1 t12 t10 t8 t5 t2$
bgord = 4${\mathbf{bgord}}=4$ (processing right to left, rounding up)
 T t7 t10 t4 t12 t9 t6 t2 t11 t8 t5 t3 t1
.
$T t7 t10 t4 t12 t9 t6 t2 t11 t8 t5 t3 t1 .$
The four construction methods described above can be further customized through the use of the input array move. To see the effect of this parameter, suppose that an array A$A$ holds the output of nag_rand_bb_make_bridge_order (g05xe) when nmove = 0${\mathbf{nmove}}=0$ (i.e., the unmodified bridge construction order as specified by bgord). Let
 B = {tj : j = move(i),i = 1,2, … ,nmove} $B = { tj : j=movei, i=1,2,…,nmove }$
be the array of all times identified by move, and let C$C$ be the array A$A$ with all the elements in B$B$ removed, i.e.,
 C = {A(i) : A(i) ≠ B(i),i = 1,2, … ,ntimes,j = 1,2, … ,nmove} . $C = { A(i) : A(i) ≠ B(i) , i=1,2,…,ntimes , j=1,2,…,nmove } .$
Then the output of nag_rand_bb_make_bridge_order (g05xe) when nmove > 0${\mathbf{nmove}}>0$ is given by
 B(1) B(2) ⋯ B(nmove) C(1) C(2) ⋯ C(ntimes − nmove)
$B(1) B(2) ⋯ B(nmove) C(1) C(2) ⋯ C(ntimes-nmove)$
When the Brownian bridge is used with quasi-random variates, this functionality can be used to allow specific sections of the bridge to be constructed using the lowest dimensions of the quasi-random points.

## References

Glasserman P (2004) Monte Carlo Methods in Financial Engineering Springer

## Parameters

### Compulsory Input Parameters

1:     bgord – int64int32nag_int scalar
The bridge construction order to use.
Constraint: bgord = 1${\mathbf{bgord}}=1$, 2$2$, 3$3$ or 4$4$.
2:     t0 – double scalar
t0${t}_{0}$, the start value of the time interval on which the Wiener process is to be constructed.
3:     tend – double scalar
T$T$, the largest time at which the Wiener process is to be constructed.
4:     intime(ntimes) – double array
ntimes, the dimension of the array, must satisfy the constraint ntimes1${\mathbf{ntimes}}\ge 1$.
The time points, t1,t2,,tN${t}_{1},{t}_{2},\dots ,{t}_{N}$, at which the Wiener process is to be constructed. Note that the final time T$T$ is not included in this array.
Constraints:
• t0 < intime(i)${\mathbf{t0}}<{\mathbf{intime}}\left(\mathit{i}\right)$ and intime(i) < intime(i + 1)${\mathbf{intime}}\left(\mathit{i}\right)<{\mathbf{intime}}\left(\mathit{i}+1\right)$, for i = 1,2,,ntimes1$\mathit{i}=1,2,\dots ,{\mathbf{ntimes}}-1$;
• ${\mathbf{intime}}\left({\mathbf{ntimes}}\right)<{\mathbf{tend}}$.
5:     move(nmove) – int64int32nag_int array
nmove, the dimension of the array, must satisfy the constraint 0nmoventimes$0\le {\mathbf{nmove}}\le {\mathbf{ntimes}}$.
The indices of the entries in intime which should be moved to the front of the times array, with move(j) = i${\mathbf{move}}\left(j\right)=i$ setting the j$j$th element of times to ti${t}_{i}$. Note that i$i$ ranges from 1$1$ to ntimes. When nmove = 0${\mathbf{nmove}}=0$, move is not referenced.
Constraint: 1move(j)ntimes$1\le {\mathbf{move}}\left(\mathit{j}\right)\le {\mathbf{ntimes}}$, for j = 1,2,,nmove$\mathit{j}=1,2,\dots ,{\mathbf{nmove}}$.
The elements of move must be unique.

### Optional Input Parameters

1:     ntimes – int64int32nag_int scalar
Default: The dimension of the array intime.
N$N$, the number of time points in the Weiner process, excluding t0${t}_{0}$ and T$T$.
Constraint: ntimes1${\mathbf{ntimes}}\ge 1$.
2:     nmove – int64int32nag_int scalar
Default: The dimension of the array move.
The number of elements in the array move.
Constraint: 0nmoventimes$0\le {\mathbf{nmove}}\le {\mathbf{ntimes}}$.

None.

### Output Parameters

1:     times(ntimes) – double array
The output bridge construction order. This should be passed to nag_rand_bb_init (g05xa) or nag_rand_bb_inc_init (g05xc).
2:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
Constraint: bgord = 1${\mathbf{bgord}}=1$, 2$2$, 3$3$ or 4$4$
ifail = 2${\mathbf{ifail}}=2$
Constraint: ntimes1${\mathbf{ntimes}}\ge 1$.
ifail = 3${\mathbf{ifail}}=3$
Constraint: 0nmoventimes$0\le {\mathbf{nmove}}\le {\mathbf{ntimes}}$.
ifail = 4${\mathbf{ifail}}=4$
Constraint: intime(1) > t0${\mathbf{intime}}\left(1\right)>{\mathbf{t0}}$.
Constraint: ${\mathbf{intime}}\left({\mathbf{ntimes}}\right)<{\mathbf{tend}}$.
Constraint: the elements in intime must be in increasing order.
ifail = 5${\mathbf{ifail}}=5$
Constraint: move(i)ntimes${\mathbf{move}}\left(i\right)\le {\mathbf{ntimes}}$ for all i$i$.
Constraint: move(i)1${\mathbf{move}}\left(i\right)\ge 1$ for all i$i$.
ifail = 6${\mathbf{ifail}}=6$
Constraint: all elements in move must be unique.

Not applicable.

None.

## Example

```function nag_rand_bb_make_bridge_order_example
% 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] = nag_rand_bb_make_bridge_order(t0, tend, intime, move, 'bgord', bgord);

% Initialize the Brownian bridge generator
[rcomm, ifail] = nag_rand_bb_init(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] = nag_rand_bb(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 = nag_file_print_matrix_real_gen('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];

% nag_rand_bb works with the Cholesky factorization of the covariance matrix c
% so perform the decomposition
[c, info] = nag_lapack_dpotrf('l', c);
if info ~= 0
error('Specified covariance matrix is not positive definite: info=%d', info);
end
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] = nag_rand_quasi_normal(xmean, std, npaths, iref);
z = z';
function [state] = initialize_prng(genid, subid, seed)
% Initialize the generator to a repeatable sequence
[state, ifail] = nag_rand_init_repeat(genid, subid, seed);

function [iref, state] = initialize_scrambled_qrng(genid,stype,idim,state)
iskip = int64(0);
nsdigits = int64(32);
[iref, state, ifail] = nag_rand_quasi_init_scrambled(genid, stype, ...
int64(idim), iskip, nsdigits, state);
```
```
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

```
```function g05xe_example
% 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];

% g05xb works with the Cholesky factorization of the covariance matrix c
% so perform the decomposition
[c, info] = f07fd('l', c);
if info ~= 0
error('Specified covariance matrix is not positive definite: info=%d', info);
end
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);
```
```
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

```