function result = acf (series,ilag)
% PURPOSE: computes the autocorrelation function of the columns of series up to lag ilag
% -----------------------------------------------------
% USAGE:     m = acf(series,ilag)
%  where:    series = a matrix of time series	
%            ilag = a integer of upbound lag
% -----------------------------------------------------
% RETURNS:  m = a correlogram of series.
% -----------------------------------------------------             
% NOTE: 
% -----------------------------------------------------
% ACKNOWLEDGEMENT: The original ARFIMA code was written by Jurgen A Doornik
% -----------------------------------------------------

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

    if  length(series) == 0 | nargin < 2
        error('The series and lag must not be empty.');
    end
    
    
    [nT n] = size(series);

    m = series - ones(nT,1)*mean(series,1);
    
    result = zeros(ilag+1,n); %initial results
    
    vec0 = diag(m'*m);
    result(1,:) = vec0';
    result(1,:) = result(1,:)./vec0';
    
    for i=2:ilag+1
        vec = diag(m'*lagmatrix(m,i-1));
        result(i,:) =   vec';
        result(i,:) = result(i,:)./vec0';
    end