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

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_tsa_multi_kalman_sqrt_invar (g13eb)

## Purpose

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

## Syntax

[a, b, c, s, k, h, u, ifail] = g13eb(transf, a, b, stq, q, c, r, s, 'n', n, 'm', m, 'l', l, 'tol', tol)
[a, b, c, s, k, h, u, ifail] = nag_tsa_multi_kalman_sqrt_invar(transf, 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=AXi+BWi, VarWi=Qi Yi=CXi+Vi, VarVi=Ri$
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 matrices $A,B$ and $C$ are time invariant.
The estimate of ${X}_{i}$ given observations ${Y}_{1}$ to ${Y}_{i-1}$ is denoted by ${\stackrel{^}{X}}_{i\mid i-1}$ with state covariance matrix $\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 ${X}_{i}$ given observations ${Y}_{1}$ to ${Y}_{i}$ is denoted by ${\stackrel{^}{X}}_{i\mid i}$ with covariance matrix $\mathrm{Var}\left({\stackrel{^}{X}}_{i\mid i}\right)={P}_{i\mid i}$. The update of the estimate, ${\stackrel{^}{X}}_{i\mid i-1}$, from time $i$ to time $\left(i+1\right)$ is computed in two stages. First, the measurement-update is given by
 $X^i∣i=X^i∣i-1+KiYi-CX^i∣i-1$ (1)
where ${K}_{i}={P}_{i\mid i}{C}^{\mathrm{T}}{\left[C{P}_{i\mid i}{C}^{\mathrm{T}}+{R}_{i}\right]}^{-1}$ is the Kalman gain matrix. The second stage is the time-update for $X$, which is given by
 $X^i+1∣i=AX^i∣i+DiUi$ (2)
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
 $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 ${K}_{i}$ and ${G}_{i}$ 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 $\left(A,C\right)$ 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 nag_tsa_multi_kalman_sqrt_invar (g13eb) will, optionally, compute them. nag_tsa_multi_kalman_sqrt_invar (g13eb) returns transformed matrices ${U}^{*}A{U}^{*T}$, ${U}^{*}B$, $C{U}^{*T}$ and ${U}^{*}A{K}_{i}$, the Cholesky factor of the updated transformed state covariance matrix ${S}_{i+1}^{*}$ (where ${U}^{*}{P}_{i+1\mid i}{U}^{*T}={S}_{i+1}^{*}{S}_{i+1}^{*T}$) and the matrix ${H}_{i}^{1/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 ${Q}_{i}$ and ${R}_{i}$ can be time-varying.

## 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:     $\mathrm{transf}$ – string (length ≥ 1)
Indicates whether to transform the input matrix pair $\left(A,C\right)$ to lower observer Hessenberg form. The transformation will only be required on the first call to nag_tsa_multi_kalman_sqrt_invar (g13eb).
${\mathbf{transf}}=\text{'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 Description.
${\mathbf{transf}}=\text{'H'}$
The matrices in arrays a, c and b should be as returned from a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with ${\mathbf{transf}}=\text{'T'}$.
Constraint: ${\mathbf{transf}}=\text{'T'}$ or $\text{'H'}$.
2:     $\mathrm{a}\left(\mathit{lds},{\mathbf{n}}\right)$ – double array
lds, the first dimension of the array, must satisfy the constraint $\mathit{lds}\ge {\mathbf{n}}$.
If ${\mathbf{transf}}=\text{'T'}$, the state transition matrix, $A$.
If ${\mathbf{transf}}=\text{'H'}$, the transformed matrix as returned by a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with ${\mathbf{transf}}=\text{'T'}$.
3:     $\mathrm{b}\left(\mathit{lds},{\mathbf{l}}\right)$ – double array
lds, the first dimension of the array, must satisfy the constraint $\mathit{lds}\ge {\mathbf{n}}$.
If ${\mathbf{transf}}=\text{'T'}$, the noise coefficient matrix $B$.
If ${\mathbf{transf}}=\text{'H'}$, the transformed matrix as returned by a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with ${\mathbf{transf}}=\text{'T'}$.
4:     $\mathrm{stq}$ – logical scalar
If ${\mathbf{stq}}=\mathit{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.
5:     $\mathrm{q}\left(\mathit{ldq},:\right)$ – double array
The first dimension, $\mathit{ldq}$, of the array q must satisfy
• if ${\mathbf{stq}}=\mathit{false}$, $\mathit{ldq}\ge {\mathbf{l}}$;
• otherwise $\mathit{ldq}\ge 1$.
The second dimension of the array q must be at least ${\mathbf{l}}$ if ${\mathbf{stq}}=\mathit{false}$ and at least $1$ if ${\mathbf{stq}}=\mathit{true}$.
If ${\mathbf{stq}}=\mathit{false}$, q must contain the lower triangular Cholesky factor of the state noise covariance matrix, ${Q}_{i}^{1/2}$. Otherwise q is not referenced.
6:     $\mathrm{c}\left(\mathit{ldm},{\mathbf{n}}\right)$ – double array
ldm, the first dimension of the array, must satisfy the constraint $\mathit{ldm}\ge {\mathbf{m}}$.
If ${\mathbf{transf}}=\text{'T'}$, the measurement coefficient matrix, $C$.
If ${\mathbf{transf}}=\text{'H'}$, the transformed matrix as returned by a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with ${\mathbf{transf}}=\text{'T'}$.
7:     $\mathrm{r}\left(\mathit{ldm},{\mathbf{m}}\right)$ – double array
ldm, the first dimension of the array, must satisfy the constraint $\mathit{ldm}\ge {\mathbf{m}}$.
The lower triangular Cholesky factor of the measurement noise covariance matrix ${R}_{i}^{1/2}$.
8:     $\mathrm{s}\left(\mathit{lds},{\mathbf{n}}\right)$ – double array
lds, the first dimension of the array, must satisfy the constraint $\mathit{lds}\ge {\mathbf{n}}$.
If ${\mathbf{transf}}=\text{'T'}$ the lower triangular Cholesky factor of the state covariance matrix, ${S}_{i}$.
If ${\mathbf{transf}}=\text{'H'}$ the lower triangular Cholesky factor of the covariance matrix of the transformed state vector ${S}_{i}^{*}$ as returned from a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with ${\mathbf{transf}}=\text{'T'}$.

### Optional Input Parameters

1:     $\mathrm{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$, the size of the state vector.
Constraint: ${\mathbf{n}}\ge 1$.
2:     $\mathrm{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$, the size of the observation vector.
Constraint: ${\mathbf{m}}\ge 1$.
3:     $\mathrm{l}$int64int32nag_int scalar
Default: the second dimension of the array b.
$l$, the dimension of the state noise.
Constraint: ${\mathbf{l}}\ge 1$.
4:     $\mathrm{tol}$ – double scalar
Default: $0.0$
The tolerance used to test for the singularity of ${H}_{i}^{1/2}$. If , then  is used instead. The inverse of the condition number of ${H}^{1/2}$ is estimated by a call to nag_lapack_dtrcon (f07tg). If this estimate is less than tol then ${H}^{1/2}$ is assumed to be singular.
Constraint: ${\mathbf{tol}}\ge 0.0$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lds},{\mathbf{n}}\right)$ – double array
If ${\mathbf{transf}}=\text{'T'}$, the transformed matrix, ${U}^{*}A{U}^{*T}$, otherwise a is unchanged.
2:     $\mathrm{b}\left(\mathit{lds},{\mathbf{l}}\right)$ – double array
If ${\mathbf{transf}}=\text{'T'}$, the transformed matrix, ${U}^{*}B$, otherwise b is unchanged.
3:     $\mathrm{c}\left(\mathit{ldm},{\mathbf{n}}\right)$ – double array
If ${\mathbf{transf}}=\text{'T'}$, the transformed matrix, $C{U}^{*T}$, otherwise c is unchanged.
4:     $\mathrm{s}\left(\mathit{lds},{\mathbf{n}}\right)$ – double array
The lower triangular Cholesky factor of the transformed state covariance matrix, ${S}_{i+1}^{*}$.
5:     $\mathrm{k}\left(\mathit{lds},{\mathbf{m}}\right)$ – double array
The Kalman gain matrix for the transformed state vector premultiplied by the state transformed transition matrix, ${U}^{*}A{K}_{i}$.
6:     $\mathrm{h}\left(\mathit{ldm},{\mathbf{m}}\right)$ – double array
The lower triangular matrix ${H}_{i}^{1/2}$.
7:     $\mathrm{u}\left(\mathit{lds},:\right)$ – double array
The first dimension of the array u will be ${\mathbf{n}}$.
The second dimension of the array u will be ${\mathbf{n}}$ if ${\mathbf{transf}}=\text{'T'}$ and $1$ otherwise.
If ${\mathbf{transf}}=\text{'T'}$ the $n$ by $n$ transformation matrix ${U}^{*}$, otherwise u is not referenced.
8:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{transf}}\ne \text{'T'}$ or $\text{'H'}$, or ${\mathbf{n}}<1$, or ${\mathbf{m}}<1$, or ${\mathbf{l}}<1$, or $\mathit{lds}<{\mathbf{n}}$, or $\mathit{ldm}<{\mathbf{m}}$, or ${\mathbf{stq}}=\mathit{true}$ and $\mathit{ldq}<1$, or ${\mathbf{stq}}=\mathit{false}$ and $\mathit{ldq}<{\mathbf{l}}$, or ${\mathbf{tol}}<0.0$.
${\mathbf{ifail}}=2$
The matrix ${H}_{i}^{1/2}$ is singular.
${\mathbf{ifail}}=-99$
An unexpected error has been triggered by this routine. Please contact NAG.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## 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-varying $A,B$ and $C$, nag_tsa_multi_kalman_sqrt_var (g13ea) can be used.
If ${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θ = κ - 12 ∑ i=1 t l n detHi - 12 ∑ i=1 t Yi - Ci Xi∣i-1 T Hi-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=AXi+Wi, VarWi=Qi Yi=CXi+Vi, VarVi=Ri$
can be specified either with b set to the identity matrix and ${\mathbf{stq}}=\mathit{false}$ and the matrix ${Q}^{1/2}$ input in q or with ${\mathbf{stq}}=\mathit{true}$ and b set to ${Q}^{1/2}$.
The algorithm requires $\frac{1}{6}{n}^{3}+{n}^{2}\left(\frac{3}{2}m+l\right)+2n{m}^{2}+\frac{2}{3}{p}^{3}$ operations and is backward stable (see Verhaegen and van Dooren (1986)). The transformation to lower observer Hessenberg form requires $\mathit{O}\left(\left(n+m\right){n}^{2}\right)$ operations.

## 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,{R}^{1/2}$ and optionally ${Q}^{1/2}$ (the Cholesky factors of the covariance matrices being input). At the first update the matrices are transformed using the ${\mathbf{transf}}=\text{'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, ${\stackrel{^}{U}}^{*}{X}_{i\mid i-1}$, and the deviance are updated. The deviance is $-2×\text{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 VARMA$\left(1,1\right)$ 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$, ${P}_{0}$, is the solution to
 $P0=AP0AT+BQBT.$
```function g13eb_example

fprintf('g13eb example results\n\n');

% Constant matrices
stq = false;
s = [ 2.8648  0.0000  0.0000  0.0000  0.0000  0.0000;
0.7191  2.7290  0.0000  0.0000  0.0000  0.0000;
0.5169  0.2194  0.7810  0.0000  0.0000  0.0000;
0.1266  0.0449  0.1899  0.0098  0.0000  0.0000;
0.0000  0.0000  0.0000  0.0000  0.0000  0.0000;
0.0000  0.0000  0.0000  0.0000  0.0000  0.0000];
ax = [ 0.000;  0.000;  0.000;  0.000;  4.404;  7.991];
a  = [ 0.607 -0.033  1.000  0.000  0.000  0.000;
0.000  0.543  0.000  1.000  0.000  0.000;
0.000  0.000  0.000  0.000  0.000  0.000;
0.000  0.000  0.000  0.000  0.000  0.000;
0.000  0.000  0.000  0.000  1.000  0.000;
0.000  0.000  0.000  0.000  0.000  1.000];
b  = [ 1.000  0.000;
0.000  1.000;
0.543  0.125;
0.134  0.026;
0.000  0.000;
0.000  0.000];
c  = [ 1.000  0.000  0.000  0.000  1.000  0.000;
0.000  1.000  0.000  0.000  0.000  1.000];
r  = [ 0.000  0.000;
0.000  0.000];
q =  [ 1.612  0.000;
0.347  2.282];
n = size(s,1);
l = size(b,2);
m = size(r,2);

% data
ym = [-1.49  7.34; -1.62  6.35;  5.20  6.96;  6.23  8.54;  6.21  6.62;
5.86  4.97;  4.09  4.55;  3.18  4.81;  2.62  4.75;  1.49  4.76;
1.17 10.88;  0.85 10.01; -0.35 11.62;  0.24 10.36;  2.44  6.40;
2.58  6.24;  2.04  7.93;  0.40  4.04;  2.26  3.73;  3.34  5.60;
5.09  5.35;  5.00  6.81;  4.78  8.27;  4.11  7.68;  3.45  6.65;
1.65  6.08;  1.29 10.25;  4.09  9.14;  6.32 17.75;  7.50 13.30;
3.89  9.63;  1.58  6.80;  5.21  4.08;  5.25  5.06;  4.93  4.94;
7.38  6.65;  5.87  7.94;  5.81 10.76;  9.68 11.89;  9.07  5.85;
7.29  9.01;  7.84  7.50;  7.55 10.02;  7.32 10.38;  7.97  8.15;
7.76  8.37;  7.00 10.73;  8.35 12.14];
ny = size(ym,1);

% Loop over data
fprintf('      Residuals\n\n');
dev = 0;
for j = 1:ny
y = ym(j,:)';

if j == 1
[a, b, c, s, k, h, u, ifail] = ...
g13eb( ...
'T', a, b, stq, q, c, r, s);
x = u*ax;
else
[a, b, c, s, k, h, ~, ifail] = ...
g13eb( ...
'H', a, b, stq, q, c, r, s);

end

% time and measurement update
y = y - c*x;
x = a*x + k*y;

fprintf('%12.4f', y);
fprintf('\n');

% Update loglikelihood
y = h\y;
dev = dev + dot(y,y);
for i = 1:m
dev = dev + 2*log(h(i,i));
end
end

% Back-transform x
x = transpose(u)*x;
us = transpose(u)*s;
p = us*us';

fprintf('\nFinal x\n\n')
disp (x');

[ifail] = x04ca(...
'Lower','N', p, 'Final Value of P');

fprintf('\nDeviance = %13.4e\n', dev);

```
```g13eb example results

Residuals

-5.8940     -0.6510
-1.4710     -1.0407
5.1658      0.0447
-1.3281      0.4580
1.3653     -1.5066
-0.2337     -2.4192
-0.8685     -1.7065
-0.4624     -1.1519
-0.7510     -1.4218
-1.3526     -1.3335
-0.6707      4.8593
-1.7389      0.4138
-1.6376      2.7549
-0.6137      0.5463
0.9067     -2.8093
-0.8255     -0.9355
-0.7494      1.0247
-2.2922     -3.8441
1.8812     -1.7085
-0.7112     -0.2849
1.6747     -1.2400
-0.6619      0.0609
0.3271      1.0074
-0.8165     -0.5325
-0.2759     -1.0489
-1.9383     -1.1186
-0.3131      3.5855
1.3726     -0.1289
1.4153      8.9545
0.3672     -0.4126
-2.3659     -1.2823
-1.0130     -1.7306
3.2472     -3.0836
-1.1501     -1.1623
0.6855     -1.2751
2.3432      0.2570
-1.6892      0.3565
1.3871      3.0138
3.3840      2.1312
-0.5118     -4.7670
0.8569      2.3741
0.9558     -1.2209
0.6778      2.1993
0.4304      1.1393
1.4987     -1.2255
0.5361      0.1237
0.2649      2.4582
2.0095      2.5623

Final x

3.6698    2.5888    0.0000    0.0000    4.4040    7.9910

Final Value of P
1            2            3            4            5
1    2.5985E+00
2    5.5936E-01   5.3279E+00
3    1.4809E+00   9.6973E-01   9.2536E-01
4    3.6275E-01   2.1348E-01   2.2366E-01   5.4159E-02
5   -4.0547E-16  -8.7283E-17  -2.3108E-16  -5.6603E-17   9.6581E-32
6    4.4742E-17   9.6312E-18   2.5499E-17   6.2458E-18  -9.3231E-33

6
1
2
3
4
5
6    1.3378E-32

Deviance =    2.2287e+02
```

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015