/* nag_heston_greeks (s30nbc) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <nag.h>
#include <nag_stdlib.h>
#include <nags.h>

int main(void)
{
#ifdef NAG_COLUMN_MAJOR
#define K(I, J)     (J-1)*pdp + I-1
#else
#define K(I, J)     (I-1)*pdp + J-1
#endif

  /* Scalars */
  Integer       exit_status = 0;
  double        corr, eta, grisk, kappa, q, r, s, sigmav, var0;
  Integer       i, j, pdp, m, n;
  /* Arrays */
  double        *charm = 0, *delta = 0, *gamma = 0, *p = 0, *rho = 0,
                *speed = 0, *t = 0, *theta = 0, *vanna = 0, *vega = 0,
                *vomma = 0, *x = 0, *zomma = 0;
  char          put[8+1];
  /* Nag types */
  Nag_OrderType order;
  Nag_CallPut   putnum;
  NagError      fail;

  INIT_FAIL(fail);

  printf("nag_heston_greeks (s30nbc) Example Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n]");
  /* Read put */
  scanf("%s%*[^\n]", put);
  /*
   * nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  putnum = (Nag_CallPut) nag_enum_name_to_value(put);
  /* Read s, r, q */
  scanf("%lf%lf%lf%*[^\n] ", &s, &r, &q);
  /* Read kappa,eta,var0,sigmav,corr,grisk */
  scanf("%lf%lf%lf%*[^\n]", &kappa, &eta, &var0);
  scanf("%lf%lf%lf%*[^\n]", &sigmav, &corr, &grisk);
  /* Read m, n */
  scanf("%ld%ld%*[^\n]", &m, &n);
  if (
    !(charm = NAG_ALLOC(m*n, double)) ||
    !(delta = NAG_ALLOC(m*n, double)) ||
    !(gamma = NAG_ALLOC(m*n, double)) ||
    !(p = NAG_ALLOC(m*n, double)) ||
    !(rho = NAG_ALLOC(m*n, double)) ||
    !(speed = NAG_ALLOC(m*n, double)) ||
    !(t = NAG_ALLOC((n), double)) ||
    !(theta = NAG_ALLOC(m*n, double)) ||
    !(vanna = NAG_ALLOC(m*n, double)) ||
    !(vega = NAG_ALLOC(m*n, double)) ||
    !(vomma = NAG_ALLOC(m*n, double)) ||
    !(x = NAG_ALLOC((m), double)) ||
    !(zomma = NAG_ALLOC(m*n, double))
    )
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
#ifdef NAG_COLUMN_MAJOR
  order = Nag_ColMajor;
  pdp = m;
#else
  order = Nag_RowMajor;
  pdp = n;
#endif

  for (i = 0; i < m; i++)
    scanf("%lf", &x[i]);
  scanf("%*[^\n] ");
  for (i = 0; i < n; i++)
    scanf("%lf", &t[i]);
  scanf("%*[^\n] ");

  /* nag_heston_greeks (s30nbc).
     Heston's model option pricing formula with Greeks
   */
  nag_heston_greeks(order, putnum, m, n, x, s, t, sigmav, kappa, corr, var0,
		    eta, grisk, r, q, p, delta, gamma, vega, theta, rho, vanna,
		    charm, speed, zomma, vomma, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_heston_greeks(s30nbc).\n%s\n",
	     fail.message);
      exit_status = 1;
      goto END;
    }

  printf("\nHeston's Stochastic volatility Model\n");
  switch (putnum)
    {
    case Nag_Call:
      printf("European Call :\n\n");
      break;
    case Nag_Put:
      printf("European Put :\n\n");
    }
  printf("  Spot                   = %10.4f\n", s);
  printf("  Volatility of vol      = %10.4f\n", sigmav);
  printf("  Mean reversion         = %10.4f\n", kappa);
  printf("  Correlation            = %10.4f\n", corr);
  printf("  Variance               = %10.4f\n", var0);
  printf("  Mean of variance       = %10.4f\n", eta);
  printf("  Risk aversion          = %10.4f\n", grisk);
  printf("  Rate                   = %10.4f\n", r);
  printf("  Dividend               = %10.4f\n\n", q);

  for (j = 1; j <= n; j++)
    {
      printf("Time to Expiry : %8.4f\n", t[j-1]);

      printf("%10s%11s%11s%11s%11s%11s%11s\n",
	     "Strike", "Price", "Delta", "Gamma", "Vega", "Theta", "Rho");
      for (i = 1; i <= m; i++)
	printf("%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n", x[i-1],
	       p[K(i, j)], delta[K(i, j)], gamma[K(i, j)], vega[K(i, j)],
	       theta[K(i, j)], rho[K(i, j)]);

      printf("%32s%11s%11s%11s%11s\n",
	     "Vanna", "Charm", "Speed", "Zomma", "Vomma");
      for (i = 1; i <= m; i++)
	printf("%21s %10.4f %10.4f %10.4f %10.4f %10.4f\n", "", vanna[K(i, j)],
	       charm[K(i, j)], speed[K(i, j)], zomma[K(i, j)], vomma[K(i, j)]);
    }
 END:
  if (charm)   NAG_FREE(charm);
  if (delta)   NAG_FREE(delta);
  if (gamma)   NAG_FREE(gamma);
  if (p)       NAG_FREE(p);
  if (rho)     NAG_FREE(rho);
  if (speed)   NAG_FREE(speed);
  if (t)       NAG_FREE(t);
  if (theta)   NAG_FREE(theta);
  if (vanna)   NAG_FREE(vanna);
  if (vega)    NAG_FREE(vega);
  if (vomma)   NAG_FREE(vomma);
  if (x)       NAG_FREE(x);
  if (zomma)   NAG_FREE(zomma);

  return exit_status;
}