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_inc_init (g05xc)

## Purpose

nag_rand_bb_inc_init (g05xc) initializes the Brownian bridge increments generator nag_rand_bb_inc (g05xd). It must be called before any calls to nag_rand_bb_inc (g05xd).

## Syntax

[rcomm, ifail] = g05xc(t0, tend, times, 'ntimes', ntimes)
[rcomm, ifail] = nag_rand_bb_inc_init(t0, tend, times, 'ntimes', ntimes)

## Description

### Brownian Bridge Algorithm

Details on the Brownian bridge algorithm and the Brownian bridge process (sometimes also called a non-free Wiener process) can be found in Section [Brownian Bridge] in the G05 Chapter Introduction. We briefly recall some notation and definitions.
Fix two times t0 < T${t}_{0} and let (ti)1iN ${\left({t}_{i}\right)}_{1\le i\le N}$ be any set of time points satisfying t0 < t1 < t2 <⋯< tN < T${t}_{0}<{t}_{1}<{t}_{2}<\cdots <{t}_{N}. Let (X ti )1iN ${\left({X}_{{t}_{i}}\right)}_{1\le i\le N}$ denote a d$d$-dimensional Wiener sample path at these time points, and let C$C$ be any d × d$d×d$ matrix such that CCT$C{C}^{\mathrm{T}}$ is the desired covariance structure for the Wiener process. Each point Xti${X}_{{t}_{i}}$ of the sample path is constructed according to the Brownian bridge interpolation algorithm (see Glasserman (2004) or Section [Brownian Bridge] in the G05 Chapter Introduction). We always start at some fixed point Xt0 = xd ${X}_{{t}_{0}}=x\in {ℝ}^{d}$. If we set XT = x + C sqrt(Tt0) Z ${X}_{T}=x+C\sqrt{T-{t}_{0}}Z$ where Z$Z$ is any d$d$-dimensional standard Normal random variable, then X$X$ will behave like a normal (free) Wiener process. However if we fix the terminal value XT = wd ${X}_{T}=w\in {ℝ}^{d}$, then X$X$ will behave like a non-free Wiener process.
The Brownian bridge increments generator uses the Brownian bridge algorithm to construct sample paths for the (free or non-free) Wiener process X$X$, and then uses this to compute the scaled Wiener increments
 ( Xt1 − Xt0 )/( t1 − t0 ) , ( Xt2 − Xt1 )/( t2 − t1 ) , … , ( XtN − XtN − 1 )/( tN − tN − 1 ) , ( XT − XtN )/( T − tN ) . $Xt1 - Xt0 t1 - t0 , Xt2 - Xt1 t2 - t1 ,…, XtN - XtN-1 tN - tN-1 , XT - XtN T - tN .$
Such increments can be useful in computing numerical solutions to stochastic differential equations driven by (free or non-free) Wiener processes.

### Implementation

Conceptually, the output of the Wiener increments generator is the same as if nag_rand_bb_init (g05xa) and nag_rand_bb (g05xb) were called first, and the scaled increments then constructed from their output. The implementation adopts a much more efficient approach whereby the scaled increments are computed directly without first constructing the Wiener sample path.
Given the start and end points of the process, the order in which successive interpolation times tj${t}_{j}$ are chosen is called the bridge construction order. The construction order is given by the array times. Further information on construction orders is given in Section [Brownian Bridge Algorithm] in the G05 Chapter Introduction. For clarity we consider here the common scenario where the Brownian bridge algorithm is used with quasi-random points. If pseudorandom numbers are used instead, these details can be ignored.
Suppose we require the increments of P$P$ Wiener sample paths each of dimension d$d$. The main input to the Brownian bridge increments generator is then an array of quasi-random points Z1,Z2 ,…, ZP${Z}^{1},{Z}^{2},\dots ,{Z}^{P}$ where each point Zp = (Z1p,Z2p,,ZDp) ${Z}^{p}=\left({Z}_{1}^{p},{Z}_{2}^{p},\dots ,{Z}_{D}^{p}\right)$ has dimension D = d(N + 1)$D=d\left(N+1\right)$ or D = dN$D=dN$ depending on whether a free or non-free Wiener process is required. When nag_rand_bb_inc (g05xd) is called, the p$p$th sample path for 1pP$1\le p\le P$ is constructed as follows: if a non-free Wiener process is required set XT${X}_{T}$ equal to the terminal value w$w$, otherwise construct XT${X}_{T}$ as
XT = Xt_0 + C sqrt(T − t0)
 [ Z1p ⋮ Zdp ]
