/* nag_rand_bb_inc (g05xdc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */
#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("%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");
}