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

NAG CL Interface Introduction
Example description
/* nag_tsa_multi_kalman_sqrt_var (g13eac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group
 *
 * Mark 29.0, 2023.
 *
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL objfun(Integer n, const double theta_phi[], double *objf,
                            double g[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

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

int main(void) {
  /* Integer scalar and array declarations */
  Integer exit_status_ex1 = 0;
  Integer exit_status_ex2 = 0;

  printf("nag_tsa_multi_kalman_sqrt_var (g13eac) Example Program Results\n\n");

  /* Run example 1 */
  exit_status_ex1 = ex1();

  /* Run example 2 */
  exit_status_ex2 = ex2();

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

/* Start of the first example ... */
#define A(I, J) a[(I)*tda + J]
#define B(I, J) b[(I)*tdb + J]
#define C(I, J) c[(I)*tdc + J]
#define K(I, J) k[(I)*tdk + J]
#define Q(I, J) q[(I)*tdq + J]
#define R(I, J) r[(I)*tdr + J]
#define S(I, J) s[(I)*tds + J]

static int ex1() {
  /* Integer scalar and array declarations */
  Integer exit_status = 0;
  Integer i, j, m, n, p, istep, tda, tdb, tdc, tdk, tdq, tdr, tds, tdh;

  /* Double scalar and array declarations */
  double tol;
  double *a = 0, *b = 0, *c = 0, *k = 0, *q = 0, *r = 0, *s = 0, *h = 0;

  /* NAG structures */
  NagError fail;

  /* Initialize the error structure */
  INIT_FAIL(fail);

  printf("Example 1\n");

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

  /* Get the problem size */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%lf", &n, &m, &p, &tol);
  if (n < 1 || m < 1 || p < 1) {
    printf("Invalid n or m or p.\n");
    exit_status = 1;
    goto END;
  }

  tda = tdc = tds = n;
  tdb = tdq = m;
  tdk = tdr = tdh = p;

  /* Allocate arrays */
  if (!(a = NAG_ALLOC(n * tda, double)) || !(b = NAG_ALLOC(n * tdb, double)) ||
      !(c = NAG_ALLOC(p * tdc, double)) || !(k = NAG_ALLOC(n * tdk, double)) ||
      !(q = NAG_ALLOC(m * tdq, double)) || !(r = NAG_ALLOC(p * tdr, double)) ||
      !(s = NAG_ALLOC(n * tds, double)) || !(h = NAG_ALLOC(n * tdh, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read data */
  for (i = 0; i < n; ++i)
    for (j = 0; j < n; ++j)
      scanf("%lf", &S(i, j));
  for (i = 0; i < n; ++i)
    for (j = 0; j < n; ++j)
      scanf("%lf", &A(i, j));
  for (i = 0; i < n; ++i)
    for (j = 0; j < m; ++j)
      scanf("%lf", &B(i, j));
  if (q) {
    for (i = 0; i < m; ++i)
      for (j = 0; j < m; ++j)
        scanf("%lf", &Q(i, j));
  }
  for (i = 0; i < p; ++i)
    for (j = 0; j < n; ++j)
      scanf("%lf", &C(i, j));
  for (i = 0; i < p; ++i)
    for (j = 0; j < p; ++j)
      scanf("%lf", &R(i, j));

  /* Perform three iterations of the Kalman filter recursion */
  for (istep = 1; istep <= 3; ++istep) {
    /* Run the kalman filter */
    nag_tsa_multi_kalman_sqrt_var(n, m, p, s, tds, a, tda, b, tdb, q, tdq, c,
                                  tdc, r, tdr, k, tdk, h, tdh, tol, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_tsa_multi_kalman_sqrt_var (g13eac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }

  /* Print the results */
  printf("\nThe square root of the state covariance matrix is\n\n");
  for (i = 0; i < n; ++i) {
    for (j = 0; j < n; ++j)
      printf("%8.4f ", S(i, j));
    printf("\n");
  }
  if (k) {
    printf("\nThe matrix AK (the product of the Kalman gain\n");
    printf("matrix with the state transition matrix) is\n\n");
    for (i = 0; i < n; ++i) {
      for (j = 0; j < p; ++j)
        printf("%8.4f ", K(i, j));
      printf("\n");
    }
  }

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(k);
  NAG_FREE(q);
  NAG_FREE(r);
  NAG_FREE(s);
  NAG_FREE(h);

  return exit_status;
}

/* ... end of the first example */

/* Start of the second example ... */
#define NY 2000

/* This illustrates the use of the kalman filter to estimate the
   parameters of an ARMA(1,1) time series model.
   Note : theta_phi[0] contains theta (moving average coefficient), and
   theta_phi[1] contains phi (autoregressive coefficient)
 */
static int ex2(void) {
  /* Integer scalar and array declarations */
  Integer exit_status = 0;
  Integer n, lr;
  Integer lstate;
  Integer *state = 0;

  /* Double scalar and array declarations */
  double a[1], sy[NY];
  double *theta_phi = 0, *g = 0, *bl = 0, *bu = 0, *r = 0;
  double objf, var;

  /* NAG structures and data types */
  Nag_Comm comm;
  Nag_E04_Opt options;
  NagError fail;
  Nag_ModeRNG mode;

  /* Choose the base generator */
  Nag_BaseRNG genid = Nag_Basic;
  Integer subid = 0;

  /* Set the seed */
  Integer seed[] = {1762543};
  Integer lseed = 1;

  /* Set the autoregressive coefficients for the (randomly generated)
     time series */
  Integer ip = 1;
  double sphi[] = {0.4};

  /* Set the moving average coefficients for the (randomly generated)
     time series */
  Integer iq = 1;
  double stheta[] = {0.9};

  /* Number of terms in the (randomly generated) time series */
  Integer nterms = NY;

  /* Mean and variance of the (randomly generated) time series */
  double mean = 0.0, vara = 1.0;

  /* Initialize the error structure */
  INIT_FAIL(fail);

  printf("\n\nExample 2\n\n");

  /* Get the length of the state array */
  lstate = -1;
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Calculate the size of the reference vector */
  lr = (ip > iq + 1) ? ip : iq + 1;
  lr += ip + iq + 6;

  /* Allocate arrays */
  if (!(state = NAG_ALLOC(lstate, Integer)) || !(r = NAG_ALLOC(lr, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize the generator to a repeatable sequence */
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Generate a time series with nterms terms */
  mode = Nag_InitializeAndGenerate;
  nag_rand_times_arma(mode, nterms, mean, ip, sphi, iq, stheta, vara, r, lr,
                      state, &var, sy, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_times_arma (g05phc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Number of parameters to estimate */
  n = ip + iq;

  /* Allocate arrays */
  if (!(theta_phi = NAG_ALLOC(n, double)) || !(g = NAG_ALLOC(n, double)) ||
      !(bl = NAG_ALLOC(n, double)) || !(bu = NAG_ALLOC(n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Make an initial guess of the parameters */
  theta_phi[0] = 0.6;
  theta_phi[1] = 0.4;

  /* Set the bounds */
  bl[0] = -1.0;
  bu[0] = 1.0;
  bl[1] = -1.0;
  bu[1] = 1.0;
  comm.user = &sy[0];

  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options);
  strcpy(options.outfile, "stdout");
  options.print_level = Nag_NoPrint;
  options.list = Nag_FALSE;

  options.obj_deriv = Nag_FALSE;
  options.con_deriv = Nag_FALSE;
  options.print_deriv = Nag_D_NoPrint;

  /* nag_opt_nlp1_solve (e04ucc).
   * Minimization with nonlinear constraints using a sequential QP method.
   */
  nag_opt_nlp1_solve(n, 0, 0, a, 1, bl, bu, objfun, NULLFN, theta_phi, &objf, g,
                     &options, &comm, &fail);
  if (fail.code != NE_NOERROR && fail.code != NW_KT_CONDITIONS) {
    printf("Error from nag_opt_nlp1_solve (e04ucc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* nag_opt_free (e04xzc).
   * Memory freeing function for use with option setting
   */
  nag_opt_free(&options, "all", &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Display the results */
  printf("The estimates are :  theta = %7.3f, phi = %7.3f \n", theta_phi[0],
         theta_phi[1]);

END:
  NAG_FREE(state);
  NAG_FREE(theta_phi);
  NAG_FREE(g);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(r);

  return exit_status;
}

/* Define  objective function used in the non-linear optimization routine ... */

static void NAG_CALL objfun(Integer n, const double theta_phi[], double *objf,
                            double g[], Nag_Comm *comm) {
  /* Routine to evaluate objective function. */

  Integer ione = 1, itwo = 2, k, m1 = 1, n1 = 2, nsteps = NY, nsum = 0, p1 = 1;
  double a[2][2], ak[2][1], b[2][1], c[1][2], h[1][1], hs, k11, logdet = 0.0;
  double phi, q[1][1], r[1][1];
  double s[2][2], ss = 0.0, temp1, temp2, theta, tol = 0.0;
  double v, xp[2], *y;
  NagError fail;

  /* Initialize the error structure */
  INIT_FAIL(fail);

  y = comm->user;
  /* The expectation of the mean of an ARMA(1,1) is 0.0 */
  xp[0] = 0.0;
  xp[1] = 0.0;
  q[0][0] = 1.0;

  c[0][0] = 1.0;
  c[0][1] = 0.0;

  /* There is no measurement noise */
  r[0][0] = 0.0;
  theta = theta_phi[0];
  phi = theta_phi[1];

  b[0][0] = 1.0;
  b[1][0] = -theta;

  a[0][0] = phi;
  a[1][0] = 0.0;
  a[0][1] = 1.0;
  a[1][1] = 0.0;

  /* set value for Cholesky factor of state covariance matrix  */
  temp1 = 1.0 + (theta * theta) - (2.0 * theta * phi);
  temp2 = 1.0 - (phi * phi);
  k11 = temp1 / temp2;
  s[0][0] = sqrt(k11);
  s[0][1] = 0.0;
  s[1][0] = -theta / s[0][0];
  s[1][1] = theta * sqrt(1.0 - (1.0 / k11));

  /* iterate kalman filter for number of observations  */
  for (k = 1; k <= nsteps; ++k) {
    /* nag_tsa_multi_kalman_sqrt_var (g13eac).
     * One iteration step of the time-varying Kalman filter
     * recursion using the square root covariance implementation
     */
    nag_tsa_multi_kalman_sqrt_var(n1, m1, p1, &s[0][0], itwo, &a[0][0], itwo,
                                  &b[0][0], ione, &q[0][0], ione, &c[0][0],
                                  itwo, &r[0][0], ione, &ak[0][0], ione,
                                  &h[0][0], ione, tol, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from  nag_tsa_multi_kalman_sqrt_var (g13eac).\n%s\n",
             fail.message);
      *objf = 0.0;
      return;
    }

    v = y[k - 1] - c[0][0] * xp[0];
    hs = h[0][0] * h[0][0];
    logdet = logdet + log(hs);
    ss = ss + (v * v / hs);
    nsum = nsum + 1;

    xp[0] = a[0][0] * xp[0] + a[0][1] * xp[1] + ak[0][0] * v;
    xp[1] = ak[1][0] * v;
  }
  *objf = nsum * log(ss / nsum) + logdet;
}

/* ... end of the objective function definition */
/* ... end of the second example */