$XT = Xt0 + C T-t0 [ Z1p ⋮ Zdp ]$
where C$C$ is the matrix described in Section [Brownian Bridge Algorithm]. The array times holds the remaining time points t1 , t2 ,… tN ${t}_{1},{t}_{2},\dots {t}_{N}$ in the order in which the bridge is to be constructed. For each j = 1 ,…, N$j=1,\dots ,N$ set r = times(j)$r={\mathbf{times}}\left(j\right)$, find
 q = max { t0, times(i) : 1 ≤ i < j , times(i) < r } $q = max { t0, timesi : 1≤i
and
 s = min {T, times(i) : 1 ≤ i < j , times(i) > r } $s = min {T, timesi : 1≤i r }$
and construct the point Xr${X}_{r}$ as
Xr = ( Xq (s − r) + Xs (r − q) )/( s − q ) + C sqrt( ( (s − r) (r − q) )/( (s − q) ) )
 [ Zjd − ad + 1p ⋮ Zjd − ad + dp ]
$Xr = Xq (s-r) + Xs (r-q) s-q + C (s-r) (r-q) (s-q) [ Zjd-ad+1p ⋮ Zjd-ad+dp ]$
where a = 0$a=0$ or a = 1$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
 ( Xt1 − Xt0 )/( t1 − t0 ) , ( Xt2 − Xt1 )/( t2 − t1 ) ,…, ( XtN − XtN − 1 )/( tN − tN − 1 ) , ( XT − XtN )/( T − tN ) $Xt1 - Xt0 t1 - t0 , Xt2 - Xt1 t2 - t1 ,…, XtN - XtN-1 tN - tN-1 , XT - XtN T - tN$
are computed.

## References

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

## Parameters

### Compulsory Input Parameters

1:     t0 – double scalar
The starting value t0${t}_{0}$ of the time interval.
2:     tend – double scalar
The end value T$T$ of the time interval.
Constraint: ${\mathbf{tend}}>{\mathbf{t0}}$.
3:     times(ntimes) – double array
ntimes, the dimension of the array, must satisfy the constraint ntimes1${\mathbf{ntimes}}\ge 1$.
The points in the time interval (t0,T)$\left({t}_{0},T\right)$ at which the Wiener process is to be constructed. The order in which points are listed in times determines the bridge construction order. The function nag_rand_bb_make_bridge_order (g05xe) can be used to create predefined bridge construction orders from a set of input times.
Constraints:
• t0 < times(i) < tend${\mathbf{t0}}<{\mathbf{times}}\left(\mathit{i}\right)<{\mathbf{tend}}$, for i = 1,2,,ntimes$\mathit{i}=1,2,\dots ,{\mathbf{ntimes}}$;
• times(i) times(j)${\mathbf{times}}\left(i\right)\ne {\mathbf{times}}\left(j\right)$, for i,j = 1,2,ntimes$i,j=1,2,\dots {\mathbf{ntimes}}$ and ij$i\ne j$.

### Optional Input Parameters

1:     ntimes – int64int32nag_int scalar
Default: The dimension of the array times.
The length of times, denoted by N$N$ in Section [Brownian Bridge Algorithm].
Constraint: ntimes1${\mathbf{ntimes}}\ge 1$.

None.

### Output Parameters

1:     rcomm(12 × (ntimes + 1)$12×\left({\mathbf{ntimes}}+1\right)$) – double array
Communication array, used to store information between calls to nag_rand_bb_inc (g05xd). This array must not be directly modified.
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: ${\mathbf{tend}}>{\mathbf{t0}}$.
ifail = 2${\mathbf{ifail}}=2$
Constraint: ntimes1${\mathbf{ntimes}}\ge 1$.
ifail = 3${\mathbf{ifail}}=3$
Constraint: t0 < times(i) < tend${\mathbf{t0}}<{\mathbf{times}}\left(i\right)<{\mathbf{tend}}$ for all i$i$.
ifail = 4${\mathbf{ifail}}=4$
Constraint: all elements of times must be unique.

## Accuracy

Not applicable.

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 N/4 $\frac{N}{4}$ or even N/2 $\frac{N}{2}$, which could be very inefficient when N$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

function nag_rand_bb_inc_init_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
[rcommb, ifail] = nag_rand_bb_init(t0, tend, times);
[rcommd, ifail] = nag_rand_bb_inc_init(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] = nag_rand_bb(npaths, start, term, z, c, rcommb, 'a', a);
[zd, bd, ifail] = nag_rand_bb_inc(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', i, ntimes+1, d);
fprintf('     Output of          Output of          Sum of \n');
fprintf('     nag_rand_bb        nag_rand_bb_inc    nag_rand_bb_inc\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));

for j = 2:ntimes
unscaled = bd(1+(j-1)*d: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(1+(j-1)*d:j*d,n), bd(1+(j-1)*d:j*d,n), cum(1:d));
end

j = ntimes+1;
unscaled = bd(1+(j-1)*d: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(1+(j-1)*d:j*d,n), bd(1+(j-1)*d: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;
ntimes = int64(10);

% We want to generate the Wiener process at these time points
intime = double(1:ntimes) + t0;
tend = t0 + double(ntimes) + 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];

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

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

z = reshape(z, idim, npaths);
function [state] = initialize_prng(genid, subid, seed)
% Initialize the generator to a repeatable sequence
[state, ifail] = nag_rand_init_repeat(genid, subid, seed)

state =

61
4826
6
0
0
816500395
0
305647340
0
460678824
0
582503647
0
529178745
0
547734508
0
0
0
1403580
0
810728
0
527612
0
0
0
1370589
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-209
0
-22853
0
4826
1
1
4826

ifail =

0

Weiner Path 0, 11 time steps, 2 dimensions
Output of          Output of          Sum of
nag_rand_bb        nag_rand_bb_inc    nag_rand_bb_inc
1  -1.3804   0.5723   -1.3804  -1.4277   -1.3804   0.5723
2  -3.4259   0.8209   -2.0456   0.2486   -3.4259   0.8209
3  -7.0455  -1.9291   -3.6196  -2.7500   -7.0455  -1.9291
4 -10.9339   0.0135   -3.8884   1.9427  -10.9339   0.0135
5  -8.4741  -0.8315    2.4598  -0.8451   -8.4741  -0.8315
6  -7.1120  -0.7627    1.3621   0.0688   -7.1120  -0.7627
7  -7.7983   0.8370   -0.6862   1.5998   -7.7983   0.8370
8  -9.9296   1.4177   -2.1314   0.5807   -9.9296   1.4177
9 -10.9003   3.7806   -0.9707   2.3629  -10.9003   3.7806
10  -8.3229   1.9173    2.5775  -1.8633   -8.3229   1.9173
11  -4.9941   0.1205    3.3288  -1.7968   -4.9941   0.1205

Weiner Path 0, 11 time steps, 2 dimensions
Output of          Output of          Sum of
nag_rand_bb        nag_rand_bb_inc    nag_rand_bb_inc
1   3.9280   2.1713    3.9280   0.1713    3.9280   2.1713
2   4.6497   1.9695    0.7217  -0.2019    4.6497   1.9695
3   9.2210   0.3321    4.5713  -1.6373    9.2210   0.3321
4  11.0763   1.3017    1.8553   0.9695   11.0763   1.3017
5  10.1419   0.8611   -0.9343  -0.4406   10.1419   0.8611
6   4.7693   4.4237   -5.3726   3.5626    4.7693   4.4237
7   4.7254   8.6772   -0.0439   4.2535    4.7254   8.6772
8   4.5599  10.2097   -0.1655   1.5325    4.5599  10.2097
9   3.5059  10.6774   -1.0540   0.4677    3.5059  10.6774
10   2.5222   9.6792   -0.9837  -0.9982    2.5222   9.6792
11   1.0553  15.1533   -1.4669   5.4741    1.0553  15.1533

function g05xc_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
[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', i, 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));

for j = 2:ntimes
unscaled = bd(1+(j-1)*d: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(1+(j-1)*d:j*d,n), bd(1+(j-1)*d:j*d,n), cum(1:d));
end

j = ntimes+1;
unscaled = bd(1+(j-1)*d: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(1+(j-1)*d:j*d,n), bd(1+(j-1)*d: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;
ntimes = int64(10);

% We want to generate the Wiener process at these time points
intime = double(1:ntimes) + t0;
tend = t0 + double(ntimes) + 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];

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

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

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

state =

61
4826
6
0
0
816500395
0
305647340
0
460678824
0
582503647
0
529178745
0
547734508
0
0
0
1403580
0
810728
0
527612
0
0
0
1370589
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-209
0
-22853
0
4826
1
1
4826

ifail =

0

Weiner Path 0, 11 time steps, 2 dimensions
Output of g05xb    Output of g05xd    Sum of g05xd
1  -1.3804   0.5723   -1.3804  -1.4277   -1.3804   0.5723
2  -3.4259   0.8209   -2.0456   0.2486   -3.4259   0.8209
3  -7.0455  -1.9291   -3.6196  -2.7500   -7.0455  -1.9291
4 -10.9339   0.0135   -3.8884   1.9427  -10.9339   0.0135
5  -8.4741  -0.8315    2.4598  -0.8451   -8.4741  -0.8315
6  -7.1120  -0.7627    1.3621   0.0688   -7.1120  -0.7627
7  -7.7983   0.8370   -0.6862   1.5998   -7.7983   0.8370
8  -9.9296   1.4177   -2.1314   0.5807   -9.9296   1.4177
9 -10.9003   3.7806   -0.9707   2.3629  -10.9003   3.7806
10  -8.3229   1.9173    2.5775  -1.8633   -8.3229   1.9173
11  -4.9941   0.1205    3.3288  -1.7968   -4.9941   0.1205

Weiner Path 0, 11 time steps, 2 dimensions
Output of g05xb    Output of g05xd    Sum of g05xd
1   3.9280   2.1713    3.9280   0.1713    3.9280   2.1713
2   4.6497   1.9695    0.7217  -0.2019    4.6497   1.9695
3   9.2210   0.3321    4.5713  -1.6373    9.2210   0.3321
4  11.0763   1.3017    1.8553   0.9695   11.0763   1.3017
5  10.1419   0.8611   -0.9343  -0.4406   10.1419   0.8611
6   4.7693   4.4237   -5.3726   3.5626    4.7693   4.4237
7   4.7254   8.6772   -0.0439   4.2535    4.7254   8.6772
8   4.5599  10.2097   -0.1655   1.5325    4.5599  10.2097
9   3.5059  10.6774   -1.0540   0.4677    3.5059  10.6774
10   2.5222   9.6792   -0.9837  -0.9982    2.5222   9.6792
11   1.0553  15.1533   -1.4669   5.4741    1.0553  15.1533