/* nag_pde_bs_1d_means (d03nec) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd03.h>

#define F(I,J) f[ns*((J)-1)+(I)-1]
#define THETA(I,J) theta[ns*((J)-1)+(I)-1]
#define DELTA(I,J) delta[ns*((J)-1)+(I)-1]
#define GAMMA(I,J) gamma[ns*((J)-1)+(I)-1]
#define LAMBDA(I,J) lambda[ns*((J)-1)+(I)-1]
#define RHO(I,J) rho[ns*((J)-1)+(I)-1]

int main(void)
{
  double ds, dt, tmat, x;
  Integer i, igreek, j, ns, nt, ntd, exit_status;
  double *delta, *f, *gamma, *lambda, q[3], ra[3], *rho, *s,
    siga[3], *t, *theta, *td, *rd, *sigd, smin, smax, tmin, tmax;
  Nag_Boolean gprnt[5]={Nag_TRUE, Nag_TRUE, Nag_TRUE, Nag_TRUE, Nag_TRUE};
  Nag_Boolean tdpar[3];
  const char *gname[5]={"Theta", "Delta", "Gamma", "Lambda", "Rho"};
  NagError fail;

  /* Skip heading in data file */

  Vscanf("%*[^\n] ");
  exit_status = 0;

  /* Read problem parameters */

  Vscanf("%lf", &x);
  Vscanf("%lf", &tmat);
  Vscanf("%ld%ld", &ns, &nt);
  Vscanf("%lf%lf", &smin, &smax);
  Vscanf("%lf%lf", &tmin, &tmax);
  Vscanf("%ld", &ntd);

  /* Allocate memory */

  if ( !(s = NAG_ALLOC(ns, double)) ||
       !(t = NAG_ALLOC(nt, double)) ||
       !(f = NAG_ALLOC(ns*nt, double)) ||
       !(theta = NAG_ALLOC(ns*nt, double)) ||
       !(delta = NAG_ALLOC(ns*nt, double)) ||
       !(gamma = NAG_ALLOC(ns*nt, double)) ||
       !(lambda = NAG_ALLOC(ns*nt, double)) ||
       !(rho = NAG_ALLOC(ns*nt, double)) ||
       !(td = NAG_ALLOC(ntd, double)) ||
       !(rd = NAG_ALLOC(ntd, double)) ||
       !(sigd = NAG_ALLOC(ntd, double)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }

  /* Read discrete times and data values for r and sigma */

  for (i=0; i<ntd; i++) Vscanf("%lf", &td[i]);
  for (i=0; i<ntd; i++) Vscanf("%lf", &rd[i]);
  for (i=0; i<ntd; i++) Vscanf("%lf", &sigd[i]);

  INIT_FAIL(fail);
  Vprintf("nag_pde_bs_1d_means (d03nec) Example Program Results\n\n");

  /* Set up input parameters for nag_pde_bs_1d_analytic (d03ndc) */

  s[0] = smin;
  s[ns-1] = smax;
  t[0] = tmin;
  t[nt-1] = tmax;
  tdpar[0] = Nag_TRUE;
  tdpar[1] = Nag_FALSE;
  tdpar[2] = Nag_TRUE;
  q[0] = 0.0;

  ds = (s[ns-1]-s[0])/(ns-1.0);
  dt = (t[nt-1]-t[0])/(nt-1.0);

  /* Loop over times */

  for (j=1; j<=nt; j++)
    {
      t[j-1] = t[0] + (j-1)*dt;

      /* Find average values of r and sigma */

      /* nag_pde_bs_1d_means (d03nec).
       * Compute average values for nag_pde_bs_1d_analytic
       * (d03ndc)
       */
      nag_pde_bs_1d_means(t[j-1], tmat, ntd, td, rd, ra, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from nag_pde_bs_1d_means (d03nec).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      /* nag_pde_bs_1d_means (d03nec), see above. */
      nag_pde_bs_1d_means(t[j-1], tmat, ntd, td, sigd, siga, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from nag_pde_bs_1d_means (d03nec).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      /* Loop over stock prices */

      for (i=1; i<=ns; i++)
        {
          s[i-1] = s[0] + (i-1)*ds;

          /* Evaluate analytic solution of Black-Scholes equation */

          /* nag_pde_bs_1d_analytic (d03ndc).
           * Analytic solution of the Black-Scholes equations
           */
          nag_pde_bs_1d_analytic(Nag_AmericanCall, x, s[i-1], t[j-1], tmat,
                                 tdpar, ra, q, siga, &F(i,j), &THETA(i,j),
                                 &DELTA(i,j), &GAMMA(i,j), &LAMBDA(i,j),
                                 &RHO(i,j), &fail);

          if (fail.code != NE_NOERROR)
            {
              Vprintf("Error from nag_pde_bs_1d_analytic (d03ndc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

        }
    }

  /* Output option values */

  Vprintf("\n");
  Vprintf("Option Values\n");
  Vprintf("-------------\n");
  Vprintf(" Stock Price  |   Time to Maturity (months)\n");
  Vprintf("              | ");
  for (i=0; i<nt; i++) Vprintf(" %11.4e", 12.0*(tmat-t[i]));
  Vprintf("\n");
  for (i=0; i<64; i++) Vprintf("-");
  Vprintf("\n");
  for (i=1; i<=ns; i++)
    {
      Vprintf(" %11.4e  | ", s[i-1]);
      for (j=1; j<=nt; j++) Vprintf(" %11.4e", F(i,j));
      Vprintf("\n");
    }

  for (igreek = 0; igreek < 5; igreek++)
    {
      if (!gprnt[igreek]) continue;

      Vprintf("\n");
      Vprintf("%s\n", gname[igreek]);
      for (i=0; i<(Integer)strlen(gname[igreek]); i++) Vprintf("-");
      Vprintf("\n");
      Vprintf(" Stock Price  |   Time to Maturity (months)\n");
      Vprintf("              | ");
      for (i=0; i<nt; i++) Vprintf(" %11.4e", 12.0*(tmat-t[i]));
      Vprintf("\n");
      for (i=0; i<64; i++) Vprintf("-");
      Vprintf("\n");

      for (i=1; i<=ns; i++)
        {
          Vprintf(" %11.4e  | ", s[i-1]);
          switch (igreek)
            {
            case 0:
              for (j=1; j<=nt; j++) Vprintf(" %11.4e", THETA(i,j));
              break;
            case 1:
              for (j=1; j<=nt; j++) Vprintf(" %11.4e", DELTA(i,j));
              break;
            case 2:
              for (j=1; j<=nt; j++) Vprintf(" %11.4e", GAMMA(i,j));
              break;
            case 3:
              for (j=1; j<=nt; j++) Vprintf(" %11.4e", LAMBDA(i,j));
              break;
            case 4:
              for (j=1; j<=nt; j++) Vprintf(" %11.4e", RHO(i,j));
              break;
            default:
              break;
            }
          Vprintf("\n");
        }
    }
 END:
  if (s) NAG_FREE(s);
  if (t) NAG_FREE(t);
  if (f) NAG_FREE(f);
  if (theta) NAG_FREE(theta);
  if (delta) NAG_FREE(delta);
  if (gamma) NAG_FREE(gamma);
  if (lambda) NAG_FREE(lambda);
  if (rho) NAG_FREE(rho);
  if (td) NAG_FREE(td);
  if (rd) NAG_FREE(rd);
  if (sigd) NAG_FREE(sigd);

  return exit_status;
}