% IDENTIFY.M

function [B0inv]=identify(A,P)

% Number of variables
n=size(P,1);

% Horizon 0: Structural impact multiplier matrix
L0=P;  

L=L0;

% Impact exclusion restriction
E1=[1 0 0;
    0 1 0];    % Restrictions for MSC expectations shock

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start ALGORITHM
% First iteration, j=1
Q0=zeros(n);
R1=[E1*L];
N0=null(R1)';
x1=randn(n,1); 
q1=N0'*(N0*x1./norm(N0*x1));

% Second iteration, j=2 
R2=[q1'];
N1=null(R2)';
x2=randn(n,1);
q2=N1'*(N1*x2./norm(N1*x2));

% Third iteration, j=3
R3=[q1';q2'];
N2=null(R3)';
x3=randn(n,1);
q3=N2'*(N2*x3./norm(N2*x3));


% Assemble results and check that Q*Q'=I
Q=[q1 q2 q3];
Q*Q';              
            

% Compute posterior draw for structural impact multiplier
B0inv=L0*Q;

% Flip signs of columns as needed
if B0inv(3,1)<0
    B0inv(:,1)=-B0inv(:,1);
end;

if B0inv(1,2)<0 & B0inv(2,2)<0 & B0inv(3,2)<0
    B0inv(:,2)=-B0inv(:,2);
end;    

if B0inv(1,3)>0 & B0inv(2,3)<0 & B0inv(3,3)<0
    B0inv(:,3)=-B0inv(:,3);
end;

  
