function [THIS,DATA,FITTED,RR,COUNT] = estimate(THIS,varargin)
% estimate  Estimate a reduced-form VAR or BVAR.
%
% Syntax with database input
% ===========================
%
%     [V,DATA,FITTED] = estimate(W,D,RANGE,...)
%     [V,DATA,FITTED] = estimate(W,D,LIST,RANGE,...)
%
% Syntax with tseries input
% ==========================
%
%     [V,DATA,FITTED] = estimate(W,X,RANGE,...)
%
% Input arguments
% ================
%
% * `V` [ VAR ] - Empty VAR object.
%
% * `D` [ struct ] - Input database.
%
% * `LIST` [ cellstr ] - List of series from the database on which the VAR
% will be estimated; `LIST` can be omitted if variable names have been
% already specified at the time of creating the VAR object.
%
% * `X` [ tseries ] - Tseries objects with input data.
%
% * `RANGE` [ numeric ] - Estimation range, including `P` pre-sample
% observations where `P` is the order of the VAR.
%
% Output arguments
% =================
%
% * `W` [ VAR ] - Estimated reduced-form VAR object.
%
% * `DATA` [ struct | tseries ] - Output data that include the endogenous
% variables and the estimated residuals; the output data is a struct
% (database) when `estimate` is called with a database input, and the
% output data is a tseries when `estimate` is called with a tseries input.
%
% * `FITTED` [ numeric ] - Periods for which fitted values can be
% calculated.
%
% Options
% ========
%
% * `'A='` [ numeric | *empty* ] - Restrictions on the individual values in
% the transition matrix, `A`.
%
% * `'BVAR='` [ numeric ] - Prior dummy observations for estimating a BVAR;
% construct the dummy observations using the one of the `BVAR` functions.
%
% * `'C='` [ numeric | *empty* ] - Restrictions on the individual values in
% the constant vector, `C`.
%
% * `'diff='` [ `true` | *`false`* ] - Difference the series before estimating
% the VAR; integrate the series back afterwards.
%
% * `'G='` [ numeric | *empty* ] - Restrictions on the individual values in
% the matrix at the co-integrating vector, `G`.
%
% * `'cointeg='` [ numeric | *empty* ] - Co-integrating vectors (in rows)
% that will be imposed on the estimated VAR.
%
% * `'constraints='` [ char ] - General linear constraints on the VAR
% parameters.
%
% * `'constant='` [ *`true`* | `false` ] - Include a constant vector in the VAR.
%
% * `'covParameters='` [ `true` | *`false`* ] - Calculate the covariance matrix
% of estimated parameters.
%
% * `'eqtnByEqtn='` [ `true` | *`false`* ] - Estimate the VAR equation by
% equation.
%
% * `'maxIter='` [ numeric | *1* ] - Maximum number of iterations when
% generalised least squares algorithm is involved.
%
% * `'mean='` [ numeric | *empty* ] - Assume a particular mean of the VAR
% process.
%
% * `'order='` [ numeric | *1* ] - Order of the VAR.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% * `'schur='` [ *`true`* | `false` ] - Calculate triangular (Schur)
% representation of the estimated VAR.
%
% * `'stdize='` [ `true` | *`false`* ] - Adjust the prior dummy observations by
% the std dev of the observations.
%
% * `'tolerance='` [ numeric | *`1e-5`* ] - Convergence tolerance when
% generalised least squares algorithm is involved.
%
% * `'yNames='` [ cellstr | function_handle | *`@(n) sprintf('y%g',n)`* ] -
% Use these names for the VAR variables.
%
% * `'eNames='` [ cellstr | function_handle | *`@(yname,n) ['res_',yname]`*
% ] - Use these names for the VAR residuals.
%
% * `'warning='` [ *`true`* | `false` ] - Display warnings produced by this
% function.
%
% Description
% ============
%
% If the VAR is estimated with a database input, `D`, the variable names,
% `LIST`, will be stored within the VAR object, and used whenever you call
% a VAR or SVAR function that returns output data, such as
% [`VAR/simulate`](VAR/simulate) or [`SVAR/srf`](SVAR/srf).
%
% If the VAR is estimated with a tseries input, `X`, the VAR functions will
% return tseries as output data.
%
% Example
% ========
%

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

% Get input data; the user range is supposed to *include* the pre-sample
% initial condition.
[flag,y,range,ynames,inputformat,varargin] = ...
    VAR.inputdata(THIS,varargin{:});
ny = size(y,1);
nper = size(y,2);
ndata = size(y,3);

% Invalid input data.
if ~flag
    utils.error('VAR','Invalid input data.');
end

% Pass and validate options.
opt = passvalopt('VAR.estimate',varargin{:});

% Set names of variables and residuals.
if length(opt.ynames) == ny
    ynames = opt.ynames;
end
if length(opt.enames) == ny
    enames = opt.enames;
else
    enames = {};
end
THIS = setnames(THIS,ynames,enames);

% Both 'const' and 'constant' allowed for bkw cmp.
if ~isempty(opt.const)
    opt.constant = opt.const;
end

if strcmpi(opt.output,'auto')
    outputformat = inputformat;
else
    outputformat = opt.output;
end

if ~isempty(opt.userdata)
    THIS = userdata(THIS,opt.userdata);
end

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

THIS.range = range;

% All data sets must have the same structure of missing observations.
if ndata > 1
    nanindex = isnan(y);
    allnan = all(nanindex,3);
    anynan = any(nanindex,3);
    if any(anynan & ~allnan)
        VAR.error(32);
    end
end

if ~isempty(opt.cointeg)
    opt.diff = true;
