function [a11iSig,a21iSig,SigmaTSigmaTp ]=my_MR_identification(MyCov,Wc, DO_ZERO)
% Uses the Mertens-Raven code for IV-based identification
% Offers option to zero out eta(1,:) to eliminate indirect effects of second shock

%the covariance between proxies and eps_st is Phi. Phi B' is the covariance of the proxies with the observed
%residuals, i.e. the variance of interest
%             Phib = [ones(T,1) nonemptyIV]\U;
%             Phib = Phib(2:end,:);

%             Phib = [nonemptyIV]\U;
%             Phib = (nonemptyIV-repmat(mean(nonemptyIV),T,1))'*U(:,1:end-1)/nu_T_F;
if ~exist('DO_ZERO', 'var')
    DO_ZERO=false;
end

Phib = MyCov;
mz=size(MyCov,1);
Phib11  = Phib(1:mz,1:mz);
Phib21  = Phib(1:mz,mz+1:size(MyCov,2));
kappa = (Phib11\Phib21)';             %kappa
%             Sigma_i=Wc;
Sig11   = Wc(1:mz,1:mz);
Sig21   = Wc(mz+1:size(MyCov,2),1:mz);
Sig22   = Wc(mz+1:size(MyCov,2),mz+1:size(MyCov,2));
ZZp     = kappa*Sig11*kappa'-(Sig21*kappa'+kappa*Sig21')+Sig22; %simply Z in Appendix
a12a12p = (Sig21- kappa*Sig11)'*(ZZp\(Sig21- kappa*Sig11));
a11a11p = Sig11-a12a12p; %ResidualVarInv=inv([eye(mz); kappa']'*(Wc\[eye(mz); kappa']))
% Wc(1:mz,1:mz)-Wc(1:mz,:)*[kappa; -eye(size(Wc,2)-mz)] *inv([kappa', -eye(size(Wc,2)-mz)]*Wc*[kappa; -eye(size(Wc,2)-mz)]) * [kappa', -eye(size(Wc,2)-mz)] * Wc(:,1:mz)
a22a22p = Sig22+kappa*(a12a12p-Sig11)*kappa';
% a12ia22 is also called eta
eta   = ((Sig21- kappa*Sig11)'+a12a12p*kappa')/(a22a22p'); %my version is correct
SigmaTSigmaTp =( eye(mz)-eta*kappa ) * a11a11p * ( eye(mz)-eta*kappa )'; %a11iSig\a11a11p/a11iSig';

if DO_ZERO
    eta(1,:) = 0;
end
a11iSig = inv( eye(mz)-eta*kappa );%eye(mz)/(eye(mz)-a12ia22*kappa);
a21iSig = kappa/(eye(mz)-eta*kappa); %kappa*a11iSig