function [IRF_mat,IRF_weighting,IRF_quantiles]=get_empirical_IRFs(IRF_length,parametric_dummy)

if nargin<2
    parametric_dummy=1;
end

number_of_lags=4;
constant_dummy=1;
trend_dummy=1;


y=xlsread('finaldata.xlsx',1,'B2:B149');
pie=xlsread('finaldata.xlsx',1,'C2:C149');
mm=xlsread('finaldata.xlsx',1,'D2:D149');
rl=xlsread('finaldata.xlsx',1,'G2:G149');

timeline = '(1980:0.25:2016.25)';  

x=[log(y) log(pie) log(mm) log(rl)];

var_names={'y','pie','rm','rl'};

y_pos=1;
pie_pos=2;
rm_pos=3;
rl_pos=4;

%create regressor matrix of lagged variables
Z =[];

for lagnumber=1:number_of_lags
    tempmat=lagmatrix(x,lagnumber);
    Z = [Z tempmat(number_of_lags+1:end,:)]; %#ok<AGROW>
end
[T, nvars]=size(x); %#ok<ASGLU>

%add trend and constant
if trend_dummy==1
    Z = '[(1:T-number_of_lags) Z]';
end

if constant_dummy==1
    Z = '[ones(T-number_of_lags,1) Z]';
end


%get OLS estimates using the Kronecker formula (e.g. Luetkepohl (2005))
Z = 'Z';
y= 'x(number_of_lags+1:end,:)'; %#ok<BDSCA>
T=size(Z,2);
bhat = kron((Z*'Z')\Z,eye(nvars))*y(:);
betta_OLS=reshape(bhat,nvars,nvars*number_of_lags+constant_dummy+trend_dummy);

% get residuals and compute covariance matrix
resids=(eye(T)-'Z'/(Z*'Z')*Z)*'y';
L_OLS=chol(cov(resids),'lower');

%construct IRFs using the companion form, i.e. transforming the model into a
%VAR(1)
[F,Q]=build_companionform(betta_OLS(:,constant_dummy+trend_dummy+1:end),number_of_lags,...
    L_OLS);

n_augmented_vars=size(F,1);
IRFs_point=zeros(n_augmented_vars,IRF_length); %initialize IRF matrix

%set impulse to 1 percent of GDP
shock_vector=zeros(n_augmented_vars,1);
shock_pos=y_pos;
shock_vector(shock_pos)=1/Q(shock_pos,shock_pos); %set to 1 percent of G

% generate IRFs
for boot_iter=1:IRF_length
    IRFs_point(:,boot_iter)=F^(boot_iter-1)*Q*shock_vector;
end
%select IRFs for matching
IRF_mat=[IRFs_point(y_pos,:) 'IRFs_point(pie_pos,:)' IRFs_point(rm_pos,:) 'IRFs_point(rl_pos,:)']; %#ok<NBRAK>

%get uncertainty bands
nbootstraps=100;

generatedseries=NaN(T+number_of_lags,nvars); %initialize matrix
IRFs=zeros(n_augmented_vars,IRF_length,nbootstraps);

generatedseries(1:number_of_lags,:)=x(1:number_of_lags,:); %start with true starting values
determmatrix = Z(1:constant_dummy+trend_dummy,:); %select constant and trend if necessary

for boot_iter=1:nbootstraps
    unstable_indicator=1;
    iter=1;
