function varout = varest(y,var_par,smpl_par)

% Computes VAR and covariance matrix of estimated parameters
%    
% Input:
%   y = txn matrix of series
%   var_par:
%     (1) nlag = number of VAR lags
%     (2) icomp = 0 Do not compute Companion Matrices
%                 1 Compute companion matrix ignoring contsant term
%                 2 Compute Companion matrix adding constant as last element of state -- loadings are zero if iconst = 0
%     (3) iconst  = 0 Do not include constant
%                   1 Include constant   		   
%     Note: icomp = 2 ... constant term is added to companion state vector (if iconst == 1)
%   smpl_par: sampling parameters
%      (1) calvec: Calendar vector
%      (2) nper:   number of periods per year
%      (3) nfirst: First observation for estimation
%      (4) nlast:  Last observation for estimation
%    		            
% Output:
%   varout
%      (1) betahat= estimated VAR coefficients (including constant if i_const == 1)
%                   each column has coefficients for a single equation
%      (2) seps = innovation covariance matrix
%      (3) resid =  VAR residuals
%      (4) coef: Companion form matrices
%           (a) Q
%           (b) M
%           (c) G
%           Companion Form Parameter, computed in icomp = 1;
%           Q, M, G for model written in SS form
%           y(t) = Q*z(t)
%           z(t) = M*z(t-1) + G*u(t)
%           var(u(t)) = I

% Extract parameters
nlag   = var_par.nlag;
icomp  = var_par.icomp;
iconst = var_par.iconst;
% % % q      = var_par.q;
% Set Up VAR
ns = size(y,2);
T = size(y,1);
x = [ones(T,1),lagmatrix(y,1:nlag)];
if iconst == 0;
 x = x(:,2:end);  % Eliminate Constant if i_const == 0;
end;

if isempty(smpl_par)==0
    [istart, iend] = smpl_HO(smpl_par);
    istart = max(istart,1);
    iend = min(iend,size(y,1));
else
    istart=1;
    iend=size(y,1);
end




trnd = (1:1:T)';
trnd = trnd(istart:iend);
y = y(istart:iend,:);
x = x(istart:iend,:);
tmp = packr([y x trnd]);
y = tmp(:,1:ns);
x = tmp(:,ns+1:end-1);
trnd = tmp(:,end);

betahat = x\y;
e = y - x*betahat;
ndf=size(x,1)-size(x,2);
seps=(e'*e)/ndf;
% seps=cov(e);
resid = NaN*zeros(T,ns);
resid(trnd,:)=e;

if icomp > 0;
   %  
   %    Transform the VAR so that it is written in Standard form as:
   %    s(t)=P1*s(t-1) + P2*s(t-2) + ... + Pvarlag*s(t-varlag) + eps(t)
   %  

   % ---- Calculate Companion Matrix ---- ;
   if iconst == 0
      b = betahat';
      const_coef = zeros(ns,1);
   elseif iconst == 1
      b=betahat(2:end,:)';
      const_coef = betahat(1,:)';      % Coefficients on Constant Term ;
   else
      error('Invalid value of i_const in VAREST');
   end;
   comp=zeros(size(b,2),size(b,2));
   comp(1:size(b,1),:)=b;
   if size(b,2) > size(b,1);
     comp(size(b,1)+1:end,1:end-size(b,1))=eye(size(b,2)-size(b,1));
   end;
   %   -- Write Model in SS Form --
   %      y(t) = Q*z(t)
   %      z(t) = M*z(t-1) + G*u(t)
   %      var(u(t)) = I
   %  
   M=comp;
   Q=zeros(ns,size(M,1));
   Q(1:ns,1:ns)=eye(ns);
   
%    [U,V,~]=svd(seps);
%    G=zeros(size(M,1),q);
%    G(1:ns,1:q)=U(:,1:q)*diag(sqrt(diag(V(1:q,1:q))));
%    [ K, MM ] = eigs(seps, q, 'LM');
% % %    if q<ns
% % %        G=zeros(size(M,1),q);
% % %        seps_eta=seps(1:q,1:q);
% % %        R=[eye(q);...
% % %            seps(q+1:end,1:q)*inv(seps_eta)];
% % %        G(1:ns,1:q)=R*(chol(seps_eta))';       
% % %    else
       G=zeros(size(M,1),ns);
       G(1:ns,1:ns)=(chol(seps))';
% % %    end
   
%    const_coef=[const_coef;zeros(size(b,2)-size(b,1),1)];
end;

if icomp == 2;  % Add Constant Term 
 G=[G ; zeros(1,size(G,2))];
 Q=[Q zeros(size(Q,1),1)];
 M=[M  zeros(size(M,1),1)];
 M=[M ; zeros(1,size(M,2))];
 M(end,end)=1;
 M(1:ns,end)=const_coef;
end;

% SAVE OUTPUT
varout.betahat = betahat;
varout.seps = seps;
varout.resid = resid;
varout.coef.Q = Q;
varout.coef.M = M;
varout.coef.G = G;
varout.coef.const_coef  = const_coef ;


end