NAG Library Function Document
nag_kalman_sqrt_filt_cov_var (g13eac)
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
the estimate of
given observations
to
is denoted by
with
.
nag_kalman_sqrt_filt_cov_var (g13eac) performs one recursion of the square root covariance filter algorithm, summarised as follows:
where
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
is
where
is the Kalman gain matrix, whilst the time-update for
is
where
represents any deterministic control used. The relationship between the Kalman gain matrix
and
is given by
The function returns the product of the matrices
and
represented as
, and the state covariance matrix
factorized as
(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:
– IntegerInput
-
On entry: the actual state dimension, , i.e., the order of the matrices and .
Constraint:
.
- 2:
– IntegerInput
-
On entry: the actual input dimension, , i.e., the order of the matrix .
Constraint:
.
- 3:
– IntegerInput
-
On entry: the actual output dimension, , i.e., the order of the matrix .
Constraint:
.
- 4:
– doubleInput/Output
-
Note: the th element of the matrix is stored in .
On entry: the leading by lower triangular part of this array must contain , the left Cholesky factor of the state covariance matrix .
On exit: the leading by lower triangular part of this array contains , the left Cholesky factor of the state covariance matrix .
- 5:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
s.
Constraint:
.
- 6:
– const doubleInput
-
Note: the th element of the matrix is stored in .
On entry: the leading by part of this array must contain , the state transition matrix of the discrete system.
- 7:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
a.
Constraint:
.
- 8:
– const doubleInput
-
Note: the th element of the matrix is stored in .
On entry: if
q is not
NULL then the leading
by
part of this array must contain the matrix
, otherwise if
q is
NULL then the leading
by
part of the array must contain the matrix
.
is the input weight matrix and
is the noise covariance matrix.
- 9:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
b.
Constraint:
.
- 10:
– const doubleInput
-
Note: the th element of the matrix is stored in .
On entry: if the noise covariance matrix is to be supplied separately from the input weight matrix then the leading
by
lower triangular part of this array must contain
, 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
then the array
q must be set to
NULL.
- 11:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
q.
Constraint:
if
q is defined.
- 12:
– const doubleInput
-
Note: the th element of the matrix is stored in .
On entry: the leading by part of this array must contain , the output weight matrix of the discrete system.
- 13:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
c.
Constraint:
.
- 14:
– const doubleInput
-
Note: the th element of the matrix is stored in .
On entry: the leading by lower triangular part of this array must contain , the left Cholesky factor of the measurement noise covariance matrix.
- 15:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
r.
Constraint:
.
- 16:
– doubleOutput
-
Note: the th element of the matrix is stored in .
On exit: if
k is not
NULL then the leading
by
part of
k contains
, the product of the Kalman filter gain matrix
with the state transition matrix
. If
is not required then
k must be set to
NULL.
- 17:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
k.
Constraint:
if
k is not
NULL,
- 18:
– doubleOutput
-
Note: the th element of the matrix is stored in .
On exit: if
k is not
NULL then the leading
by
lower triangular part of this array contains
. If
k is
NULL then
h is not referenced and may be set to
NULL.
- 19:
– IntegerInput
-
On entry: the stride separating matrix column elements in the array
h.
Constraint:
if both
k and
h are not
NULL,
- 20:
– doubleInput
-
On entry: if both
k and
h are not
NULL then
tol is used to test for near singularity of the matrix
. If you set
tol to be less than
then the tolerance is taken as
, where
is the
machine precision. Otherwise,
tol need not be set by you.
- 21:
– NagError *Input/Output
-
The NAG error argument (see
Section 3.7 in How to Use the NAG Library and its Documentation).
6
Error Indicators and Warnings
- NE_2_INT_ARG_LT
-
On entry while . These arguments must satisfy .
On entry while . These arguments must satisfy .
On entry while . These arguments must satisfy .
On entry while . These arguments must satisfy .
On entry while . These arguments must satisfy .
On entry while . These arguments must satisfy .
On entry while . These arguments must satisfy .
On entry, while . These arguments must satisfy .
- NE_ALLOC_FAIL
-
Dynamic memory allocation failed.
- NE_INT_ARG_LT
-
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
- NE_MAT_SINGULAR
-
The matrix sqrt() 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
nag_kalman_sqrt_filt_cov_var (g13eac) is not threaded in any implementation.
The algorithm requires
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 . 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
and
) are generated using the function
nag_rand_arma (g05phc). The Kalman filter and optimization function
nag_opt_nlp (e04ucc) are then used to find the maximum likelihood estimate for the time series arguments
and
. The ARMA(1,1) time series is defined by
This has the following state space representation (
Harvey and Phillips (1979))
where the state vector
and
is uncorrelated white noise with zero mean and variance
, i.e.,
Since
we arrive at the following Kalman Filter matrices
The initial estimates for the state vector,
, and state covariance matrix,
, are:
Since
(
Wei (1990))
Using
gives an initial Cholesky ‘square root’ of
10.1
Program Text
Program Text (g13eace.c)
10.2
Program Data
None.
10.3
Program Results
Program Results (g13eace.r)