function [C,Q] = acf(this,varargin)
% acf  Autocovariance and autocorrelation functions for VAR variables.
%
% Syntax
% =======
%
%     [C,R] = acf(V,...)
%
% Input arguments
% ================
%
% * `V` [ VAR ] - VAR object for which the ACF will be computed.
%
% Output arguments
% =================
%
% * `C` [ numeric ] - Auto/cross-covariance matrices.
%
% * `R` [ numeric ] - Auto/cross-correlation matrices.
%
% Options
% ========
%
% * `'applyTo='` [ logical | *`Inf`* ] - Logical index of variables to
% which the `'filter='` will be applied; the default Inf means all
% variables.
%
% * `'filter='` [ char  | *empty* ] - Linear filter that is applied to
% variables specified by 'applyto'.
%
% * `'nfreq='` [ numeric | *`256`* ] - Number of equally spaced frequencies
% over which the 'filter' is numerically integrated.
%
% * `'order='` [ numeric | *`0`* ] - Order up to which ACF will be
% computed.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

[ny,p,nalt] = sizeof(this);

opt = passvalopt('VAR.acf',varargin{:});

%**************************************************************************

[isfilter,filter,freq,applyto] = ...
    freqdom.applyfilteropt(this,opt,[],ny,this.ynames);

C = nan(ny,ny,opt.order+1,nalt);
% find explosive parameterisations
explosive = isexplosive(this);

if opt.progress
    pbar = progressbar('IRIS VAR.acf progress');
end

for ialt = find(~explosive)
    [T,R,K,Z,H,D,U,Omega] = sspace(this,ialt); %#ok<ASGLU>
    if isfilter
        S = freqdom.xsfvar(this.A(:,:,ialt),Omega,freq,filter,applyto);
        C(:,:,:,ialt) = freqdom.xsf2acf(S,freq,opt.order);
    else
        % Compute contemporaneous ACF for its first-order state space form.
        % This gives us autocovariances up to order p-1.
        c = covfun.acovf(T,R,[],[],[],[],U,Omega,this.eigval(1,:,ialt),0);
        if p > 1
            c0 = c;
            c = reshape(c0(1:ny,:),ny,ny,p);
        end
        if p == 0
            c(:,:,end+1:opt.order+1) = 0;
        elseif opt.order > p - 1
            % Compute higher-order acfs using Yule-Walker equations.
            c = xxacovfyw(this.A(:,:,ialt),c,opt.order);
        else
            c = c(:,:,1:1+opt.order);
        end
        C(:,:,:,ialt) = c;
    end
    % Update the progress bar.
    if opt.progress
        update(pbar,ialt/sum(~explosive));
    end
end

if any(explosive)
    % Explosive parameterisations.
    VAR.warning(10,sprintf(' #%g',find(explosive)));
end

% Fix entries with negative variances.
C = timedom.fixcov(C);

% Autocorrelation function.
if nargout > 1
    % Convert covariances to correlations.
    Q = covfun.cov2corr(C,'acf');
end

end

% Subfunctions.

%**************************************************************************
function C = xxacovfyw(A,C,order)
[ny,aux] = size(A);
p = aux/ny;

% residuals included or not in ACF
ne = size(C,1) - ny;

A = reshape(A(:,1:ny*p),[ny,ny,p]);
C = C(:,:,1+(0:p-1));
for i = p : order
    aux = zeros([ny,ny+ne]);
    for j = 1 : size(A,3)
        aux = aux + A(:,:,j)*C(1:ny,:,end-j+1);
    end
    C(1:ny,:,1+i) = aux;
end

end
% xxacovfyw().
