/* nag_tsa_arma_filter (g13bac) Example Program.
 *
 * Copyright 2002 Numerical Algorithms Group.
 *
 * Mark 7, 2002.
 */

#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, nmr, 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;

  /* Initialise the options structure used by g13bjc */
  g13bxc(&options);

  Vprintf("g13bac Example Program Results\n");

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

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

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

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

      Vscanf("%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)))
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

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

          /* Read model by which to filter series */
          for (i = 1; i <= 7; ++i)
            Vscanf("%ld", &mr[i-1]);
          Vscanf("%*[^\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)))
                {
                  Vprintf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }
              for (i = 1; i <= npar; ++i)
                Vscanf("%lf", &par[i-1]);
              Vscanf("%*[^\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)))
                {
                  Vprintf("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)) )
                    {
                      Vprintf("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)) )
                  {
                    Vprintf("Allocation failure\n");
                    exit_status = -1;
                    goto END;
                  }

                  g13byc(nser, &transfj, &fail);
                  if (fail.code != NE_NOERROR)
                    {
                      Vprintf("Error from 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 g13bjc not to print parameters on entry */
                  options.list = FALSE;

                  g13bjc(&arimaj, nser, &transfj, parx, npara, nx, iqxd, x,
                          tdxxy, rms, mrxx, tdmrx, parxx, ldparx, tdparx,
                          fva, fsd, &options, &fail);
                  if (fail.code != NE_NOERROR)
                    {
                      Vprintf("Error from 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 g13bac */
              nmr = 14;
              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)) )
                {
                  Vprintf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }

              /* Filter series by call to g13bac */
              g13bac(y, ny, &arimaf, &arimav, par, npar, cy, b,
                      nb, &fail);
              if (fail.code != NE_NOERROR)
                {
                  Vprintf("Error from g13bac.\n%s\n", fail.message);
                  exit_status = 1;
                  goto END;
                }

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

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

 END:

  /* Free the options structure used by g13bjc */
  g13xzc(&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;
}