/* nag_pde_bs_1d_means (d03nec) Example Program.
 *
 * Copyright 2014 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 = 0;
  double      *delta = 0, *f = 0, *gamma = 0, *lambda = 0, q[3], ra[3],
              *rho = 0, *s = 0;
  double      siga[3], *t = 0, *theta = 0, *td = 0, *rd = 0, *sigd = 0,
              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;

  INIT_FAIL(fail);

  printf("nag_pde_bs_1d_means (d03nec) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Read problem parameters */
  scanf("%lf", &x);
  scanf("%lf", &tmat);
  scanf("%ld%ld", &ns, &nt);
  scanf("%lf%lf", &smin, &smax);
  scanf("%lf%lf", &tmin, &tmax);
  scanf("%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)))
    {
      printf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }

  /* Read discrete times and data values for r and sigma */
  for (i = 0; i < ntd; i++) scanf("%lf", &td[i]);
  for (i = 0; i < ntd; i++) scanf("%lf", &rd[i]);
  for (i = 0; i < ntd; i++) scanf("%lf", &sigd[i]);

  /* 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)
        {
          printf("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)
        {
          printf("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)
            {
              printf(
                      "Error from nag_pde_bs_1d_analytic (d03ndc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        }
    }

  /* Output option values */

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

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

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

      for (i = 1; i <= ns; i++)
        {
          printf(" %13.4e  | ", s[i-1]);
          switch (igreek)
            {
            case 0:
              for (j = 1; j <= nt; j++) printf(" %13.4e", THETA(i, j));
              break;
            case 1:
              for (j = 1; j <= nt; j++) printf(" %13.4e", DELTA(i, j));
              break;
            case 2:
              for (j = 1; j <= nt; j++) printf(" %13.4e", GAMMA(i, j));
              break;
            case 3:
              for (j = 1; j <= nt; j++) printf(" %13.4e", LAMBDA(i, j));
              break;
            case 4:
              for (j = 1; j <= nt; j++) printf(" %13.4e", RHO(i, j));
              break;
            default:
              break;
            }
          printf("\n");
        }
    }

 END:
  NAG_FREE(s);
  NAG_FREE(t);
  NAG_FREE(f);
  NAG_FREE(theta);
  NAG_FREE(delta);
  NAG_FREE(gamma);
  NAG_FREE(lambda);
  NAG_FREE(rho);
  NAG_FREE(td);
  NAG_FREE(rd);
  NAG_FREE(sigd);

  return exit_status;
}