# Kalman Filter Likelihood

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