function [B,COVRES,R2] = regress(THIS,LHS,RHS,varargin)
% regress  Centred population regression for selected model variables.
%
% Syntax
% =======
%
%     [B,COVRES,R2] = regress(M,LHS,RHS,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Model on whose covariance matrices the popolation
% regression will be based.
%
% * `LHS` [ char | cellstr ] - LHS variables in the regression; each of the
% variables must be part of the state-space vector.
%
% * `RHS` [ char | cellstr ] - RHS variables in the regression; each of the
% variables must be part of the state-space vector, or must refer to a
% larger lag of a transition variable present in the state-space vector.
%
% Output arguments
% =================
%
% * `B` [ namedmat | numeric ] - Population regression coefficients.
%
% * `COVRES` [ namedmat | numeric ] - Covariance matrix of residuals from
% the population regression.
%
% * `R2` [ numeric ] - Coefficient of determination (R-squared).
%
% Options
% ========
%
% * `'output='` [ *`'namedmat'`* | `'numeric'` ] - Output matrices will be
% either namedmat objects or plain numeric arrays.
%
% Description
% ============
%
% Population regressions calculated by this function are always centred.
% This means the regressions are always calculated as if estimated on
% observations with their uncondional means (the steady-state levels)
% removed from them.
%
% The LHS and RHS variables that are log-variables must include
% `log(...)` explicitly in their names. For instance, if `X` is declared
% to be a log variable, then you must refer to `log(X)` or `log(X{-1})`.
%
% Example
% ========
%
%     [B,C] = regress('log(R)',{'log(R{-1})','log(dP)'});
%

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

opt = passvalopt('model.regress',varargin{:});
isnamedmat = strcmpi(opt.output,'namedmat');

if ischar(LHS)
    LHS = regexp(LHS,'[\w\(\)\{\}\d+\-]+','match');
end

if ischar(RHS)
    RHS = regexp(RHS,'[\w\(\)\{\}\d+\-]+','match');
end

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

% Remove blank spaces from user names.
LHS = regexprep(LHS(:).','\s+','');
RHS = regexprep(RHS(:).','\s+','');

nalt = size(THIS.Assign,3);
% `lhspos` is a vector of positions in [solutionid{1:2}];
[~,~,lhspos] = myfindsspacepos(THIS,LHS,'error');
% `rhspos` is a vector of either positions in `[solutionid{1:2}]` (when
% imag is zero) or the position of the max lag of that variables in
% `[solutionid{1:2}]` with imag being the distance to the lag requested.
% For example, if `x` enters `[solutionid{1:2}]` as `x{-3}` at maximum, and
% requested is x{-5}, `rhspos` is the position of `x{-3}` and imag is -2.
% The minimum imag also determines the order up to which ACF needs to be
% calculated.
[~,~,rhspos] = myfindsspacepos(THIS,RHS,'error');
nlhs = length(lhspos);
nrhs = length(rhspos);

p = -min([imag(rhspos),imag(lhspos)]);
C = acf(THIS,opt.acf{:},'order=',p,'output=','numeric');
nc = size(C,1);

% Convert `lhspos` and `rhspos` to positions in
% `[solutionid{1:2},solutionid{1:2}{-1},...]`.
lhspos = real(lhspos) - nc*imag(lhspos);
rhspos = real(rhspos) - nc*imag(rhspos);

B = nan(nlhs,nrhs,nalt);
COVRES = nan(nlhs,nlhs,nalt);
R2 = nan(nlhs,nalt);
for ialt = 1 : nalt
    row = reshape(C(:,:,:,ialt),[nc,nc*(p+1)]);
    % CC := [C0,C1,C2,...;C1.',C0,C1,...;C2.',C1.',C0,...;...].
    CC = row;
    for k = 1 : p
        row = [C(:,:,k+1).',row(:,1:end-nc)];
        CC = [CC;row]; %#ok<AGROW>
    end
    % Y = B*X + RES.
    % YY is the cov matrix of LHS variables.
    % YX is the cov matrix of LHS versus RHS variables.
    % XX is the cov matrix of RHS variables.
    YY = CC(lhspos,lhspos);
    YX = CC(lhspos,rhspos);
    XX = CC(rhspos,rhspos);
    [B(:,:,ialt),COVRES(:,:,ialt)] = covfun.popregress(YY,YX,XX);
    R2(:,ialt) = 1 - diag(COVRES(:,:,ialt))./diag(YY);
end

if isnamedmat
    B = namedmat(B,LHS,RHS);
    COVRES = namedmat(COVRES,LHS,LHS);
end

end