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
