JMSLTM Numerical Library 4.0

com.imsl.stat
Class KalmanFilter

java.lang.Object
  extended bycom.imsl.stat.KalmanFilter
All Implemented Interfaces:
Cloneable, Serializable

public class KalmanFilter
extends Object
implements Serializable, Cloneable

Performs Kalman filtering and evaluates the likelihood function for the state-space model.

Class KalmanFilter is based on a recursive algorithm given by Kalman (1960), which has come to be known as the Kalman filter. The underlying model is known as the state-space model. The model is specified stage by stage where the stages generally correspond to time points at which the observations become available. KalmanFilter avoids many of the computations and storage requirements that would be necessary if one were to process all the data at the end of each stage in order to estimate the state vector. This is accomplished by using previous computations and retaining in storage only those items essential for processing of future observations.

The notation used here follows that of Sallas and Harville (1981). Let y_k (input in y using method update) be the n_k times 1 vector of observations that become available at time k. The subscript k is used here rather than t, which is more customary in time series, to emphasize that the model is expressed in stages k = 1, 2, ldots and that these stages need not correspond to equally spaced time points. In fact, they need not correspond to time points of any kind. The observation equation for the state-space model is

y_k = Z_kb_k + e_k,,,,,k = 1, 2, 
  ldots

Here, Z_k (input in z using method update) is an n_k times q known matrix and b_k is the q times 1 state vector. The state vector b_k is allowed to change with time in accordance with the state equation

b_{k+1} = T_{k+1}b_k + w_{k+1} ,,,,, 
  k = 1, 2, ldots

starting with b_1 = mu_1 + w_1.

The change in the state vector from time k to k + 1 is explained in part by the transition matrix T_{k+1} (the identity matrix by default, or optionally using method setTransitionMatrix), which is assumed known. It is assumed that the q-dimensional w_ks (k = 1, 2,ldots) are independently distributed multivariate normal with mean vector 0 and variance-covariance matrix sigma^2Q_{k}, that the n_k-dimensional e_ks (k = 1, 2,ldots) are independently distributed multivariate normal with mean vector 0 and variance-covariance matrix sigma^2R_{k}, and that the w_ks and e_ks are independent of each other. Here, mu_1is the mean of b_1 and is assumed known, sigma^2 is an unknown positive scalar. Q_{k+1} (input in Q) and R_k (input in R) are assumed known.

Denote the estimator of the realization of the state vector b_k given the observations y_1, y_2, ldots, y_j by

hat beta _{k|j}

By definition, the mean squared error matrix for

hat beta _{k|j}

is

sigma ^2 C_{left. k right|j}  = E(hat 
  beta _{left. k right|j}  - b_k ) (hat beta _{left. k right|j}  
  - b_k )^T

At the time of the k-th invocation, we have

hat beta _{k|k-1}

and

C_{kleft| {k - 1} right.}, which were computed from the k-1-st invocation, input in b and covb, respectively. During the k-th invocation, KalmanFilter computes the filtered estimate

hat beta _{k|k}

along with C_{kleft| k right.}. These quantities are given by the update equations:

hat beta _{k|k}  = hat beta _{k|k-1}  + 
  C_{k|k-1}Z_k^{T}H_k^{-1}v_k

C_{left. k right|k}  = C_{left. k 
  right|k - 1}  - C_{left. k right|k - 1} Z_k^T H_k^{ - 1} Z_k C_{left. k 
  right|k - 1}

where

v_k  = y_k  - Z_k hat beta _{left. k 
  right|k - 1}

and where

H_k  = R_k  + Z_k C_{left. k right|k - 1} 
  Z_k^T

Here, v_k (stored in getPredictionError) is the one-step-ahead prediction error, and sigma^2{H_k} is the variance-covariance matrix for v_k. H_k is obtained from method getCovV. The "start-up values" needed on the first invocation of KalmanFilter are

hat beta _{1left| 0 right.}  = 
  mu _{_1 }

and C_{1left| 0 right.} = Q{}_1 input via b and covb, respectively. Computations for the k-th invocation are completed by KalmanFilter computing the one-step-ahead estimate

hat beta _{k+1|k}

along with C_{k + 1left| k right.} given by the prediction equations:

hat beta _{k + left. 1 right|k}  = 
  T_{k + 1} hat beta _{left. k right|k}

C_{k+1|k} = T_{k+1}C_{k|k}T_{k+1}^{T} + 
  Q_{k+1}

If both the filtered estimates and one-step-ahead estimates are needed by the user at each time point, KalmanFilter can be used twice for each time point-first without methods SetTransitionMatrix and setQ to produce

hat beta _{left. k right|k}

and C_{kleft| k right.}, and second without method update to produce

hat beta _{k + left. 1 right|k}

