/* nag_tsa_noise_spectrum_bivar (g13cgc) Example Program.

 *

 * Copyright 1996 Numerical Algorithms Group.

 *

 * Mark 4, 1996.

 * Mark 8 revised, 2004.

 *

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <naga02.h>

#include <nagg13.h>



#define LMAX 80

#define KC 8*L

#define NGMAX KC

#define L LMAX

#define NXYMAX 300



int main(void)

{



  Complex *xyg;

  Integer exit_status=0, i, is, j, kc=KC, l=L, mw, ng, nxy;

  NagError fail;

  double *er=0, erlw, erup, pw, pxy, *rf=0, rfse, *stats=0, *x=0, *xg, *y=0;

  double *yg;



  INIT_FAIL(fail);



  Vprintf("nag_tsa_noise_spectrum_bivar (g13cgc) Example Program Results\n");



  /* Skip heading in data file */

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



  Vscanf("%ld ", &nxy);

  if (nxy > 0 && nxy <= NXYMAX)

    {

      if ( !( stats = NAG_ALLOC(4, double)) ||

           !( x = NAG_ALLOC(KC, double)) ||

           !( y = NAG_ALLOC(KC, double)) ||

           !( er = NAG_ALLOC(NGMAX, double)) ||

           !( rf = NAG_ALLOC(LMAX, double)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

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

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

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

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



      /* Set parameters for call to nag_tsa_spectrum_univar (g13cbc) and g13cdc

       *  with mean correction and 10 percent taper

       */

      pxy = 0.1;

      /* Window shape parameter and zero covariance at lag 16 */

      pw = 0.5;

      mw = 16;

      /* Alignment shift of 3 */

      is = 3;



      /* Obtain univariate spectrum for the x and the y series */

      /* nag_tsa_spectrum_univar (g13cbc).

       * Univariate time series, smoothed sample spectrum using

       * spectral smoothing by the trapezium frequency (Daniell)

       * window

       */

      nag_tsa_spectrum_univar(nxy, Nag_Mean, pxy, mw, pw, l, kc, Nag_Unlogged,

                              x, &xg, &ng, stats, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_tsa_spectrum_univar (g13cbc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }

      /* nag_tsa_spectrum_univar (g13cbc), see above. */

      nag_tsa_spectrum_univar(nxy, Nag_Mean, pxy, mw, pw, l, kc, Nag_Unlogged,

                              y, &yg, &ng, stats, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_tsa_spectrum_univar (g13cbc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      /* Obtain cross spectrum of the bivariate series */

      /* nag_tsa_spectrum_bivar (g13cdc).

       * Multivariate time series, smoothed sample cross spectrum

       * using spectral smoothing by the trapezium frequency

       * (Daniell) window

       */

      nag_tsa_spectrum_bivar(nxy, Nag_Mean, pxy, mw, is, pw, l, kc, x, y, &xyg,

                             &ng, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_tsa_spectrum_bivar (g13cdc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      /* nag_tsa_noise_spectrum_bivar (g13cgc).

       * Multivariate time series, noise spectrum, bounds, impulse

       * response function and its standard error

       */

      nag_tsa_noise_spectrum_bivar(xg, yg, xyg, ng, stats, l, nxy, er, &erlw,

                                   &erup, rf, &rfse, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_tsa_noise_spectrum_bivar (g13cgc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      Vprintf("\n");

      Vprintf("           Noise spectrum\n\n");

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

        Vprintf("%6ld%16.4f\n",j - 1, er[j - 1]);



      Vprintf("\nNoise spectrum bounds multiplying factors\n\n");

      Vprintf("Lower =%10.4f", erlw);

      Vprintf("     Upper =%10.4f\n\n", erup);

      Vprintf("Impulse response function\n\n");

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

        Vprintf("%6ld%16.4f\n",j - 1, rf[j - 1]);

      Vprintf("\nImpulse response function standard error =%10.4f\n", rfse);

    }

  NAG_FREE(xg);

  NAG_FREE(yg);

  NAG_FREE(xyg);

 END:

  if (stats) NAG_FREE(stats);

  if (x) NAG_FREE(x);

  if (y) NAG_FREE(y);

  if (er) NAG_FREE(er);

  if (rf) NAG_FREE(rf);

  return exit_status;

}