function a = glsq(a,opt)
% glsq  [Not a public function] Generalised least squares estimator for reduced-form VARs.
%
% Backend IRIS function.
% No help provided.

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

y0 = a.y0;
k0 = a.k0;
y1 = a.y1;
g1 = a.g1;
ci = a.ci;
bvary0 = a.bvary0;
bvark0 = a.bvark0;
bvary1 = a.bvary1;
bvarg1 = a.bvarg1;
Rr = a.Rr;

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

ny = size(y0,1);
nk = size(k0,1);
ng = size(g1,1);

% Number of lags included in regression; needs to be decreased by one for
% difference VARs or VECs.
p = opt.order;
if opt.diff
    p = p - 1;
end

% Find effective estimation range and exclude NaNs.
fitted = all(~isnan([y0;k0;y1;g1]),1);
nobs = sum(fitted);
y0 = y0(:,fitted);
k0 = k0(:,fitted);
y1 = y1(:,fitted);
g1 = g1(:,fitted);

if ~isempty(opt.mean)
    ymean = opt.mean;
    y0 = y0 - ymean(:,ones(1,nobs));
    y1 = y1 - repmat(ymean(:,ones(1,nobs)),p,1);
end

% Number of dummy observations.
nb = size(bvary0,2);

if opt.stdize && nb > 0
    % Create a matrix of observations (that will be possibly demeaned)
    % including pre-sample initial condition.
    yd = [reshape(y1(:,1),ny,p),y0];
    if nk > 0
        % Demean the observations if the constant is included in the
        % regression.
        ymean = mean(yd,2);
        yd = yd - ymean(:,ones(1,size(yd,2)));
    end
    % Calculate the std dev on the demeaned observations, and adjust the
    % prior dummy observations. This is equivalent to standardizing the
    % observations with given dummies.
    ystd = std(yd,1,2);
    bvary0 = bvary0 .* ystd(:,ones(1,nb));
    bvary1 = bvary1 .* repmat(ystd(:,ones(1,nb)),p,1);
end

% Matrix of explanatory observations.
X = [k0;y1;g1];

% Add prior dummy observations to the LHS and RHS data matrices.
if nb > 0
    y0 = [bvary0,y0];
    X = [[bvark0;bvary1;bvarg1],X];
end

if ~isempty(Rr)
    R = Rr(:,1:end-1);
    r = Rr(:,end);
else
    R = [];
    r = [];
end

% `Omg0` is covariance of residuals based on unrestricted non-bayesian VAR.
% It is used to compute covariance of parameters.
Omg0 = [];
count = 0;
if ~isempty(R) && opt.eqtnbyeqtn
    % Estimate equation by equation with parameter restrictions. No check
    % for cross-equation restrictions is performed. This is the user's
    % responsibility.
    pos = (1:ny).';
    pos = pos(:,ones(1,nk+ny*p+ng));
    pos = pos(:);
    M = X*X.';
    beta = nan(ny*(nk+ny*p+ng),1);
    realsmall = getrealsmall();
    for i = 1 : ny
        % Get restrictions for equation i.
        betaindex = pos == i;
        Ri = R(betaindex,:);
        gammaindex = any(abs(Ri) > realsmall,1);
        Ri = Ri(:,gammaindex);
        ri = r(betaindex);
        % Estimate free hyperparameters.
        c = y0(i,:).' - X.'*ri;
        gammai = (Ri.'*M*Ri) \ (Ri.'*X*c);
        beta(betaindex) = Ri*gammai + ri;
    end
    beta = reshape(beta,[ny,ny*p+nk+ng]);
    e = y0 - beta*X;
    e = e(:,nb+1:end);
    Omg = e*e.' / nobs;
    count = count + 1;
else
    % Test for empty(r) not empty(R). This is because if all parameters are
    % fixed to a number, R is empty but we still need to run LSQ with
    % restrictions.
    if isempty(r)
        % Ordinary least squares for unrestricted VAR or BVAR.
        beta = y0 / X;
        e = y0 - beta*X;
        e = e(:,nb+1:end);
        Omg = e*e.' / nobs;
        Omg0 = Omg;
        count = count + 1;
    else
        % Generalized least squares for parameter restrictions.
        Omg = eye(ny);
        Omgi = eye(ny);
        beta = Inf;
        M = X*X.';
        maxdiff = Inf;
        while maxdiff > opt.tolerance && count <= opt.maxiter
            lastbeta = beta;
            c = y0(:) - kron(X.',eye(ny))*r;
            % Estimate free hyperparameters.
            gamma = (R.'*kron(M,Omgi)*R) \ (R.'*kron(X,Omgi)*c);
            % Compute parameters.
            beta = reshape(R*gamma + r,[ny,ny*p+nk+ng]);
            e = y0 - beta*X;
            e = e(:,nb+1:end);
            Omg = e*e.' / nobs;
            Omgi = inv(Omg);
            maxdiff = max(abs(beta(:) - lastbeta(:)));
            count = count + 1;
        end
    end
end

% Covariance of parameter estimates, not available for VECM and diff VARs.
Sgm = [];
if opt.covparameters && ~opt.diff
    dosigma();
end

% Constant vector.
if nk > 0
    K = beta(:,1:nk);
    beta(:,1:nk) = [];
else
    K = zeros(ny,1);
end

% Transition matrices.
A = beta(:,1:ny*p);
beta(:,1:ny*p) = [];

% Convert VEC to co-integrated VAR.
if opt.diff
    % Coefficients on co-integrating vectors.
    G = beta(:,1:ng);
    K = K + G*ci(:,1);
    A = reshape(A,ny,ny,p);
    A = poly.polyprod(A,cat(3,eye(ny),-eye(ny)));
    A = poly.polysum(A,eye(ny)+G*ci(:,2:end));
    p = p + 1;
    A = reshape(A,ny,ny*p);
else
    G = zeros(ny,0);
end

% Add mean to the VAR process.
if ~isempty(opt.mean)
    K = K + (eye(ny) - sum(reshape(A,ny,ny,p),3))*ymean;
end

a.A = A;
a.K = K;
a.G = G;
a.Omg = Omg;
a.Sgm = Sgm;
a.resid = nan(ny,length(fitted));
a.resid(:,fitted) = e;
a.count = count;
a.fitted = fitted;

% Nested functions.

    function dosigma()
        % Asymptotic covariance of parameters is based on covariance of residuals
        % computed from a non-restricted, non-bayesian VAR. The risk exists that we
        % bump into singularity or near-singularity.
        if isempty(Omg0)
            if ~isempty(X)
                beta0 = y0 / X;
                e0 = y0 - beta0*X;
                e0 = e0(:,nb+1:end);
                Omg0 = e0*e0.' / nobs;
            else
                Omg0 = nan(ny);
            end
        end
        if isempty(r)
            % Unrestricted parameters, `M` may not be available.
            if ~isempty(X)
                M = X*X.';
                Sgm = kron(inv(M),Omg0);
            else
                Sgm = nan(size(X,1)*ny);
            end
        elseif ~isempty(R)
            % If `R` is empty, all parameters are fixed, and we do not have to
            % calculate `Sgm`. If not, then `M` and `Omegainv` are guaranteed
            % to exist.
            if ~isempty(X)
                Sgm = R*((R.'*kron(M,inv(Omg0))*R) \ R.');
            else
                Sgm = nan(size(X,1)*ny);
            end
        end
    end

end