Maximum Likelihood Estimation for PDF via Kalman Filter

I'm currently trying to do a MLE for the parameter estimation of a linear gaussian model. (Schwartz Smith Model for spot prices of commodities.)

My model consists of an Ornstein Uhlenbeck process and a brownian motion. The sum of both of the two resulting states represents the observable spot price of the commodity. Schwartz and Smith show in their paper, that the state variables are multivariat normal.

I managed to implement the Kalman Filter for my model and to compute Likelihoods, but here is my problem:

My observable variable is given directly by the sum of the two state variables and therefore I don't want to include any measurement noise in the Filter. But if the measurement noise approaches 0, the log-likelihood approaches minus infinity. If i include any measurement noise, the log likelihoods don't even get close to reasonable values. I got the best results for a very small measurement noise, but the log likelihoods reach values up tp 10^(-10), which might cause some numerical issues in the MLE? Any advice, what to do?

The Code to produce the simulation data is:

bm_obj = bm(t_drift,t_sd_bm);
bm_data = bm_obj.simulate(n);

ou_obj = sdemrd(t_reversion,t_mean,0,t_sd_ou);
ou_data = ou_obj.simulate(n);

data = bm_data + ou_data;

And the Kalman Filter:

    %System Matrices 
    A = [exp(-rev) 0; 0 1];

    Ex = [(1-exp(-2*rev))*(sd_ou^2/(2*rev)) 0; 0 sd_bm^2];

    Ez = sd_measurement^2;

    B = [0;1];

    C = [1 1];

    %Start

    results = zeros(n,8);

    X = [0 ; data(1)];

    P = Ex;
    X_est = X;

    for i = 2:n+1

        X_est = A * X_est + B * drift;

        data_est = C * X_est;

        P = A * P * A' + Ex;

        K = P * C' * inv(C * P * C' + Ez);

        X_est = X_est + K * (data(i) - data_est);

        P =  (eye(2)-K*C)*P;

        results(i-1,1) = X_est(1);
        results(i-1,2) = X_est(2);
        results(i-1,3) = data(i) - data_est;
        results(i-1,4) = P(1,1);
        results(i-1,5) = P(2,2);
        results(i-1,6) = C * P * C' + sd_measurement^2;
        results(i-1,7) = K(1);
        results(i-1,8) = K(2);

    end

    Sum = 0;

    for i = 1:n
        Sum = Sum + log(results(i,6))+ ((results(i,3)^2) / results(i,6));
        LogLikelihood = - n / 2 * log(2 * pi) - 1/2 * Sum;
    end

L = LogLikelihood;

I'm grateful for any help!

Answer 1

[Shangw Song]

Hi,

Have you settle this problem?