Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_tsa_multi_kalman_sqrt_var (g13ea)

Purpose

nag_tsa_multi_kalman_sqrt_var (g13ea) performs a combined measurement and time update of one iteration of the time-varying Kalman filter using a square root covariance filter.

Syntax

[s, k, h, ifail] = g13ea(a, b, stq, q, c, r, s, 'n', n, 'm', m, 'l', l, 'tol', tol)
[s, k, h, ifail] = nag_tsa_multi_kalman_sqrt_var(a, b, stq, q, c, r, s, 'n', n, 'm', m, 'l', l, 'tol', tol)

Description

The Kalman filter arises from the state space model given by:
 Xi + 1 = AiXi + BiWi, Var(Wi) = Qi Yi = CiXi + Vi, Var(Vi) = Ri
$Xi+1=AiXi+BiWi, Var(Wi)=Qi Yi=CiXi+Vi, Var(Vi)=Ri$
where Xi${X}_{i}$ is the state vector of length n$n$ at time i$i$, Yi${Y}_{i}$ is the observation vector of length m$m$ at time i$i$, and Wi${W}_{i}$ of length l$l$ and Vi${V}_{i}$ of length m$m$ are the independent state noise and measurement noise respectively.
The estimate of Xi${X}_{i}$ given observations Y1${Y}_{1}$ to Yi1${Y}_{i-1}$ is denoted by ii1${\stackrel{^}{X}}_{i\mid i-1}$ with state covariance matrix Var(ii1) = Pii1 = SiSiT$\mathrm{Var}\left({\stackrel{^}{X}}_{i\mid i-1}\right)={P}_{i\mid i-1}={S}_{i}{S}_{i}^{\mathrm{T}}$, while the estimate of Xi${X}_{i}$ given observations Y1${Y}_{1}$ to Yi${Y}_{i}$ is denoted by ii${\stackrel{^}{X}}_{i\mid i}$ with covariance matrix Var(ii) = Pii$\mathrm{Var}\left({\stackrel{^}{X}}_{i\mid i}\right)={P}_{i\mid i}$. The update of the estimate, ii1${\stackrel{^}{X}}_{i\mid i-1}$, from time i$i$ to time (i + 1)$\left(i+1\right)$, is computed in two stages. First, the measurement-update is given by
 X̂i ∣ i = X̂i ∣ i − 1 + Ki[Yi − CiX̂i ∣ i − 1] $X^i∣i=X^i∣i-1+Ki[Yi-CiX^i∣i-1]$ (1)
and
 Pi ∣ i = [I − KiCi]Pi ∣ i − 1 $Pi∣i=[I-KiCi]Pi∣i-1$ (2)
where Ki = Pii1 CiT [CiPii1CiT + Ri]1${K}_{i}={P}_{i\mid i-1}{C}_{i}^{\mathrm{T}}{\left[{C}_{i}{P}_{i\mid i-1}{C}_{i}^{\mathrm{T}}+{R}_{i}\right]}^{-1}$ is the Kalman gain matrix. The second stage is the time-update for X$X$ which is given by
 X̂i + 1 ∣ i = AiX̂i ∣ i + DiUi $X^i+1∣i=AiX^i∣i+DiUi$ (3)
and
 Pi + 1 ∣ i = AiPi ∣ i AiT + BiQi BiT $Pi+1∣i=AiPi∣i AiT +BiQi BiT$ (4)
