KalmanFilter Example 2

Example 2: Kalman Filter Maximum Likelihood Estimate

KalmanFilter is used with the MinUnconMultivar class to find a maximum likelihood estimate of the paramether $\theta$ in an MA(1) time series represented by $y_k=\epsilon_k-\theta \epsilon_{k-1}$ . The input is 200 observations from an MA(1) time series with $\theta=0.5$ and $\sigma^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):

$y_k=\begin{pmatrix}1 & 0\end{pmatrix}b_k$

$b_k=\begin{pmatrix}0 & 1\\ 0 & 0\end{pmatrix}b_{k-1}+w_k$

where the two-dimensional s are independently distributed bivariate normal with mean 0 and variance $\sigma^2 Q_k$ and

$Q_1=\begin{pmatrix}1+\theta^2 & -\theta\\ -\theta & \theta^2\end{pmatrix}$

$\begin{matrix}Q_k=\begin{pmatrix} 1 & -\theta\\ -\theta & \theta^2\end{pmatrix} & k=2,3,...,200\end{matrix}$

using System;
using Imsl.Math;
using Imsl.Stat;

public class KalmanFilterEx2 : MinUnconMultiVar.IFunction
{

    public double F(double[] theta)
    {
        int nobs = 200;
        int rank = 0;
        double ss = 0.0;
        double logDeterminant = 0.0;
        double res;
        double[,] covb = new double[2, 2];
        double[,] q = new double[2, 2];
        double[,] r = { { 0.0 } };
        double[] b = { 0.0, 0.0 };
        double[,] z = { { 1.0, 0.0 } };
        double[,] t = {{0.0, 1.0},
                        {0.0, 0.0}};
        double[] y = new double[1];
        double[] ydata = {0.057466,-0.459067,1.236446,-1.864825,0.951507,
                        -1.489367,-0.864401,0.302984,-0.376017,0.610208,
                        -1.701583,2.271734,-0.917110,0.582576,1.018681,
                        0.107443,-0.557858,0.502818,0.858002,-1.417659,
                        0.839981,0.047021,0.404448,0.383749,-0.383227,
                        1.179607,-1.154339,1.927541,0.996344,1.019415,
                        -0.816815,-0.100467,-0.125334,-0.068065,-1.891685,
                        0.657241,-1.070823,1.023510,2.672653,-2.141434,
                        -1.232266,0.925311,-0.201665,-1.325580,-0.086926,
                        -0.647157,1.314749,-0.085708,1.430485,0.304040,
                        0.305559,1.669671,1.004800,-1.678350,0.631133,
                        0.502284,0.247378,-1.345484,0.994982,1.145546,
                        -1.248822,0.616784,-2.127822,2.264872,-1.590343,
                        0.365785,-1.056652,0.969750,1.028274,0.332050,
                        0.430686,0.364553,0.482446,-0.303248,1.581931,
                        -0.140942,-0.265280,-0.939284,0.464963,-0.778145,
                        0.583486,-0.113080,-0.009839,1.580189,-1.116377,
                        1.744513,-0.298106,0.332944,-0.228859,1.101747,
                        0.772369,-1.608111,2.671822,-0.504800,-1.647797,
                        -0.596313,0.845472,0.507869,0.833377,-0.460099,
                        0.416891,-1.139069,0.159028,-0.193971,-0.154656,
                        1.720997,-0.403189,0.400026,0.285921,-1.914338,
                        1.296864,1.426898,-0.426181,0.255961,-1.790193,
                        0.721048,1.150173,-0.980386,0.940958,0.313898,
                        0.505735,-1.058126,0.111918,0.185493,-0.296146,
                        -0.104457,1.151733,0.683025,-0.714269,-0.787972,
                        -1.277062,1.378816,-0.658115,-0.259860,-2.614008,
                        2.251646,-0.006773,-0.738467,-0.260685,1.896505,
                        -0.094919,0.089954,-1.627679,-1.675018,0.896835,
                        0.498690,-0.368775,0.131849,-2.060292,0.272666,
                        2.115804,-1.323451,0.557106,-0.602031,1.424524,
                        -0.107996,-1.580671,0.672012,-1.668931,2.474710,
                        -1.471518,1.780317,-0.419588,2.008474,-1.246855,
                        0.231161,0.706104,-0.474320,-0.705431,0.599358,
                        -2.469494,2.024374,0.849572,-2.410618,2.812321,
                        -1.066520,-0.539768,-0.067784,1.978078,0.592879,
                        -0.184623,-1.403912,-0.995537,1.727320,-0.313251,
                        0.472437,-0.241800,-0.875680,-0.159557,0.508238,
                        -0.116888,-0.981759,-0.472451,0.847273,-1.713030,
                        2.010192,-0.981891,1.190901,0.453348,-0.743333};

        if (Math.Abs(theta[0]) > 1.0)
        {
            res = 1.0e10;
        }
        else
        {
            q[0, 0] = 1.0;
            q[0, 1] = -theta[0];
            q[1, 0] = -theta[0];
            q[1, 1] = theta[0] * theta[0];

            covb[0, 0] = 1.0 + theta[0] * theta[0];
            covb[0, 1] = -theta[0];
            covb[1, 0] = -theta[0];
            covb[1, 1] = theta[0] * theta[0];

            KalmanFilter kalman = new KalmanFilter(b, covb, rank, ss, logDeterminant);

            for (int i = 0; i < nobs; i++)
            {
                y[0] = ydata[i];
                kalman.Update(y, z, r);
                kalman.SetTransitionMatrix(t);
                kalman.SetQ(q);
                kalman.Filter();
            }
            ss = kalman.SumOfSquares;
            logDeterminant = kalman.LogDeterminant;
            rank = kalman.Rank;
            res = rank * Math.Log(ss / rank) + logDeterminant;
        }
        return (res);
    }

    public static void Main(String[] args)
    {
        MinUnconMultiVar solver = new MinUnconMultiVar(1);
        double[] x = solver.ComputeMin(new KalmanFilterEx2());
        Console.Out.WriteLine("Maximum likelihood estimate, THETA = {0}", x[0]);
    }
}

Output

Maximum likelihood estimate, THETA = 0.452940213528753

Link to C# source.