/* 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. * Mark 8 revised, 2004. * */ #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 int ex1(void); static int ex2(void); int main(void) { Integer exit_status_ex1=0; Integer exit_status_ex2=0; exit_status_ex1 = ex1(); exit_status_ex2 = ex2(); return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1; } #define NY 2000 static void objfun(Integer n, 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; NagError fail; 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]; /* There is no measurement noise */ double s[2][2], ss=0.0, temp1, temp2, theta, tol=0.0; double v, xp[2], *y; INIT_FAIL(fail); 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) { /* nag_kalman_sqrt_filt_cov_var (g13eac). * One iteration step of the time-varying Kalman filter * recursion using the square root covariance implementation */ nag_kalman_sqrt_filt_cov_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, 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 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] #define H(I,J) h[(I)*tdh + J] static int ex1(void) { double *a=0, *b=0, *c=0, *k=0, *q=0, *r=0, *s=0, *h=0; Integer exit_status=0; Integer i, j; Integer m, n, p; Integer istep; double tol; Integer tda, tdb, tdc, tdk, tdq, tdr, tds, tdh; NagError fail; INIT_FAIL(fail); Vprintf("nag_kalman_sqrt_filt_cov_var (g13eac) Example Program Results\n\n"); Vprintf("Example 1\n"); /* Skip the heading in the data file */ Vscanf("%*[^\n]"); Vscanf("%ld%ld%ld%lf",&n,&m,&p,&tol); if (n>=1 && m>=1 && p>=1) { if ( !( a = NAG_ALLOC(n*n, double)) || !( b = NAG_ALLOC(n*m, double)) || !( c = NAG_ALLOC(p*n, double)) || !( k = NAG_ALLOC(n*p, double)) || !( q = NAG_ALLOC(m*m, double)) || !( r = NAG_ALLOC(p*p, double)) || !( s = NAG_ALLOC(n*n, double)) || !( h = NAG_ALLOC(n*p, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdb = m; tdc = n; tdk = p; tdq = m; tdr = p; tds = n; tdh = p; } else { Vprintf("Invalid n or m or p.\n"); exit_status = 1; return exit_status; } /* Read data */ for (i=0; i