and C_{k + 1left| k right.} (Without methods SetTransitionMatrix and setQ, the prediction equations are skipped. Without method update, the update equations are skipped.).

Often, one desires the estimate of the state vector more than one-step-ahead, i.e., an estimate of

hat beta _{k|j}

is needed where k gt j + 1. At time j, KalmanFilter is invoked with method update to compute

hat beta _{j + 1left| j right.}

Subsequent invocations of KalmanFilter without method update can compute

hat beta _{j+2|j}, , hat beta _{j+3|j}, , 
  dots , , hat beta _{k|j}

Computations for

hat beta _{left. k right|j}

and C_{k left| j right.} assume the variance-covariance matrices of the errors in the observation equation and state equation are known up to an unknown positive scalar multiplier, sigma^2. The maximum likelihood estimate of sigma^2 based on the observations y_1, y_2, ldots, y_m, is given by

hat sigma^2  = SS/N

where

N = sumnolimits _{k = 1}^m n_k ,,{rm{and}},,SS = 
  sumnolimits _{k = 1}^m v_k^T H_k^{ - 1} v_k

N and SS are input arguments rank and SumofSquares. Updated values are obtained from methods getRank and getSumofSquares

If sigma^2 is known, the R_ks and Q_ks can be input as the variance-covariance matrices exactly. The earlier discussion is then simplified by letting sigma^2 = 1.

In practice, the matrices T_k, Q_k, and R_k are generally not completely known. They may be known functions of an unknown parameter vector theta. In this case, KalmanFilter can be used in conjunction with an optimization class (see MinUnconMultiVar, JMSL Math package), to obtain a maximum likelihood estimate of theta. The natural logarithm of the likelihood function for y_1, y_2, ldots, y_m differs by no more than an additive constant from

L(theta ,sigma ^2 ;y_1 ,y_2 ,; ldots 
  ,;y_m ) = - frac{1}{2}N,{rm{ln}}, sigma ^{rm{2}}  - 
  frac{1}{2}sumlimits_{k = 1}^m {{rm{ln}}[{rm{det}}(H_k )] - 
  frac{1}{2}sigma ^{ - 2} sumlimits_{k = 1}^m {v_k^T H_k^{ - 1} 
  v_k } }

(Harvey 1981, page 14, equation 2.21).

Here,

sumnolimits_{k = 1}^m 
  {{rm{ln}}[{rm{det}}(H_k )]}

(input in logDeterminant, updated by getLogDeterminant) is the natural logarithm of the determinant of V where sigma^2V is the variance-covariance matrix of the observations.

Minimization of -2L(theta, sigma^2; y_1, y_2, ldots, y_m) over all theta and sigma^2 produces maximum likelihood estimates. Equivalently, minimization of -2L_c(theta; y_1, y_2, ldots, y_m) where

L_c (theta ;y_1 ,y_2 ,; ldots ,;y_m ) = 
  - frac{1}{2}N,{rm{ln}}left( {frac{{SS}}{N}} right) - 
  frac{1}{2}sumlimits_{k = 1}^m {{rm{ln}}[{rm{det}}(H_k )]}

produces maximum likelihood estimates

hat theta , {rm{ and }} , hat sigma ^2  
  = SS/N

Minimization of -2L_c(theta; y_1, y_2, ldots, y_m) instead of -2L(theta, sigma^2; y_1, y_2, ldots, y_m), reduces the dimension of the minimization problem by one. The two optimization problems are equivalent since

hat sigma ^2 (theta ) = SS(theta )/N

minimizes -2L(theta, sigma^2; y_1, y_2, ldots, y_m) for all theta, consequently,

hat sigma ^{^2 } left( theta  right)

can be substituted for sigma^2 in L(theta, sigma^2; y_1, y_2, ldots, y_m) to give a function that differs by no more than an additive constant from L_c(theta; y_1, y_2, ldots, y_m).

The earlier discussion assumed H_k to be nonsingular. If H_k is singular, a modification for singular distributions described by Rao (1973, pages 527-528) is used. The necessary changes in the preceding discussion are as follows:

  1. Replace H_k^{ - 1} by a generalized inverse.
  2. Replace det(H_k) by the product of the nonzero eigenvalues of H_k.
  3. Replace N by sumnolimits_{k = 1}^m {{rm{rank}}left( 
  {H_k } right)}

Maximum likelihood estimation of parameters in the Kalman filter is discussed by Sallas and Harville (1988) and Harvey (1981, pages 111-113).

See Also:
Example, Serialized Form

Constructor Summary
KalmanFilter(double[] b, double[] covb, int rank, double sumOfSquaress, double logDeterminant)
          Constructor for KalmanFilter.
 
