/* nag_rand_bb_inc (g05xdc) Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.2, 2021.
*/
#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(Integer npaths, Integer ntimes, double *St,
double *analytic);
#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) {
Integer exit_status = 0;
NagError fail;
/* Scalars */
double t0, tend, r, S0, sigma;
Integer a, d, pdb, pdc, pdz, nmove, npaths, ntimesteps, i, p;
/* Arrays */
double *b = 0, c[1], *t = 0, *rcomm = 0, *diff = 0, *times = 0, *z = 0,
*analytic = 0, *St = 0;
Integer *move = 0;
INIT_FAIL(fail);
/* We wish to solve the stochastic differential equation (SDE)
* dSt = r * St * dt + sigma * St * dXt
* where X is a one dimensional Wiener process. This means we have
* a = 0
* d = 1
* c = 1
* We now set the other parameters of the SDE and the Euler-Maruyama scheme
*
* Initial value of the process */
S0 = 1.0;
r = 0.05;
sigma = 0.12;
/* Number of paths to simulate */
npaths = 3;
/* The time interval [t0,T] on which to solve the SDE */
t0 = 0.0;
tend = 1.0;
/* The time steps in the discretization of [t0,T] */
ntimesteps = 20;
/* Other bridge parameters */
c[0] = 1.0;
a = 0;
nmove = 0;
d = 1;
pdz = d * (ntimesteps + 1 - a);
pdb = d * (ntimesteps + 1);
pdc = d;
/* Allocate memory */
if (!(t = NAG_ALLOC((ntimesteps), double)) ||
!(times = NAG_ALLOC((ntimesteps), double)) ||
!(rcomm = NAG_ALLOC((12 * (ntimesteps + 1)), double)) ||
!(diff = NAG_ALLOC(d, double)) ||
!(z = NAG_ALLOC(pdz * npaths, double)) ||
!(b = NAG_ALLOC(pdb * npaths, double)) ||
!(move = NAG_ALLOC(nmove, Integer)) ||
!(St = NAG_ALLOC(npaths * (ntimesteps + 1), double)) ||
!(analytic = NAG_ALLOC(npaths, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Fix the time points at which the bridge is required */
for (i = 0; i < ntimesteps; i++) {
t[i] = t0 + (i + 1) * (tend - t0) / (double)(ntimesteps + 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, ntimesteps, t, nmove,
move, times, &fail);
CHECK_FAIL("nag_rand_bb_make_bridge_order", fail);
/* 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, ntimesteps, rcomm, &fail);
CHECK_FAIL("nag_rand_bb_inc_init", fail);
/* Generate the random numbers */
if (get_z(Nag_RowMajor, ntimesteps, d, a, npaths, z, pdz) != 0) {
printf("Error generating random numbers\n");
exit_status = -1;
goto END;
}
/* nag_rand_bb_inc (g05xdc). Backs out the increments from sample paths
* generated by a Brownian bridge algorithm */
nag_rand_bb_inc(Nag_RowMajor, npaths, d, a, diff, z, pdz, c, pdc, b, pdb,
rcomm, &fail);
CHECK_FAIL("nag_rand_bb_inc", fail);
/* Do the Euler-Maruyama time stepping for SDE
* dSt = r * St * dt + sigma * St * dXt */
// Definitions consistent with Nag_RowMajor
#define B(I, J) b[(I - 1) * pdb + J - 1]
#define ST(I, J) St[(I - 1) * (ntimesteps + 1) + J - 1]
for (p = 1; p <= npaths; p++) {
ST(p, 1) = S0 + (t[0] - t0) * (r * S0 + sigma * S0 * B(p, 1));
}
for (i = 2; i <= ntimesteps; i++) {
for (p = 1; p <= npaths; p++) {
ST(p, i) = ST(p, i - 1) +
(t[i - 1] - t[i - 2]) *
(r * ST(p, i - 1) + sigma * ST(p, i - 1) * B(p, i));
}
}
for (p = 1; p <= npaths; p++) {
ST(p, i) =
ST(p, i - 1) + (tend - t[ntimesteps - 1]) *
(r * ST(p, i - 1) + sigma * ST(p, i - 1) * B(p, i));
}
/* Compute the analytic solution:
* ST = S0*exp( (r-sigma*sigma/2)*T + sigma*WT )
* where WT =law sqrt(T)Z is the Wiener process at time T.
*
* The first quasi-random point in z is always used to compute
* the final value of WT (equivalently BT, the final value of the
* Brownian bridge). Hence we have that
* WT = sqrt(tend-t0)*z[0]
*/
for (p = 0; p < npaths; p++) {
analytic[p] = S0 * exp((r - 0.5 * sigma * sigma) * (tend - t0) +
sigma * sqrt(tend - t0) * z[p * (ntimesteps + 1)]);
}
/* Display the results */
display_results(npaths, ntimesteps, St, analytic);
END:;
NAG_FREE(b);
NAG_FREE(t);
NAG_FREE(rcomm);
NAG_FREE(analytic);
NAG_FREE(diff);
NAG_FREE(St);
NAG_FREE(times);
NAG_FREE(z);
NAG_FREE(move);
return exit_status;
#undef C
}
int get_z(Nag_OrderType order, Integer ntimes, Integer d, Integer a,
Integer npaths, double *z, Integer pdz) {
NagError fail;
Integer lseed, lstate, seed[1], idim, liref, *iref = 0, state[80], i;
Integer exit_status = 0;
double *xmean = 0, *stdev = 0;
lstate = 80;
lseed = 1;
INIT_FAIL(fail);
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;
/* 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);
/* g05ync. Initializes a scrambled quasi-random generator */
nag_rand_quasi_init_scrambled(Nag_QuasiRandom_Sobol, Nag_FaureTezuka, idim,
iref, liref, 0, 32, state, &fail);
CHECK_FAIL("nag_rand_quasi_init_scrambled", fail);
for (i = 0; i < idim; i++) {
xmean[i] = 0.0;
stdev[i] = 1.0;
}
/* g05skc. Generates a Normal quasi-random number sequence */
nag_rand_quasi_normal(order, xmean, stdev, npaths, z, pdz, iref, &fail);
CHECK_FAIL("nag_rand_quasi_normal", fail);
END:
NAG_FREE(iref);
NAG_FREE(xmean);
NAG_FREE(stdev);
return exit_status;
}
void display_results(Integer npaths, Integer ntimesteps, double *St,
double *analytic) {
Integer i, p;
printf("nag_rand_bb_inc (g05xdc) Example Program Results\n\n");
printf("Euler-Maruyama solution for Geometric Brownian motion\n");
printf(" ");
for (p = 1; p <= npaths; p++) {
printf("Path");
printf("%2" NAG_IFMT " ", p);
}
printf("\n");
for (i = 1; i <= ntimesteps + 1; i++) {
printf("%2" NAG_IFMT " ", i);
for (p = 1; p <= npaths; p++)
printf("%10.4f", ST(p, i));
printf("\n");
}
printf("\nAnalytic solution at final time step\n");
printf(" ");
for (p = 1; p <= npaths; p++) {
printf("%7s", "Path");
printf("%2" NAG_IFMT " ", p);
}
printf("\n ");
for (p = 0; p < npaths; p++)
printf("%10.4f", analytic[p]);
printf("\n");
}