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

Routine KALMN 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. The routine KALMN 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 yk (input in Y) be the nk × 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, … 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

Here, Zk (input in Z) is an nk × q known matrix and bk is the q × 1 state vector. The state vector bk is allowed to change with time in accordance with the state equation

starting with b1 = μ1 + w1.

The change in the state vector from time k to k + 1 is explained in part by the transition matrix Tk+1(input in T), which is assumed known. It is assumed that the q‑dimensional wk’s (k = 1, 2, ).are independently distributed multivariate normal with mean vector 0 and variance‑covariance matrix σ2Qk, that the nk‑dimensional ek’s (k = 1, 2,).are independently distributed multivariate normal with mean vector 0 and variance‑covariance matrix σ2Rk, and that the wk’s and ek’s are independent of each other. Here, μ1 is the mean of b1 and is assumed known, σ2 is an unknown positive scalar. Qk+1 (input in Q) and Rk (input in R) are assumed known.

Denote the estimator of the realization of the state vector bk given the observations y1, y2, …, yj by

By definition, the mean squared error matrix for

is

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

and Ck∣k−1, which were computed from the (k‑1)-st invocation, input in B and COVB, respectively. During the k‑th invocation, routine KALMN computes the filtered estimate

along with Ck∣k. These quantities are given by the update equations:

where

and where

Here, vk (stored in V) is the one‑step‑ahead prediction error, and σ2Hk is the variance‑covariance matrix for vk. Hk is stored in COVV. The “start‑up values” needed on the first invocation of KALMN are

and C1∣0 = Q1 input via B and COVB, respectively. Computations for the k‑th invocation are completed by KALMN computing the one‑step‑ahead estimate

along with Ck+1∣k given by the prediction equations:

If both the filtered estimates and one‑step‑ahead estimates are needed by the user at each time point, KALMN can be invoked twice for each time point—first with IT = 1 and IQ = 1 to produce

and Ck∣k, and second with NY = 0 to produce

and Ck+1∣k (With IT = 1 and IQ = 1, the prediction equations are skipped. With NY = 0, the update equations are skipped.)

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

is needed where k > j + 1. At time j, KALMN is invoked to compute

Subsequent invocations of KALMN with NY = 0 can compute

Computations for

and Ck∣j assume the variance‑covariance matrices of the errors in the observation equation and state equation are known up to an unknown positive scalar multiplier, σ2. The maximum likelihood estimate of σ2 based on the observations y1, y2, …, ym, is given by

where

If σ2 is known, the Rk’s and Qk’s can be input as the variance‑covariance matrices exactly. The earlier discussion is then simplified by letting σ2 = 1.

In practice, the matrices Tk, Qk, and Rk are generally not completely known. They may be known functions of an unknown parameter vector θ. In this case, KALMN can be used in conjunction with an optimization program (see routine UMINF, (IMSL MATH/LIBRARY)) to obtain a maximum likelihood estimate of θ. The natural logarithm of the likelihood function for y1, y2, …, ym differs by no more than an additive constant from

(Harvey 1981, page 14, equation 2.21). Here,

(stored in ALNDET) is the natural logarithm of the determinant of V where σ2V is the variance‑covariance matrix of the observations.

Minimization of ‑2L(θ, σ2; y1, y2, …, ym) over all θ and σ2 produces maximum likelihood estimates. Equivalently, minimization of ‑2Lc(θ; y1, y2, …, ym) where

produces maximum likelihood estimates

The minimization of ‑2Lc(θ; y1, y2, …, ym) instead of ‑2L(θ, σ2; y1, y2, …, ym), reduces the dimension of the minimization problem by one. The two optimization problems are equivalent since

minimizes ‑2L(θ, σ2; y1, y2, …, ym) for all θ, consequently,

can be substituted for σ2 in L(θ, σ2; y1, y2, …, ym) to give a function that differs by no more than an additive constant from Lc(θ; y1, y2, …, ym).

The earlier discussion assumed Hk to be nonsingular. If Hk 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:

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

Routine KALMN is used to compute the filtered estimates and one‑step‑ahead estimates for a scalar problem discussed by Harvey (1981, pages 116–117). The observation equation and state equation are given by

where the ek’s are identically and independently distributed normal with mean 0 and variance σ2, the wk’s are identically and independently distributed normal with mean 0 and variance 4σ2, and b1is distributed normal with mean 4 and variance 16σ2. Two invocations of KALMN are needed for each time point in order to compute the filtered estimate and the one‑step‑ahead estimate. The first invocation uses Default: IQ = 1 and IT = 1 so that the prediction equations are skipped in the computations. The second invocation uses NY = 0 so that the update equations are skipped in the computations.

This example also computes the one‑step‑ahead prediction errors. Harvey (1981, page 117) contains a misprint for the value v4 that he gives as 1.197. The correct value of v4 = 1.003 is computed by KALMN.

Routine KALMN is used with routine UMINF (IMSL MATH/LIBRARY) to find a maximum likelihood estimate of the parameter θ in a MA(1) time series represented by yk = ɛk ‑ θɛk−1. Routine RNARM (see Chapter 18, “Random Number Generation”) is used to generate 200 random observations from an MA(1) time series with θ = 0.5 and σ2 = 1.

The MA(1) time series is cast as a state‑space model of the following form (see Harvey 1981, pages 103–104, 112):

where the two‑dimensional wk’s are independently distributed bivariate normal with mean 0 and variance σ2 Qk and

The warning error that is printed as part of the output is not serious and indicates that UMINF is generally used for multi‑parameter minimization.