G13EBF (PDF version)
G13 Chapter Contents
G13 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

G13EBF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

G13EBF performs a combined measurement and time update of one iteration of the time-invariant Kalman filter using a square root covariance filter.

2  Specification

SUBROUTINE G13EBF ( TRANSF, N, M, L, A, LDS, B, STQ, Q, LDQ, C, LDM, R, S, K, H, U, TOL, IWK, WK, IFAIL)
INTEGER  N, M, L, LDS, LDQ, LDM, IWK(M), IFAIL
REAL (KIND=nag_wp)  A(LDS,N), B(LDS,L), Q(LDQ,*), C(LDM,N), R(LDM,M), S(LDS,N), K(LDS,M), H(LDM,M), U(LDS,*), TOL, WK((N+M)*(N+M+L))
LOGICAL  STQ
CHARACTER(1)  TRANSF

3  Description

The Kalman filter arises from the state space model given by
Xi+1=AXi+BWi, VarWi=Qi Yi=CXi+Vi, VarVi=Ri
where Xi is the state vector of length n at time i, Yi is the observation vector of length m at time i and Wi of length l and Vi of length m are the independent state noise and measurement noise respectively. The matrices A,B and C are time invariant.
The estimate of Xi given observations Y1 to Yi-1 is denoted by X^ii-1 with state covariance matrix VarX^ii-1=Pii-1=SiSiT while the estimate of Xi given observations Y1 to Yi is denoted by X^ii with covariance matrix VarX^ii=Pii. The update of the estimate, X^ii-1, from time i to time i+1 is computed in two stages. First, the measurement-update is given by
X^ii=X^ii-1+KiYi-CX^ii-1 (1)
where Ki=PiiCTCPiiCT+Ri -1 is the Kalman gain matrix. The second stage is the time-update for X, which is given by
X^i+1i=AX^ii+DiUi (2)
where DiUi 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
Ri1/2 0 CSi 0 BQi1/2 ASi U= Hi1/2 0 0 Gi Si+1 0
where U is an orthogonal transformation triangularizing the left-hand pre-array to produce the right-hand post-array. The triangularization is carried out via Householder transformations exploiting the zero pattern of the pre-array. The relationship between the Kalman gain matrix Ki and Gi is given by
AKi=Gi Hi1/2 -1.
In order to exploit the invariant parts of the model to simplify the computation of U the results for the transformed state space U*X are computed where U* is the transformation that reduces the matrix pair A,C to lower observer Hessenberg form. That is, the matrix U* is computed such that the compound matrix
CU*T U*AU*T
is a lower trapezoidal matrix. Further the matrix B is transformed to U*B. These transformations need only be computed once at the start of a series, and G13EBF will, optionally, compute them. G13EBF returns transformed matrices U*AU*T, U*B, CU*T and U*AKi, the Cholesky factor of the updated transformed state covariance matrix Si+1* (where U*Pi+1iU*T=Si+1*Si+1 *T) and the matrix Hi1/2, valid for both transformed and original models, which is used in the computation of the likelihood for the model. Note that the covariance matrices Qi and Ri can be time-varying.

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  Parameters

1:     TRANSF – CHARACTER(1)Input
On entry: indicates whether to transform the input matrix pair A,C to lower observer Hessenberg form. The transformation will only be required on the first call to G13EBF.
TRANSF='T'
The matrices in arrays A and C are transformed to lower observer Hessenberg form and the matrices in B and S are transformed as described in Section 3.
TRANSF='H'
The matrices in arrays A, C and B should be as returned from a previous call to G13EBF with TRANSF='T'.
Constraint: TRANSF='T' or 'H'.
2:     N – INTEGERInput
On entry: n, the size of the state vector.
Constraint: N1.
3:     M – INTEGERInput
On entry: m, the size of the observation vector.
Constraint: M1.
4:     L – INTEGERInput
On entry: l, the dimension of the state noise.
Constraint: L1.
5:     A(LDS,N) – REAL (KIND=nag_wp) arrayInput/Output
On entry: if TRANSF='T', the state transition matrix, A.
If TRANSF='H', the transformed matrix as returned by a previous call to G13EBF with TRANSF='T'.
On exit: if TRANSF='T', the transformed matrix, U*AU*T, otherwise A is unchanged.
6:     LDS – INTEGERInput
On entry: the first dimension of the arrays A, B, S, K and U as declared in the (sub)program from which G13EBF is called.
Constraint: LDSN.
7:     B(LDS,L) – REAL (KIND=nag_wp) arrayInput/Output
On entry: if TRANSF='T', the noise coefficient matrix B.
If TRANSF='H', the transformed matrix as returned by a previous call to G13EBF with TRANSF='T'.
On exit: if TRANSF='T', the transformed matrix, U*B, otherwise B is unchanged.
8:     STQ – LOGICALInput
On entry: if STQ=.TRUE., the state noise covariance matrix Qi is assumed to be the identity matrix. Otherwise the lower triangular Cholesky factor, Qi1/2, must be provided in Q.
9:     Q(LDQ,*) – REAL (KIND=nag_wp) arrayInput
Note: the second dimension of the array Q must be at least L if STQ=.FALSE. and at least 1 if STQ=.TRUE..
On entry: if STQ=.FALSE., Q must contain the lower triangular Cholesky factor of the state noise covariance matrix, Qi1/2. Otherwise Q is not referenced.
10:   LDQ – INTEGERInput
On entry: the first dimension of the array Q as declared in the (sub)program from which G13EBF is called.
Constraints:
  • if STQ=.FALSE., LDQL;
  • otherwise LDQ1.
