1
$\begingroup$

I am trying to implement the exact maximum likelihood estimation of ARMA(p,q) models using the Kalman filter. I wrote the following code in MATLAB and seems right to me, however when i compare the estimates that it returns with those obtained with other software such as gretl or matlab they are different. The code is the following

   function [Estimate] = ExactMLmissing_ARMA(y,ARorder,MAorder)

   %% Compute the exact maximum likelihood estimate of an ARMA model
   % the model is placed in the state space form
   %
   %   y(t) = Z*a(t) + e(t)
   %
   %   a(t+1) = T*a(t) + R*v(t)
   %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % number of observations
   n = length(y);

   % define the order of the model to estimate;
   p = ARorder;
   q = MAorder;
   m = max([p,q+1]);


   % initialize the hyperparameters values
   % the vector is partitioned as follows
   %[p AR parameters, q MA parameters, error variance]
   Psi0 = repmat(0.3,1,(p+q+1));


   % initialize v that will contain the one step ahead prediction errors
   % computed by the filter
   v = zeros(n,1);

   function dz = PEdecomposition(Psi)

   % initialize all the system matrices

   % initialize matrix T
   I = eye(m-1);
   Il = zeros(1,m-1);
   arpar = zeros(m,1);
   arpar(1:p) = Psi(1:p);
   T = [arpar,[I;Il]];


   % vector Z
   Z = [1,zeros(1,m-1)];

   % vector R
   R = zeros(m,1);
   R(1) = 1;
   if q > 0
   R(2:end) = Psi(p+1:end-1);
   end


   % matrix Q
   Q = Psi(end);

   % since the process is stationary a(1) can be placed equal to 0
   % (Reference: pp 126-137, Koopman Durbin)
   a_t = zeros(m,1);

   %exact initialization of state covariance matrix (Reference: pp 138, 
   Koopman Durbin)
   LL = (R*Q*R');
   LL = LL(:);
   vecP_t = (eye(size(kron(T,T))) - kron(T,T))\LL;
   P_t = reshape(vecP_t,m,m);


    % initialize llh
    llh1 = 0;
    llh2 = 0;


    for i = 1:n

        % Check whether the unit is missing or not
        if ~isnan(y(i))

        % one step ahead prediction error
        v_t = y(i) - Z*a_t;
        %store the one step ahead prediction error
        v(i) = v_t;

        % variance of the one step ahead prediction error
        F_t = Z*P_t*Z';

        %% updating equations
        % estimate state vector at time t given the information at time t
        a_tt = a_t + P_t*Z'*inv(F_t)*v_t;

        % variance of state vector at time t given information at time t
        P_tt = P_t - P_t*Z'*inv(F_t)*Z*P_t;

        % compute the Kalman Gain
        K = T*P_t*Z'*inv(F_t);

        %% prediction equations
        a_t = T*a_tt +K*v_t;
        P_t = T*P_tt*(T-K*Z)'+R*Q*R';



        else
            % compute Kalman filter step for missing values
            % Z equal to 0 and computing the standard filter
            % recursions

        % one step ahead prediction error
        v_t = y(i);
        v(i) = v_t;

        % variance of the one step ahead prediction error
        F_t = 0;

        %% updating equations
        % estimate state vector at time t given the information at time t
        a_tt = a_t ;

        % variance of state vector at time t given information at time t
        P_tt = P_t ;


        %% prediction equations
        a_t = T*a_tt;
        P_t = T*P_tt*T'+R*Q*R';

        llhincrement = 0;

        end

         %% log-likelihood quantities
        llh1 = llh1 +log(F_t);
        llh2 = llh2 + v_t'*inv(F_t)*v_t;

        negllh = n/2*log(Psi(end))+ 0.5*llh1+1/(2*Psi(end))*llh2;

    end

    dz = negllh;


   end
   %  tt = toc;
   %  fprintf('Total estimated time to complete LMS: %5.2f seconds \n', tt);

   options = optimset('Display','none');
   [thetaHat,fval] = fminsearch(@PEdecomposition,Psi0,options);
   Estimate.AR = thetaHat(1:p);
   Estimate.MA = thetaHat(p+1:p+q);
   Estimate.s2 = thetaHat(end);
   Estimate.v = v;
   Estimate.llh = -fval;
   end

Is there anyone that can help me. Thank you

$\endgroup$

1 Answer 1

1
$\begingroup$

Perhaps a different constant term? One easy way to check is to compute yours and the competing likelihood for different values of the parameters and see if they are proportional (or differ only by a fixed additive ammount, if you use log-likelihood).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.