function [X,Xc,e,inx,nz,IIE]=reduce_sys(A);
% Reduce transition matrix A of system y=A*y(-1)+M*epsilon to core (linearly
% independent) set of variables using QR decomposition 
% INPUTS: Square Matrix A
% OUTPUTS: Reduced (core) matrix X
%            Ordering matrix E (new-variables=old-variables(e))
%          Index of core variables
%          vector of non-zero columns (some variables might be place
%          holders of lagged innovations)
% NOTE THAT  E'*A*E=E'Q*R; so E'y=E'Q*R+E'M*epsilon

[Q,R,e]=qr(A,'vector'); % 
II=eye(size(R));
IIE=II(:,e);% readjust variables to get them in order A.
inx=1:rank(R);

nz=find(sum(abs(R),1)>1e-10);

Qtild=IIE'*Q;
R11=R(nz,nz);
Qtild11=Qtild(nz,nz);
X=Qtild11*R11;
Xc=Qtild(nz(end)+1:end,nz)*R11;

