/* nag_tsa_arma_filter (g13bac) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg13.h>

int main(void)
{
  /* Scalars */
  double a1, a2, cx, cy;
  Integer i, idd, ii, ij, iqxd,
         j, k, n, nb, ni, npar, nparx, nx, ny,
         nser, npara, tdxxy, tdmrx, ldparx, tdparx;
  Integer exit_status = 0;

  /* Arrays */
  double *b = 0, *fsd = 0, *fva = 0, *par = 0, *parx = 0,
         *x = 0, *y = 0, *rms = 0, *parxx = 0;
  Integer mr[14], mrx[7], *mrxx = 0;
  Nag_ArimaOrder arimaj, arimaf, arimav;
  Nag_TransfOrder transfj;
  Nag_G13_Opt options;
  NagError fail;

  INIT_FAIL(fail);

  exit_status = 0;

  /* Initialize the options structure used by nag_tsa_multi_inp_model_forecast
   * (g13bjc) */
  /* nag_tsa_options_init (g13bxc).
   * Initialization function for option setting
   */
  nag_tsa_options_init(&options);

  printf("nag_tsa_arma_filter (g13bac) Example Program Results\n");

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

  if (nx > 0) {
    /* Allocate array x */
    if (!(x = NAG_ALLOC(nx + 2, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    for (i = 1; i <= nx; ++i)
      scanf("%lf", &x[i - 1]);
    scanf("%*[^\n] ");

    /* Read univariate Arima for series */
    for (i = 1; i <= 7; ++i)
      scanf("%" NAG_IFMT "", &mrx[i - 1]);
    scanf("%*[^\n] ");
    scanf("%lf%*[^\n] ", &cx);

    nparx = mrx[0] + mrx[2] + mrx[3] + mrx[5];

    arimaj.p = mrx[0];
    arimaj.d = mrx[1];
    arimaj.q = mrx[2];
    arimaj.bigp = mrx[3];
    arimaj.bigd = mrx[4];
    arimaj.bigq = mrx[5];
    arimaj.s = mrx[6];

    nser = 1;

    if (nparx > 0) {
      /* Allocate array parx */
      if (!(parx = NAG_ALLOC(nparx + nser, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      for (i = 1; i <= nparx; ++i)
        scanf("%lf", &parx[i - 1]);
      scanf("%*[^\n] ");

      /* Read model by which to filter series */
      for (i = 1; i <= 7; ++i)
        scanf("%" NAG_IFMT "", &mr[i - 1]);
      scanf("%*[^\n] ");

      arimaf.p = mr[0];
      arimaf.d = mr[1];
      arimaf.q = mr[2];
      arimaf.bigp = mr[3];
      arimaf.bigd = mr[4];
      arimaf.bigq = mr[5];
      arimaf.s = mr[6];

      npar = mr[0] + mr[2] + mr[3] + mr[5];
      if (npar > 0) {
        /* Allocate array par */
        if (!(par = NAG_ALLOC(npar + nparx, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
        for (i = 1; i <= npar; ++i)
          scanf("%lf", &par[i - 1]);
        scanf("%*[^\n] ");

        /* Initially backforecast QY values */
        /* (1) Reverse series in situ */
        n = nx / 2;
        ni = nx;
        for (i = 1; i <= n; ++i) {
          a1 = x[i - 1];
          a2 = x[ni - 1];
          x[i - 1] = a2;
          x[ni - 1] = a1;
          --ni;
        }
        idd = mrx[1] + mrx[4];
        /* (2) Possible sign reversal for ARIMA constant */
        if (idd % 2 != 0)
          cx = -cx;

        /* (3) Calculate number of backforecasts required */
        iqxd = mrx[2] + mrx[5] * mrx[6];

        /* Calculate series length */
        ny = nx + iqxd;

        /* Allocate array y */
        if (!(y = NAG_ALLOC(ny, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }

        if (iqxd != 0) {
          /* Allocate arrays fsd, fva and st. */
          if (!(fsd = NAG_ALLOC(iqxd, double)) ||
              !(fva = NAG_ALLOC(iqxd, double)))
          {
            printf("Allocation failure\n");
            exit_status = -1;
            goto END;
          }
          /* (4) Set up parameter list for call to forecast
           * routine g13bjc
           */
          npara = nparx + nser;
          parx[npara - 1] = cx;
          tdxxy = nser;
          tdmrx = nser - 1;
          ldparx = nser - 1;
          tdparx = nser - 1;
          if (!(rms = NAG_ALLOC(nser, double)) ||
              !(parxx = NAG_ALLOC(nser, double)) ||
              !(mrxx = NAG_ALLOC(7 * nser, Integer)))
          {
            printf("Allocation failure\n");
            exit_status = -1;
            goto END;
          }

          /* nag_tsa_transf_orders (g13byc).
           * Allocates memory to transfer function model orders
           */
          nag_tsa_transf_orders(nser, &transfj, &fail);
          if (fail.code != NE_NOERROR) {
            printf("Error from nag_tsa_transf_orders (g13byc)"
                   ".\n%s\n", fail.message);
            exit_status = 1;
            goto END;
          }

          rms[0] = 0;
          transfj.nag_b = 0;
          transfj.nag_q = 0;
          transfj.nag_p = 0;
          transfj.nag_r = 1;
          for (i = 1; i <= 7; ++i)
            mrxx[i - 1] = 0;
          parxx[0] = 0;

          /* Tell nag_tsa_multi_inp_model_forecast (g13bjc) not to
           * print parameters on entry */
          options.list = Nag_FALSE;

          /* nag_tsa_multi_inp_model_forecast (g13bjc).
           * Forecasting function
           */
          nag_tsa_multi_inp_model_forecast(&arimaj, nser, &transfj,
                                           parx, npara, nx, iqxd, x,
                                           tdxxy, rms, mrxx, tdmrx,
                                           parxx, ldparx, tdparx,
                                           fva, fsd, &options, &fail);
          if (fail.code != NE_NOERROR) {
            printf("Error from nag_tsa_multi_inp_model_forecast "
                   "(g13bjc).\n%s\n", fail.message);
            exit_status = 1;
            goto END;
          }

          j = iqxd;
          for (i = 1; i <= iqxd; ++i) {
            y[i - 1] = fva[j - 1];
            --j;
          }

          /* Move series into y */
          j = iqxd + 1;
          k = nx;
          for (i = 1; i <= nx; ++i) {
            if (j > 305)
              goto END;
            y[j - 1] = x[k - 1];
            ++j;
            --k;
          }
        }

        /* Move ARIMA for series into mr */
        for (i = 1; i <= 7; ++i)
          mr[i + 6] = mrx[i - 1];

        arimav.p = mr[7];
        arimav.d = mr[8];
        arimav.q = mr[9];
        arimav.bigp = mr[10];
        arimav.bigd = mr[11];
        arimav.bigq = mr[12];
        arimav.s = mr[13];

        /* Move parameters of ARIMA for y into par */
        for (i = 1; i <= nparx; ++i)
          par[npar + i - 1] = parx[i - 1];
        npar += nparx;

        /* Move constant and reset sign reversal */
        cy = cx;
        if (idd % 2 != 0)
          cy = -cy;

        /* Set parameters for call to filter routine
         * nag_tsa_arma_filter (g13bac) */
        nb = ny + MAX(mr[2] + mr[5] * mr[6],
                      mr[0] + mr[1] + (mr[3] + mr[4]) * mr[6]);

        /* Allocate array b */
        if (!(b = NAG_ALLOC(nb, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }

        /* Filter series by call to nag_tsa_arma_filter (g13bac) */
        /* nag_tsa_arma_filter (g13bac).
         * Multivariate time series, filtering (pre-whitening) by an
         * ARIMA model
         */
        nag_tsa_arma_filter(y, ny, &arimaf, &arimav, par, npar, cy, b,
                            nb, &fail);
        if (fail.code != NE_NOERROR) {
          printf("Error from nag_tsa_arma_filter (g13bac).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

        printf("\n");
        printf("                 Original        Filtered\n");
        printf("Backforecasts    y-series         series\n");
        if (iqxd != 0) {
          ij = -iqxd;
          for (i = 1; i <= iqxd; ++i) {
            printf("%8" NAG_IFMT "%17.4f%15.4f\n", ij, y[i - 1], b[i - 1]);
            ++ij;
          }
        }

        printf("\n");
        printf("       Filtered        Filtered        "
               "Filtered        Filtered\n");
        printf("        series          series  "
               "        series          series\n");
        for (i = iqxd + 1; i <= ny; i += 4) {
          for (ii = i; ii <= MIN(ny, i + 3); ++ii) {
            printf("%5" NAG_IFMT "", ii - iqxd);
            printf("%9.4f  ", b[ii - 1]);
          }
          printf("\n");
        }
      }
    }
  }

END:

  /* Free the options structure used by nag_tsa_multi_inp_model_forecast
   * (g13bjc) */
  /* nag_tsa_free (g13xzc).
   * Freeing function for use with g13 option setting
   */
  nag_tsa_free(&options);

  NAG_FREE(b);
  NAG_FREE(fsd);
  NAG_FREE(fva);
  NAG_FREE(par);
  NAG_FREE(parx);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(rms);
  NAG_FREE(parxx);
  NAG_FREE(mrxx);

  return exit_status;
}