function acvf = armaacvf(phi,theta,n)
% ARMAACVF(PHI,THETA,N) computes acvf out to lag N of ARMA model with
%   given coefficients in (PHI,THETA) assuming sigma_a^2 = 1.
% phi, theta must be column vectors (px1 and qx1, respectively).
% returns acvf = [ gamma_0 gamma_1 ... gamma_n ]'
[p m] = size(phi);
[q m] = size(theta);
phi1 = [1 ; -phi ];
if q>p
   phi1 = [phi1 ; zeros(q-p,1) ];
end
theta1 = [1 ; -theta];
if p>q
   theta1 = [ theta1 ; zeros(p-q,1) ];
end
m = 1+ max(p,q);
% find gamma_0 , ... , gamma_{m-1} by solving linear system
% setting up the matrix
R = zeros(m,m);
T = toeplitz(1:m);
for i=1:m
   for j=1:m
      R(i,T(i,j)) = R(i,T(i,j))+phi1(j);
   end
end
% setting up r.h.s.
psi1 = [1];
for i=2:m
   psi1 = [ psi1 ; theta1(i) - psi1(1:(i-1))'*phi1(i:-1:2) ];
end
rhs = zeros(m,1);
for i=1:m
   rhs(i) = theta1(i:m)'*psi1(1:(m-i+1));
end
acvf = R\rhs;
if m > n
   acvf = acvf(1:(n+1));
else
   for i=m:n
      temp = phi'*acvf(i:-1:(i-p+1));
      acvf = [ acvf ; temp ];
   end
end

