nag_kalman_sqrt_filt_cov_var (g13eac) (PDF version)
g13 Chapter Contents
g13 Chapter Introduction
NAG Library Manual

NAG Library Function Document

nag_kalman_sqrt_filt_cov_var (g13eac)

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

nag_kalman_sqrt_filt_cov_var (g13eac) performs a combined measurement and time update of one iteration of the time-varying Kalman filter. The method employed for this update is the square root covariance filter with the system matrices in their original form.

2  Specification

#include <nag.h>
#include <nagg13.h>
void  nag_kalman_sqrt_filt_cov_var (Integer n, Integer m, Integer p, double s[], Integer tds, const double a[], Integer tda, const double b[], Integer tdb, const double q[], Integer tdq, const double c[], Integer tdc, const double r[], Integer tdr, double k[], Integer tdk, double h[], Integer tdh, double tol, NagError *fail)

3  Description

For the state space system defined by
X i+1 = A i X i + B i W i var W i = Q i Y i = C i X i + V i var V i = R i
the estimate of X i  given observations Y 1  to Y i-1  is denoted by X ^ ii - 1  with var X ^ ii - 1 = P ii - 1 = S i SiT . nag_kalman_sqrt_filt_cov_var (g13eac) performs one recursion of the square root covariance filter algorithm, summarised as follows:
R i 1/2 C i S i 0 0 A i S i B i Q i 1/2 U = H i 1/2 0 0 G i S i+1 0 Pre-array Post-array
where U  is an orthogonal transformation triangularizing the pre-array. The triangularisation is carried out via Householder transformations exploiting the zero pattern in the pre-array. The measurement-update for the estimated state vector X  is
X ^ ii = X ^ ii - 1 - K i C i X ^ ii - 1 - Y i (1)
where K i  is the Kalman gain matrix, whilst the time-update for X  is
X ^ i + 1i = A i X ^ ii + D i U i (2)
where D i U i  represents any deterministic control used. The relationship between the Kalman gain matrix K i  and G i  is given by
A i K i = G i H i 1/2 -1
The function returns the product of the matrices A i  and K i  represented as AK i , and the state covariance matrix P ii - 1  factorized as P ii - 1 = S i SiT  (see the g13 Chapter Introduction for more information concerning the covariance filter).

4  References

Anderson B D O and Moore J B (1979) Optimal Filtering Prentice–Hall
Harvey A C and Phillips G D A (1979) Maximum likelihood estimation of regression models with autoregressive — moving average disturbances Biometrika 66 49–58
Vanbegin M, van Dooren P and Verhaegen M H G (1989) Algorithm 675: FORTRAN subroutines for computing the square root covariance filter and square root information filter in dense or Hessenberg forms ACM Trans. Math. Software 15 243–256
Verhaegen M H G and van Dooren P (1986) Numerical aspects of different Kalman filter implementations IEEE Trans. Auto. Contr. AC-31 907–917
Wei W W S (1990) Time Series Analysis: Univariate and Multivariate Methods Addison–Wesley

5  Arguments

1:     nIntegerInput
On entry: the actual state dimension, n , i.e., the order of the matrices S i  and A i .
Constraint: n1 .
2:     mIntegerInput
On entry: the actual input dimension, m , i.e., the order of the matrix Q i 1/2 .
Constraint: m1 .
3:     pIntegerInput
On entry: the actual output dimension, p , i.e., the order of the matrix R i 1/2 .
Constraint: p1 .
4:     s[n×tds]doubleInput/Output
Note: the i,jth element of the matrix S is stored in s[i-1×tds+j-1].
On entry: the leading n  by n  lower triangular part of this array must contain S i , the left Cholesky factor of the state covariance matrix P ii - 1 .
On exit: the leading n  by n  lower triangular part of this array contains S i+1 , the left Cholesky factor of the state covariance matrix P i + 1i .
5:     tdsIntegerInput
On entry: the stride separating matrix column elements in the array s.
Constraint: tdsn .
6:     a[n×tda]const doubleInput
Note: the i,jth element of the matrix A is stored in a[i-1×tda+j-1].
On entry: the leading n  by n  part of this array must contain A i , the state transition matrix of the discrete system.
7:     tdaIntegerInput
On entry: the stride separating matrix column elements in the array a.
Constraint: tdan .
8:     b[n×tdb]const doubleInput
Note: the i,jth element of the matrix B is stored in b[i-1×tdb+j-1].
On entry: if q is not NULL then the leading n  by m  part of this array must contain the matrix B i , otherwise if q is NULL then the leading n  by m  part of the array must contain the matrix B i Q i 1/2 . B i  is the input weight matrix and Q i  is the noise covariance matrix.
9:     tdbIntegerInput
On entry: the stride separating matrix column elements in the array b.
Constraint: tdbm .
10:   q[m×tdq]const doubleInput
Note: the i,jth element of the matrix Q is stored in q[i-1×tdq+j-1].
On entry: if the noise covariance matrix is to be supplied separately from the input weight matrix then the leading m  by m  lower triangular part of this array must contain Q i 1/2 , the left Cholesky factor of the input process noise covariance matrix. If the noise covariance matrix is to be input with the weight matrix as B i Q i 1/2  then the array q must be set to NULL.
11:   tdqIntegerInput
On entry: the stride separating matrix column elements in the array q.
Constraint: tdqm  if q is defined.
12:   c[p×tdc]const doubleInput
Note: the i,jth element of the matrix C is stored in c[i-1×tdc+j-1].
On entry: the leading p  by n  part of this array must contain C i , the output weight matrix of the discrete system.
13:   tdcIntegerInput
On entry: the stride separating matrix column elements in the array c.
Constraint: tdcn .
14:   r[p×tdr]const doubleInput
Note: the i,jth element of the matrix R is stored in r[i-1×tdr+j-1].
On entry: the leading p  by p  lower triangular part of this array must contain R i 1/2 , the left Cholesky factor of the measurement noise covariance matrix.
15:   tdrIntegerInput
On entry: the stride separating matrix column elements in the array r.
Constraint: tdrp .
16:   k[n×tdk]doubleOutput
Note: the i,jth element of the matrix K is stored in k[i-1×tdk+j-1].
On exit: if k is not NULL then the leading n  by p  part of k contains AK i , the product of the Kalman filter gain matrix K i  with the state transition matrix A i . If AKi is not required then k must be set to NULL.
17:   tdkIntegerInput
On entry: the stride separating matrix column elements in the array k.
Constraint: if k is not NULL, tdkp
18:   h[p×tdh]doubleOutput
Note: the i,jth element of the matrix H is stored in h[i-1×tdh+j-1].
On exit: if k is not NULL then the leading p  by p  lower triangular part of this array contains H i 1/2 . If k is NULL then h is not referenced and may be set to NULL.
19:   tdhIntegerInput
On entry: the stride separating matrix column elements in the array h.
Constraint: if both k and h are not NULL, tdhp
20:   toldoubleInput
On entry: if both k and h are not NULL then tol is used to test for near singularity of the matrix H i 1/2 . If you set tol to be less than p 2 ε  then the tolerance is taken as p 2 ε , where ε  is the machine precision. Otherwise, tol need not be set by you.
21:   failNagError *Input/Output
The NAG error argument (see Section 3.6 in the Essential Introduction).