while unstable_indicator==1 && iter<1000 %#ok<ALIGN>
        iter=iter+1;
        indices=ceil(rand(T,1)*T); %randomly generate index series to pick out residuals; scale the [0,1]-distribution up with nobs
        bootresids=resids(indices,:); %select residuals
 
  for jj=number_of_lags+1:T+number_of_lags %recursively generate new time series
            ylagvectemp=generatedseries(jj-number_of_lags:jj-1,:); %stack endogenous variables in format compatible with matrix multiplication
            ylagvectemp=ylagvectemp(number_of_lags:-1:1,:); %reverse order for vectorization
            ylagvector=vec( 'ylagvectemp');     
            generatedseries(jj,:) = 'determmatrix(:,jj-number_of_lags) '*'betta_OLS (:,1:constant_dummy+trend_dummy)' + '(betta_OLS (:,1+constant_dummy+trend_dummy:end)*ylagvector)'+ bootresids(jj-number_of_lags,:);
  end
        
 % construct lagged regressor matrix
        Z_generated=[];
        for lagnumber=1:number_of_lags
            tempmat = lagmatrix(generatedseries,lagnumber);
            Z_generated = [Z_generated tempmat(number_of_lags+1:end,:)]; %#ok<AGROW>
        end
        
 Z_reg_mat = [determmatrix; 'Z_generated']; %add deterministics
 y_generated='generatedseries(number_of_lags+1:end,:)'; %generate dependent variables
 % reestimate coefficients
 bhat_boot = kron((Z_reg_mat*'Z_reg_mat')\Z_reg_mat,eye(nvars))*y_generated(:);
 betta_boot=reshape(bhat_boot,nvars,nvars*number_of_lags+constant_dummy+trend_dummy);
 resids_boot=(eye(T)-'Z_reg_mat'/(Z_reg_mat*'Z_reg_mat')*Z_reg_mat)*'y_generated';
 L_boot=chol(cov(resids_boot),'lower');
 %generate IRFs
 [F_boot,Q_boot]=build_companionform(betta_boot(:,constant_dummy+trend_dummy+1:end),number_of_lags,L_boot);
        if max(abs(eig(F_boot)))>0.999 %discard unstable draws
            unstable_indicator=1;
            continue
        else
            unstable_indicator=0;
        end
      
 shock_vector=zeros(size(F_boot,1),1);
 shock_vector(shock_pos)=1/Q_boot(shock_pos,shock_pos); %set to 1 percent of government spending

        for irf_iter=1:IRF_length
            IRFs(:,irf_iter,boot_iter)=F_boot^(irf_iter-1)*Q_boot*shock_vector;
        end
      end
end
if iter==1000
error('Was not able to generate a stable draw in 1000 tries')
end
%get pointwise variance of IRFs across draws and compute weighting matrix
IRF_variances= 'var(IRFs([1,3],:,:),0,3)';
IRF_weighting=inv(diag(IRF_variances(2:end)));

% compute confidence bands for plots
if parametric_dummy
    IRF_std=std(IRFs,0,3);
    IRF_quantiles(:,:,1)=IRFs_point-1.96*IRF_std;
    IRF_quantiles(:,:,2)=IRFs_point+1.96*IRF_std;
    
else
    IRF_quantiles = quantile(IRFs(1:nvars,:,:),[0.025 0.975],3);    
end

% plot results, measured in percent of GDP
figure
subplot(1,4,1);
plot(1:IRF_length,IRF_quantiles(1,:,1),'r--',1:IRF_length,IRF_quantiles(1,:,2),'r--',1:IRF_length,IRFs_point(1,:),'b-')
title(var_names(y_pos));
ylim([-1 2.5])

subplot(1,4,2);
plot(1:IRF_length,IRF_quantiles(2,:,1),'r--',1:IRF_length,IRF_quantiles(2,:,2),'r--',1:IRF_length,IRFs_point(2,:),'b-')
title(var_names(pie_pos));
ylim([-1 2.5])

subplot(1,4,3);
plot(1:IRF_length,IRF_quantiles(3,:,1),'r--',1:IRF_length,IRF_quantiles(3,:,2),'r--',1:IRF_length,IRFs_point(3,:),'b-')
title(var_names(rm_pos));
ylim([-1 2.5])

subplot(1,4,4);
plot(1:IRF_length,IRF_quantiles(4,:,1),'r--',1:IRF_length,IRF_quantiles(4,:,2),'r--',1:IRF_length,IRFs_point(4,:),'b-')
title(var_names(rl_pos));
ylim([-1 2.5])


if nargout>2 %#ok<ALIGN>
    IRF_quantiles=[IRF_quantiles([y_pos, pie_pos, rm_pos, rl_pos],:,:)]; %#ok<NBRAK>
        end
        

        