

function [SED, SEDalt, low, high] = boot_second(A,U,q,p,T,V,h,hh,IRF,B,BIAS,X,RR,pp);

N = RR;
EE = U(1:q,:);
[bootstat11,bootsam11] = bootstrp(N,'corrcoef',EE(1,:)); 
[n,m] = size(bootsam11);
iter = 0;
m = RR;
for i = 1:m    
    for j = 1:(T-p)    
        ee(:,j) = EE(:,bootsam11(j,i));
    end
    
    Y = [];
    for j = 1:p
        Y = [X(j,:)';Y];
    end
       
    for j = 1:T-p
        ysim(:,j) = V(1:q,1) + A(1:q,:)*Y(1:p*q) + ee(:,j);
        Y = [ysim(:,j);Y];
    end
    
    z = ysim';    
    [AA,SIGMAA,UU,VV]=olsvarc(z,p); 
%     if max(abs(eig(AA)))<0.999
        iter = iter + 1;
        BB = chol(SIGMAA(1:q,1:q));
        BB = BB';

        AA = AA + BIAS;
        IRFF = IRFVAR(AA,BB,p,hh);

        for j = 1:hh
            for k = 1:q^2
                IRFFc(k,j,iter) = IRFF(k,j);          
            end
        end
        
        % Performs the max share calculation---oprion I
        V9boot(:,:,iter) = IRFFc((q+1:q:(q-1)*q+1),1:hh,iter)';
        WWWboot = zeros(q-1);
        for j = 1:h         
            WWWboot = WWWboot + (h+1-j)*V9boot(j,:,iter)'*V9boot(j,:,iter);
        end
        
        [Qboot,QQboot] = eig(WWWboot);       % Decompose WWW
        RQboot = diag(QQboot);           % Main diagonal has eigenvalues
        [maxvboot,maxiboot] = max(RQboot);   % Maximum is first eigenvalue
        xboot = [0 Qboot(:,maxiboot)'];      % Take the 1st eigenvector (PC estimate)

        %newshockboot = zeros(q,hh,iter);
        for k = 1:q       
            newshockboot(k,1:hh,iter) = IRFFc((0:q:(q-1)*q)+k,1:hh,iter)'*xboot';           
        end
        
        if newshockboot(1,h,iter)<0
            newshockboot = -1*newshockboot;   
            xboot= -xboot;
        end
%     end
end

for j = 1:hh       
    for k = 1:q
        SED(k,j) = std(newshockboot(k,j,:)); % 1 std. deviation of boot newshock: This is what Barsky Sims show in their paper          
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Performs the max share calculation---option II
for j = 1:hh
    for k = 1:q^2   
        SEDalt(k,j) = std(IRFFc(k,j,:));
    end
end

for j = 1:hh       
    for k = 1:q^2           
        low(k,j) = prctile(IRFFc(k,j,:),pp);
    end
end
% %
% for j = 1:hh
%     for k = 1:q^2
%         SED(k,j) = std(IRFFc(k,j,:));
%     end
% end
% 
% for j = 1:hh
%     for k = 1:q^2
%         VAR(k,j) = var(IRFFc(k,j,:));
%     end
% end
% 
for j = 1:hh
    for k = 1:q^2   
        high(k,j) = prctile(IRFFc(k,j,:),100-pp);   
    end
end

