NAG Library Routine Document

g05xcf (bb_inc_init)

1
Purpose

g05xcf initializes the Brownian bridge increments generator g05xdf. It must be called before any calls to g05xdf.

2
Specification

Fortran Interface
Subroutine g05xcf ( t0, tend, times, ntimes, rcomm, ifail)
Integer, Intent (In):: ntimes
Integer, Intent (Inout):: ifail
Real (Kind=nag_wp), Intent (In):: t0, tend, times(ntimes)
Real (Kind=nag_wp), Intent (Out):: rcomm(12*(ntimes+1))
C Header Interface
#include <nagmk26.h>
void  g05xcf_ (const double *t0, const double *tend, const double times[], const Integer *ntimes, double rcomm[], Integer *ifail)

3
Description

3.1
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 2.6 in the G05 Chapter Introduction. We briefly recall some notation and definitions.
Fix two times t0<T and let ti 1iN  be any set of time points satisfying t0<t1<t2<⋯<tN<T. Let Xti 1iN  denote a d-dimensional Wiener sample path at these time points, and let C be any d by d matrix such that CCT is the desired covariance structure for the Wiener process. Each point Xti of the sample path is constructed according to the Brownian bridge interpolation algorithm (see Glasserman (2004) or Section 2.6 in the G05 Chapter Introduction). We always start at some fixed point Xt0 = xd . If we set XT =x+ C T-t0 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 XT = wd , then 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, 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 .  
Such increments can be useful in computing numerical solutions to stochastic differential equations driven by (free or non-free) Wiener processes.

3.2
Implementation

Conceptually, the output of the Wiener increments generator is the same as if g05xaf and g05xbf 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 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 2.6.2 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 Wiener sample paths each of dimension d. The main input to the Brownian bridge increments generator is then an array of quasi-random points Z1,Z2,…,ZP where each point Zp = Z1p,Z2p,,ZDp  has dimension D=dN+1 or D=dN depending on whether a free or non-free Wiener process is required. When g05xdf is called, the pth sample path for 1pP is constructed as follows: if a non-free Wiener process is required set XT equal to the terminal value w, otherwise construct XT as
XT = Xt0 + C T-t0 Z1p Zdp  
where C is the matrix described in Section 3.1. The array times holds the remaining time points t1 , t2 ,… tN  in the order in which the bridge is to be constructed. For each j=1,…,N set r=timesj, find
q = max t0, timesi : 1i<j , timesi < r  
and
s = min T, timesi : 1i<j , timesi > r  
and construct the point Xr as
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 or a=1 depending on whether a free or non-free Wiener process is required. The routine g05xef 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  
are computed.

4
References

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

5
Arguments

1:     t0 – Real (Kind=nag_wp)Input
On entry: the starting value t0 of the time interval.
2:     tend – Real (Kind=nag_wp)Input
On entry: the end value T of the time interval.
Constraint: tend>t0.
3:     timesntimes – Real (Kind=nag_wp) arrayInput
On entry: the points in the time interval t0,T at which the Wiener process is to be constructed. The order in which points are listed in times determines the bridge construction order. The routine g05xef can be used to create predefined bridge construction orders from a set of input times.
Constraints:
  • t0<timesi<tend, for i=1,2,,ntimes;
  • timesi timesj, for i,j=1,2,ntimes and ij.
4:     ntimes – IntegerInput
On entry: the length of times, denoted by N in Section 3.1.
Constraint: ntimes1.
5:     rcomm12×ntimes+1 – Real (Kind=nag_wp) arrayCommunication Array
On exit: communication array, used to store information between calls to g05xdf. This array must not be directly modified.
6:     ifail – IntegerInput/Output
On entry: ifail must be set to 0, -1 or 1. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1 or 1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this argument, the recommended value is 0. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6
Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, tend=value and t0=value.
Constraint: tend>t0.
ifail=2
On entry, ntimes=value.
Constraint: ntimes1.
ifail=3
On entry, timesvalue=value, t0=value and tend=value.
Constraint: t0<timesi<tend for all i.
ifail=4
On entry, timesi=timesj=value, for i=value and j=value.
Constraint: all elements of times must be unique.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7
Accuracy

Not applicable.

8
Parallelism and Performance

g05xcf is not threaded in any implementation.

9
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 N 4  or even N 2 , which could be very inefficient when N is large. g05xcf 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.

10
Example

The following example program calls g05xaf and g05xbf to generate two sample paths from a two-dimensional free Wiener process. It then calls g05xcf and g05xdf with the same input arguments to obtain the scaled increments of the Wiener sample paths. Lastly, the program prints the Wiener sample paths from g05xbf, the scaled increments from g05xdf, 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 g05xbf.
Please see Section 10 in g05xdf for additional examples.

10.1
Program Text

Program Text (g05xcfe.f90)

10.2
Program Data

None.

10.3
Program Results

Program Results (g05xcfe.r)