/* nag_rand_bb_inc_init (g05xcc) Example Program.
*
* Copyright 2020 Numerical Algorithms Group.
*
* Mark 27.1, 2020.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
int get_z(Integer nelements, double *z);
void display_results(Nag_OrderType order, Integer npaths, double t0,
double tend, Integer ntimes, double intime[], Integer d,
double start[], double bb[], double bd[], Integer pdb);
#define CHECK_FAIL(name, fail) \
if (fail.code != NE_NOERROR) { \
printf("Error calling %s\n%s\n", name, fail.message); \
exit_status = -1; \
goto END; \
}
int main(void) {
#define C(I, J) c[(J - 1) * d + I - 1]
Integer exit_status = 0;
NagError fail;
/* Scalars */
double t0, tend;
Integer a, d, pdb, pdc, pdz, nmove, npaths, ntimes, i;
/* Arrays */
double *bb = 0, *c = 0, *intime = 0, *rcommb = 0, *start = 0, *term = 0,
*times = 0, *zb = 0, *rcommd = 0, *zd = 0, *diff = 0, *bd = 0;
Integer *move = 0;
INIT_FAIL(fail);
/* Parameters which determine the bridge */
ntimes = 10;
t0 = 0.0;
npaths = 2;
a = 0;
nmove = 0;
d = 2;
pdz = npaths;
pdb = npaths;
pdc = d;
/* Allocate memory */
if (!(intime = NAG_ALLOC((ntimes), double)) ||
!(times = NAG_ALLOC((ntimes), double)) ||
!(rcommb = NAG_ALLOC((12 * (ntimes + 1)), double)) ||
!(rcommd = NAG_ALLOC((12 * (ntimes + 1)), double)) ||
!(start = NAG_ALLOC(d, double)) || !(term = NAG_ALLOC(d, double)) ||
!(diff = NAG_ALLOC(d, double)) || !(c = NAG_ALLOC(pdc * d, double)) ||
!(zb = NAG_ALLOC(pdz * d * (ntimes + 1 - a), double)) ||
!(zd = NAG_ALLOC(pdz * d * (ntimes + 1 - a), double)) ||
!(bb = NAG_ALLOC(pdb * d * (ntimes + 1), double)) ||
!(bd = NAG_ALLOC(pdb * d * (ntimes + 1), 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);
/* g05xec. Creates a Brownian bridge construction order out of a set of */
/* input times */
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);
/* g05xac. Initializes the Brownian bridge generator */
nag_rand_bb_init(t0, tend, times, ntimes, rcommb, &fail);
CHECK_FAIL("nag_rand_bb_init", fail);
start[0] = 0.0, start[1] = 2.0;
/* We want the following covariance matrix */
C(1, 1) = 6.0;
C(2, 1) = -1.0;
C(1, 2) = -1.0;
C(2, 2) = 5.0;
/* nag_rand_bb uses Cholesky factorization of the covariance matrix C */
/* f07fdc. Cholesky factorization of real positive definite matrix */
nag_lapacklin_dpotrf(Nag_ColMajor, Nag_Lower, d, c, d, &fail);
CHECK_FAIL("nag_lapacklin_dpotrf", fail);
/* Generate the random numbers */
if (get_z(npaths * d * (ntimes + 1 - a), zb) != 0) {
printf("Error generating random numbers\n");
exit_status = -1;
goto END;
}
/* Copy the random numbers for call to g05xdc */
for (i = 0; i < npaths * d * (ntimes + 1 - a); i++)
zd[i] = zb[i];
/* g05xbc. Generate paths for a free or non-free Wiener process using the */
/* Brownian bridge algorithm */
nag_rand_bb(Nag_ColMajor, npaths, d, start, a, term, zb, pdz, c, pdc, bb, pdb,
rcommb, &fail);
CHECK_FAIL("nag_rand_bb", fail);
/* nag_rand_bb_inc_init (g05xcc). Initializes the generator which backs out
* the increments of sample paths generated by a Brownian bridge algorithm */
nag_rand_bb_inc_init(t0, tend, times, ntimes, rcommd, &fail);
CHECK_FAIL("nag_rand_bb_inc_init", fail);
/* g05xdc. Backs out the increments from sample paths generated by a */
/* Brownian bridge algorithm */
nag_rand_bb_inc(Nag_ColMajor, npaths, d, a, diff, zd, pdz, c, pdc, bd, pdb,
rcommd, &fail);
CHECK_FAIL("nag_rand_bb_inc", fail);
/* Display the results */
display_results(Nag_ColMajor, npaths, t0, tend, ntimes, intime, d, start, bb,
bd, pdb);
END:;
NAG_FREE(bb);
NAG_FREE(bd);
NAG_FREE(c);
NAG_FREE(intime);
NAG_FREE(rcommb);
NAG_FREE(rcommd);
NAG_FREE(start);
NAG_FREE(term);
NAG_FREE(diff);
NAG_FREE(times);
NAG_FREE(zb);
NAG_FREE(zd);
NAG_FREE(move);
return exit_status;
#undef C
}
int get_z(Integer nelements, double *z) {
NagError fail;
Integer lseed, lstate, exit_status = 0;
/* Arrays */
Integer seed[1];
Integer state[80];
lstate = 80;
lseed = 1;
INIT_FAIL(fail);
/* We now need to generate the input pseudorandom numbers */
seed[0] = 1023401;
/* g05kfc. Initializes a pseudorandom number generator */
/* to give a repeatable sequence */
nag_rand_init_repeat(Nag_MRG32k3a, 0, seed, lseed, state, &lstate, &fail);
CHECK_FAIL("nag_rand_init_repeat", fail);
/* nag_rand_dist_normal. Generates a vector of pseudorandom numbers from */
/* a Normal distribution */
nag_rand_dist_normal(nelements, 0.0, 1.0, state, z, &fail);
CHECK_FAIL("nag_rand_dist_normal", fail);
END:;
return exit_status;
}
void display_results(Nag_OrderType order, Integer npaths, double t0,
double tend, Integer ntimes, double intime[], Integer d,
double start[], double bb[], double bd[], Integer pdb) {
// Order consistend with Nag_RowMajor
#define BB(I, J) bb[(I - 1) * pdb + J - 1]
#define BD(I, J) bd[(I - 1) * pdb + J - 1]
#define CUM(I) cum[I - 1]
Integer i, p, k;
double *cum = 0;
if (!(cum = NAG_ALLOC(d, double))) {
printf("Error allocating memory in display_results\n");
return;
}
printf("nag_rand_bb_inc_init (g05xcc) Example Program Results\n\n");
for (p = 1; p <= npaths; p++) {
printf("Weiner Path ");
printf("%3" NAG_IFMT " ", p);
printf(", ");
printf("%3" NAG_IFMT " ", ntimes + 1);
printf(" time steps, ");
printf("%3" NAG_IFMT " ", d);
printf(" dimensions \n");
printf(" Output of g05xbc Output of g05xdc Sum of g05xdc \n");
for (k = 1; k <= d; k++)
CUM(k) = start[k - 1];
/* Print first point */
i = 0;
printf("%2" NAG_IFMT " ", i + 1);
if (order == Nag_RowMajor) {
for (k = 1; k <= d; k++)
CUM(k) += BD(p, k) * (intime[i] - t0);
for (k = 1; k <= d; k++)
printf(" %8.4f", BB(p, k));
for (k = 1; k <= d; k++)
printf(" %8.4f", BD(p, k));
for (k = 1; k <= d; k++)
printf(" %8.4f", CUM(k));
} else {
for (k = 1; k <= d; k++)
CUM(k) += BD(k, p) * (intime[i] - t0);
for (k = 1; k <= d; k++)
printf(" %8.4f", BB(k, p));
for (k = 1; k <= d; k++)
printf(" %8.4f", BD(k, p));
for (k = 1; k <= d; k++)
printf(" %8.4f", CUM(k));
}
printf("\n");
/* Print intermediate points */
for (i = 1; i < ntimes; i++) {
printf("%2" NAG_IFMT " ", i + 1);
if (order == Nag_RowMajor) {
for (k = 1; k <= d; k++)
CUM(k) += BD(p, i * d + k) * (intime[i] - intime[i - 1]);
for (k = 1; k <= d; k++)
printf(" %8.4f", BB(p, i * d + k));
for (k = 1; k <= d; k++)
printf(" %8.4f", BD(p, i * d + k));
for (k = 1; k <= d; k++)
printf(" %8.4f", CUM(k));
} else {
for (k = 1; k <= d; k++)
CUM(k) += BD(i * d + k, p) * (intime[i] - intime[i - 1]);
for (k = 1; k <= d; k++)
printf(" %8.4f", BB(i * d + k, p));
for (k = 1; k <= d; k++)
printf(" %8.4f", BD(i * d + k, p));
for (k = 1; k <= d; k++)
printf(" %8.4f", CUM(k));
}
printf("\n");
}
/* Print final point */
i = ntimes;
printf("%2" NAG_IFMT " ", i + 1);
if (order == Nag_RowMajor) {
for (k = 1; k <= d; k++)
CUM(k) += BD(p, i * d + k) * (tend - intime[i - 1]);
for (k = 1; k <= d; k++)
printf(" %8.4f", BB(p, i * d + k));
for (k = 1; k <= d; k++)
printf(" %8.4f", BD(p, i * d + k));
for (k = 1; k <= d; k++)
printf(" %8.4f", CUM(k));
} else {
for (k = 1; k <= d; k++)
CUM(k) += BD(i * d + k, p) * (tend - intime[i - 1]);
for (k = 1; k <= d; k++)
printf(" %8.4f", BB(i * d + k, p));
for (k = 1; k <= d; k++)
printf(" %8.4f", BD(i * d + k, p));
for (k = 1; k <= d; k++)
printf(" %8.4f", CUM(k));
}
printf("\n\n");
}
NAG_FREE(cum);
}