arma_ACV(a, c, v, n)

% gamma = arma_ACV(a, c, v, n)
% Compute ARMA autocovariance sequence
%
% Parameters:
%   a         = Autoregressive parameters, sans leading 1 [p x 1]
%   c         = Moving average parameters, sans leading 1 [q x 1]
%   v         = Innovation variance [1 x 1]
%   n         = Number of autocovariances to generate, starting from zero [1 x 1]
%
% Returns:
%   gamma     = Autocovariances 0 to (n-1) [n x 1]
%
% Reference
%
%  (1) Computation of the Exact Information Matrix of Gaussian Time Series with Stationary Random Components
%      B. Porat and B. Friedlander
%      IEEE Transactions on Acoustics, Speech and Signal Processing, Vol. 34, No. 1, 1986
%
% (c) Copyright Daniel F. Schmidt 2008

function gamma = arma_ACV(a, c, v, n)

% Augment the MA component with the leading 1
a = a';
c = [1,c'];

p = length(a);
q = length(c);

% First p elements of AR autocovariance
A1 = triu(toeplitz([1 a]));
A2 = tril(toeplitz([1 a]*fliplr(eye(p+1))))*fliplr(eye(p+1));

r = inv(A1+A2) * [zeros(1, p), 1]';
r = fliplr(r');
r(1) = r(1)*2;

% Next values of gamma computed recursively 
for i = p+2:n+max([p, q])+1
    r(i) = - a*fliplr(r(i-p:i-1))';
end

% Now compute the ARMA covariances
gamma = zeros(n,1);
for k = 1:n
    gamma(k) = 0;
    for i = 0:q-1
        for j = 0:q-1
            gamma(k) = gamma(k) + c(i+1)*c(j+1) * r(abs((k-1) - i + j) + 1) * v;
        end
    end
end