/* nag_tsa_transf_filter (g13bbc) Example Program.
 *
 * Copyright 2002 Numerical Algorithms Group.
 *
 * Mark 7, 2002.
 * Mark 7b revised, 2004.
 */

#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;
  Integer         nser, npara, tdxxy, tdmrx, ldparx, tdparx;
  Integer         exit_status = 0;

  /* Arrays */
  double          *b = 0, *fsd = 0, *fva = 0, *par = 0, *parx = 0;
  double          *x = 0, *y = 0, *rms = 0, *parxx = 0;
  Integer         mr[10], mrx[7], *mrxx = 0;

  Nag_TransfOrder transfj, transfv;
  Nag_ArimaOrder  arimaj, arimas;
  Nag_G13_Opt     options;
  NagError        fail;

  INIT_FAIL(fail);

  exit_status = 0;

  /* Initialise 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_transf_filter (g13bbc) Example Program Results\n");

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

  scanf("%ld%*[^\n] ", &nx);

  printf("\n");
  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("%ld", &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 <= 3; ++i)
            scanf("%ld", &mr[i-1]);
          scanf("%*[^\n] ");

          transfv.nag_b = mr[0];
          transfv.nag_q = mr[1];
          transfv.nag_p = mr[2];

          npar = mr[1] + mr[2] + 1;
          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];
              if (iqxd != 0)
                {
                  if (!(fsd = NAG_ALLOC(iqxd, double)) ||
                      !(fva = NAG_ALLOC(iqxd, double)))
                    {
                      printf("Allocation failure\n");
                      exit_status = -1;
                      goto END;
                    }
                  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;
                    }
                }

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

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

              /* Move backforecasts to start of y array */
              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 > 215)
                    goto END;
                  y[j-1] = x[k-1];
                  ++j;
                  --k;
                }
            }

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

          arimas.p = mr[3];
          arimas.d = mr[4];
          arimas.q = mr[5];
          arimas.bigp = mr[6];
          arimas.bigd = mr[7];
          arimas.bigq = mr[8];
          arimas.s = mr[9];

          /* 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_transf_filter (g13bbc) */
          nb = ny + MAX(mr[0] + mr[1], mr[2]);

          /* 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_transf_filter (g13bbc) */
          /* nag_tsa_transf_filter (g13bbc).
           * Multivariate time series, filtering by a transfer
           * function model
           */
          nag_tsa_transf_filter(y, ny, &transfv, &arimas, par, npar, cy, b, nb,
                                &fail);
          if (fail.code != NE_NOERROR)
            {
              printf(
                      "Error from nag_tsa_transf_filter (g13bbc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          printf("                  Original        Filtered\n");
          printf(" Backforecasts    y-series         series\n");
          if (iqxd != 0)
            {
              ij = -iqxd;
              for (i = 1; i <= iqxd; ++i)
                {
                  printf("%8ld%17.1f%16.1f\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("%5ld", ii-iqxd);
                      printf("%10.1f", 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);

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

  return exit_status;
}