/* nag_specfun_opt_heston_more_greeks (s30ndc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 */

#include <nag.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, 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, *q = 0, *r = 0, *dp_dq = 0, *dp_deta = 0, *dp_dkappa = 0,
         *dp_dsigmav = 0, *dp_dcorr = 0, *dp_dx = 0, *dp_dgrisk = 0;

  char put[8 + 1];
  /* Nag types */
  Nag_OrderType order;
  Nag_CallPut putnum;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_specfun_opt_heston_more_greeks (s30ndc) Example Program Results\
\n\n");
  /* Skip heading in data file */
  scanf("%*[^\n]");
  /* Read put */
  scanf("%8s%*[^\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%*[^\n] ", &s);
  /* 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("%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &m, &n);
  if (!(x          = NAG_ALLOC((m), double))   ||
      !(t          = NAG_ALLOC((n), double))   ||
      !(r          = NAG_ALLOC((n), double))   ||
      !(q          = NAG_ALLOC((n), double))   ||      
      !(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)) || 
      !(theta      = NAG_ALLOC(m * n, double)) ||
      !(vanna      = NAG_ALLOC(m * n, double)) ||
      !(vega       = NAG_ALLOC(m * n, double)) ||
      !(vomma      = NAG_ALLOC(m * n, double)) ||
      !(zomma      = NAG_ALLOC(m * n, double)) ||
      !(dp_dq      = NAG_ALLOC(m * n, double)) ||
      !(dp_deta    = NAG_ALLOC(m * n, double)) ||
      !(dp_dkappa  = NAG_ALLOC(m * n, double)) ||
      !(dp_dsigmav = NAG_ALLOC(m * n, double)) ||
      !(dp_dcorr   = NAG_ALLOC(m * n, double)) ||
      !(dp_dx      = NAG_ALLOC(m * n, double)) ||
      !(dp_dgrisk  = 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] ");

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

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

  /* nag_specfun_opt_heston_more_greeks (s30ndc).
   *  Heston's model option pricing formula with More_Greeks
   */

  nag_specfun_opt_heston_more_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, dp_dx, dp_dq,
                                     dp_deta, dp_dkappa, dp_dsigmav, dp_dcorr,
                                     dp_dgrisk, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_specfun_opt_heston_more_greeks (s30ndc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
 
  printf(" Heston's Stochastic volatility Model\n");
  switch (putnum) {
  case Nag_Call:
    printf(" European Call :\n");
    break;
  case Nag_Put:
    printf(" European Put :\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\n", grisk);

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

    for (i = 1; i <= m; i++){
      printf("     Strike     Price      Rate       Dividend\n");
      printf(" %10.4f %10.4f %10.4f %10.4f\n", x[i - 1],
             p[K(i, j)], r[j - 1], q[j - 1]);

      printf("     Delta      Gamma      Vega       Theta      Rho\n");
      printf(" %10.4f %10.4f %10.4f %10.4f %10.4f\n", delta[K(i, j)],
             gamma[K(i, j)], vega[K(i, j)], theta[K(i, j)], rho[K(i, j)]);

      printf("     Vanna      Charm      Speed      Zomma      Vomma\n");
      printf(" %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)]);

      printf("     dp_dx      dp_dq      dp_deta    dp_dkappa  dp_dsigmav \
dp_dcorr   dp_dgrisk\n");
      printf(" %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n",
             dp_dx[K(i, j)], dp_dq[K(i, j)], dp_deta[K(i, j)],
             dp_dkappa[K(i, j)], dp_dsigmav[K(i, j)], dp_dcorr[K(i, j)],
             dp_dgrisk[K(i, j)]);
      
    }
  }
END:
  NAG_FREE(charm);
  NAG_FREE(delta);
  NAG_FREE(gamma);
  NAG_FREE(p);
  NAG_FREE(rho);
  NAG_FREE(speed);
  NAG_FREE(t);
  NAG_FREE(theta);
  NAG_FREE(vanna);
  NAG_FREE(vega);
  NAG_FREE(vomma);
  NAG_FREE(x);
  NAG_FREE(zomma);
  NAG_FREE(dp_dx);
  NAG_FREE(dp_dq);
  NAG_FREE(dp_deta);
  NAG_FREE(dp_dkappa);
  NAG_FREE(dp_dsigmav);
  NAG_FREE(dp_dcorr);
  NAG_FREE(dp_dgrisk);

  return exit_status;
}