function [CC,RR,list] = acf(this,varargin)
% acf  Autocovariance and autocorrelation functions for model variables.
%
% Syntax
% =======
%
%     [C,R,LIST] = acf(M,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Solved model object for which the ACF will be computed.
%
% Output arguments
% =================
%
% * `C` [ namedmat | numeric ] - Auto/cross-covariance matrices.
%
% * `R` [ namedmat | numeric ] - Auto/cross-correlation matrices.
%
% * `LIST` [ cellstr ] - List of variables in rows and columns of `C` and
% `R`.
%
% Options
% ========
%
% * `'applyTo='` [ cellstr | char | *`Inf`* ] - List of variables to which the
% `'filter='` will be applied; `Inf` means all variables.
%
% * `'contributions='` [ `true` | *`false`* ] - If `true` the contributions
% of individual shocks to ACFs will be computed and stored in the 5th
% dimension of the `C` and `R` matrices.
%
% * `'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 in the option `'filter='` is numerically integrated.
%
% * `'order='` [ numeric | *`0`* ] - Order up to which ACF will be computed.
%
% * `'output='` [ *`'namedmat'`* | `'numeric'` ] - Output matrices `C` and `R`
% will be either namedmat objects or plain numeric arrays; if the option
% `'select='` is used, `'output='` is always a namedmat object.
%
% * `'select='` [ cellstr | *`Inf`* ] - Return the ACF matrices for selected
% variables only; `Inf` means all variables.
%
% Description
% ============
%
% ACF with linear filters 
% ------------------------
%
% You can use the option `'filter='` to get the ACF for variables as though
% they were filtered through a linear filter. You can specify the filter in
% both the time domain (such as first-difference filter, or
% Hodrick-Prescott) and the frequncy domain (such as a band of certain
% frequncies or periodicities). The filter is a text string in which you
% can use the following references:
%
% * `'L='`, the lag operator, which will be replaced with `exp(-1i*freq)`;
% * `'per='`, the periodicity,
% * `'freq='`, the frequency.
% 
% Example 1
% ==========
%
% A first-difference filter (i.e. computes the ACF for the first
% differences of the respective variables):
%
%     [C,R] = acf(m,'filter=','1-L')
%
% Example 2
% ==========
%
% The cyclical component of the Hodrick-Prescott filter with the
% smoothing parameter, $lambda$, 1,600. The formula for the filter follows
% from the classical Wiener-Kolmogorov signal extraction theory,
%
% $$w(L) = \frac{\lambda}{\lambda + \frac{1}{ | (1-L)(1-L) | ^2}}$$
%
%     [C,R] = acf(m,'filter','1600/(1600 + 1/abs((1-L)^2)^2)')
%
% Example 3
% ==========
%
% A band-pass filter with user-specified lower and upper bands. The
% band-pass filters can be defined either in frequencies or periodicities;
% the latter is usually more convenient. The following is a filter which
% retains periodicities between 4 and 40 periods (if the model has been
% estimated on quarterly data, this would correspond to periodicities
% between 1 and 10 years),
%
%     [C,R] = acf(m,'filter','per >= 4 & per <= 40')

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

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

if ischar(opt.select)
    opt.select = regexp(opt.select, ...
        '[a-zA-Z][\w\(\)\{\}\+\-]*','match');
elseif isempty(opt.select)
    opt.select = {''};
end

isselect = iscellstr(opt.select);
isnamedmat = strcmpi(opt.output,'namedmat') || isselect;

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

ny = length(this.solutionid{1});
nx = length(this.solutionid{2});
ne = length(this.solutionid{3});
nalt = size(this.Assign,3);

if opt.contributions
    ncont = ne;
else
    ncont = 1;
end
CC = nan(ny+nx,ny+nx,opt.order+1,ncont,nalt);

[isfilter,filter,freq,applyto] = freqdom.applyfilteropt( ...
    this,opt,[],ny+nx,this.name(this.nametype <= 2));

% Solution not available for some parameterisations.
[flag,nansolution] = isnan(this,'solution');
if flag
    utils.warning('model', ...
        '#Solution_not_available',strfun.alt2str(nansolution));
end

% Autocovariance function.
for ialt = 1 : nalt
    if nansolution(ialt)
        continue
    end
    [T,R,K,Z,H,D,U,Omega] = mysspace(this,ialt,false);
    for icont = 1 : ncont
        if opt.contributions
            index = false([1,ne]);
            index(icont) = true;
            if Omega(index,index) == 0
                CC(:,:,:,icont,ialt) = 0;
                continue
            end
        else
            index = true(1,ne);
        end
        if isfilter
            S = freqdom.xsf( ...
                T,R(:,index),K,Z,H(:,index),D,U,Omega(index,index), ...
                freq,filter,applyto);
            CC(:,:,:,icont,ialt) = freqdom.xsf2acf(S,freq,opt.order);
        else
            CC(:,:,:,icont,ialt) = covfun.acovf( ...
                T,R(:,index),K,Z,H(:,index),D,U,Omega(index,index), ...
                this.eigval(1,:,ialt),opt.order);
        end
    end
end

% Squeeze the covariance matrices if ~contributions.
if ~opt.contributions
    CC = reshape(CC,[ny+nx,ny+nx,opt.order+1,nalt]);
end

% Fix negative variances.
tmpsize = size(CC);
if length(tmpsize) < 4
    tmpsize(end+1:4) = 1;
end
CC0 = reshape(CC(:,:,1,:),[tmpsize(1:2),tmpsize(4)]);
CC0 = timedom.fixcov(CC0);
CC(:,:,1,:) = reshape(CC0,[tmpsize(1:2),1,tmpsize(4)]);

list = [this.solutionvector{1:2}];

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

% Convert double arrays to namedmat objects.
if isnamedmat
    CC = namedmat(CC,list,list);
    if exist('RR','var')
        RR = namedmat(RR,list,list);
    end
end

% Select variables.
if isselect
    [CC,index] = select(CC,opt.select);
    if exist('RR','var')
        RR = RR(index{1},index{2},:,:,:);
    end
end

end