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 w_ks 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}


import com.imsl.stat.*;
import com.imsl.math.*;

public class KalmanFilterEx2 implements MinUnconMultiVar.Function {

    public double f(double[] theta) {
        int nobs = 200, 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.getSumOfSquares();
            logDeterminant = kalman.getLogDeterminant();
            rank = kalman.getRank();
            res = rank * Math.log(ss / rank) + logDeterminant;
        }
        return (res);
    }

    public static void main(String args[]) throws com.imsl.IMSLException {
        MinUnconMultiVar solver = new MinUnconMultiVar(1);
        double[] x = solver.computeMin(new KalmanFilterEx2());
        System.out.println("Maximum likelihood estimate, THETA = " + x[0]);
    }
}

Output

Maximum likelihood estimate, THETA = 0.4529402135287534
Link to Java source.