end

if ~isempty(opt.mean)
    if length(opt.mean) == 1
        opt.mean = opt.mean(ones([ny,1]));
    else
        opt.mean = opt.mean(:);
    end
end

if ~isempty(opt.mean)
    opt.constant = false;
end

% Arrange data on LHS and RHS.
[y0,k0,y1,g1,ci] = VAR.stackdata(y,opt);
ng = size(g1,1);
nk = size(k0,1);

% Read parameter restrictions, and set up their hyperparameter form.
% They are organised as follows:
% * Rr = [R,r],
% * beta = R*gamma + r.
THIS.Rr = VAR.restrict(ny,nk,ng,opt);

% Get the number of hyperparameters.
if isempty(THIS.Rr)
    % Unrestricted VAR.
    if ~opt.diff
        % Level VAR.
        THIS.nhyper = ny*(nk+opt.order*ny+ng);
    else
        % Difference VAR or VEC.
        THIS.nhyper = ny*(nk+(opt.order-1)*ny+ng);
    end
else
    % Parameter restrictions in the hyperparameter form:
    % beta = R*gamma + r;
    % The number of hyperparams is given by the number of columns of R.
    % The Rr matrix is [R,r], so we need to subtract 1.
    THIS.nhyper = size(THIS.Rr,2) - 1;
end

% Number of priors.
nprior = size(opt.bvar,3);

% Total number of alternatives.
nloop = max([ndata,nprior]);

% Estimate reduced-form VAR parameters. The size of coefficient matrices
% will always be determined by p whether this is a~level VAR or
% a~difference VAR.
p = opt.order;
resid = nan(ny,nper,nloop);
COUNT = zeros(1,nloop);
doprealloc();

% Do not pass options.bvar into glsq().
bvar = opt.bvar;
opt = rmfield(opt,'bvar');

% Create command-window progress bar.
if opt.progress
    progress = progressbar('IRIS VAR.estimate progress');
end

% Main loop.
use = struct();
use.Rr = THIS.Rr;
use.k0 = k0;
use.ci = ci;
for iloop = 1 : nloop
    if iloop <= ndata
        use.y0 = y0(:,:,iloop);
        use.y1 = y1(:,:,iloop);
        use.g1 = g1(:,:,iloop);
    end
    if isempty(bvar)
        use.bvary0 = [];
        use.bvark0 = [];
        use.bvary1 = [];
        use.bvarg1 = [];
    elseif iloop <= nprior
        temp = bvar(:,:,iloop);
        use.bvary0 = temp(1:ny,:);
        temp(1:ny,:) = [];
        use.bvark0 = temp(1:nk,:);
        temp(1,:) = [];
        if ~opt.diff
            use.bvary1 = temp(1:ny*p,:);
        else
            use.bvary1 = temp(1:ny*(p-1),:);
        end
        temp(1:ny*p,:) = [];
        use.bvarg1 = temp(1:ng,:);
    end
    % Run generalised least squares. We assign the individual properties
    % computed within `glsq` in a separate set of assignments to help trace
    % down errors from dimension mismatches.
    use = VAR.glsq(use,opt);
    THIS.A(:,:,iloop) = use.A;
    THIS.K(1:end,iloop) = use.K;
    THIS.G(:,:,iloop) = use.G;
    THIS.Omega(:,:,iloop) = use.Omg;
    THIS.Sigma(:,:,iloop) = use.Sgm;
    resid(:,p+1:end,iloop) = use.resid;
    COUNT(iloop) = use.count;
    THIS.fitted(1,p+1:end,iloop) = use.fitted;
    if opt.progress
        update(progress,iloop/nloop);
    end
end
% End of main loop.

% Triangular representation.
if opt.schur
    THIS = myschur(THIS);
end

% Information criteria AIC and SBC.
THIS = myinfocrit(THIS);

if size(y,3) == 1 && nloop > 1
    y = y(:,:,ones(1,nloop));
end

% Report observations that could not be fitted.
dochkobsnotfitted();

% Range fitted.
FITTED = cell(1,nloop);
for iloop = 1 : nloop
    FITTED{iloop} = THIS.range(THIS.fitted(:,:,iloop));
end

if nargout > 1
    % Convert data to output time series.
    DATA = VAR.outputdata(THIS,outputformat,THIS.range,[y;resid],[], ...
        [THIS.ynames,THIS.enames]);
end

if nargout > 2
    RR = THIS.Rr;
end

% Nested functions.

%**************************************************************************
    function doprealloc()
        THIS.K = nan(ny,nloop);
        THIS.A = nan(ny,ny*p,nloop);
        THIS.G = nan(ny,ng,nloop);
        THIS.Omega = nan(ny,ny,nloop);
        THIS.T = nan(ny*p,ny*p,nloop);
        THIS.U = nan(ny*p,ny*p,nloop);
        THIS.eigval = nan(1,ny*p,nloop);
        THIS.Sigma = [];
        THIS.fitted = false(1,nper,nloop);
        THIS.aic = nan(1,nloop);
        THIS.sbc = nan(1,nloop);
    end
% doprealloc().

%**************************************************************************
    function dochkobsnotfitted()
        allfitted = all(THIS.fitted,3);
        if opt.warning && any(~allfitted(p+1:end))
            missing = THIS.range(p+1:end);
            missing = missing(~allfitted(p+1:end));
            [~,consec] = datconsecutive(missing);
            utils.warning('VAR', ...
                ['The following period(s) not fitted ', ...
                'because of missing observations: %s.'], ...
                consec{:});
        end
    end
    % dochkobsnotfitted().

end