/* nag_rand_bb_inc (g05xdc) 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(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_repeatable(Nag_MRG32k3a,0,seed, lseed, state, &lstate, &fail);
CHECK_FAIL("nag_rand_init_repeatable",fail);
/* g05ync. Initializes a scrambled quasi-random generator */
nag_quasi_init_scrambled(Nag_QuasiRandom_Sobol, Nag_FaureTezuka, idim,
iref, liref, 0, 32, state, &fail);
CHECK_FAIL("nag_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_quasi_rand_normal(order, xmean, stdev, npaths, z, pdz, iref, &fail);
CHECK_FAIL("nag_quasi_rand_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("%2ld ", p);
}
printf("\n");
for(i=1; i<=ntimesteps+1; i++) {
printf("%2ld ", 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("%2ld ", p);
}
printf("\n ");
for(p=0; p<npaths; p++) printf("%10.4f", analytic[p]);
printf("\n");
}