/* nag_tsa_multi_auto_corr_part (g13dbc) Example Program.

 *

 * Copyright 2002 Numerical Algorithms Group.

 *

 * Mark 7, 2002.

 * Mark 7b revised, 2004.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagg13.h>



int main(void)

{

  /* Scalars */

  double v0;

  Integer exit_status, i1, i, j, j1, k, nk, nl, ns, nvp,

    pdc0, pddb;

  NagError fail;



  /* Arrays */

  double *c0 = 0, *c = 0, *d = 0, *db = 0, *p = 0, *v = 0, *w = 0,

    *wb = 0;



#define C(I,J,K) c[((K-1)*ns + (J-1))*ns + I - 1]

#define D(I,J,K) d[((K-1)*ns + (J-1))*ns + I - 1]

#define W(I,J,K) w[((K-1)*ns + (J-1))*ns + I - 1]

#define WB(I,J,K) wb[((K-1)*ns + (J-1))*ns + I - 1]



#ifdef NAG_COLUMN_MAJOR

#define C0(I,J) c0[(J-1)*pdc0 + I - 1]

#define DB(I,J) db[(J-1)*pddb + I - 1]

#else

#define C0(I,J) c0[(I-1)*pdc0 + J - 1]

#define DB(I,J) db[(I-1)*pddb + J - 1]

#endif



  INIT_FAIL(fail);

  exit_status = 0;



  Vprintf("nag_tsa_multi_auto_corr_part (g13dbc) Example Program Results\n");



  /* Skip heading in data file */

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



  /* Read series length, and numbers of lags */

  Vscanf("%ld%ld%ld%*[^\n] ", &ns, &nl, &nk);



  if (ns > 0 && nl > 0 && nk > 0)

    {

      /* Allocate arrays */

      if ( !(c0 = NAG_ALLOC(ns * ns, double)) ||

           !(c = NAG_ALLOC(ns * ns * nl, double)) ||

           !(d = NAG_ALLOC(ns * ns * nk, double)) ||

           !(db = NAG_ALLOC(ns * ns, double)) ||

           !(p = NAG_ALLOC(nk, double)) ||

           !(v = NAG_ALLOC(nk, double)) ||

           !(w = NAG_ALLOC(ns * ns * nk, double)) ||

           !(wb = NAG_ALLOC(ns * ns * nk, double)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }



      pdc0 = ns;

      pddb = ns;



      /* Read autocovariances */

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

        {

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

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

        }

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



      for (k = 1; k <= nl; ++k)

        {

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

            {

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

                Vscanf("%lf", &C(i,j,k));

            }

        }

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



      /* Call routine to calculate multivariate partial

         autocorrelation function */



      /* nag_tsa_multi_auto_corr_part (g13dbc).

       * Multivariate time series, multiple squared partial

       * autocorrelations

       */

      nag_tsa_multi_auto_corr_part(c0, c, ns, nl, nk, p, &v0, v, d, db, w, wb,

                                   &nvp, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_tsa_multi_auto_corr_part (g13dbc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      if (fail.code == NE_NOERROR || fail.code == NE_POS_DEF)

        {

          Vprintf("\n");

          Vprintf("Number of valid parameters =%10ld\n", nvp);



          Vprintf("\n");

          Vprintf("Multivariate partial autocorrelations\n");



          for (i1 = 1; i1 <= nk; ++i1)

            {

              Vprintf("%13.5f", p[i1-1]);

              if (i1 % 5 == 0 || i1 == nk)

                Vprintf("\n");

            }



          Vprintf("\n");

          Vprintf("Zero lag predictor error variance determinant\n");

          Vprintf("followed by error variance ratios\n");

          Vprintf("%12.5f", v0);



          for (i1 = 1; i1 <= nk; ++i1)

            {

              Vprintf("%13.5f", v[i1-1]);

              if (i1 % 5 == 0 || i1 == nk)

                Vprintf("\n");

            }



          Vprintf("\n");

          Vprintf("Prediction error variances\n");

          Vprintf("\n");



          for (k = 1; k <= nk; ++k)

            {

              Vprintf("Lag =%5ld\n", k);

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

                {

                  for (j1 = 1; j1 <= ns; ++j1)

                    {

                      Vprintf("%13.5f", D(i,j1,k));

                      if (j1 % 5 == 0 || j1 == ns)

                        Vprintf("\n");

                    }

                }

              Vprintf("\n");

            }



          Vprintf("Last backward prediction error variances\n");

          Vprintf("\n");

          Vprintf("Lag =%5ld\n", nvp);



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

            {

              for (j1 = 1; j1 <= ns; ++j1)

                {

                  Vprintf("%13.5f", DB(i,j1));

                  if (j1 % 5 == 0 || j1 == ns)

                    Vprintf("\n");

                }

            }



          Vprintf("\n");

          Vprintf("Prediction coefficients\n");

          Vprintf("\n");



          for (k = 1; k <= nk; ++k)

            {

              Vprintf("Lag =%5ld\n", k);

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

                {

                  for (j1 = 1; j1 <= ns; ++j1)

                    {

                      Vprintf("%13.5f", W(i,j1,k));

                      if (j1 % 5 == 0 || j1 == ns)

                        Vprintf("\n");

                    }

                }

              Vprintf("\n");

            }

          Vprintf("Backward prediction coefficients\n");

          Vprintf("\n");



          for (k = 1; k <= nk; ++k)

            {

              Vprintf("Lag =%5ld\n", k);

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

                {

                  for (j1 = 1; j1 <= ns; ++j1)

                    {

                      Vprintf("%13.5f", WB(i,j1, k));

                      if (j1 % 5 == 0 || j1 == ns)

                        Vprintf("\n");

                    }

                }

              Vprintf("\n");

            }

        }

    }



 END:

  if (c0) NAG_FREE(c0);

  if (c) NAG_FREE(c);

  if (d) NAG_FREE(d);

  if (db) NAG_FREE(db);

  if (p) NAG_FREE(p);

  if (v) NAG_FREE(v);

  if (w) NAG_FREE(w);

  if (wb) NAG_FREE(wb);



  return exit_status;

}