6  Error Indicators and Warnings

NE_2_INT_ARG_LT
On entry tda=value  while n=value . These arguments must satisfy tdan .
On entry tdb=value  while m=value . These arguments must satisfy tdbm .
On entry tdc=value  while n=value . These arguments must satisfy tdcn .
On entry tdh=value  while p=value . These arguments must satisfy tdhp .
On entry tdk=value  while p=value . These arguments must satisfy tdkp .
On entry tdq=value  while m=value . These arguments must satisfy tdqm .
On entry tdr=value  while p=value . These arguments must satisfy tdrp .
On entry, tds=value  while n=value . These arguments must satisfy tdsn .
NE_ALLOC_FAIL
Dynamic memory allocation failed.
NE_INT_ARG_LT
On entry, m=value.
Constraint: m1.
On entry, n=value.
Constraint: n1.
On entry, p=value.
Constraint: p1.
NE_MAT_SINGULAR
The matrix sqrt(H) is singular.
NE_NULL_ARRAY
Array h has null address.

7  Accuracy

The use of the square root algorithm improves the stability of the computations.

8  Parallelism and Performance

Not applicable.

9  Further Comments

The algorithm requires 7 6 n 3 + n 2 5 2 p + m + n 1 2 m 2 + p 2  operations and is backward stable (see Vanbegin et al. (1989)).

10  Example

For this function two examples are presented. There is a single example program for nag_kalman_sqrt_filt_cov_var (g13eac), with a main program and the code to solve the two example problems is given in the functions ex1 and ex2.
Example 1 (ex1)
To apply three iterations of the Kalman filter (in square root covariance form) to the system A i , B i , C i . The same data is used for all three iterative steps.
Example 2 (ex2)
In the second example 2000 terms of an ARMA(1,1) time series (with σ 2 = 1.0 , θ = 0.9  and ϕ=0.4 ) are generated using the function nag_rand_arma (g05phc). The Kalman filter and optimization function nag_opt_nlp_solve (e04wdc) are then used to find the maximum likelihood estimate for the time series arguments θ  and ϕ . The ARMA(1,1) time series is defined by
y k = ϕ y k-1 + ε k - θ ε k-1
This has the following state space representation (Harvey and Phillips (1979))
x k = ϕ 1 0 0 x k-1 + 1 -θ ε k y k = 1 0 x k
where the state vector x k = y k -θ ε k  and ε k  is uncorrelated white noise with zero mean and variance σ 2 , i.e.,
E ε k = 0 , E ε k ε k = σ 2 , E y k ε k = σ 2 ​ and ​ E ε k ε k-1 = 0 .
Since σ 2 = 1  we arrive at the following Kalman Filter matrices
A k = ϕ 1 0 0 , B k = 1 -θ C k = 1 0 , Q k = 0 ​ and ​ R k = 1 .
The initial estimates for the state vector, x 10 , and state covariance matrix, P 10 , are:
x 10 = E x k = 0 ​ and ​ P 10 = E x k xkT = E y k y k -θ E y k ε k -θ E y k ε k θ 2 E ε k ε k .
Since E y k y k = γ = 1 + θ 2 - 2 ϕ θ σ 2 1 - ϕ 2  (Wei (1990))
P 10 = γ -θ -θ θ 2 .
Using P 10 = S 10 S10T  gives an initial Cholesky ‘square root’ of
S 10 = γ 0 -θ γ θ γ - 1 γ .

10.1  Program Text

Program Text (g13eace.c)

10.2  Program Data

None.

10.3  Program Results

Program Results (g13eace.r)


nag_kalman_sqrt_filt_cov_var (g13eac) (PDF version)
g13 Chapter Contents
g13 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2014