11:   C(LDM,N) – REAL (KIND=nag_wp) arrayInput/Output
On entry: if TRANSF='T', the measurement coefficient matrix, C.
If TRANSF='H', the transformed matrix as returned by a previous call to G13EBF with TRANSF='T'.
On exit: if TRANSF='T', the transformed matrix, CU*T, otherwise C is unchanged.
12:   LDM – INTEGERInput
On entry: the first dimension of the arrays C, R and H as declared in the (sub)program from which G13EBF is called.
Constraint: LDMM.
13:   R(LDM,M) – REAL (KIND=nag_wp) arrayInput
On entry: the lower triangular Cholesky factor of the measurement noise covariance matrix Ri1/2.
14:   S(LDS,N) – REAL (KIND=nag_wp) arrayInput/Output
On entry: if TRANSF='T' the lower triangular Cholesky factor of the state covariance matrix, Si.
If TRANSF='H' the lower triangular Cholesky factor of the covariance matrix of the transformed state vector Si* as returned from a previous call to G13EBF with TRANSF='T'.
On exit: the lower triangular Cholesky factor of the transformed state covariance matrix, Si+1*.
15:   K(LDS,M) – REAL (KIND=nag_wp) arrayOutput
On exit: the Kalman gain matrix for the transformed state vector premultiplied by the state transformed transition matrix, U*AKi.
16:   H(LDM,M) – REAL (KIND=nag_wp) arrayOutput
On exit: the lower triangular matrix Hi1/2.
17:   U(LDS,*) – REAL (KIND=nag_wp) arrayOutput
Note: the second dimension of the array U must be at least N if TRANSF='T', and at least 1 otherwise.
On exit: if TRANSF='T' the n by n transformation matrix U*, otherwise U is not referenced.
18:   TOL – REAL (KIND=nag_wp)Input
On entry: the tolerance used to test for the singularity of Hi1/2. If 0.0TOL<m2×machine precision, then m2×machine precision is used instead. The inverse of the condition number of H1/2 is estimated by a call to F07TGF (DTRCON). If this estimate is less than TOL then H1/2 is assumed to be singular.
Suggested value: TOL=0.0.
Constraint: TOL0.0.
19:   IWK(M) – INTEGER arrayWorkspace
20:   WK(N+M×N+M+L) – REAL (KIND=nag_wp) arrayWorkspace
21:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to 0, -1​ or ​1. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1​ or ​1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is 0. When the value -1​ or ​1 is used it is essential to test the value of IFAIL on exit.
On exit: IFAIL=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6  Error Indicators and Warnings

If on entry IFAIL=0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF).
Errors or warnings detected by the routine:
IFAIL=1
On entry,TRANSF'T' or 'H',
orN<1,
orM<1,
orL<1,
orLDS<N,
orLDM<M,
orSTQ=.TRUE. and LDQ<1,
orSTQ=.FALSE. and LDQ<L,
orTOL<0.0.
IFAIL=2
The matrix Hi1/2 is singular.

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  Further Comments

For models with time-varying A,B and C, G13EAF can be used.
The initial estimate of the transformed state vector can be computed from the estimate of the original state vector X^10, say, by premultiplying it by U* as returned by G13EBF with TRANSF='T'; that is, X^10*=U*X^10. The estimate of the transformed state vector X^i+1i* can be computed from the previous value X^ii-1* by
X^i+1i*=U*AU*TX^ii-1*+U*AKiri
where
ri=Yi-CU*TX^ii- 1*
are the independent one-step prediction residuals for both the transformed and original model. The estimate of the original state vector can be computed from the transformed state vector as U*TX^1+1i*. The required matrix-vector multiplications can be performed by F06PAF (DGEMV).
If Wi and Vi are independent multivariate Normal variates then the log-likelihood for observations i=1,2,,t is given by
lθ = κ - 12 i=1 t l n detHi - 12 i=1 t Yi - Ci Xii-1 T Hi-1 Yi - Ci X ii-1
where κ is a constant.
The Cholesky factors of the covariance matrices can be computed using F07FDF (DPOTRF).
Note that the model
Xi+1=AXi+Wi, VarWi=Qi Yi=CXi+Vi, VarVi=Ri
can be specified either with B set to the identity matrix and STQ=.FALSE. and the matrix Q1/2 input in Q or with STQ=.TRUE. and B set to Q1/2.
The algorithm requires 16n3+n232m+l+2nm2+23p3 operations and is backward stable (see Verhaegen and van Dooren (1986)). The transformation to lower observer Hessenberg form requires On+mn2 operations.

9  Example

This example first inputs the number of updates to be computed and the problem sizes. The initial state vector and the Cholesky factor of the state covariance matrix are input followed by the model matrices A,B,C,R1/2 and optionally Q1/2 (the Cholesky factors of the covariance matrices being input). At the first update the matrices are transformed using the TRANSF='T' option and the initial value of the state vector is transformed. At each update the observed values are input and the residuals are computed and printed and the estimate of the transformed state vector, U^*Xii-1, and the deviance are updated. The deviance is -2×log-likelihood ignoring the constant. After the final update the estimate of the state vector is computed from the transformed state vector and the state covariance matrix is computed from S and these are printed along with the value of the deviance.
The data is for a two-dimensional time series to which a VARMA1,1 has been fitted. For the specification of a VARMA model as a state space model see the G13 Chapter Introduction. The means of the two series are included as additional states that do not change over time. The initial value of P, P0, is the solution to
P0=AP0AT+BQBT.

9.1  Program Text

Program Text (g13ebfe.f90)

9.2  Program Data

Program Data (g13ebfe.d)

9.3  Program Results

Program Results (g13ebfe.r)


G13EBF (PDF version)
G13 Chapter Contents
G13 Chapter Introduction
NAG Library Manual

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