g13eb:: Time Series Analysis (NAG Toolbox)

Description

The Kalman filter arises from the state space model given by

\begin{array}{l} X_{i + 1} = A X_{i} + B W_{i}, & Var (W_{i}) = Q_{i} \\ Y_{i} = C 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 matrices

A, B

and

C

are time invariant.

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

(1)

where

K_{i} = P_{i ∣ i} C^{T} {[C P_{i ∣ i} C^{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 {\hat{X}}_{i ∣ i} + D_{i} U_{i}

(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

(\begin{matrix} R_{i}^{1 / 2} & 0 & C S_{i} \\ 0 & B Q_{i}^{1 / 2} & A S_{i} \end{matrix}) U = (\begin{matrix} H_{i}^{1 / 2} & 0 & 0 \\ G_{i} & S_{i + 1} & 0 \end{matrix})

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

A K_{i} = G_{i} {(H_{i}^{1 / 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

[\begin{matrix} C U^{* T} \\ U^{*} A U^{* T} \end{matrix}]

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 ∣ 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

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

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

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,

{\hat{U}}^{*} X_{i ∣ i - 1}

, and the deviance are updated. The deviance is

- 2 \times 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.

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

On entry,	$transf \neq'T'$ or $'H'$ ,
or	$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$ .

NAG Toolbox: nag_tsa_multi_kalman_sqrt_invar (g13eb)

▸▿ Contents

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example