function VAR = doProxySVARx(VAR)
 
tmp=packr([VAR.m1 VAR.res]);
VAR.T = size(tmp,1);
VAR.n=size(VAR.res,2);
VAR.k=size(VAR.m1,2);
VAR.m=tmp(:,1:VAR.k);
VAR.res1=tmp(:,VAR.k+1:end);
VAR.T1=size(VAR.res1,1);
VAR.Sigma = cov(VAR.res);%(VAR.res'*VAR.res)/(VAR.T-VAR.n*VAR.p-1);

      
% Identification (See also appendix in the paper)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Phib = [ones(length(VAR.m),1) VAR.m]\VAR.res;
Phib = Phib(2:end,:);
VAR.Phib=Phib;
Phib11  = Phib(1:VAR.k,1:VAR.k);
Phib21  = Phib(1:VAR.k,VAR.k+1:VAR.n);
b21ib11 = (Phib11\Phib21)';
Sig11   = VAR.Sigma(1:VAR.k,1:VAR.k);
Sig21   = VAR.Sigma(VAR.k+1:VAR.n,1:VAR.k);
Sig22   = VAR.Sigma(VAR.k+1:VAR.n,VAR.k+1:VAR.n);
ZZp     = b21ib11*Sig11*b21ib11'-(Sig21*b21ib11'+b21ib11*Sig21')+Sig22;
b12b12p = (Sig21- b21ib11*Sig11)'*(ZZp\(Sig21- b21ib11*Sig11));
b11b11p = Sig11-b12b12p;
b22b22p = Sig22+b21ib11*(b12b12p-Sig11)*b21ib11';
b12ib22   = ((Sig21- b21ib11*Sig11)'+b12b12p*b21ib11')/(b22b22p');
b11iSig = eye(VAR.k)/(eye(VAR.k)-b12ib22*b21ib11);
b21iSig = b21ib11*b11iSig;

SigmaTSigmaTp =b11iSig\b11b11p/b11iSig';
VAR.S1=SigmaTSigmaTp;
% Reliability
Sigmm   = VAR.m'*VAR.m/VAR.T;
ED      = eye(VAR.k)*sum(sum(VAR.m,2)~=0)/VAR.T;
mu1     = VAR.m'*VAR.res(:,1:VAR.k)/VAR.T;
PhiPhip = mu1*inv(b11b11p)*mu1';
VAR.RM  = inv(Sigmm)*PhiPhip*inv(ED);
VAR.RMeigs = sort(eig(VAR.RM));

% Impulse Responses
%%%%%%%%%%%%%%%%%%%%%%
for i=1:VAR.k   
    
 if VAR.ord(1)==1;
 s1 = sqrt(SigmaTSigmaTp(1,1));
 a  = SigmaTSigmaTp(2,1)/s1;
 s2 = sqrt(SigmaTSigmaTp(2,2)-a^2);
 SigmaT = [s1 0 ; a s2];
 elseif VAR.ord(1)==2;
 s2 = sqrt(SigmaTSigmaTp(2,2));
 b  = SigmaTSigmaTp(1,2)/s2;
 s1 = sqrt(SigmaTSigmaTp(1,1)-b^2);
 SigmaT = [s1 b; 0 s2];
 end
 %[b11iSig;b21iSig]*SigmaT
 VAR.b1(:,:,i) = [b11iSig;b21iSig]*SigmaT;
 %calculate shocks
 Sigma=VAR.Sigma;
 iSigma=invpd(Sigma);
 B1= VAR.b1(:,:,1);
 A1=iSigma*B1;
 VAR.a1=A1;
 shocks=A1'*VAR.res';
 VAR.shocks=shocks';
        
end



