function [this,d,c,f,u,e] = filter(this,data,range,varargin)
% filter  Re-estimate the factors by Kalman filtering the data taking FAVAR coefficients as given.
%
% Syntax
% =======
%
%     [a,d,cc,f,u,e] = filter(a,d,range,...)
%
% Input arguments
% ================
%
% * `f` [ FAVAR ] - Estimated FAVAR object.
%
% * `d` [ struct | tseries ] - Input database or tseries object with the
% FAVAR observables.
%
% * `range` [ numeric ] - Filter date range.
%
% Output arguments
% =================
%
% * `a` [ FAVAR ] - FAVAR object.
%
% * `d` [ struct ] - Output database or tseries object with the FAVAR
% observables.
%
% * `cc` [ struct | tseries ] - Re-estimated common components in the
% observables.
%
% * `f` [ tseries ] - Re-estimated common factors.
%
% * `u` [ tseries ] - Re-estimated idiosyncratic residuals.
%
% * `e` [ tseries ] - Re-estimated structural residuals.
%
% Options
% ========
%
% * `'cross='` [ *true* | false | numeric ] - Run the filter with the
% off-diagonal elements in the covariance matrix of idiosyncratic
% residuals; if false all cross-covariances are reset to zero; if a number
% between zero and one, all cross-covariances are multiplied by that
% number.
%
% * `'invFunc='` [ *'auto'* | function_handle ] - Inversion method for the
% FMSE matrices.
%
% * `'meanOnly='` [ `true` | *`false`* ] - Return only mean data, i.e. point
% estimates.
%
% * `'persist='` [ `true` | *`false`* ] - If `filter` or `forecast` is used with
% `'persist='` set to true for the first time, the forecast MSE matrices and
% their inverses will be stored; subsequent calls of the `filter` or
% `forecast` functions will re-use these matrices until `filter` or
% `forecast` is called.
%
% * `'output='` [ *'auto'* | 'dbase' | 'tseries' ] - Format of output data.
%
% * `'tolerance='` [ numeric | *0* ] - Numerical tolerance under which two
% FMSE matrices computed in two consecutive periods will be treated as
% equal and their inversions will be re-used, not re-computed.
%
% Description
% ============
%
% It is the user's responsibility to make sure that `filter` and `forecast`
% called with `'persist='` set to true are valid, i.e. that the
% previously computed FMSE matrices can be really re-used in the current
% run.
%
% Example
% ========
% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

% Parse input arguments.
P = inputParser();
P.addRequired('a',@(x) isa(x,'FAVAR'));
P.addRequired('d',@(x) isstruct(x) || isa(x,'tseries'));
P.addRequired('range',@isnumeric);
P.parse(this,data,range);

% Parse options.
options = passvalopt('FAVAR.filter',varargin{:});

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

nx = size(this.C,2);
p = size(this.A,2)/nx;
range = range(1) : range(end);

% Retrieve and standardise input data.
[outputformat,range,y] = VAR.datarequest(this,data,range,options);
[this,y] = standardise(this,y);
ndata = size(y,3);
nper = size(y,2);

% Initialise Kalman filter.
x0 = zeros([p*nx,1]);
R = this.U(1:nx,:)'*this.B;
P0 = covfun.acovf(this.T,R,[],[],[],[],this.U,1,[],0);

Sigma = this.Sigma;
% Reduce or zero off-diagonal elements in the cov matrix of idiosyncratic
% residuals if requested.
options.cross = double(options.cross);
if options.cross < 1
    index = logical(eye(size(Sigma)));
    Sigma(~index) = options.cross*Sigma(~index);
end

% Inversion method for the FMSE matrix. It is safe to use `inv` if
% cross-correlations are pulled down because then the idiosyncratic cov
% matrix is non-singular.
if isequal(options.invfunc,'auto')
    if this.cross == 1 && options.cross == 1
        invfunc = @pinv;
    else
        invfunc = @inv;
    end
else
    invfunc = options.invfunc;
end

% Run Kalman smoother to re-estimate the common factors taking the
% coefficient matrices as given. If `allobserved` is true then the VAR
% smoother makes an assumption that the factors are uniquely determined
% whenever all observables are available; this is only true if the
% idiosyncratic covariance matrix is not scaled down.
allobserved = this.cross == 1 && options.cross == 1;
[x,Px,e,u,y,Py,yindex] = timedom.varsmoother( ...
    this.A,this.B,[],this.C,[],1,Sigma,y,[],x0,P0,invfunc,allobserved, ...
    options.tolerance,options.persist);

if options.meanonly
    Px = Px(:,:,[],[]);
    Py = [];
end

if nargout > 1
    ynames = get(this,'ynames');
    [y,Py] = FAVAR.destandardise(this.Mean,this.Std,y,Py);
    d = VAR.outputdata(this,outputformat,range,y,Py,ynames);
end

if nargout > 2
    % Common components.
    [c,Pc] = FAVAR.cc(this.C,x(1:nx,:,:),Px(1:nx,1:nx,:,:));
    [c,Pc] = FAVAR.destandardise(this.Mean,this.Std,c,Pc);
    c = VAR.outputdata(this,outputformat,range,c,Pc,ynames);
end

if nargout > 3
    % Factors.
    f = VAR.outputdata(this,'tseries',range,x(1:nx,:,:),Px(1:nx,1:nx,:,:));
    if ~options.meanonly
        % Means and MSEs that can be used as initial condition.
        lastobs = find(any(yindex ~= 0,1),1,'last');
        if isempty(lastobs)
            index = 1 : nper;
        else
            index = lastobs : nper;
        end
        f.init = {x(:,index,:),Px(:,:,index,:),range(index)};
    end
end

if nargout > 4
    u = FAVAR.destandardise(0,this.Std,u);
    u = VAR.outputdata(this,outputformat,range,u,NaN,ynames);
end

if nargout > 5
    template = tseries();
    e = replace(template,permute(e,[2,1,3]),range(1));
    if ~options.meanonly
        e = struct( ...
            'mean',e, ...
            'std',replace(template,zeros([0,size(e,1),size(e,3)]),range(1)) ...
            );
    end
end

end
