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:

\begin{array}{l} 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} \end{array}

where

X_{i}

is the state vector of length

n

at time

i

Y_{i}

is the observation vector of length

m

at time

i

, and

W_{i}

of length

l

and

V_{i}

of length

m

are the independent state noise and measurement noise respectively.

The estimate of

X_{i}

given observations

Y_{1}

Y_{i - 1}

is denoted by

{\hat{X}}_{i ∣ i - 1}

with state covariance matrix

Var ({\hat{X}}_{i ∣ i - 1}) = P_{i ∣ i - 1} = S_{i} S_{i}^{T}

, while the estimate of

X_{i}

given observations

Y_{1}

Y_{i}

is denoted by

{\hat{X}}_{i ∣ i}

with covariance matrix

Var ({\hat{X}}_{i ∣ i}) = P_{i ∣ i}

. The update of the estimate,

{\hat{X}}_{i ∣ i - 1}

, from time

i

to time

(i + 1)

, is computed in two stages. First, the measurement-update is given by

{\hat{X}}_{i ∣ i} = {\hat{X}}_{i ∣ i - 1} + K_{i} [Y_{i} - C_{i} {\hat{X}}_{i ∣ i - 1}]

(1)

and

P_{i ∣ i} = [I - K_{i} C_{i}] P_{i ∣ i - 1}

(2)

where

K_{i} = P_{i ∣ i - 1} C_{i}^{T} {[C_{i} P_{i ∣ i - 1} C_{i}^{T} + R_{i}]}^{- 1}

is the Kalman gain matrix. The second stage is the time-update for

X

which is given by

{\hat{X}}_{i + 1 ∣ i} = A_{i} {\hat{X}}_{i ∣ i} + D_{i} U_{i}

(3)

and

P_{i + 1 ∣ i} = A_{i} P_{i ∣ i} A_{i}^{T} + B_{i} Q_{i} B_{i}^{T}

(4)

where

D_{i} U_{i}

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

(\begin{matrix} R_{i}^{1 / 2} & C_{i} S_{i} & 0 \\ 0 & A_{i} S_{i} & B_{i} Q_{i}^{1 / 2} \end{matrix}) U = (\begin{matrix} H_{i}^{1 / 2} & 0 & 0 \\ G_{i} & S_{i + 1} & 0 \end{matrix})

(5)

where

U

is an orthogonal transformation triangularizing the left-hand pre-array to produce the right-hand post-array. 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} .

g13eaf requires the input of the lower triangular Cholesky factors of the noise covariance matrices

R_{i}^{1 / 2}

and, optionally,

Q_{i}^{1 / 2}

and the lower triangular Cholesky factor of the current state covariance matrix,

S_{i}

, and returns the product of the matrices

A_{i}

and

K_{i}

A_{i} K_{i}

, the Cholesky factor of the updated state covariance matrix

S_{i + 1}

and the matrix

H_{i}^{1 / 2}

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: $n$ – IntegerInput

On entry:

n

, the size of the state vector.

Constraint:

n \geq 1

2: $m$ – IntegerInput

On entry:

m

, the size of the observation vector.

Constraint:

m \geq 1

3: $l$ – IntegerInput

On entry:

l

, the dimension of the state noise.

Constraint:

l \geq 1

4: $a (lds, n)$ – Real (Kind=nag_wp) arrayInput

On entry: the state transition matrix,

A_{i}

5: $lds$ – 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:

lds \geq n

6: $b (lds, l)$ – Real (Kind=nag_wp) arrayInput

On entry: the noise coefficient matrix

B_{i}

7: $stq$ – LogicalInput

On entry: if

stq = .TRUE.

, the state noise covariance matrix

Q_{i}

is assumed to be the identity matrix. Otherwise the lower triangular Cholesky factor,

Q_{i}^{1 / 2}

, must be provided in q.

8: $q (ldq *)$ – Real (Kind=nag_wp) arrayInput

Note: the second dimension of the array q must be at least

l

stq = .FALSE.

On entry: if

stq = .FALSE.

, q must contain the lower triangular Cholesky factor of the state noise covariance matrix,

Q_{i}^{1 / 2}

. Otherwise q is not referenced.

9: $ldq$ – IntegerInput

On entry: the first dimension of the array q as declared in the (sub)program from which g13eaf is called.

Constraints:

if $stq = .FALSE.$ , $ldq \geq l$ ;
otherwise $ldq \geq 1$ .

10: $c (ldm, n)$ – Real (Kind=nag_wp) arrayInput

On entry: the measurement coefficient matrix,

C_{i}

11: $ldm$ – 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:

ldm \geq m