where DiUi${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
 Ri1 / 2 CiSi 0 0 AiSi BiQi1 / 2
U =
 Hi1 / 2 0 0 Gi Si + 1 0
$Ri1/2 CiSi 0 0 AiSi BiQi1/2 U= Hi1/2 0 0 Gi Si+1 0$
(5)
where U$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, Ki${K}_{i}$, and Gi${G}_{i}$ is given by
 AiKi = Gi (Hi1 / 2) − 1. $AiKi=Gi (Hi1/2) -1.$
nag_tsa_multi_kalman_sqrt_var (g13ea) requires the input of the lower triangular Cholesky factors of the noise covariance matrices Ri1 / 2${R}_{i}^{1/2}$ and, optionally, Qi1 / 2${Q}_{i}^{1/2}$ and the lower triangular Cholesky factor of the current state covariance matrix, Si${S}_{i}$, and returns the product of the matrices Ai${A}_{i}$ and Ki${K}_{i}$, AiKi${A}_{i}{K}_{i}$, the Cholesky factor of the updated state covariance matrix Si + 1${S}_{i+1}$ and the matrix Hi1 / 2${H}_{i}^{1/2}$ used in the computation of the likelihood for the model.

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

Parameters

Compulsory Input Parameters

1:     a(lds,n) – double array
lds, the first dimension of the array, must satisfy the constraint ldsn$\mathit{lds}\ge {\mathbf{n}}$.
The state transition matrix, Ai${A}_{i}$.
2:     b(lds,l) – double array
lds, the first dimension of the array, must satisfy the constraint ldsn$\mathit{lds}\ge {\mathbf{n}}$.
The noise coefficient matrix Bi${B}_{i}$.
3:     stq – logical scalar
If stq = true${\mathbf{stq}}=\mathbf{true}$, the state noise covariance matrix Qi${Q}_{i}$ is assumed to be the identity matrix. Otherwise the lower triangular Cholesky factor, Qi1 / 2${Q}_{i}^{1/2}$, must be provided in q.
4:     q(ldq, : $:$) – double array
The first dimension, ldq, of the array q must satisfy
• if stq = false${\mathbf{stq}}=\mathbf{false}$, ldql$\mathit{ldq}\ge {\mathbf{l}}$;
• otherwise ldq1$\mathit{ldq}\ge 1$.
The second dimension of the array must be at least l${\mathbf{l}}$ if stq = false${\mathbf{stq}}=\mathbf{false}$ and at least 1$1$ if stq = true${\mathbf{stq}}=\mathbf{true}$
If stq = false${\mathbf{stq}}=\mathbf{false}$, q must contain the lower triangular Cholesky factor of the state noise covariance matrix, Qi1 / 2${Q}_{i}^{1/2}$. Otherwise q is not referenced.
5:     c(ldm,n) – double array
ldm, the first dimension of the array, must satisfy the constraint ldmm$\mathit{ldm}\ge {\mathbf{m}}$.
The measurement coefficient matrix, Ci${C}_{i}$.
6:     r(ldm,m) – double array
ldm, the first dimension of the array, must satisfy the constraint ldmm$\mathit{ldm}\ge {\mathbf{m}}$.
The lower triangular Cholesky factor of the measurement noise covariance matrix Ri1 / 2${R}_{i}^{1/2}$.
7:     s(lds,n) – double array
lds, the first dimension of the array, must satisfy the constraint ldsn$\mathit{lds}\ge {\mathbf{n}}$.
The lower triangular Cholesky factor of the state covariance matrix, Si${S}_{i}$.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the arrays a, b, s and the second dimension of the arrays a, c, s. (An error is raised if these dimensions are not equal.)
n$n$, the size of the state vector.
Constraint: n1${\mathbf{n}}\ge 1$.
2:     m – int64int32nag_int scalar
Default: The first dimension of the arrays c, r and the second dimension of the array r. (An error is raised if these dimensions are not equal.)
m$m$, the size of the observation vector.
Constraint: m1${\mathbf{m}}\ge 1$.
3:     l – int64int32nag_int scalar
Default: The second dimension of the array b.
l$l$, the dimension of the state noise.
Constraint: l1${\mathbf{l}}\ge 1$.
4:     tol – double scalar
The tolerance used to test for the singularity of Hi1 / 2${H}_{i}^{1/2}$. If 0.0tol < m2 × machine precision, then m2 × machine precision is used instead. The inverse of the condition number of H1 / 2${H}^{1/2}$ is estimated by a call to nag_lapack_dtrcon (f07tg). If this estimate is less than tol then H1 / 2${H}^{1/2}$ is assumed to be singular.
Default: 0.0$0.0$
Constraint: tol0.0${\mathbf{tol}}\ge 0.0$.

Input Parameters Omitted from the MATLAB Interface

lds ldq ldm iwk wk

Output Parameters

1:     s(lds,n) – double array
ldsn$\mathit{lds}\ge {\mathbf{n}}$.
The lower triangular Cholesky factor of the state covariance matrix, Si + 1${S}_{i+1}$.
2:     k(lds,m) – double array
ldsn$\mathit{lds}\ge {\mathbf{n}}$.
The Kalman gain matrix, Ki${K}_{i}$, premultiplied by the state transition matrix, Ai${A}_{i}$, AiKi${A}_{i}{K}_{i}$.
3:     h(ldm,m) – double array
ldmm$\mathit{ldm}\ge {\mathbf{m}}$.
The lower triangular matrix Hi1 / 2${H}_{i}^{1/2}$.
4:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, n < 1${\mathbf{n}}<1$, or m < 1${\mathbf{m}}<1$, or l < 1${\mathbf{l}}<1$, or lds < n$\mathit{lds}<{\mathbf{n}}$, or ldm < m$\mathit{ldm}<{\mathbf{m}}$, or stq = true${\mathbf{stq}}=\mathbf{true}$ and ldq < 1$\mathit{ldq}<1$, or stq = false${\mathbf{stq}}=\mathbf{false}$ and ldq < l$\mathit{ldq}<{\mathbf{l}}$, or tol < 0.0${\mathbf{tol}}<0.0$.
ifail = 2${\mathbf{ifail}}=2$
The matrix Hi1 / 2${H}_{i}^{1/2}$ is singular.

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.

For models with time-invariant A,B$A,B$ and C$C$, nag_tsa_multi_kalman_sqrt_invar (g13eb) can be used.
If Wi${W}_{i}$ and Vi${V}_{i}$ are independent multivariate Normal variates then the log-likelihood for observations i = 1,2,,t$i=1,2,\dots ,t$ is given by
 t t l(θ) = κ − (1/2) ∑ l n (det(Hi)) − (1/2) ∑ (Yi − CiXi ∣ i − 1)THi − 1(Yi − CiXi ∣ i − 1) i = 1 i = 1
$l(θ) = κ - 12 ∑ i=1 t l n ( det(Hi) ) - 12 ∑ i=1 t ( Yi - Ci X i∣i-1 )T H i -1 ( Yi - Ci X i∣i-1 )$
where κ$\kappa$ is a constant.
The Cholesky factors of the covariance matrices can be computed using nag_lapack_dpotrf (f07fd).
Note that the model
 Xi + 1 = AiXi + Wi, Var(Wi) = Qi Yi = CiXi + Vi, Var(Vi) = Ri
$Xi+1=AiXi+Wi, Var(Wi)=Qi Yi=CiXi+Vi, Var(Vi)=Ri$
can be specified either with b set to the identity matrix and stq = false${\mathbf{stq}}=\mathbf{false}$ and the matrix Q1 / 2${Q}^{1/2}$ input in q or with stq = true${\mathbf{stq}}=\mathbf{true}$ and b set to Q1 / 2${Q}^{1/2}$.
The algorithm requires (7/6)n3 + n2 ((5/2)m + l) + n((1/2)l2 + m2)$\frac{7}{6}{n}^{3}+{n}^{2}\left(\frac{5}{2}m+l\right)+n\left(\frac{1}{2}{l}^{2}+{m}^{2}\right)$ operations and is backward stable (see Verhaegen and van Dooren (1986)).

Example

```function nag_tsa_multi_kalman_sqrt_var_example
a = [0.607, -0.033, 1, 0;
0, 0.543, 0, 1;
0, 0, 0, 0;
0, 0, 0, 0];
b = [1, 0;
0, 1;
0.543, 0.125;
0.134, 0.026];
stq = false;
q = [1.611831256676703, 0.56;
0.3474309098302363, 2.282387294675146];
c = [1, 0, 0, 0;
0, 1, 0, 0];
r = [0, 0;
0, 0];
s = [2.864751298105998, 2.0599, 1.4807, 0.3627;
0.7190502021456042, 2.729004728246978, 0.9703, 0.2136;
0.516868602513227, 0.2193640489823183, 0.7810417797724439, 0.2236;
0.1266078490791838, 0.04491109863526066, 0.1898854854878202, 0.009846226280729009];
[sOut, k, h, ifail] = nag_tsa_multi_kalman_sqrt_var(a, b, stq, q, c, r, s)
```
```

sOut =

-1.7911    2.0599    1.4807    0.3627
-0.3955   -2.2825    0.9703    0.2136
-0.8267   -0.2819    0.4030    0.2236
-0.2025   -0.0585    0.0986   -0.0000

k =

0.7672    0.0474
0.0401    0.5595
0         0
0         0

h =

2.8648         0
0.7191    2.7290

ifail =

0

```
```function g13ea_example
a = [0.607, -0.033, 1, 0;
0, 0.543, 0, 1;
0, 0, 0, 0;
0, 0, 0, 0];
b = [1, 0;
0, 1;
0.543, 0.125;
0.134, 0.026];
stq = false;
q = [1.611831256676703, 0.56;
0.3474309098302363, 2.282387294675146];
c = [1, 0, 0, 0;
0, 1, 0, 0];
r = [0, 0;
0, 0];
s = [2.864751298105998, 2.0599, 1.4807, 0.3627;
0.7190502021456042, 2.729004728246978, 0.9703, 0.2136;
0.516868602513227, 0.2193640489823183, 0.7810417797724439, 0.2236;
0.1266078490791838, 0.04491109863526066, 0.1898854854878202, 0.009846226280729009];
[sOut, k, h, ifail] = g13ea(a, b, stq, q, c, r, s)
```
```

sOut =

-1.7911    2.0599    1.4807    0.3627
-0.3955   -2.2825    0.9703    0.2136
-0.8267   -0.2819    0.4030    0.2236
-0.2025   -0.0585    0.0986   -0.0000

k =

0.7672    0.0474
0.0401    0.5595
0         0
0         0

h =

2.8648         0
0.7191    2.7290

ifail =

0

```