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

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

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

int main(void)
{
  Integer exit_status_ex1=0;
  Integer exit_status_ex2=0;

  /* 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("nag_tsa_multi_inp_model_forecast (g13bjc) Example Program "
          "Results\n");
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();

  return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1;
}

#define PARX(I,J) parx[(I)*tdparx + J]
#define XXY(I,J) xxy[(I)*tdxxy + J]
#define MRX(I,J) mrx[(I)*tdmrx + J]
static int ex1(void)
{
  Integer exit_status=0, i, inser, j, ldparx, *mrx=0, n, nev, nfv, npara;
  Integer nseries, tdmrx, tdparx, tdxxy;
  NagError fail;
  Nag_ArimaOrder arimav;
  Nag_TransfOrder transfv;
  double *fsd=0, *fva=0, *para=0, *parx=0, *rmsxy=0, *xxy=0;

  INIT_FAIL(fail);
  Vprintf("\nnag_tsa_multi_inp_model_forecast (g13bjc) 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 >=1 && nev>0 && nfv>0)
    {
      /*
       * 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(nseries, &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%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;
      ldparx = 8;
      if (npara>=1)
        {
          if ( !( fsd = NAG_ALLOC(nfv, double)) ||
               !( fva = NAG_ALLOC(nfv, double)) ||
               !( para = NAG_ALLOC(npara, double)) ||
               !( parx = NAG_ALLOC(ldparx*(nseries-1), double)) ||
               !( rmsxy = NAG_ALLOC(nseries, double)) ||
               !( xxy = NAG_ALLOC((nev+nfv)*nseries, double))  ||
               !( mrx = NAG_ALLOC(7*(nseries-1), Integer)) )
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          tdmrx =  nseries-1;
          tdparx = nseries-1;
          tdxxy =  nseries;
          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));

          /* nag_tsa_multi_inp_model_forecast (g13bjc).
           * Forecasting function
           */
          nag_tsa_multi_inp_model_forecast(&arimav, nseries, &transfv, para,
                                           npara, nev, nfv, xxy, tdxxy, rmsxy,
                                           mrx, tdmrx, parx, nev+nfv, tdparx,
                                           fva, fsd, G13_DEFAULT, &fail);
          if (fail.code != NE_NOERROR)
            {
              Vprintf("Error from nag_tsa_multi_inp_model_forecast (g13bjc)."
                      "\n%s\n", fail.message);
              exit_status = 1;
              goto END;
            }

          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
        {
          Vprintf("npara is out of range: npara = %-3ld\n", npara);
          /* nag_tsa_trans_free (g13bzc).
           * Freeing function for the structure holding the transfer
           * function model orders
           */
          nag_tsa_trans_free(&transfv);
          exit_status = 1;
          goto END;
        }
    }
  else
    {
      Vprintf("One or more of nseries, nev and nfv are out of "
              "range: nseries = %-3ld,  nev = %-3ld while  "
              "nfv = %-3ld\n", nseries, nev, nfv);
      exit_status = 1;
      goto END;
    }
  /* nag_tsa_trans_free (g13bzc), see above. */
  nag_tsa_trans_free(&transfv);
 END:
  if (fsd) NAG_FREE(fsd);
  if (fva) NAG_FREE(fva);
  if (para) NAG_FREE(para);
  if (parx) NAG_FREE(parx);
  if (rmsxy) NAG_FREE(rmsxy);
  if (xxy) NAG_FREE(xxy);
  if (mrx) NAG_FREE(mrx);
  return exit_status;
}

static int ex2(void)
{
  Integer exit_status=0, i, inser, j, ldparx, *mrx=0, n, nev, nfv, npara;
  Integer nseries, tdmrx, tdparx, tdxxy;
  NagError fail;
  Nag_ArimaOrder arimav;
  Nag_G13_Opt options;
  Nag_TransfOrder transfv;
  double *fsd=0, *fva=0, *para=0, *parx=0, *rmsxy=0, *xxy=0;

  INIT_FAIL(fail);


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

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

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

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

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

      options.cfixed = Nag_TRUE;
      /*
       * Allocate memory to the arrays in structure transfv containing
       * the transfer function model orders of the input series.
       */
      /* nag_tsa_transf_orders (g13byc), see above. */
      nag_tsa_transf_orders(nseries, &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%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;
      ldparx = 8;
      if (npara>=1)
        {

          if ( !( fsd = NAG_ALLOC(nfv, double)) ||
               !( fva = NAG_ALLOC(nfv, double)) ||
               !( para = NAG_ALLOC(npara, double)) ||
               !( parx = NAG_ALLOC(ldparx*(nseries-1), double)) ||
               !( rmsxy = NAG_ALLOC(nseries, double)) ||
               !( xxy = NAG_ALLOC((nev+nfv)*(nseries), double))  ||
               !( mrx = NAG_ALLOC(7*(nseries-1), Integer)) )
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          tdmrx = nseries-1;
          tdparx = nseries-1;
          tdxxy = nseries;

          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));

          /* nag_tsa_multi_inp_model_forecast (g13bjc), see above. */
          nag_tsa_multi_inp_model_forecast(&arimav, nseries, &transfv, para,
                                           npara, nev, nfv, xxy, tdxxy, rmsxy,
                                           mrx, tdmrx, parx, ldparx,
                                           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
            {
              Vprintf("Error from nag_tsa_multi_inp_model_forecast (g13bjc)."
                      "\n%s\n", fail.message);
              exit_status = 1;
              goto END;
            }
        }
      else
        {
          Vprintf("npara is out of range: npara = %-3ld\n", npara);
          /* nag_tsa_free (g13xzc).
           * Freeing function for use with g13 option setting
           */
          nag_tsa_free(&options);
          /* nag_tsa_trans_free (g13bzc), see above. */
          nag_tsa_trans_free(&transfv);
          exit_status = 1;
          goto END;

        }
    }
  else
    {
      Vprintf("One or more of nseries, nev and nfv are out of range:"
              " nseries = %-3ld,  nev = %-3ld while  "
              "nfv = %-3ld\n", nseries, nev, nfv);
      exit_status = 1;
      goto END;

    }
  /* nag_tsa_free (g13xzc), see above. */
  nag_tsa_free(&options);
  /* nag_tsa_trans_free (g13bzc), see above. */
  nag_tsa_trans_free(&transfv);
 END:
  if (fsd) NAG_FREE(fsd);
  if (fva) NAG_FREE(fva);
  if (para) NAG_FREE(para);
  if (parx) NAG_FREE(parx);
  if (rmsxy) NAG_FREE(rmsxy);
  if (xxy) NAG_FREE(xxy);
  if (mrx) NAG_FREE(mrx);
  return exit_status;
}