function [arrImpResp,arrmatA0]= IRFs_oil_largeVAR_signs(vecBeta,VCV,KEEP,intN,p,Sign)


h_IRF               = Sign.h_IRF;
SIGNS               = Sign.SIGNS;
per_SIGNS           = Sign.per_SIGNS;
max_check_periods   = Sign.max_check_periods;
arrImpResp          = NaN(intN,intN,h_IRF,KEEP);
arrmatA0            = NaN(intN,intN,KEEP);
MaxRandQ            = 1;

% counters
nKEEP   = 0;        % number of IRFs to keep
nRandQ=1;
MaxPermutations     = 1;
keepTrackPermutationsIndex = zeros(1,intN);
while nKEEP<KEEP %&& nRandQ<=MaxRandQ
    pos     = SIGNS>-9999; %seleziona le variabili per le quali uno shock è identificato
    Nshocks = size(SIGNS,2);
    [P,D]   = eig(VCV);
    KK      = normrnd(0,1,intN,intN,MaxRandQ);
    K=KK(:,:,nRandQ);
    [Q,R]=qr(K);
    for i=1:intN;
        if R(i,i)<0
            Q(:,i)=-Q(:,i);
        end
    end
    aa0=P*D.^0.5*Q';
    nPerm=1;
        %disp(strcat('nRandQ=',num2str(nRandQ)));
    while nKEEP<KEEP && nPerm<=MaxPermutations
        % disp(strcat('nRandQ=',num2str(nRandQ),' and nPerm=',num2str(nPerm)))
        % draw permutations and make sure they do not repeat;
        AcceptPermutation=0;
        while AcceptPermutation==0
            PermutationsIndex = randperm(intN);
            if  ismember(PermutationsIndex,keepTrackPermutationsIndex,'rows')==1
                AcceptPermutation=0;
            else
                AcceptPermutation=1;
                keepTrackPermutationsIndex = [keepTrackPermutationsIndex;PermutationsIndex];
            end
        end
        
        A0=aa0(:,PermutationsIndex);
        %% compute IRFS
        % now simulate the path of the t-varying coefficients and compute
        % IRFs
        % Impact IRF is given by the impact matrix
        B0t             = A0(:,:);
        
        if max_check_periods>1
            IRF_cand =cumsum(get_choleskIRF(intN,p,B0t,vecBeta',max_check_periods),3);
        else
            IRF_cand(:,:,1)=B0t;
        end
        [Index,A0_cand]=CheckSignsA0_2(B0t,IRF_cand,SIGNS,per_SIGNS,Nshocks,pos,intN); % This is a subroutine checking the signs; it returns 1 if all the restrictions
        %         % are satisfied, 0 if at least one restriction is not satisfied
        %         % Now check dynamic sign restriction
        %         if Index==1
        %             %   OIL PRICES after an
        %             % oil supply or global demand shock 4 quarters later cannot be lower
        %             % than the baseline level before the shock (in cumulated sums, cannot be lower than 0)
        %             % compute IRFs with the right sign (some columns could have changed sign)
        %             IRF_oilsup= NaN(intN,intN,4);
        %             for jj=1:5
        %                 IRF_oilsup(:,:,jj)=w(:,:,jj)*A0_cand;
        %             end
        %             check_dynamic1 = cumsum(squeeze(IRF_oilsup(1,1,:)));
        %             check_dynamic2 = cumsum(squeeze(IRF_oilsup(1,2,:)));
        %             check_all       = check_dynamic2(end).*check_dynamic1(end); % if this is positive ok, else Index =0
        %             Index = Index-(check_all(end)<0);
        %         end
        %         % check elasticity bounds
        %         % supply price elasticity bound: 0.1, the upper end of the distribution
        %         % in Kilian and Murphy
        %         if Index==1
        %             %   DeltaQ/DeltaP after a demand shock lower than
        %             %   .1 (slope of the supply curve)
        %             check_bound = A0_cand(end,2)./A0_cand(1,2);
        %             Index       = Index-(check_bound>0.1);
        %         end
        %         % demand price elasticity bound: -0.8, the lower end of the distribution
        %         % in Kilian and Murphy, conditional on eta_supply being .1
        if Index==1
            %             %   DeltaQ/DeltaP after a supply shock higher than -0.8
            %             %   (slope of the demand curve)
            check_bound = A0_cand(end,1)./A0_cand(1,1);
            Index = Index-(check_bound<-.8);
        end
        %        % bounds on the elasticity of world GDP to TFP shocks
        %         if Index==1
        %             check_bound = A0_cand(3,2);
        %             Index = Index-(check_bound<0.002);
        %         end
        %         % are satisfied, 0 if at least one restriction is not satisfied
        %         zz=zz+1;
        %     end
        if Index==1
            IRF_right =get_choleskIRF(intN,p,A0_cand,vecBeta',h_IRF);
            %            IRF_accept(:,:,:)    = IRF_right;
            nKEEP   = nKEEP+1;
            % Store I-R and output for this draw
            arrImpResp(:, :, :, nKEEP) = IRF_right;
            arrmatA0(:, :, nKEEP) = A0_cand;
            if rem(nKEEP ,1)==0
                disp(['found ',num2str(nKEEP),' IRFs out of ',num2str(KEEP)]);
            end
            % when you accept get another draw
            KK      = normrnd(0,1,intN,intN,MaxRandQ);            
            K       = KK(:,:,nRandQ);
            [Q,R]=qr(K);
            for i=1:intN;
                if R(i,i)<0
                    Q(:,i)=-Q(:,i);
                end
            end
            aa0=P*D.^0.5*Q';
            nPerm=0;            
        end
        nPerm=nPerm+1;
    end
%    nRandQ=nRandQ+1;
end
if nKEEP<KEEP
    error('you do not have enough IRFs, try increasing MaxRandQ');
end