/* nag_specfun_opt_heston_greeks (s30nbc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.1, 2024.
 */

#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, 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_specfun_opt_heston_greeks (s30nbc) Example Program Results\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%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("%" NAG_IFMT "%" NAG_IFMT "%*[^\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_specfun_opt_heston_greeks (s30nbc).
     Heston's model option pricing formula with Greeks
   */
  nag_specfun_opt_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_specfun_opt_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:
  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);

  return exit_status;
}