% QUICK_VAR: Compute VAR coefficients, including controls if desired.
%
% usage
%
% [beta,beta1,beta2 c, mu] = quick_var(Y,nlag, const, ctrl, nl)
%
% where
%
% Y     = vector of observables in VAR
% nlag  = number of lags in VAR (1,2,...)
% const = dummy if a constant is desired
% ctrl  = a vector of controls
% nl    = number of lags (0,1,2,..) for controls
%
% beta  = the full set of coefficients
% beta1 = the set of coeffs related to the VAR
% beta2 = the set of coeffs related to the controls


function [beta, beta1, beta2, c, mu ] = quick_var(Y,nlag, const, ctrl, nl)



%Compute max lag
if nargin>3
    mlag = max(nlag,nl);
else
    mlag = nlag;
end

%Number of Variables
nvar = size(Y,2);
T = length(Y)-mlag;


%*****************
% Arrange Data
%*****************

%Stack lagged Y
X = zeros(T, nlag*nvar);
for j =1:nlag
    X(:, (j-1)*nvar+1:j*nvar) = Y(mlag+1-j:end-j,:);
end

if nargin > 3 
    %Stack lagged controls
    nc = size(ctrl,2);
    CT = zeros(T, nl*nc);
    for j = 0:nl
        CT(:,j*nc+1:(j+1)*nc) = ctrl(mlag+1-j:end-j,:);
    end

    %Include a constant in VAR
    X = [X, CT, ones(const*T,1)];
else
    %No controls
    X = [X, ones(const*T,1)];
end

%Exclude lags from DY
YY = Y(mlag+1:end,:);


%*********************************************
% Regress, get covariance matrix and cholesky
%*********************************************

%VAR coefficients, stderrs, etc
beta = (X'*X)\(X'*YY);
mu = YY - X*beta;
omega = mu'*mu/(T);%-nvar*nlag - nrank - 1);
c = chol(omega)';
beta = beta(1:end-const,:);
beta1 = beta(1:nlag*nvar,:);
beta2 = beta(nlag*nvar+1:end,:);
