function DATA = bn(THIS,DATA,RANGE,varargin)
% bn  Beveridge-Nelson trends.
%
% Syntax
% =======
%
%     D = bn(M,D,RANGE,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Solved model object.
%
% * `DATA` [ struct | cell ] - Input data on which the BN trends will be
% computed.
%
% * `RANGE` [ numeric ] - Date range on which the BN trends will be
% computed.
%
% Output arguments
% =================
%
% * `D` [ struct | cell ] - Output data with the BN trends.
%
% Options
% ========
%
% * `'deviations='` [ `true` | *`false`* ] - Input and output data are
% deviations from balanced-growth paths.
%
% * `'dtrends='` [ *`'auto'`* | `true` | `false` ] - Measurement variables
% in input and output data include deterministic trends specified in
% [`!dtrends`](modellang/dtrends) equations.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright 2007-2012 Jaromir Benes.

% Parse required input arguments.
P = inputParser();
P.addRequired('m',@ismodel);
P.addRequired('data',@(x) isstruct(x) || iscell(x));
P.addRequired('range',@(x) isnumeric(x));
P.parse(THIS,DATA,RANGE);

options = passvalopt('model.bn',varargin{:});

% Auto set the 'dtrends' option.
if ischar(options.dtrends)
    options.dtrends = ~options.deviation;
end

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

tol = getrealsmall();
ny = length(THIS.solutionid{1});
ne = length(THIS.solutionid{2});
nx = length(THIS.solutionid{3});
nb = size(THIS.solution{1},2);
nf = nx - nb;
nalt = size(THIS.Assign,3);
nper = numel(RANGE);

a = datarequest('alpha',THIS,DATA,RANGE);
ndata = size(a,3);

nloop = max([ndata,nalt]);

% Pre-allocate output data.
DATA = { ...
    nan(ny,nper,nloop), ...
    nan(nx,nper,nloop), ...
    zeros(ne,nper,nloop), ...
    RANGE(1) : RANGE(end), ...
    };

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

repeat = ones([1,nper]);
nondiffstationary = [];

for iloop = 1 : nloop
    
    if nansolution(min([iloop,nalt]))
        continue
    end
    
    if iloop <= nalt
        Tfi = THIS.solution{1}(1:nf,:,iloop);
        Tai = THIS.solution{1}(nf+1:end,:,iloop);
        abseig = abs(ordeig(Tai));
        nunit = sum(abs(abseig - 1) <= tol);
        if ~iseye(Tai(1:nunit,1:nunit))
            nondiffstationary(end+1) = iloop; %#ok<AGROW>
        end
        Zi = THIS.solution{4}(:,:,iloop);
        % non-stationary xf and y
        %xfdiffuse = any(abs(Tfi(:,1:nunit)) > tol,2);
        %ydiffuse = any(abs(Zi(:,1:nunit)) > tol,2);
        % steady state for [xf;alpha]
        if ~options.deviation
            Kai = THIS.solution{3}(nf+1:end,1,iloop);
            abar = zeros([nb,1]);
            abar(nunit+1:end) = (eye(nb-nunit) - Tai(nunit+1:end,nunit+1:end)) \ Kai(nunit+1:end);
            abar = abar(:,repeat);
            Kfi = THIS.solution{3}(1:nf,1,iloop);
            Di = THIS.solution{6};
            Kfi = Kfi(:,repeat);
            Di = Di(:,repeat);
        end
        if options.dtrends
            Wi = mydtrendsrequest(THIS,'range',RANGE,iloop);
        end
    end
    
    tmpindex = min([iloop,ndata]);
    ai = a(:,:,tmpindex);
    if ~options.deviation
        ai = ai - abar;
    end
    
    %**********************************************************************
    % Forward cumsum of stable alpha.
    
    acum = ...
        (eye(nb-nunit) - Tai(nunit+1:end,nunit+1:end)) \ ai(nunit+1:end,:);
    
    %**********************************************************************
    % Beveridge Nelson for non-stationary variables.
    
    ai(1:nunit,:) = ai(1:nunit,:) + ...
        Tai(1:nunit,nunit+1:end)*acum;
    if options.deviation
        ai(nunit+1:end,:) = 0;
    else
        ai(nunit+1:end,:) = abar(nunit+1:end,:);
    end
    
    xf = Tfi*ai;
    y = Zi*ai;
    if ~options.deviation
        xf = xf + Kfi;
        y = y + Di;
    end
    if options.dtrends
        y = y + Wi;
    end
    
    DATA{1}(:,:,iloop) = y;
    DATA{2}(:,:,iloop) = [xf;ai];
    
end

% Parameterisations that are not difference-stationary.
if ~isempty(nondiffstationary)
    utils.error('model', ...
        ['Cannot run BN on models or parameterisations ', ...
        'that are not difference-stationary:%s'], ...
        strfun.alt2str(nondiffstationary));
end

% Keep only current dates so that no pre-sample observations are added.
imagid = imag(THIS.solutionid{2});
DATA{2}(imagid ~= 0,:,:) = NaN;
DATA = myalpha2xb(THIS,DATA);

DATA = dp2db(THIS,DATA);

end