function [IRF_IV, IRF_IV_noDir, IRF_IV_noInd, IRF_IV_Neither]=my_IRF_IV(MyCov,Wc,IV, MySwitch)
% calls Mertens-Raven or my Choleski-based IV identification, depending on value on MySwitch.
% Returns vector of impact IRF

switch MySwitch
    case 1 
       [b11iSig,b21iSig,SigmaTSigmaTp ]=my_MR_identification(MyCov,Wc);
       SigmaT=chol(SigmaTSigmaTp,'lower');

       IRF_IV= [b11iSig;b21iSig]*SigmaT;
       noDir_SigmaT = SigmaT;
       noDir_SigmaT(1,2) = 0;
       IRF_IV_noDir= [b11iSig;b21iSig]*noDir_SigmaT;
       [noInd_b11iSig,noInd_b21iSig, ~ ]=my_MR_identification(MyCov,Wc, true);
       IRF_IV_noInd = [noInd_b11iSig;noInd_b21iSig]*SigmaT;
       IRF_IV_Neither = [noInd_b11iSig;noInd_b21iSig]*noDir_SigmaT;

    case 2
        [b11iSig,b21iSig,SigmaTSigmaTp ]=my_MR_identification(MyCov,Wc);
        s1 = sqrt(SigmaTSigmaTp(1,1));
        rho_a  = SigmaTSigmaTp(2,1)/s1;
        s2 = sqrt(SigmaTSigmaTp(2,2)-rho_a^2);
        SigmaT = [s1 0 ; rho_a s2];

       IRF_IV= [b11iSig;b21iSig]*SigmaT;
       noDir_SigmaT = SigmaT;
       noDir_SigmaT(1,2) = 0;
       IRF_IV_noDir= [b11iSig;b21iSig]*noDir_SigmaT;
       [noInd_b11iSig,noInd_b21iSig, ~ ]=my_MR_identification(MyCov,Wc, true);
       IRF_IV_noInd = [noInd_b11iSig;noInd_b21iSig]*SigmaT;
       IRF_IV_Neither = [noInd_b11iSig;noInd_b21iSig]*noDir_SigmaT;
    case 3
        [b11iSig,b21iSig,SigmaTSigmaTp ]=my_MR_identification(MyCov,Wc);
        s2 = sqrt(SigmaTSigmaTp(2,2));
        rho2  = SigmaTSigmaTp(1,2)/s2;
        s1 = sqrt(SigmaTSigmaTp(1,1)-rho2^2);
        SigmaT = [s1 rho2; 0 s2];

       IRF_IV= [b11iSig;b21iSig]*SigmaT;
       noDir_SigmaT = SigmaT;
       noDir_SigmaT(1,2) = 0;
       IRF_IV_noDir= [b11iSig;b21iSig]*noDir_SigmaT;
       [noInd_b11iSig,noInd_b21iSig, ~ ]=my_MR_identification(MyCov,Wc, true);
       IRF_IV_noInd = [noInd_b11iSig;noInd_b21iSig]*SigmaT;
       IRF_IV_Neither = [noInd_b11iSig;noInd_b21iSig]*noDir_SigmaT;
    case -3
        pv1=IV;
        pz=size(MyCov,1);
        pv2=size(MyCov,2)-pz-pv1;
        CovV2_Z=MyCov(:,pv1+1:end-pz)';
        CovV1_Z=MyCov(:,1:pv1)';
        VarV1=Wc(1:pv1,1:pv1);
        VarV2=Wc(pv1+1:end,pv1+1:end);
        CovV1_V2=Wc(1:pv1,pv1+1:end);
        VarZ=MyCov(:,end-(pz-1):end);
        kappa=CovV2_Z * inv(VarZ) * CovV1_Z' * inv( CovV1_Z * inv(VarZ) * CovV1_Z') ;
        BB=VarV2-CovV1_V2'*kappa';
        v2rTv2r = VarV2 - (CovV1_V2'*kappa' + kappa * CovV1_V2) + kappa*VarV1*kappa';
        eta= inv(BB * inv(v2rTv2r ) * BB') * BB * inv(v2rTv2r) * (CovV1_V2'-kappa*VarV1) ; 
        
        DD = eta' * inv(eye(pv2)-kappa*eta') * kappa;
        vS1S1T = (eye(pv1^2) + kron(eye(pv1), DD) + kron(DD, eye(pv1)) ) \ reshape(VarV1 - eta' * VarV2 * eta, pv1^2,1);
        SigmaTSigmaTp = reshape(vS1S1T, pv1, pv1) ;

        s2 = sqrt(SigmaTSigmaTp(2,2));
        rho2  = SigmaTSigmaTp(1,2)/s2;
        s1 = sqrt(SigmaTSigmaTp(1,1)-rho2^2);
        SigmaT = [s1 rho2; 0 s2];

        IRF_IV= [inv(eye(pv1)-eta'*kappa)*SigmaT;inv(eye(pv2)-kappa*eta')*kappa*SigmaT];        
        noDir_SigmaT = SigmaT;
        noDir_SigmaT(1,2) = 0;
        IRF_IV_noDir= [inv(eye(pv1)-eta'*kappa)*noDir_SigmaT;inv(eye(pv2)-kappa*eta')*kappa*noDir_SigmaT];   
        eta(1,:)=0;
        IRF_IV_noInd = [inv(eye(pv1)-eta'*kappa)*SigmaT;inv(eye(pv2)-kappa*eta')*kappa*SigmaT];        
        IRF_IV_Neither = [inv(eye(pv1)-eta'*kappa)*noDir_SigmaT;inv(eye(pv2)-kappa*eta')*kappa*noDir_SigmaT];        
    case -1
        pv1=IV;
        pz=size(MyCov,1);
        pv2=size(MyCov,2)-pz-pv1;
        CovV2_Z=MyCov(:,pv1+1:end-pz)';
        CovV1_Z=MyCov(:,1:pv1)';
        VarV1=Wc(1:pv1,1:pv1);
        VarV2=Wc(pv1+1:end,pv1+1:end);
        CovV1_V2=Wc(1:pv1,pv1+1:end);
        VarZ=MyCov(:,end-(pz-1):end);
        kappa=CovV2_Z * inv(VarZ) * CovV1_Z' * inv( CovV1_Z * inv(VarZ) * CovV1_Z') ;
        BB=VarV2-CovV1_V2'*kappa';
        v2rTv2r = VarV2 - (CovV1_V2'*kappa' + kappa * CovV1_V2) + kappa*VarV1*kappa';
        eta= inv(BB * inv(v2rTv2r ) * BB') * BB * inv(v2rTv2r) * (CovV1_V2'-kappa*VarV1) ; 
        
        DD = eta' * inv(eye(pv2)-kappa*eta') * kappa;
        vS1S1T = (eye(pv1^2) + kron(eye(pv1), DD) + kron(DD, eye(pv1)) ) \ reshape(VarV1 - eta' * VarV2 * eta, pv1^2,1);
        SigmaTSigmaTp = reshape(vS1S1T, pv1, pv1) ;

        SigmaT=chol(SigmaTSigmaTp,'lower');

        IRF_IV= [inv(eye(pv1)-eta'*kappa)*SigmaT;inv(eye(pv2)-kappa*eta')*kappa*SigmaT];           
        noDir_SigmaT = SigmaT;
        noDir_SigmaT(1,2) = 0;
        IRF_IV_noDir= [inv(eye(pv1)-eta'*kappa)*noDir_SigmaT;inv(eye(pv2)-kappa*eta')*kappa*noDir_SigmaT];   
        eta(1,:)=0;
        IRF_IV_noInd = [inv(eye(pv1)-eta'*kappa)*SigmaT;inv(eye(pv2)-kappa*eta')*kappa*SigmaT];        
        IRF_IV_Neither = [inv(eye(pv1)-eta'*kappa)*noDir_SigmaT;inv(eye(pv2)-kappa*eta')*kappa*noDir_SigmaT];        
    case 99 
        % is this where we take the first variable as independent and use
        % the second one as the residual
       [ind_b11iSig,ind_b21iSig,ind_SigmaTSigmaTp ]=my_MR_identification(MyCov(end,:),Wc);
       ind_SigmaT=chol(ind_SigmaTSigmaTp,'lower');

%        ind_IRF_IV= [ind_b11iSig;ind_b21iSig]*ind_SigmaT;
    
       [b11iSig,b21iSig,SigmaTSigmaTp ]=my_MR_identification(MyCov,Wc);
       q2=(b11iSig)\([ind_b11iSig; ind_b21iSig(1)]*ind_SigmaT);
       [aa,bb]=eig(SigmaTSigmaTp-q2*q2');
       bb=diag(bb);
       q1=sqrt(max(bb))*aa(:,bb==max(bb));
%        SigmaT=chol(SigmaTSigmaTp,'lower');

       IRF_IV= [b11iSig;b21iSig]*[q1 q2];
       noDir_q2 = q2;
       noDir_q2(1) = 0;
       IRF_IV_noDir= [b11iSig;b21iSig]*[q1 noDir_q2];
       [noInd_b11iSig,noInd_b21iSig, ~ ]=my_MR_identification(MyCov,Wc, true);
       IRF_IV_noInd = [noInd_b11iSig;noInd_b21iSig]*[q1 q2];
       IRF_IV_Neither = [noInd_b11iSig;noInd_b21iSig]*[q1 noDir_q2];
    case -99
       [ind_b11iSig,ind_b21iSig,ind_SigmaTSigmaTp ]=my_MR_identification(MyCov(1,:),Wc);
       ind_SigmaT=chol(ind_SigmaTSigmaTp,'lower');

%        ind_IRF_IV= [ind_b11iSig;ind_b21iSig]*ind_SigmaT;
    
       [b11iSig,b21iSig,SigmaTSigmaTp ]=my_MR_identification(MyCov,Wc);
       q2=(b11iSig)\([ind_b11iSig; ind_b21iSig(1)]*ind_SigmaT);
       [aa,bb]=eig(SigmaTSigmaTp-q2*q2');
       bb=diag(bb);
       q1=sqrt(max(bb))*aa(:,bb==max(bb));
%        SigmaT=chol(SigmaTSigmaTp,'lower');

       IRF_IV= [b11iSig;b21iSig]*[q2 q1];
       noDir_q1 = q1;
       noDir_q1(1) = 0;
       IRF_IV_noDir= [b11iSig;b21iSig]*[q2 noDir_q1];
       [noInd_b11iSig,noInd_b21iSig, ~ ]=my_MR_identification(MyCov,Wc, true);
       IRF_IV_noInd = [noInd_b11iSig;noInd_b21iSig]*[q2 q1];
       IRF_IV_Neither = [noInd_b11iSig;noInd_b21iSig]*[q2 noDir_q1];
    otherwise
        [a11,kappa_p]=my_IV_Identification(Wc,MyCov,size(MyCov,1)); %adjustment for obs/degrees of freedom in covariance doesn't matter
        IRF_IV_comp=[eye(size(MyCov,1)); kappa_p']*a11;
        IRF_IV= IRF_IV_comp;
end
