/* nag_rand_bb (g05xbc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.4, 2022.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
int get_z(Nag_OrderType order, Integer ntimes, Integer d, Integer a,
Integer npaths, double *z, Integer pdz);
void display_results(Nag_OrderType order, Integer npaths, Integer ntimes,
Integer d, double *b, Integer pdb);
#define CHECK_FAIL(name, fail) \
if (fail.code != NE_NOERROR) { \
printf("Error from %s.\n%s\n", name, fail.message); \
exit_status = -1; \
goto END; \
}
int main(void) {
/* Scalars */
Integer exit_status = 0;
double t0, tend;
Integer a, d, pdb, pdc, pdz, nmove, npaths, ntimes, i;
/* Arrays */
double *b = 0, *c = 0, *intime = 0, *rcomm = 0, *start = 0, *term = 0,
*times = 0, *z = 0;
Integer *move = 0;
/* Nag Types */
NagError fail;
Nag_OrderType order;
INIT_FAIL(fail);
/* Parameters which determine the bridge. */
ntimes = 10;
t0 = 0.0;
npaths = 2;
/* Create a non-free bridge. */
a = 1;
nmove = 0;
d = 3;
#ifdef NAG_COLUMN_MAJOR
order = Nag_ColMajor;
pdz = npaths;
pdb = npaths;
#else
order = Nag_RowMajor;
pdz = d * (ntimes + 1 - a);
pdb = d * (ntimes + 1);
#endif
pdc = d;
#define C(I, J) c[(J - 1) * pdc + I - 1]
/* Allocate memory */
if (!(intime = NAG_ALLOC((ntimes), double)) ||
!(times = NAG_ALLOC((ntimes), double)) ||
!(rcomm = NAG_ALLOC((12 * (ntimes + 1)), double)) ||
!(start = NAG_ALLOC(d, double)) || !(term = NAG_ALLOC(d, double)) ||
!(c = NAG_ALLOC(d * pdc, double)) ||
!(z = NAG_ALLOC(d * (ntimes + 1 - a) * npaths, double)) ||
!(b = NAG_ALLOC(d * (ntimes + 1) * npaths, double)) ||
!(move = NAG_ALLOC(nmove, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Fix the time points at which the bridge is required */
for (i = 0; i < ntimes; i++)
intime[i] = t0 + (double)(i + 1);
tend = t0 + (double)(ntimes + 1);
/* Create a Brownian bridge construction order out of a set of input times
* using nag_rand_bb_make_bridge_order (g05xec).
*/
nag_rand_bb_make_bridge_order(Nag_RLRoundDown, t0, tend, ntimes, intime,
nmove, move, times, &fail);
CHECK_FAIL("nag_rand_bb_make_bridge_order", fail);
/* Initialize the Brownian bridge generator using
* nag_rand_bb_init (g05xac).
*/
nag_rand_bb_init(t0, tend, times, ntimes, rcomm, &fail);
CHECK_FAIL("nag_rand_bb_init (g05xac)", fail);
/* We want the following covariance matrix ... */
C(1, 1) = 6.0;
C(2, 1) = C(1, 2) = 1.0;
C(3, 1) = C(1, 3) = -0.2;
C(2, 2) = 5.0;
C(3, 2) = C(2, 3) = 0.3;
C(3, 3) = 4.0;
/* Cholesky factorize the covariance matrix C, as required by
* nag_rand_bb (g05xbc), using nag_lapacklin_dpotrf (f07fdc).
*/
nag_lapacklin_dpotrf(Nag_ColMajor, Nag_Lower, d, c, pdc, &fail);
CHECK_FAIL("nag_lapacklin_dpotrf", fail);
/* Generate the random numbers z. */
if (get_z(order, ntimes, d, a, npaths, z, pdz) != 0) {
printf("Error generating random numbers\n");
exit_status = -1;
goto END;
}
/* Give start and terminal values of pinned bridge */
start[0] = start[1] = start[2] = 0.0;
term[0] = 1.0;
term[1] = 0.5;
term[2] = 0.0;
/* Generate paths for a free or non-free Wiener process using the
* Brownian bridge algorithm: nag_rand_bb (g05xbc).
*/
nag_rand_bb(order, npaths, d, start, a, term, z, pdz, c, pdc, b, pdb, rcomm,
&fail);
CHECK_FAIL("nag_rand_bb", fail);
/* Display the results */
display_results(order, npaths, ntimes, d, b, pdb);
END:
NAG_FREE(b);
NAG_FREE(c);
NAG_FREE(intime);
NAG_FREE(rcomm);
NAG_FREE(start);
NAG_FREE(term);
NAG_FREE(times);
NAG_FREE(z);
NAG_FREE(move);
return exit_status;
}
int get_z(Nag_OrderType order, Integer ntimes, Integer d, Integer a,
Integer npaths, double *z, Integer pdz) {
/* Scalars */
Integer exit_status = 0;
Integer lseed, lstate, idim, liref, i;
/* Arrays */
Integer seed[1], *iref = 0, state[80];
double *xmean = 0, *stdev = 0;
/* Nag Types */
NagError fail;
INIT_FAIL(fail);
lstate = 80;
lseed = 1;
idim = d * (ntimes + 1 - a);
liref = 32 * idim + 7;
if (!(iref = NAG_ALLOC((liref), Integer)) ||
!(xmean = NAG_ALLOC((idim), double)) ||
!(stdev = NAG_ALLOC((idim), double))) {
printf("Allocation failure in get_z\n");
exit_status = -1;
goto END;
}
/* We now need to generate the input pseudorandom numbers. */
seed[0] = 1023401;
/* Initialize a pseudorandom number generator to give a repeatable sequence
* using nag_rand_init_repeat (g05kfc).
*/
nag_rand_init_repeat(Nag_MRG32k3a, 0, seed, lseed, state, &lstate, &fail);
CHECK_FAIL("nag_rand_init_repeat (g05kfc)", fail);
/* Initialize a scrambled quasi-random number generator using
* nag_rand_quasi_init_scrambled (g05ync).
*/
nag_rand_quasi_init_scrambled(Nag_QuasiRandom_Sobol, Nag_FaureTezuka, idim,
iref, liref, 0, 32, state, &fail);
CHECK_FAIL("nag_rand_quasi_init_scrambled (g05ync)", fail);
for (i = 0; i < idim; i++) {
xmean[i] = 0.0;
stdev[i] = 1.0;
}
/* Generate a (repeatable) Normal quasi-random number sequence using
* nag_rand_quasi_normal (g05yjc).
*/
nag_rand_quasi_normal(order, xmean, stdev, npaths, z, pdz, iref, &fail);
CHECK_FAIL("nag_rand_quasi_normal (g05yjc)", fail);
END:
NAG_FREE(iref);
NAG_FREE(xmean);
NAG_FREE(stdev);
return exit_status;
}
void display_results(Nag_OrderType order, Integer npaths, Integer ntimes,
Integer d, double *b, Integer pdb) {
#define B(I, J) \
(order == Nag_RowMajor ? b[(I - 1) * pdb + J - 1] : b[(J - 1) * pdb + I - 1])
Integer i, p, k;
printf("nag_rand_bb (g05xbc) Example Program Results\n\n");
for (p = 1; p <= npaths; p++) {
printf("Wiener Path %1" NAG_IFMT ", %1" NAG_IFMT "", p, ntimes + 1);
printf(" time steps, %1" NAG_IFMT " dimensions\n", d);
for (k = 1; k <= d; k++)
printf("%10" NAG_IFMT " ", k);
printf("\n");
for (i = 0; i < ntimes + 1; i++) {
printf("%2" NAG_IFMT " ", i + 1);
for (k = 1; k <= d; k++)
printf("%10.4f", B(p, k + i * d));
printf("\n");
}
printf("\n");
}
}