/* nag_tsa_multi_inp_model_forecast(g13bjc) Example Program
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 */

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

static void ex1(void);
static void ex2(void);

#define  NSERMX  6
#define  NPMAX  10
#define  LDPARX  8
#define  NFVMAX 10
#define  NEVMAX 40
#define  LDXXY   NEVMAX + NFVMAX

#define  TDMRX    NSERMX
#define  TDPARX   NSERMX
#define  TDXXY    NSERMX

int main(void)
{
  /* Two examples are called, ex1() which uses the
   * default settings to solve the problem and
   * ex2() which solves the same problem with
   * some optional parameters set by the user.
   */

  Vprintf("g13bjc Example Program Results\n");
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  ex1();
  ex2();
  return EXIT_SUCCESS;
}

static void ex1(void)
{
  Integer  i, j, n, nev, nfv, npara, nseries, inser;
  double  fsd[NFVMAX], fva[NFVMAX], para[NFVMAX], parx[LDPARX][NSERMX],
  rmsxy[NSERMX], xxy[LDXXY][NSERMX];
  Integer mrx[7][NSERMX];
  Nag_ArimaOrder  arimav;
  Nag_TransfOrder transfv;
  static NagError fail;

  Vprintf("\ng13bjc example 1: no option setting.\n\n");

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

  Vscanf("%ld%ld%ld", &nev, &nfv, &nseries);
  if (nseries>0 && nseries<=NSERMX && nev>0 && nev<=NEVMAX &&
      nfv>0 && nfv<=NFVMAX)
    {
      /*
       * Allocate memory to the arrays in structure transfv containing
       * the transfer function model orders of the input series.
       */
      g13byc(nseries, &transfv, NAGERR_DEFAULT);

      /*
       * Read the orders vector of the ARIMA model for the output noise
       * component into structure arimav.
       */
      Vscanf("%ld%ld%ld%ld%ld%ld%ld", &arimav.p, &arimav.d, &arimav.q,
             &arimav.bigp, &arimav.bigd, &arimav.bigq, &arimav.s);
      /*
       * Read the transfer function model orders of the input series into
       * structure transfv.
       */
      inser = nseries - 1;

      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.b[j]);
      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.q[j]);
      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.p[j]);
      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.r[j]);

      npara = 0;
      for (i=0; i<inser; ++i)
        npara = npara + transfv.q[i] + transfv.p[i];
      npara = npara + arimav.p + arimav.q + arimav.bigp + arimav.bigq
        + nseries;
      if (npara<=NPMAX)
        {
          for (i=0; i<npara; ++i)
            Vscanf("%lf", &para[i]);
          n = nev + nfv;
          for (i=0; i<n; ++i)
            for (j=0; j<nseries; ++j)
              Vscanf("%lf", &xxy[i][j]);
          for (i=0; i<nseries; ++i)
            Vscanf("%lf", &rmsxy[i]);
          for (i=0; i<7; ++i)
            for (j=0; j<inser; ++j)
              Vscanf("%ld", &mrx[i][j]);
          for (i=0; i<5; ++i)
            for (j=0; j<inser; ++j)
              Vscanf("%lf", &parx[i][j]);

          fail.print = TRUE;

          g13bjc(&arimav, nseries, &transfv, para, npara, nev, nfv,
                 &xxy[0][0], (Integer)TDXXY, rmsxy, &mrx[0][0],
                 (Integer)TDMRX, &parx[0][0], (Integer)LDPARX,
                 (Integer)TDPARX, fva, fsd, G13_DEFAULT, &fail);

          if (fail.code==NE_NOERROR || fail.code==NE_SOLUTION_FAIL_CONV ||
              fail.code==NE_MAT_NOT_POS_DEF)
            {
              Vprintf("\nThe residual mean square for the output\n");
              Vprintf("series is also derived and its value is %10.4f\n\n",
                      rmsxy[nseries-1]);
              Vprintf("The forecast values and their standard errors are\n\n");
              Vprintf("\n   i      fva      fsd\n\n");
              for (i=0; i<nfv; ++i)
                Vprintf("%4ld%10.3f%10.4f\n", i+1, fva[i], fsd[i]);
            }
        }
      else
        {
          Vfprintf(stderr, "npara is out of range: npara = %-3ld\n", npara);
          g13bzc(&transfv);
          exit(EXIT_FAILURE);
        }
    }
  else
    {
      Vfprintf(stderr, "One or more of nseries, nev and nfv are out of "
               "range: nseries = %-3ld,  nev = %-3ld while  nfv = %-3ld\n", nseries, nev, nfv);

      exit(EXIT_FAILURE);
    }
  g13bzc(&transfv);
  if (fail.code!=NE_NOERROR)
    exit(EXIT_FAILURE);
}

