/* nag_tsa_multi_inp_update (g13bgc) Example Program.

 *

 * Copyright 2004 Numerical Algorithms Group.

 *

 * Mark 8, 2004.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagg13.h>



int main (void)

{

  double df, objf, rss;

  Integer exit_status, i, inser, j, kzef, nnv, npara, nser, nxxy, tdxxy, tdxxyn;

  double *para = 0, *sd = 0, *xxy = 0, *xxyn = 0;

  /* Nag types */

  Nag_ArimaOrder arimav;

  Nag_TransfOrder transfv;

  Nag_G13_Opt options;

  NagError fail;



#define CM(I,J) options.cm[(I - 1) * options.tdcm + (J - 1)]

#define XXY(I,J) xxy[(I - 1) * tdxxy + J - 1]

#define XXYN(I,J) xxyn[(I - 1) * tdxxyn + J - 1]

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



  Vprintf("nag_tsa_multi_inp_update (g13bgc) Example Program Results\n");



  INIT_FAIL(fail);

  exit_status = 0;



  /* Skip heading in data file */

  Vscanf("%*[^\n]");

  /* Initialise the option structure */

  /* nag_tsa_options_init (g13bxc).

   * Initialization function for option setting

   */

  nag_tsa_options_init(&options);

  Vscanf("%ld%ld%ld%"NAG_IFMT, &nxxy, &nser,

         &options.max_iter, &nnv);



  if (nxxy > 0 && nser > 0)

    {

      /* 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.

       */

      /* nag_tsa_transf_orders (g13byc).

       * Allocates memory to transfer function model orders

       */

      nag_tsa_transf_orders(nser, &transfv, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_tsa_transf_orders (g13byc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      /*

       * Read the orders vector of the ARIMA model for the output noise

       * component into structure arimav.

       */

      Vscanf("%ld%ld%ld%ld%ld"

             "%ld%"NAG_IFMT, &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 = nser - 1;

      for (j = 1; j <= inser; ++j)

        {

          Vscanf("%"NAG_IFMT, &transfv.b[j-1]);

        }

      for (j = 1; j <= inser; ++j)

        {

          Vscanf("%"NAG_IFMT, &transfv.q[j-1]);

        }

      for (j = 1; j <= inser; ++j)

        {

          Vscanf("%"NAG_IFMT, &transfv.p[j-1]);

        }

      for (j = 1; j <= inser; ++j)

        {

          Vscanf("%"NAG_IFMT, &transfv.r[j-1]);

        }

      Vscanf("%*[^\n]");



      npara = 0;

      for (i = 1; i <= inser; ++i)

        {

          npara += transfv.q[i-1] + transfv.p[i-1];

        }

      npara += arimav.p + arimav.q + arimav.bigp + arimav.bigq + nser;



      tdxxy = nser;

      tdxxyn = nser;

      /* Memory allocation */

      if ( !(para = NAG_ALLOC(npara, double)) ||

           !(sd = NAG_ALLOC(npara, double)) ||

           !(xxy = NAG_ALLOC(nxxy * tdxxy, double)) ||

           !(xxyn = NAG_ALLOC(nnv * tdxxyn, double)) )

        {

          Vprintf("Memory allocation failure\n");

          exit_status = -1;

          goto END;

        }



      for (i = 1; i <= npara; ++i)

        {

          Vscanf("%lf", &para[i-1]);

        }

      Vscanf("%*[^\n]");



      for (i = 1; i <= nxxy; ++i)

        {

          for (j = 1; j <= nser; ++j)

            {

              Vscanf("%lf", &XXY(i,j));

            }

        }

      Vscanf("%*[^\n]");



      for (i = 1; i <= nnv; ++i)

        {

          for (j = 1; j <= nser; ++j)

            {

              Vscanf("%lf", &XXYN(i,j));

            }

        }

      Vscanf("%*[^\n]");



      options.print_level = Nag_NoPrint;

      /* nag_tsa_multi_inp_model_estim (g13bec).

       * Estimation for time series models

       */

      nag_tsa_multi_inp_model_estim(&arimav, nser, &transfv, para, npara, nxxy,

                                    xxy, tdxxy, sd, &rss, &objf, &df, &options,

                                    &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("\nError from nag_tsa_multi_inp_model_estim (g13bec).\n%s\n.",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      /* Calculate update */

      kzef = 1;

      /* nag_tsa_multi_inp_update (g13bgc).

       * Multivariate time series, update state set for

       * forecasting from multi-input model

       */

      nag_tsa_multi_inp_update(&arimav, nser, &transfv, para, npara, nnv, xxyn,

                               tdxxyn, kzef, &options, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("\nError from nag_tsa_multi_inp_update (g13bgc).\n%s\n.",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      Vprintf("\n\nThe residuals (after differencing)\n");

      for (i = 1; i <= nnv; i++)

        {

          Vprintf("%ld\t%6.4f\n", i, options.res[i-1]);

        }

      Vprintf("\n\n");



      Vprintf("\nThe values of z(t) and n(t)\n");

      for (i = 1; i <= nnv; i++)

        {

          Vprintf("%"NAG_IFMT,i);

          for (j = 1; j <= nser; j++)

            {

              Vprintf(" %6.4f", XXYN(i,j));

            }

          Vprintf("\n");

        }

      Vprintf("\n");

    }

  else

    {

      if (nxxy <= 0 || nser <= 0)

        {

          Vfprintf(stderr, "One or both of nxxy and nser are out of range: "

                   "nxxy = %-3ld while nser = %-3ld\n",

                   nxxy, nser);

          exit_status = -1;

          goto END;

        }

    }



 END:

  /* nag_tsa_trans_free (g13bzc).

   * Freeing function for the structure holding the transfer

   * function model orders

   */

  nag_tsa_trans_free(&transfv);

  /* nag_tsa_free (g13xzc).

   * Freeing function for use with g13 option setting

   */

  nag_tsa_free(&options);

  if (para) NAG_FREE (para);

  if (sd) NAG_FREE (sd);

  if (xxy) NAG_FREE (xxy);



  return exit_status;

}