12: $r (ldm, m)$ – Real (Kind=nag_wp) arrayInput

On entry: the lower triangular Cholesky factor of the measurement noise covariance matrix

R_{i}^{1 / 2}

13: $s (lds, n)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: the lower triangular Cholesky factor of the state covariance matrix,

S_{i}

On exit: the lower triangular Cholesky factor of the state covariance matrix,

S_{i + 1}

14: $k (lds, m)$ – Real (Kind=nag_wp) arrayOutput

On exit: the Kalman gain matrix,

K_{i}

, premultiplied by the state transition matrix,

A_{i}

A_{i} K_{i}

15: $h (ldm, m)$ – Real (Kind=nag_wp) arrayOutput

On exit: the lower triangular matrix

H_{i}^{1 / 2}

16: $tol$ – Real (Kind=nag_wp)Input

On entry: the tolerance used to test for the singularity of

H_{i}^{1 / 2}

. If

0.0 \leq tol < m^{2} \times machine precision

, then

m^{2} \times machine precision

is used instead. The inverse of the condition number of

H^{1 / 2}

is estimated by a call to f07tgf (dtrcon). If this estimate is less than tol then

H^{1 / 2}

is assumed to be singular.

Suggested value:

tol = 0.0

Constraint:

tol \geq 0.0

17: $iwk (m)$ – Integer arrayWorkspace

18: $wk ((n + m) \times (n + m + l))$ – Real (Kind=nag_wp) arrayWorkspace

19: $ifail$ – IntegerInput/Output

On entry: ifail must be set to

0

- 1 ​ or ​ 1

. 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

- 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 argument, 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

- 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,	$n < 1$ ,
or	$m < 1$ ,
or	$l < 1$ ,
or	$lds < n$ ,
or	$ldm < m$ ,
or	$stq = .TRUE.$ and $ldq < 1$ ,
or	$stq = .FALSE.$ and $ldq < l$ ,
or	$tol < 0.0$ .

$ifail = 2$: The matrix $H_{i}^{1 / 2}$ is singular.

$ifail = - 99$: 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.

$ifail = - 399$: 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.

$ifail = - 999$: 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.

9

Further Comments

For models with time-invariant

A, B

and

C

, g13ebf can be used.

The estimate of the state vector

{\hat{X}}_{i + 1 ∣ i}

can be computed from

{\hat{X}}_{i ∣ i - 1}

{\hat{X}}_{i + 1 ∣ i} = A_{i} {\hat{X}}_{i ∣ i - 1} + A K_{i} r_{i}

where

r_{i} = Y_{i} - C_{i} {\hat{X}}_{i ∣ i - 1}

are the independent one step prediction residuals. The required matrix-vector multiplications can be performed by f06paf (dgemv).

W_{i}

and

V_{i}

are independent multivariate Normal variates then the log-likelihood for observations

i = 1, 2, \dots, t

is given by

l (θ) = κ - \frac{1}{2} \sum_{i = 1}^{t} l n (\det (H_{i})) - \frac{1}{2} \sum_{i = 1}^{t} {(Y_{i} - C_{i} X_{i ∣ i - 1})}^{T} H_{i}^{- 1} (Y_{i} - C_{i} X_{i ∣ i - 1})

where

κ

is a constant.

The Cholesky factors of the covariance matrices can be computed using f07fdf (dpotrf).

Note that the model

\begin{array}{l} X_{i + 1} = A_{i} X_{i} + W_{i}, & Var (W_{i}) = Q_{i} \\ Y_{i} = C_{i} X_{i} + V_{i}, & Var (V_{i}) = R_{i} \end{array}

can be specified either with b set to the identity matrix and

stq = .FALSE.

and the matrix

Q^{1 / 2}

input in q or with

stq = .TRUE.

and b set to

Q^{1 / 2}

The algorithm requires

\frac{7}{6} n^{3} + n^{2} (\frac{5}{2} m + l) + n (\frac{1}{2} l^{2} + m^{2})

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

A_{i}, B_{i}, C_{i}, R_{i}

and optionally

Q_{i}

. 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,

{\hat{X}}_{i ∣ i - 1}

, and the deviance are updated. The deviance is

- 2 \times

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

(1, 1)

has been fitted. For the specification of a VARMA model as a state space model see the G13 Chapter Introduction. The initial value of

P

P_{0}

, is the solution to

P_{0} = A_{1} P_{0} A_{1}^{T} + B_{1} Q_{1} B_{1}^{T} .

For convenience, the mean of each series is input before the first update and subtracted from the observations before the measurement update is computed.

NAG Library Routine Document

g13eaf (multi_kalman_sqrt_var)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results