NAG Library Manual, Mark 30
Interfaces:  FL   CL   CPP   AD 

NAG CL Interface Introduction
Example description
/* nag_tsa_multi_filter_transf (g13bbc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.0, 2024.
 */

#include <nag.h>
#include <stdio.h>

int main(void) {
  /* Scalars */
  double a1, a2, cx, cy;
  Integer i, ii, ij, iqxd, j, k, n, nb, ni, npar, nparx, nx;
  Integer nser, npara, tdxxy, tdmrx, ldparx, tdparx;
  Integer exit_status = 0, idd = 0, ny = 0;

  /* Arrays */
  double *b = 0, *fsd = 0, *fva = 0, *par = 0, *parx = 0;
  double *x = 0, *y = 0, *rms = 0, *parxx = 0;
  Integer mr[10], mrx[7], *mrxx = 0;

  Nag_TransfOrder transfj, transfv;
  Nag_ArimaOrder arimaj, arimas;
  Nag_G13_Opt options;
  NagError fail;

  INIT_FAIL(fail);

  exit_status = 0;

  /* Initialize the options structure used by nag_tsa_multi_inputmod_forecast
   * (g13bjc) */
  /* nag_tsa_options_init (g13bxc).
   * Initialization function for option setting
   */
  nag_tsa_options_init(&options);

  printf("nag_tsa_multi_filter_transf (g13bbc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%*[^\n] ", &nx);

  printf("\n");
  if (nx > 0) {
    /* Allocate array x */
    if (!(x = NAG_ALLOC(nx + 2, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    for (i = 1; i <= nx; ++i)
      scanf("%lf", &x[i - 1]);
    scanf("%*[^\n] ");

    /* Read univariate ARIMA for series */
    for (i = 1; i <= 7; ++i)
      scanf("%" NAG_IFMT "", &mrx[i - 1]);
    scanf("%*[^\n] ");
    scanf("%lf%*[^\n] ", &cx);

    nparx = mrx[0] + mrx[2] + mrx[3] + mrx[5];

    arimaj.p = mrx[0];
    arimaj.d = mrx[1];
    arimaj.q = mrx[2];
    arimaj.bigp = mrx[3];
    arimaj.bigd = mrx[4];
    arimaj.bigq = mrx[5];
    arimaj.s = mrx[6];

    nser = 1;

    if (nparx > 0) {
      /* Allocate array parx */
      if (!(parx = NAG_ALLOC(nparx + nser, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
      for (i = 1; i <= nparx; ++i)
        scanf("%lf", &parx[i - 1]);
      scanf("%*[^\n] ");

      /* Read model by which to filter series */
      for (i = 1; i <= 3; ++i)
        scanf("%" NAG_IFMT "", &mr[i - 1]);
      scanf("%*[^\n] ");

      transfv.nag_b = mr[0];
      transfv.nag_q = mr[1];
      transfv.nag_p = mr[2];

      npar = mr[1] + mr[2] + 1;
      if (npar > 0) {
        /* Allocate array par */
        if (!(par = NAG_ALLOC(npar + nparx, double))) {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
        for (i = 1; i <= npar; ++i)
          scanf("%lf", &par[i - 1]);
        scanf("%*[^\n] ");

        /* Initially backforecast QY values */
        /* (1) Reverse series in situ */
        n = nx / 2;
        ni = nx;
        for (i = 1; i <= n; ++i) {
          a1 = x[i - 1];
          a2 = x[ni - 1];
          x[i - 1] = a2;
          x[ni - 1] = a1;
          --ni;
        }
        idd = mrx[1] + mrx[4];
        /* (2) Possible sign reversal for ARIMA constant */
        if (idd % 2 != 0)
          cx = -cx;

        /* (3) Calculate number of backforecasts required */
        iqxd = mrx[2] + mrx[5] * mrx[6];
        if (iqxd != 0) {
          if (!(fsd = NAG_ALLOC(iqxd, double)) ||
              !(fva = NAG_ALLOC(iqxd, double))) {
            printf("Allocation failure\n");
            exit_status = -1;
            goto END;
          }
          npara = nparx + nser;
          parx[npara - 1] = cx;
          tdxxy = nser;
          tdmrx = nser - 1;
          ldparx = nser - 1;
          tdparx = nser - 1;
          if (!(rms = NAG_ALLOC(nser, double)) ||
              !(parxx = NAG_ALLOC(nser, double)) ||
              !(mrxx = NAG_ALLOC(7 * nser, Integer))) {
            printf("Allocation failure\n");
            exit_status = -1;
            goto END;
          }

          /* nag_tsa_transf_orders (g13byc).
           * Allocates memory to transfer function model orders
           */
          nag_tsa_transf_orders(nser, &transfj, &fail);
          if (fail.code != NE_NOERROR) {
            printf("Error from nag_tsa_transf_orders (g13byc)"
                   ".\n%s\n",
                   fail.message);
            exit_status = 1;
            goto END;
          }

          rms[0] = 0;
          transfj.nag_b = 0;
          transfj.nag_q = 0;
          transfj.nag_p = 0;
          transfj.nag_r = 1;
          for (i = 1; i <= 7; ++i)
            mrxx[i - 1] = 0;
          parxx[0] = 0;

          /* Tell nag_tsa_multi_inputmod_forecast (g13bjc) not to
           * print parameters on entry */
          options.list = Nag_FALSE;

          /* nag_tsa_multi_inputmod_forecast (g13bjc).
           * Forecasting function
           */
          nag_tsa_multi_inputmod_forecast(
              &arimaj, nser, &transfj, parx, npara, nx, iqxd, x, tdxxy, rms,
              mrxx, tdmrx, parxx, ldparx, tdparx, fva, fsd, &options, &fail);
          if (fail.code != NE_NOERROR) {
            printf("Error from nag_tsa_multi_inputmod_forecast "
                   "(g13bjc).\n%s\n",
                   fail.message);
            exit_status = 1;
            goto END;
          }
        }

        /* Calculate series length */
        ny = nx + iqxd;

        /* Allocate array y */
        if (!(y = NAG_ALLOC(ny, double))) {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }

        /* Move backforecasts to start of y array */
        j = iqxd;
        for (i = 1; i <= iqxd; ++i) {
          y[i - 1] = fva[j - 1];
          --j;
        }

        /* Move series into y */
        j = iqxd + 1;
        k = nx;
        for (i = 1; i <= nx; ++i) {
          if (j > 215)
            goto END;
          y[j - 1] = x[k - 1];
          ++j;
          --k;
        }
      }

      /* Move ARIMA for series into mr */
      for (i = 1; i <= 7; ++i)
        mr[i + 2] = mrx[i - 1];

      arimas.p = mr[3];
      arimas.d = mr[4];
      arimas.q = mr[5];
      arimas.bigp = mr[6];
      arimas.bigd = mr[7];
      arimas.bigq = mr[8];
      arimas.s = mr[9];

      /* Move parameters of ARIMA for y into par */
      for (i = 1; i <= nparx; ++i)
        par[npar + i - 1] = parx[i - 1];
      npar += nparx;

      /* Move constant and reset sign reversal */
      cy = cx;
      if (idd % 2 != 0)
        cy = -cy;

      /* Set parameters for call to filter routine
       * nag_tsa_multi_filter_transf (g13bbc) */
      nb = ny + MAX(mr[0] + mr[1], mr[2]);

      /* Allocate array b */
      if (!(b = NAG_ALLOC(nb, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      /* Filter series by call to nag_tsa_multi_filter_transf (g13bbc) */
      /* nag_tsa_multi_filter_transf (g13bbc).
       * Multivariate time series, filtering by a transfer
       * function model
       */
      nag_tsa_multi_filter_transf(y, ny, &transfv, &arimas, par, npar, cy, b,
                                  nb, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_tsa_multi_filter_transf (g13bbc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }

      printf("                  Original        Filtered\n");
      printf(" Backforecasts    y-series         series\n");
      if (iqxd != 0) {
        ij = -iqxd;
        for (i = 1; i <= iqxd; ++i) {
          printf("%8" NAG_IFMT "%17.1f%16.1f\n", ij, y[i - 1], b[i - 1]);
          ++ij;
        }

        printf("\n");
        printf("        Filtered       Filtered"
               "       Filtered       Filtered\n");
        printf("         series         series"
               "         series         series\n");
        for (i = iqxd + 1; i <= ny; i += 4) {
          for (ii = i; ii <= MIN(ny, i + 3); ++ii) {
            printf("%5" NAG_IFMT "", ii - iqxd);
            printf("%10.1f", b[ii - 1]);
          }
          printf("\n");
        }
      }
    }
  }

END:

  /* Free the options structure used by nag_tsa_multi_inputmod_forecast
   * (g13bjc) */
  /* nag_tsa_free (g13xzc).
   * Freeing function for use with g13 option setting
   */
  nag_tsa_free(&options);

  NAG_FREE(b);
  NAG_FREE(fsd);
  NAG_FREE(fva);
  NAG_FREE(par);
  NAG_FREE(parx);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(rms);
  NAG_FREE(parxx);
  NAG_FREE(mrxx);

  return exit_status;
}