/* nag_tsa_multi_inp_model_estim(g13bec) 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  2
#define  NPMAX  10
#define  NXXYMX 50
#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("g13bec Example Program Results\n");
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  ex1();
  ex2();
  return EXIT_SUCCESS;
}

static void ex1(void)
{
  double   df, objf, rss;
  Integer  i, j, npara, nseries, nxxy, inser;
  double   para[NPMAX], sd[NPMAX], xxy[NXXYMX][NSERMX];
  Nag_ArimaOrder  arimav;
  Nag_TransfOrder transfv;
  static NagError fail;

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

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

  Vscanf("%ld%ld", &nxxy, &nseries);
  if (nxxy>0 && nxxy<=NXXYMX && nseries>0 && nseries<=NSERMX)
    {
      /*
       * 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]);
          for (i=0; i<nxxy; ++i)
            for (j=0; j<nseries; ++j)
              Vscanf("%lf", &xxy[i][j]);

          fail.print = TRUE;

          g13bec(&arimav, nseries, &transfv, para, npara, nxxy, &xxy[0][0],
                 (Integer)TDXXY, sd, &rss, &objf, &df, G13_DEFAULT,
                 &fail);

        }
      else
        {
          Vfprintf(stderr, "npara is out of range: npara = %-3ld\n", npara);
          g13bzc(&transfv);
          exit(EXIT_FAILURE);
        }
    }
  else
    {
      Vfprintf(stderr, "One or both of nxxy and nseries are out of range:"
               " nxxy = %-3ld  while  nseries = %-3ld\n", nxxy, nseries);
      exit(EXIT_FAILURE);
    }
  g13bzc(&transfv);
  if (fail.code!=NE_NOERROR) exit(EXIT_FAILURE);
}


static void ex2(void)
{
  double   df, objf, rss;
  Integer  i, j, npara, nseries, inser, nxxy;
  double   para[NPMAX], sd[NPMAX], xxy[NXXYMX][NSERMX];
  Nag_ArimaOrder  arimav;
  Nag_TransfOrder transfv;
  Nag_G13_Opt  options;
  static NagError fail;

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

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

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

  /*
   * Initialise the option structure.
   */
  g13bxc(&options);

  Vscanf("%ld%ld%ld", &nxxy, &nseries, &options.max_iter);
  if (nxxy>0 && nxxy<=NXXYMX && nseries>0 && nseries<=NSERMX)
    {
      /*
       * Set some specific option variables to the desired values.
       */

      options.criteria = Nag_Marginal;
      options.print_level = Nag_Soln_Iter_Full;
      /*
       * 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]);
          for (i=0; i<nxxy; ++i)
            for (j=0; j<nseries; ++j)
              Vscanf("%lf", &xxy[i][j]);

          fail.print = TRUE;

          g13bec(&arimav, nseries, &transfv, para, npara, nxxy, &xxy[0][0],
                 (Integer)TDXXY, sd, &rss, &objf, &df, &options, &fail);

          Vprintf("\nThe correlation matrix is \n\n");
          for (i=0; i<npara; ++i)
            for (j=0; j<npara; ++j)
              Vprintf("%10.4f%c", CM(i,j), (j%5==4) ? '\n' : ' ');
          Vprintf("\nThe residuals and the z and n values are\n\n");
          Vprintf("   i         res[i]          z(t)        noise(t)\n\n");
          for (i=0; i<nxxy; ++i)
            {
              if (i+1<=options.lenres)
                {
                  Vprintf("%4ld%15.3f", i+1, options.res[i]);
                  for (j=0; j<nseries-1; ++j)
                    Vprintf("%15.3f ", ZT(i,j));
                  Vprintf("%15.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 both of nxxy and nseries are out of range:"
               " nxxy = %-3ld  while  nseries = %-3ld\n", nxxy, nseries);
      exit(EXIT_FAILURE);
    }
  g13bzc(&transfv);
  g13xzc(&options);
  if (fail.code!=NE_NOERROR) exit(EXIT_FAILURE);
}