arma_KalmanLikelihood(A, C, v, y, eps);

% [L, v_est, v_var, yp] = max_KalmanLikelihood(A, C, v, y, eps)
% Efficiently evaluation Autoregressive-Moving Average Model Likelihood 
% via a Kalman Filter
%
% The model is specified in the form
%
%  y_n + A_1 y_(n-1) + A_2 y_(n-2) ... = v_n + C_1 v_(n-1) + C_2 v_(n-2) + ...
%
% where v_n ~ N(0, v) are the innovations
%
%
% Parameters:
%   A       = Autoregressive component, sans the leading 1 [p x 1]
%   C       = Moving average component, sans the leading 1 [q x 1]
%   v       = Innovation variance [1 x 1]
%   y       = Vector of observations (the data) [n x 1]
%   eps     = (Optional): Once (v_var(i) - v) < eps, the algorithm uses
%             simple model inversion to compute the remaining innovations
%             If this parameter is not specified, the Kalman Filter is used
%             to compute all innovations; eps should at least be < 0.05 [1 x 1]
% 
% Returns:
%   L       = Negative log-likelihood of y ~ ARMA(A, C, v) [1 x 1]
%   v_est   = Estimated values of the innovations [n x 1]
%   v_var   = Estimated values of the innovation variances [n x 1]
%   yp      = Predictions; y = yp + v_est [n x 1]
%
%
% References:
%
%   (1) Algorithm AS154: An Algorithm for Exact Maximum Likelihood Estimation of Autoregressive-Moving Average Models by Means of Kalman Filtering
%       G.Gardner, A.C.Harvey and G.D.A.Phillips
%       Applied Statistics, Vol. 29, No. 3, 1980, pp. 311-322
% 
% (c) Copyright Daniel F. Schmidt 2008
%
function [L, v_est, v_var, yp] = arma_KalmanLikelihood(A, C, v, y, eps);

%
A=A';
C=C';

%
if (size(y,2)==1)
    y=y';
end

%
N = length(y);

% If eps was not passed, set it to negative (thus ensuring Kalman Filtering
% is used for all N measurements)
if (exist('eps') ~= 1)
    eps = -0.01;
end

% If eps >= 0.05
if (eps >= 0.05)
    warning('Setting eps to greater than 0.05 may degrade the approximation badly');
end

% First of all, produce the SS representation of this model
[F, R, H, Q] = arma_ConvertToSS([A], [C], v);

% Seed the initial state with zeros and the initial state covariance with 
% the model autocovariances
m = size(F,1);

% Initialise the state and covariance X and P0: For the algorithm see
X = zeros(m, 1);

P = zeros(m, m);
S = sparse(eye(m^2) - kron(F,F));       
P(:) = inv(S)*Q(:);

% Sanity check - if NaN, just return very large likelihood
if (isnan(det(P)))
    L = inf;
    v_est = zeros(1,N);
    v_var = ones(1,N);
    yp = zeros(1,N);
    
    return;
end

% Evaluate the likelihood
L = 0;

yp = zeros(N,1);
v_est = zeros(N,1);
v_var = zeros(N,1);

for i = 1:N
    % Prediction step
    X_pred = F*X;
    y_pred = H*X_pred;
    P_pred = F*P*F' + Q;
    
    % Costing step
    v_est(i) = y(i) - y_pred;
    v_var(i) = H*P_pred*H';
    yp(i) = y_pred;
    L = L + (1/2)*log(v_var(i)) + (1/2/v_var(i))*v_est(i)^2;
        
    % Use inversion?
    if ( (v_var(i) - v) < eps && i > (m+1))
        break;
    end
    
    % Update step
    K = P_pred*H'*inv(v_var(i));
    X = X_pred + K*v_est(i);
    P = P_pred - K*H*P_pred;
end

% Estimate any remaining innovations via model inversion
if (i < N)
    % Estimate the innovations
    Ic = filtic([1, A], [1, C], fliplr(v_est(i-m-1:i)), fliplr(y(i-m-1:i)));
    v_est(i+1:N) = filter([1, A], [1, C], y((i+1):N), Ic);
    
    % Update the negative log-likelihood
    Ni = (N-i);
    v_var(i+1:N) = ones(1,Ni)*v;
    L = L + (Ni/2)*log(v) + (1/2/v)*v_est(i+1:N)*v_est(i+1:N)';
end

% Add in the normalisation constant 
L = L + (N/2)*log(2*pi);

return;