Method Summary
 void filter()
          Performs Kalman filtering and evaluates the likelihood function for the state-space model.
 double[] getCovB()
          Returns the mean squared error matrix for b divided by sigma squared.
 double[][] getCovV()
          Returns the variance-covariance matrix of v divided by sigma squared.
 double getLogDeterminant()
          Returns the natural log of the product of the nonzero eigenvalues of P where P * sigma2 is the variance-covariance matrix of the observations.
 double[] getPredictionError()
          Returns the one-step-ahead prediction error.
 int getRank()
          Returns the rank of the variance-covariance matrix for all the observations.
 double[] getStateVector()
          Returns the estimated state vector at time k + 1 given the observations through time k.
 double getSumOfSquares()
          Returns the generalized sum of squares.
 void setQ(double[][] q)
          Sets the Q matrix.
 void setTolerance(double tolerance)
          Sets the tolerance used in determining linear dependence.
 void setTransitionMatrix(double[][] t)
          Sets the transition matrix.
 void update(double[] y, double[][] z, double[][] r)
          Performs computation of the update equations.
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

Constructor Detail

KalmanFilter

public KalmanFilter(double[] b,
                    double[] covb,
                    int rank,
                    double sumOfSquaress,
                    double logDeterminant)
Constructor for KalmanFilter.

Parameters:
b - A double array containing the estimated state vector. b is the estimated state vector at time k given the observations through time k-1.
covb - A double array of size b.length by b.length such that covb * sigma^2 is the mean squared error matrix for b.
rank - An int scalar containing the rank of the variance-covariance matrix for all the observations.
sumOfSquaress - A double scalar containing the generalized sum of squares.
logDeterminant - A double scalar containing the natural log of the product of the nonzero eigenvalues of P where P * sigma^2 is the variance-covariance matrix of the observations.
Throws:
IllegalArgumentException - is thrown if the dimensions of b, and covb are not consistent.
Method Detail

filter

public final void filter()
Performs Kalman filtering and evaluates the likelihood function for the state-space model.


getCovB

public double[] getCovB()
Returns the mean squared error matrix for b divided by sigma squared.

Returns:
a double array of size b.length by b.length such that covb * sigma^2 is the mean squared error matrix for b.

getCovV

public double[][] getCovV()
Returns the variance-covariance matrix of v divided by sigma squared.

Returns:
a double matrix containing a y.length by y.length matrix such that covv * sigma^2 is the variance-covariance matrix of the one-step-ahead prediction error, getPredictionError.

getLogDeterminant

public double getLogDeterminant()
Returns the natural log of the product of the nonzero eigenvalues of P where P * sigma2 is the variance-covariance matrix of the observations.

Returns:
a double scalar containing the natural log of the product of the nonzero eigenvalues of P where P * sigma^2 is the variance-covariance matrix of the observations. In the usual case when P is nonsingular, logDeterminant is the natural log of the determinant of P.

getPredictionError

public double[] getPredictionError()
Returns the one-step-ahead prediction error.

Returns:
a double array of size y.length containing the one-step-ahead prediction error.

getRank

public int getRank()
Returns the rank of the variance-covariance matrix for all the observations.

Returns:
An int scalar containing the rank of the variance-covariance matrix for all the observations.

getStateVector

public double[] getStateVector()
Returns the estimated state vector at time k + 1 given the observations through time k.

Returns:
a double array containing the estimated state vector at time k + 1 given the observations through time k.

getSumOfSquares

public double getSumOfSquares()
Returns the generalized sum of squares.

Returns:
a double scalar containing the generalized sum of squares. The estimate of sigma^2 is given by sumOfSquares / rank.

setQ

public void setQ(double[][] q)
Sets the Q matrix.

Parameters:
q - A double matrix containing the b.length by b.length matrix such that q * sigma^2 is the variance-covariance matrix of the error vector in the state equation. Default: There is no error term in the state equation.

setTolerance

public void setTolerance(double tolerance)
Sets the tolerance used in determining linear dependence.

Parameters:
tolerance - A double scalar containing the tolerance used in determining linear dependence. Default: tolerance = 100.0*2.2204460492503131e-16.

setTransitionMatrix

public void setTransitionMatrix(double[][] t)
Sets the transition matrix.

Parameters:
t - A double matrix containing the b.length by b.length transition matrix in the state equation. Default: t = identity matrix

update

public void update(double[] y,
                   double[][] z,
                   double[][] r)
Performs computation of the update equations.

Parameters:
y - A double array containing the observations.
z - A double matrix containing the y.length by b.length matrix relating the observations to the state vector in the observation equation.
r - A double matrix containing the y.length by y.length matrix such that r * sigma^2 is the variance-covariance matrix of errors in the observation equation. sigma^2 is a positive unknown scalar. Only elements in the upper triangle of r are referenced.

JMSLTM Numerical Library 4.0

Copyright 1970-2006 Visual Numerics, Inc.
Built June 1 2006.