function VAR = doProxySVARx_sign(VAR)
%code from Piffer and Podowski (Economic Journal)
maxtrys=100; 
reps=100;
tmp=packr([VAR.m1 VAR.res]);
T = size(tmp,1);
VAR.n=size(VAR.res,2);
VAR.k=size(VAR.m1,2);
proxy=tmp(:,1:VAR.k);
innovations=tmp(:,VAR.k+1:end);
sigma = cov(innovations);%(VAR.res'*VAR.res)/(VAR.T-VAR.n*VAR.p-1);



% Identify B and A structural forms of the model
    E_1_loop = innovations'*proxy(:,1)/T;
    E_2_loop = innovations'*proxy(:,2)/T;
    
    
    E   = [E_1_loop, E_2_loop];
    E11 = E(1:2,:);
    E21 = E(3:end,:);

    Mu = E21*E11^(-1);
    Sigma11 = sigma(1:2,1:2);
    Sigma12 = sigma(1:2,3:end);
    Sigma21 = sigma(3:end,1:2);
    Sigma22 = sigma(3:end,3:end);


    Gamma = Sigma22 + Mu*Sigma11*Mu' - Sigma21*Mu' - Mu*Sigma21';
    B12B12prime = (Sigma21-Mu*Sigma11)'*Gamma^(-1)*(Sigma21-Mu*Sigma11);
    B11B11prime = Sigma11 - B12B12prime;
    B11_chol = cholx(B11B11prime)';
    B1_save=nan(reps,VAR.n,2);
    parfor j=1:reps
    trys=0;
    chck=-1;
    while chck<0 && trys<maxtrys
     step = randn(2,2);
    [Q, ~] = qr(step);  
    assert(max(max(Q*Q' - eye(2))) < 0.001)   
    assert(max(max(Q'*Q - eye(2))) < 0.001)   
   
    % B model 
    B11_cand = B11_chol*Q;
    B21_cand = Mu*B11_cand;
    B_firsttwo_cand = [B11_cand; B21_cand];
     step = sign( [B_firsttwo_cand(1,1), B_firsttwo_cand(1,2)]); % ensures that the first shock of positive sign increases the VXO and the second shock of positive sign decreases S&P500
    B11_cand = B11_cand*diag(step);
    B21_cand = B21_cand*diag(step);
    B_firsttwo_cand = [B11_cand; B21_cand];
    
    % A model  
    Lambda12  = Sigma12 - B11_cand*B21_cand';
    Lambda22  = Sigma22 - B21_cand*B21_cand';
    B12B22inv = Lambda12*Lambda22^(-1);
    A_firsttwo_cand              = (B11_cand-B12B22inv*B21_cand)^(-1)*[eye(2), -B12B22inv];
    shocks_estim_loop      = A_firsttwo_cand*innovations'; 
    shocks=shocks_estim_loop';
    c11=corrcoef([shocks(:,1) proxy(:,1)]);c11=c11(2,1);
    c21=corrcoef([shocks(:,1) proxy(:,2)]);c21=c21(2,1);
    c22=corrcoef([shocks(:,2) proxy(:,2)]);c22=c22(2,1);
    c12=corrcoef([shocks(:,2) proxy(:,1)]);c12=c12(2,1);
    cc=[c11 c12;
        c21 c22];
    e1=abs(cc(1,1))>0;
    e2=abs(cc(2,2))>0;
    e3=(abs(cc(1,1))-abs(cc(2,1)))>0.1;
    e4=(abs(cc(2,2))-abs(cc(1,2)))>0.1;
    
    ee=e1+e2+e3+e4;
    if ee==4
        chck=1;
    else
        trys=trys+1;
    end
    end
    if chck>0
    B1=B_firsttwo_cand;
    else
    B1=nan(rows(B_firsttwo_cand),cols(B_firsttwo_cand));
    end
    B1_save(j,:,:)=B1;
    end
    B1_median=nanmedian(B1_save);
    B1_diff=(B1_save-B1_median).^2;
    %B1_diff(:,:,VAR.shock)
    B1_diff=sum(B1_diff(:,:),2);
    [~,id]=min(B1_diff);
    
    VAR.b1(:,:,1)=squeeze(B1_save(id,:,:));
    