NAG Library Routine Document
g13eaf (multi_kalman_sqrt_var)
1
Purpose
g13eaf performs a combined measurement and time update of one iteration of the time-varying Kalman filter using a square root covariance filter.
2
Specification
Fortran Interface
Subroutine g13eaf ( |
n, m, l, a, lds, b, stq, q, ldq, c, ldm, r, s, k, h, tol, iwk, wk, ifail) |
Integer, Intent (In) | :: | n, m, l, lds, ldq, ldm | Integer, Intent (Inout) | :: | ifail | Integer, Intent (Out) | :: | iwk(m) | Real (Kind=nag_wp), Intent (In) | :: | a(lds,n), b(lds,l), q(ldq,*), c(ldm,n), r(ldm,m), tol | Real (Kind=nag_wp), Intent (Inout) | :: | s(lds,n), k(lds,m), h(ldm,m) | Real (Kind=nag_wp), Intent (Out) | :: | wk((n+m)*(n+m+l)) | Logical, Intent (In) | :: | stq |
|
C Header Interface
#include <nagmk26.h>
void |
g13eaf_ (const Integer *n, const Integer *m, const Integer *l, const double a[], const Integer *lds, const double b[], const logical *stq, const double q[], const Integer *ldq, const double c[], const Integer *ldm, const double r[], double s[], double k[], double h[], const double *tol, Integer iwk[], double wk[], Integer *ifail) |
|
3
Description
The Kalman filter arises from the state space model given by:
where
is the state vector of length
at time
,
is the observation vector of length
at time
, and
of length
and
of length
are the independent state noise and measurement noise respectively.
The estimate of
given observations
to
is denoted by
with state covariance matrix
, while the estimate of
given observations
to
is denoted by
with covariance matrix
. The update of the estimate,
, from time
to time
, is computed in two stages. First, the measurement-update is given by
and
where
is the Kalman gain matrix. The second stage is the time-update for
which is given by
and
where
represents any deterministic control used.
The square root covariance filter algorithm provides a stable method for computing the Kalman gain matrix and the state covariance matrix. The algorithm can be summarised as
where
is an orthogonal transformation triangularizing the left-hand pre-array to produce the right-hand post-array. The relationship between the Kalman gain matrix,
, and
is given by
g13eaf requires the input of the lower triangular Cholesky factors of the noise covariance matrices
and, optionally,
and the lower triangular Cholesky factor of the current state covariance matrix,
, and returns the product of the matrices
and
,
, the Cholesky factor of the updated state covariance matrix
and the matrix
used in the computation of the likelihood for the model.
4
References
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
5
Arguments
- 1: – IntegerInput
-
On entry: , the size of the state vector.
Constraint:
.
- 2: – IntegerInput
-
On entry: , the size of the observation vector.
Constraint:
.
- 3: – IntegerInput
-
On entry: , the dimension of the state noise.
Constraint:
.
- 4: – Real (Kind=nag_wp) arrayInput
-
On entry: the state transition matrix, .
- 5: – IntegerInput
-
On entry: the first dimension of the arrays
a,
b,
s and
k as declared in the (sub)program from which
g13eaf is called.
Constraint:
.
- 6: – Real (Kind=nag_wp) arrayInput
-
On entry: the noise coefficient matrix .
- 7: – LogicalInput
-
On entry: if
, the state noise covariance matrix
is assumed to be the identity matrix. Otherwise the lower triangular Cholesky factor,
, must be provided in
q.
- 8: – Real (Kind=nag_wp) arrayInput
-
Note: the second dimension of the array
q
must be at least
if
.
On entry: if
,
q must contain the lower triangular Cholesky factor of the state noise covariance matrix,
. Otherwise
q is not referenced.
- 9: – IntegerInput
-
On entry: the first dimension of the array
q as declared in the (sub)program from which
g13eaf is called.
Constraints:
- if , ;
- otherwise .
- 10: – Real (Kind=nag_wp) arrayInput
-
On entry: the measurement coefficient matrix, .
- 11: – IntegerInput
-
On entry: the first dimension of the arrays
c,
r and
h as declared in the (sub)program from which
g13eaf is called.
Constraint:
.
- 12: – Real (Kind=nag_wp) arrayInput
-
On entry: the lower triangular Cholesky factor of the measurement noise covariance matrix .
- 13: – Real (Kind=nag_wp) arrayInput/Output
-
On entry: the lower triangular Cholesky factor of the state covariance matrix, .
On exit: the lower triangular Cholesky factor of the state covariance matrix, .
- 14: – Real (Kind=nag_wp) arrayOutput
-
On exit: the Kalman gain matrix, , premultiplied by the state transition matrix, , .
- 15: – Real (Kind=nag_wp) arrayOutput
-
On exit: the lower triangular matrix .
- 16: – Real (Kind=nag_wp)Input
-
On entry: the tolerance used to test for the singularity of
. If
, then
is used instead. The inverse of the condition number of
is estimated by a call to
f07tgf (dtrcon). If this estimate is less than
tol then
is assumed to be singular.
Suggested value:
.
Constraint:
.
- 17: – Integer arrayWorkspace
- 18: – Real (Kind=nag_wp) arrayWorkspace
-
- 19: – IntegerInput/Output
-
On entry:
ifail must be set to
,
. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
.
When the value is used it is essential to test the value of ifail on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
-
On entry, .
Constraint: .
On entry, and .
Constraint: .
On entry, .
Constraint: .
On entry, and .
Constraint: .
On entry, and .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
-
The matrix is singular.
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
The use of the square root algorithm improves the stability of the computations as compared with the direct coding of the Kalman filter. The accuracy will depend on the model.
8
Parallelism and Performance
g13eaf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g13eaf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
For models with time-invariant
and
,
g13ebf can be used.
The estimate of the state vector
can be computed from
by
where
are the independent one step prediction residuals. The required matrix-vector multiplications can be performed by
f06paf (dgemv).
If
and
are independent multivariate
Normal variates then the log-likelihood for observations
is given by
where
is a constant.
The Cholesky factors of the covariance matrices can be computed using
f07fdf (dpotrf).
Note that the model
can be specified either with
b set to the identity matrix and
and the matrix
input in
q or with
and
b set to
.
The algorithm requires
operations and is backward stable (see
Verhaegen and van Dooren (1986)).
10
Example
This example first inputs the number of updates to be computed and the problem sizes. The initial state vector and state covariance matrix are input followed by the model matrices
and optionally
. The Cholesky factors of the covariance matrices can be computed if required. The model matrices can be input at each update or only once at the first step. At each update the observed values are input and the residuals are computed and printed and the estimate of the state vector,
, and the deviance are updated. The deviance is
log-likelihood ignoring the constant. After the final update the state covariance matrix is computed from
s and printed along with final estimate of the state vector and the value of the deviance.
The data is for a two-dimensional time series to which a VARMA
has been fitted. For the specification of a VARMA model as a state space model see the
G13 Chapter Introduction. The initial value of
,
, is the solution to
For convenience, the mean of each series is input before the first update and subtracted from the observations before the measurement update is computed.
10.1
Program Text
Program Text (g13eafe.f90)
10.2
Program Data
Program Data (g13eafe.d)
10.3
Program Results
Program Results (g13eafe.r)