function cc = cc_boot(var1,var2,k,n,conf)

%cc_boot    Calculates cross-correlations with bootstrapped error bands
%   cc_boot(VAR1,VAR2,K,N,CONF) calculates correlations between VAR1_t and 
%     VAR2_t-k, k = -K,-K+1,...,K, for N bootstrapped samples, and outputs 
%     median values along with error band at significance CONF.
%
%   VAR1, VAR2 = Tx1 data vectors. 
%   K = number of leads and lags.
%   N = number of boostrap replications.
%   CONF = significance of reported error band.  Must be between 0 and 1.
%
%   CC = (2*K+1)-by-3 matrix of calculated statistics
%     Row i gives correlation between VAR1_t and VAR2_t+(K+1-i).
%     Col 1 gives lower boundary of error band.
%     Col 2 gives median values across the bootstrapped samples.
%     Col 3 gives upper boundary of error band.

if size(var1) ~= size(var2)
    error('Data vectors must be same size')
end
t = size(var1,1);
var1 = var1-mean(var1);
var2 = var2-mean(var2);
var1_mat = repmat(var1,1,n);
var2_mat = repmat(var2,1,n);

boot_repl = zeros(2*k+1,n);
for i = 0:k
     draw_ld = randi(t-i,t-i,n);
     draw_lg = randi(t-i,t-i,n);
     var1_ld = var1_mat(draw_ld);
     var2_ld = var2_mat(i+draw_ld);
     var1_lg = var1_mat(draw_lg+i);
     var2_lg = var2_mat(draw_lg);
     boot_repl(k+1-i,:) = diag(corr(var1_ld,var2_ld));
     boot_repl(k+1+i,:) = diag(corr(var1_lg,var2_lg));
end
boot_repl = sort(boot_repl,2);
low = round(n*(1-conf));   
high = round(n*conf);
cc = zeros(2*k+1,3);
cc(:,1) = boot_repl(:,low);
cc(:,2) = median(boot_repl,2);
cc(:,3) = boot_repl(:,high);

%                                               G. Ramey - UCSD - 4/11/2019
