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

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

public class KalmanFilterEx2 implements MinUnconMultiVar.Function {

    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.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.