function c = acv(vParameters,nAR,nMA,nT)
% PURPOSE: computes the autocorrelation function of the arfima model
% -----------------------------------------------------
% USAGE:     c = acv(vParameters,nAR,nMA,nT)
%  where:    series = a time series matrix
%            vParameters = a matrix with all coefficients, in order:d, AR, MA, beta, gamma
%            nAR, nMA = the order number of AR, MA
% -----------------------------------------------------
% RETURNS:   c = the autocovariance coefficient of the AutoRegressive model.
% -----------------------------------------------------             
% NOTE:      the function will be used by the acfilter function.
% -----------------------------------------------------

% The matlab port written by:
% Ying Ma
% yingma@email.com

d = vParameters(1);
phi= vParameters(2:1+nAR); 
phi=[1; phi];
theta= vParameters(2+nAR:1+nAR+nMA);
theta=[1; theta];
% Get the polynomial psi from MA 
psi= conv(theta,flipud(theta));

% Autocovariance coefficients of the ARFIMA model
if (nAR ~= 0)
  % Calculate the roots of the polynomial psi
  rho=roots(phi);   

  % Check to ensure none of roots is close to or outstanding the unit circle
  if (max(abs(rho)) > 0.99)
    %error('Error:fracsigma failed because at least one root is too close to the unit circle');
    %return 1.0;  % make sure none of such root exists
    c = zeros(1,nMA+nT+1);
  end

  % Compute  the xi_i coefficients:
  xi =zeros(nAR,1);
  for j=1:nAR,
    xi(j)= 1/(prod(1 - rho(1:j-1)/rho(j))*prod(1-rho(j+1:nAR)/rho(j))*prod(1-rho*rho(j))); 
  end
end


T= nMA+nT+1;
coef= 1;
% Case 1: non fractional models d is 0: use standard evaluation 
if (abs(d) <= 1e-9) 
   c= zeros(1,T);
   if (nAR ~= 0)
      for j=1:nAR,
         c= c + xi(j)*rho(j).^(0:(T-1));
      end
      c= real(c);
   else
      c(1)= 1;
   end;
   c= c(:);
else
   %Case 2: fractional models 
   c= zeros(T,1);
   j=0;
   c(1)=1;
   for m=2:T+1,
      c(m) = c(m-1) * (d+j) / (1-d+j);
      j= j+1;
   end

   cinv= flipud(c(2:T+1));
   coef= exp(gammaln(1-2*d)-2*gammaln(1-d));
   
   if (nAR ~= 0)
      % Compute the sum
      csum= zeros(T,1);
      for j=1:nAR
         U = zeros(T,1);
         U0= hypergm(d,(1-d),rho(j),1e+5);
         U(1)= U0;
         U= filter(1, [1 -rho(j)],[0;c(2:T)],U0);
         V0= cinv(1)*hypergm(d+T+1,2-d+T,rho(j),1e+5);
         V= filter(1, [1 -rho(j)],cinv,V0);
         csum= csum + xi(j)*(U(1:T)+rho(j)*flipud(V(1:T)));
      end  
      c= csum;
  end
end
c= [flipud(c(2:nMA+1)); c];

% Convolve psi,coef and c
c= conv(psi,c);
c= coef*c(2*nMA+1:2*nMA+nT);
c= real(c);
