Example: Kilman Filter

KalmanFilter 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

k = 1, 2, 3, 4

where the s are identically and independently distributed normal with mean 0 and variance , the s are identically and independently distributed normal with mean 0 and variance , and is distributed normal with mean 4 and variance . Two KalmanFilter objects are needed for each time point in order to compute the filtered estimate and the one-step-ahead estimate. The first object does not use the methods SetTransitionMatrix and SetQ so that the prediction equations are skipped in the computations. The update equations are skipped in the computations in the second object.

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

using System;
using Imsl.Stat;

public class KalmanFilterEx1
{
	private static readonly String format =
		"{0}/{1}\t{2:0.000}\t{3:0.000}\t{4}\t{5:0.000}" + 
        "\t{6:0.000}\t{7:0.000}\t{8:0.000}";
	
	public static void  Main(String[] args)
	{
		int nobs = 4;
		int rank = 0;
		double logDeterminant = 0.0;
		double ss = 0.0;
		double[] b = new double[]{4};
		double[,] covb = new double[1,1]{{16}};
		double[,] q = {{4}};
		double[,] r = {{1}};
		double[,] t = {{1}};
		double[,] z = {{1}};
		double[] ydata = new double[]{4.4, 4.0, 3.5, 4.6};
		
		System.Object[] argFormat = 
			new System.Object[]{"k", "j", "b", "cov(b)", "rank", "ss",
								   "ln(det)", "v", "cov(v)"};
		Console.Out.WriteLine(String.Format(format, argFormat));
		
		KalmanFilter kalman =
			new KalmanFilter(b, covb, rank, ss, logDeterminant);

		for (int i = 0; i < nobs; i++)
		{
			double[] y = new double[]{ydata[i]};
			kalman.Update(y, z, r);
			kalman.Filter();
			double[] v = kalman.GetPredictionError();
			double[,] covv = kalman.GetCovV();
			argFormat[0] = i;
			argFormat[1] = i;
			argFormat[2] = kalman.GetStateVector()[0];
			argFormat[3] = kalman.GetCovB()[0,0];
			argFormat[4] = kalman.Rank;
			argFormat[5] = kalman.SumOfSquares;
			argFormat[6] = kalman.LogDeterminant;
			argFormat[7] = v[0];
			argFormat[8] = covv[0,0];
			Console.Out.WriteLine(String.Format(format, argFormat));
			kalman.ResetUpdate();
			
			kalman.SetTransitionMatrix(t);
			kalman.SetQ(q);
			kalman.Filter();
			argFormat[0] = i + 1;
			argFormat[1] = i;
			argFormat[2] = kalman.GetStateVector()[0];
			argFormat[3] = kalman.GetCovB()[0,0];
			argFormat[4] = kalman.Rank;
			argFormat[5] = kalman.SumOfSquares;
			argFormat[6] = kalman.LogDeterminant;
			argFormat[7] = v[0];
			argFormat[8] = covv[0,0];
			Console.Out.WriteLine(String.Format(format, argFormat));
			kalman.ResetTransitionMatrix();
			kalman.ResetQ();
		}
	}
}

Output

k/j	b	cov(b)	rank	ss	ln(det)	v	cov(v)
0/0	4.376	0.941	1	0.009	2.833	0.400	17.000
1/0	4.376	4.941	1	0.009	2.833	0.400	17.000
1/1	4.063	0.832	2	0.033	4.615	-0.376	5.941
2/1	4.063	4.832	2	0.033	4.615	-0.376	5.941
2/2	3.597	0.829	3	0.088	6.378	-0.563	5.832
3/2	3.597	4.829	3	0.088	6.378	-0.563	5.832
3/3	4.428	0.828	4	0.260	8.141	1.003	5.829
4/3	4.428	4.828	4	0.260	8.141	1.003	5.829