function V = VAR(THIS,SELECT,RANGE,varargin)
% VAR  Population VAR for selected model variables.
%
% Syntax
% =======
%
%     V = VAR(M,SELECT,RANGE,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Solved model object.
%
% * `SELECT` [ cellstr | char ] - List of variables selected for the VAR.
%
% * `RANGE` [ numeric ] - Hypothetical range, including pre-sample initial
% condition, on which the VAR would be estimated.
%
% Output arguments
% =================
%
% * `V` [ VAR ] - Asymptotic reduced-form VAR for selected model variables.
%
% Options
% ========
%
% * `'order='` [ numeric | *1* ] - Order of the VAR.
%
% * `'constant='` [ *`true`* | `false` ] - Include in the VAR a constant vector
% derived from the steady state of the selected variables.
%
% Description
% ============
%
% Example
% ========
%

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

% Parse required arguments.
P = inputParser();
P.addRequired('m',@ismodel);
P.addRequired('list',@(x) ischar(x) || iscellstr(x));
P.addRequired('range',@isnumeric);
P.parse(THIS,SELECT,RANGE);

% Parse options.
opt = passvalopt('model.VAR',varargin{:});

% Convert char list to cellstr.
if ischar(SELECT)
    SELECT = regexp(SELECT, ...
        '[a-zA-Z][\w\(\)\{\}\+\-]*','match');
end

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

nalt = size(THIS.Assign,3);
nz = length(SELECT);
p = opt.order;
C = acf(THIS,opt.acf{:},'order=',p,'output=','numeric');
RANGE = RANGE(1) : RANGE(end);
nper = length(RANGE);
nk = double(opt.constant);

% Find the position of selected variables in the sspace vector and in the
% model names.
[sspacepos,namepos] = myfindsspacepos(THIS,SELECT,'error');

C = C(sspacepos,sspacepos,:,:);
zbar = permute(THIS.Assign(1,namepos,:),[2,3,1]);
isLog = THIS.log(1,namepos);
zbar(isLog) = log(zbar(isLog));

% TODO: Calculate Sigma.
V = VAR();
V.A = nan(nz,nz*p,nalt);
V.K = zeros(nz,nalt);
V.Omega = nan(nz,nz,nalt);
V.Sigma = [];
V.G = nan(nz,0,nalt);
V.range = RANGE;
V.fitted = true(1,nper);
V.fitted(1:p) = false;
V.nhyper = nz*(nk+p*nz);

for ialt = 1 : nalt
    Ci = C(:,:,:,ialt);
    zbari = zbar(:,ialt);
    
    % Put together moment matrices.
    % M1 := [C1, C2, ...]
    M1 = reshape(Ci(:,:,2:end),nz,nz*p);
    
    % M0 := [C0, C1, ...; C1', C0, ... ]
    % First, produce M0' : = [C0, C1', ...; C1, C0, ...].
    M0t = [];
    for i = 0 : p-1
        M0t = [M0t; ...
            nan(nz,nz*i),reshape(Ci(:,:,1:p-i),nz,nz*(p-i)) ...
            ]; %#ok<AGROW>
    end
    M0 = M0t.';
    nanindex = isnan(M0t);
    M0t(nanindex) = M0(nanindex); %#ok<AGROW>
    % Then, tranpose M0' to get M0.
    M0 = M0t.';
    
    % Compute transition matrix.
    Ai = M1 / M0;
    
    % Estimate cov matrix of residuals.
    Omgi = Ci(:,:,1) - M1*Ai.' - Ai*M1.' + Ai*M0*Ai.';
    
    % Calculate constant vector.
    if opt.constant
        Ki = sum(poly.var2poly(Ai),3)*zbari;
    end
    
    % Populate VAR properties.
    V.A(:,:,ialt) = Ai;
    V.K(:,ialt) = Ki;
    V.Omega(:,:,ialt) = Omgi;
end

V = myschur(V);
V = myinfocrit(V);

end