/* nag_rand_bb_inc_init (g05xcc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 24, 2013.
*/
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg05.h>
#include <nagf07.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_dpotrf(Nag_ColMajor, Nag_Lower, d, c, d, &fail);
CHECK_FAIL("nag_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). Initialises 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_repeatable(Nag_MRG32k3a, 0, seed, lseed, state, &lstate, &fail);
CHECK_FAIL("nag_rand_init_repeatable",fail);
/* nag_rand_normal. Generates a vector of pseudorandom numbers from */
/* a Normal distribution */
nag_rand_normal(nelements, 0.0, 1.0, state, z, &fail);
CHECK_FAIL("nag_rand_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("%3ld ", p);
printf(", ");
printf("%3ld ", ntimes + 1);
printf(" time steps, ");
printf("%3ld ",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("%2ld ", 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("%2ld ", 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("%2ld ", 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);
}