/* nag_kalman_sqrt_filt_cov_var(g13eac) Example Program * * Copyright 1996 Numerical Algorithms Group * * Mark 4, 1996. * Mark 5 revised, 1998. * Mark 6 revised, 2000. * Mark 7 revised, 2001. * */ #include #include #include #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void objfun(Integer n, double theta_phi[], double *objf, double g[], Nag_Comm *comm); #ifdef __cplusplus } #endif static void ex1(void); static void ex2(void); int main(void) { ex1(); ex2(); return EXIT_SUCCESS; } #define NY 2000 static void objfun(Integer n, double theta_phi[], double *objf, double g[], Nag_Comm *comm) /* Routine to evaluate objective function. */ { double ak[2][1], h[1][1]; double tol = 0.0; double xp[2]; double q[1][1]; double c[1][2]; double r[1][1]; /* There is no measurement noise */ double a[2][2]; double s[2][2]; double b[2][1]; double hs, v, ss = 0.0, logdet = 0.0; Integer k, nsum = 0, ione = 1, itwo = 2; Integer nsteps = NY, n1 = 2, m1 = 1, p1 = 1; double phi, theta, temp1, temp2, k11; double *y; y = comm->user; xp[0]= 0.0; /* The expectation of the mean of an * ARMA(1,1) is 0.0 */ xp[1] = 0.0; q[0][0] = 1.0; c[0][0] = 1.0; c[0][1] = 0.0; r[0][0] = 0.0; /* There is no measurement noise */ 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) { g13eac(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, NAGERR_DEFAULT); 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; } /* objfun */ #define NMAX 20 #define MMAX 20 #define PMAX 20 #define TRADIM 20 static void ex1(void) { double a[NMAX][TRADIM], b[NMAX][TRADIM], c[PMAX][TRADIM], k[NMAX][TRADIM], q[MMAX][TRADIM], r[PMAX][TRADIM], s[NMAX][TRADIM], h[NMAX][TRADIM]; Integer i, j; Integer m, n, p; Integer istep; double tol; Integer nmax, mmax, pmax, tradim; Vprintf("g13eac Example Program Results\n\n"); Vprintf("Example 1\n"); /* Skip the heading in the data file */ Vscanf("%*[^\n]"); nmax = NMAX; mmax = MMAX; pmax = PMAX; tradim = TRADIM; Vscanf("%ld%ld%ld%lf",&n,&m,&p,&tol); if (n<=0 || m<=0 || p<=0 || n>nmax || m>mmax || p>pmax) { Vfprintf(stderr, "One of n m or p is out of range " "n = %ld, m = %ld, p = %ld\n", n, m, p); exit(EXIT_FAILURE); } /* Read data */ for (i=0; i