static void ex2(void)
{
  Integer  i, j, n, nev, nfv, npara, nseries, inser;
  double  fsd[NFVMAX], fva[NFVMAX], para[NFVMAX], parx[LDPARX][TDPARX],
  rmsxy[NSERMX], xxy[LDXXY][TDXXY];
  Integer mrx[7][NSERMX];
  Nag_ArimaOrder  arimav;
  Nag_TransfOrder transfv;
  Nag_G13_Opt  options;
  static NagError fail;

#define ZT(I,J)       options.zt[(J)+(I) * options.tdzt]

  Vprintf("\n\ng13bjc example 2: using option setting.\n\n");

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

  /*
   * Initialise the option-setting function.
   */
  g13bxc(&options);

  Vscanf("%ld%ld%ld", &nev, &nfv, &nseries);
  if (nseries>0 && nseries<=NSERMX && nev>0 && nev<=NEVMAX &&
      nfv>0 && nfv<=NFVMAX)
    {
      /*
       * Set option variable to the desired value.
       */

      options.cfixed = TRUE;
      /*
       * Allocate memory to the arrays in structure transfv containing
       * the transfer function model orders of the input series.
       */
      g13byc(nseries, &transfv, NAGERR_DEFAULT);

      /*
       * Read the orders vector of the ARIMA model for the output noise
       * component into structure arimav.
       */
      Vscanf("%ld%ld%ld%ld%ld%ld%ld", &arimav.p, &arimav.d, &arimav.q,
             &arimav.bigp, &arimav.bigd, &arimav.bigq, &arimav.s);
      /*
       * Read the transfer function model orders of the input series into
       * structure transfv.
       */
      inser = nseries - 1;

      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.b[j]);
      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.q[j]);
      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.p[j]);
      for (j=0; j<inser; ++j)
        Vscanf("%ld", &transfv.r[j]);

      npara = 0;
      for (i=0; i<inser; ++i)
        npara = npara + transfv.q[i] + transfv.p[i];
      npara = npara + arimav.p + arimav.q + arimav.bigp + arimav.bigq
        + nseries;
      if (npara<=NPMAX)
        {
          for (i=0; i<npara; ++i)
            Vscanf("%lf", &para[i]);
          n = nev + nfv;
          for (i=0; i<n; ++i)
            for (j=0; j<nseries; ++j)
              Vscanf("%lf", &xxy[i][j]);
          for (i=0; i<nseries; ++i)
            Vscanf("%lf", &rmsxy[i]);
          for (i=0; i<7; ++i)
            for (j=0; j<inser; ++j)
              Vscanf("%ld", &mrx[i][j]);
          for (i=0; i<5; ++i)
            for (j=0; j<inser; ++j)
              Vscanf("%lf", &parx[i][j]);

          fail.print = TRUE;

          g13bjc(&arimav, nseries, &transfv, para, npara, nev, nfv,
                 &xxy[0][0], (Integer)TDXXY, rmsxy, &mrx[0][0],
                 (Integer)TDMRX, &parx[0][0], (Integer)LDPARX,
                 (Integer)TDPARX, fva, fsd, &options, &fail);

          if (fail.code==NE_NOERROR || fail.code==NE_SOLUTION_FAIL_CONV ||
              fail.code==NE_MAT_NOT_POS_DEF)
            {
              Vprintf("%1ld sets of observations were processed.\n",nev);
              Vprintf("\nThe residual mean square for the output ");
              Vprintf("series is %10.4f\n\n", rmsxy[nseries-1]);
              Vprintf("The forecast values and their standard errors are\n\n");
              Vprintf("\n   i      fva      fsd\n\n");
              for (i=0; i<nfv; ++i)
                Vprintf("%4ld%10.3f%10.4f\n", i+1, fva[i], fsd[i]);
              Vprintf("\nThe values of z(t) and noise(t) are\n\n");
              Vprintf("   i       z1         z2         z3         z4"
                      "         z5       noise\n\n");
              for (i=0; i<n; ++i)
                {
                  Vprintf("%4ld", i+1);
                  for (j=0; j<nseries-1; ++j)
                    Vprintf("%10.3f ", ZT(i,j));
                  Vprintf("%10.3f\n", options.noise[i]);
                }
            }
        }
      else
        {
          Vfprintf(stderr, "npara is out of range: npara = %-3ld\n", npara);
          g13xzc(&options);
          g13bzc(&transfv);
          exit(EXIT_FAILURE);
        }
    }
  else
    {
      Vfprintf(stderr, "One or more of nseries, nev and nfv are out of range:"
               " nseries = %-3ld,  nev = %-3ld while  nfv = %-3ld\n",
               nseries, nev, nfv);
      exit(EXIT_FAILURE);
    }
  g13xzc(&options);
  g13bzc(&transfv);
  if (fail.code!=NE_NOERROR)
    exit(EXIT_FAILURE);
}