function [intervals]=bootstrapProxy_LP_rotation_b(st_f_ccvar,vardata,hor,c,cc,nboot,conf,p,n,Z,startsample,endsample,iv_index)

% Function to perform bootstrap method
% Step 1: estimate reduced-form VAR and obtain reduced-form residuals
[pi_hat,~,~,Y_initial,~,err]=VAR_LP(vardata,p,c,st_f_ccvar);
T=length(vardata)-size(Y_initial,1); % Data after losing lags
if nargin < 14 || isempty(iv_index)
    iv_index = (startsample+1):(size(err,1)-endsample);
end

if c==1
    cons=pi_hat(1,:)';
    PI=pi_hat(2:end,:)';
else, cons=0;
    PI=pi_hat';
end

Data_new=zeros(T+p,n,nboot);
pi_hat_c=zeros(size(pi_hat,1),size(pi_hat,2));
err_c=zeros(T,n,nboot);

for j=1:nboot
    
% STEP 2: generate a new sample with bootstrap

Y_new=zeros(T,n);
Y_in= reshape(flipud(Y_initial)',1,[]);                  % 1*(n*p) row vector of [y_0 ... y_-p+1]

  for i=1:T
    Y_new(i,:)= cons' + Y_in*PI' + err(randi(T),:);  
    Y_in=[Y_new(i,:) Y_in(1:n*(p-1))];
  end

Data_new(:,:,j)=[Y_initial ; Y_new];                        % Add the p initial lags to recreate the sample

% STEP 3 Correct for the biases:

[pi_hat_c(:,:,j),~,~,~,~,err_c(:,:,j)]=VAR(Data_new(:,:,j),p,c);

end

%% 

bias=mean(pi_hat_c,3)-pi_hat;
new_pi_hat=pi_hat-bias;

if c==1
    cons=new_pi_hat(1,:)';
    PI=new_pi_hat(2:end,:)';
else, cons=0;
    PI=new_pi_hat';
    
end

nboot_2 = 2000;

ExtraBigC=zeros(hor,n,nboot_2);

for j=1:nboot_2

% STEP 2: generate a new sample with bootstrap

% Wild bootstrap based on simple distribution (~Rademacher)
rr = 1-2*(rand(T,1)>0.5);
u=zeros(length(err),size(err,2));

for k=1:n
    u(:,k) = err(:,k).*rr;
end

	Znew = Z.*repmat(rr,1,size(Z,2));
	Znew(any(isnan(Znew),2),:)=[];

Y_new=zeros(T,n);
Y_in= reshape(flipud(Y_initial)',1,[]);                  % 1*(n*p) row vector of [y_0 ... y_-p+1]
  for i=1:T
    Y_new(i,:)= cons' + Y_in*PI' + u(i,:);  
    Y_in=[Y_new(i,:) Y_in(1:n*(p-1))];
  end
    Data_new=[Y_initial ; Y_new];                        % Add the p initial lags to recreate the sample
  
% STEP 3: Esimate VAR(p) with new sample and obtain IRF:

[pi_hat_c,~,~,~,~,err_c]=VAR(Data_new,p,c);


	eps_p=err_c(iv_index,1);
	eps_q=err_c(iv_index,2:end);
	
	if size(Znew,1) ~= numel(iv_index)
	    error('The bootstrapped proxy and residual samples are not aligned.');
	end

ols_instr=OLS(eps_p,Znew,cc);

if cc==1
	u_hat_p=[ones(size(Znew,1),1),Znew]*ols_instr;
else, u_hat_p=Znew*ols_instr;
end

sq_sp=(u_hat_p'*u_hat_p)\u_hat_p'*eps_q;

s=[1,sq_sp]';

if c==1
BigA_c=[pi_hat_c(2:end,:)'; eye(n*p-n) zeros(n*p-n,n)]; % BigA companion form, np+1 x np+1 matrix
else, BigA_c=[pi_hat_c'; eye(n*p-n) zeros(n*p-n,n)]; % BigA companion form, np x np matrix
end

C_c=zeros(n,n,hor);
H=zeros(n,1,hor);

for l=1:hor
    BigC_c=BigA_c^(l-1);
    C_c(:,:,l)=BigC_c(1:n,1:n); 
    H(:,:,l)=C_c(:,:,l)*s; % Impulse response functions of the IV Wold representation (Bootstrap)
end

C_c_1= reshape(permute(H,[3 2 1]),hor,n,[]);

ExtraBigC(:,:,j)=C_c_1;

national_single_irf(:,j) = ExtraBigC(:,end,j);

for jjj = 1:size(ExtraBigC,2)-1
difference_national_ea(:,jjj,j) = national_single_irf(:,j) - ExtraBigC(:,jjj,j);
end

end

ExtraBigC_cum=cumsum(ExtraBigC,1);
national_single_irf_cum = cumsum(national_single_irf,1);

for j = 1:nboot_2
    for jjj = 1:size(ExtraBigC,2)-1
        difference_national_ea_cum(:,jjj,j) = national_single_irf_cum(:,j) - ExtraBigC_cum(:,jjj,j);
    end
end

% Create bands:
LowH=zeros(hor,n);
HighH=zeros(hor,n);
LowH_cum=zeros(hor,n);
HighH_cum=zeros(hor,n);

for k=1:n

        Cmin = prctile(ExtraBigC(:,k,:),(100-conf)/2,3); %16th percentile
        LowH(:,k) = Cmin; %lower band
        
        Cmin_cum = prctile(ExtraBigC_cum(:,k,:),(100-conf)/2,3); %16th percentile
        LowH_cum(:,k) = Cmin_cum; %lower band
        
        Cmax = prctile(ExtraBigC(:,k,:),(100+conf)/2,3); %84th percentile
        HighH(:,k) = Cmax; %upper band
        
        Cmax_cum = prctile(ExtraBigC_cum(:,k,:),(100+conf)/2,3); %84th percentile
        HighH_cum(:,k) = Cmax_cum; %upper band

        % Difference with EA

        if k < n

        Cmin_diff = prctile(difference_national_ea(:,k,:),(100-conf)/2,3); %16th percentile
        LowH_diff(:,k) = Cmin_diff; %lower band
        
        Cmin_diff_cum = prctile(difference_national_ea_cum(:,k,:),(100-conf)/2,3); %16th percentile
        LowH_diff_cum(:,k) = Cmin_diff_cum; %lower band
        
        Cmax_diff = prctile(difference_national_ea(:,k,:),(100+conf)/2,3); %84th percentile
        HighH_diff(:,k) = Cmax_diff; %upper band
        
        Cmax_diff_cum = prctile(difference_national_ea_cum(:,k,:),(100+conf)/2,3); %84th percentile
        HighH_diff_cum(:,k) = Cmax_diff_cum; %upper band
        
        end

        intervals.LowH = LowH;
        intervals.LowH_cum = LowH_cum;
        intervals.HighH = HighH;
        intervals.HighH_cum = HighH_cum;

        intervals.LowH_diff = LowH_diff;
        intervals.LowH_diff_cum = LowH_diff_cum;
        intervals.HighH_diff = HighH_diff;
        intervals.HighH_diff_cum = HighH_diff_cum;
        
end

end

