NAG C Library Function Document

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

\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}

the estimate of

X_{i}

given observations

Y_{1}

Y_{i - 1}

is denoted by

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

with

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

. nag_kalman_sqrt_filt_cov_var (g13eac) performs one recursion of the square root covariance filter algorithm, summarised as follows:

\begin{matrix} (\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}) \\ (Pre-array) & (Post-array) \end{matrix}

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

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

(1)

where

K_{i}

is the Kalman gain matrix, whilst the time-update for

X

{\hat{X}}_{i + 1 ∣ i} = A_{i} {\hat{X}}_{i ∣ i} + 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

{A K}_{i}

, and the state covariance matrix

P_{i ∣ i - 1}

factorized as

P_{i ∣ i - 1} = S_{i} S_{i}^{T}

(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: $n$ – IntegerInput: On entry: the actual state dimension, $n$ , i.e., the order of the matrices $S_{i}$ and $A_{i}$ .

Constraint: $n \geq 1$ .
2: $m$ – IntegerInput: On entry: the actual input dimension, $m$ , i.e., the order of the matrix $Q_{i}^{1 / 2}$ .

Constraint: $m \geq 1$ .
3: $p$ – IntegerInput: On entry: the actual output dimension, $p$ , i.e., the order of the matrix $R_{i}^{1 / 2}$ .

Constraint: $p \geq 1$ .
4: $s [n \times tds]$ – doubleInput/Output: Note: the $(i, j)$ th element of the matrix $S$ is stored in $s [(i - 1) \times 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_{i ∣ i - 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 + 1 ∣ i}$ .
5: $tds$ – IntegerInput: On entry: the stride separating matrix column elements in the array s.

Constraint: $tds \geq n$ .
6: $a [n \times tda]$ – const doubleInput: Note: the $(i, j)$ th element of the matrix $A$ is stored in $a [(i - 1) \times 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: $tda$ – IntegerInput: On entry: the stride separating matrix column elements in the array a.

Constraint: $tda \geq n$ .
8: $b [n \times tdb]$ – const doubleInput: Note: the $(i, j)$ th element of the matrix $B$ is stored in $b [(i - 1) \times 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: $tdb$ – IntegerInput: On entry: the stride separating matrix column elements in the array b.

Constraint: $tdb \geq m$ .
10: $q [m \times tdq]$ – const doubleInput: Note: the $(i, j)$ th element of the matrix $Q$ is stored in $q [(i - 1) \times 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: $tdq$ – IntegerInput: On entry: the stride separating matrix column elements in the array q.

Constraint: $tdq \geq m$ if q is defined.
12: $c [p \times tdc]$ – const doubleInput: Note: the $(i, j)$ th element of the matrix $C$ is stored in $c [(i - 1) \times 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: $tdc$ – IntegerInput: On entry: the stride separating matrix column elements in the array c.

Constraint: $tdc \geq n$ .
14: $r [p \times tdr]$ – const doubleInput: Note: the $(i, j)$ th element of the matrix $R$ is stored in $r [(i - 1) \times 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: $tdr$ – IntegerInput: On entry: the stride separating matrix column elements in the array r.

Constraint: $tdr \geq p$ .
16: $k [n \times tdk]$ – doubleOutput: Note: the $(i, j)$ th element of the matrix $K$ is stored in $k [(i - 1) \times tdk + j - 1]$ .

On exit: if k is not NULL then the leading $n$ by $p$ part of k contains ${A K}_{i}$ , the product of the Kalman filter gain matrix $K_{i}$ with the state transition matrix $A_{i}$ . If $A K_{i}$ is not required then k must be set to NULL.
17: $tdk$ – IntegerInput: On entry: the stride separating matrix column elements in the array k.

Constraint: if k is not NULL, $tdk \geq p$
18: $h [p \times tdh]$ – doubleOutput: Note: the $(i, j)$ th element of the matrix $H$ is stored in $h [(i - 1) \times 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: $tdh$ – IntegerInput: On entry: the stride separating matrix column elements in the array h.

Constraint: if both k and h are not NULL, $tdh \geq p$
20: $tol$ – doubleInput: 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: $fail$ – 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 $tda = 〈value〉$ while $n = 〈value〉$ . These arguments must satisfy $tda \geq n$ .

On entry $tdb = 〈value〉$ while $m = 〈value〉$ . These arguments must satisfy $tdb \geq m$ .

On entry $tdc = 〈value〉$ while $n = 〈value〉$ . These arguments must satisfy $tdc \geq n$ .

On entry $tdh = 〈value〉$ while $p = 〈value〉$ . These arguments must satisfy $tdh \geq p$ .

On entry $tdk = 〈value〉$ while $p = 〈value〉$ . These arguments must satisfy $tdk \geq p$ .

On entry $tdq = 〈value〉$ while $m = 〈value〉$ . These arguments must satisfy $tdq \geq m$ .

On entry $tdr = 〈value〉$ while $p = 〈value〉$ . These arguments must satisfy $tdr \geq p$ .

On entry, $tds = 〈value〉$ while $n = 〈value〉$ . These arguments must satisfy $tds \geq n$ .
NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_INT_ARG_LT: On entry, $m = 〈value〉$ .
Constraint: $m \geq 1$ .

On entry, $n = 〈value〉$ .
Constraint: $n \geq 1$ .

On entry, $p = 〈value〉$ .
Constraint: $p \geq 1$ .
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

nag_kalman_sqrt_filt_cov_var (g13eac) is not threaded in any implementation.

9

Further Comments

The algorithm requires

\frac{7}{6} n^{3} + n^{2} (\frac{5}{2} p + m) + n (\frac{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 (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

y_{k} = ϕ y_{k - 1} + ε_{k} - θ ε_{k - 1}

This has the following state space representation (Harvey and Phillips (1979))

\begin{array}{l} x_{k} = (\begin{matrix} ϕ & 1 \\ 0 & 0 \end{matrix}) x_{k - 1} + (\begin{matrix} 1 \\ - θ \end{matrix}) ε_{k} \\ y_{k} = (\begin{matrix} 1 & 0 \end{matrix}) x_{k} \end{array}

where the state vector

x_{k} = (\begin{matrix} y_{k} \\ - θ ε_{k} \end{matrix})

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

\begin{array}{l} A_{k} = (\begin{matrix} ϕ & 1 \\ 0 & 0 \end{matrix}), B_{k} = (\begin{matrix} 1 \\ - θ \end{matrix}) \\ C_{k} = (\begin{matrix} 1 & 0 \end{matrix}), Q_{k} = 0 ​ and ​ R_{k} = 1 . \end{array}

The initial estimates for the state vector,

x_{1 ∣ 0}

, and state covariance matrix,

P_{1 ∣ 0}

, are:

x_{1 ∣ 0} = E [x_{k}] = 0 ​ and ​ P_{1 ∣ 0} = E [x_{k} x_{k}^{T}] = (\begin{matrix} E [y_{k} y_{k}] & - θ E [y_{k} ε_{k}] \\ - θ E [y_{k} ε_{k}] & θ^{2} E [ε_{k} ε_{k}] \end{matrix}) .

Since

E [y_{k} y_{k}] = γ_{\circ} = \frac{(1 + θ^{2} - 2 ϕ θ) σ^{2}}{(1 - ϕ^{2})}

(Wei (1990))

P_{1 ∣ 0} = (\begin{matrix} γ_{\circ} & - θ \\ - θ & θ^{2} \end{matrix}) .

Using

P_{1 ∣ 0} = S_{1 ∣ 0} S_{1 ∣ 0}^{T}

gives an initial Cholesky ‘square root’ of

S_{1 ∣ 0} = (\begin{matrix} \sqrt{γ_{\circ}} & 0 \\ \frac{- θ}{\sqrt{γ_{\circ}}} & θ \sqrt{\frac{γ_{\circ} - 1}{γ_{\circ}}} \end{matrix}) .

10.1

Program Text

Program Text (g13eace.c)

10.2

Program Data

None.

10.3

Program Results

Program Results (g13eace.r)

NAG C Library Manual, Mark 26

NAG Library Manual, Mark 26

NAG AD Library Manual, Mark 26

g13 Chapter Contents

g13 Chapter Introduction

NAG C Library Function Document

nag_kalman_sqrt_filt_cov_var (g13eac)

▸▿ 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