%%% This program estiamtes bootstrapped estiamtes of the VAR dynamics that
%%% adjust for small sample bias. VAR(1)

%%% Written by Jonathan Hambur
%%%% Based on code from Bauer 2014

function [Gamma_hat_boot, alpha_hat_boot, V_hat_boot]=Bootstrap(Y,nboot)

[Gamma_hat_var,~,~, ~, ~]= olsvarc(Y,1);

[T,k] = size(Y);
Z = [ones(T-1,1), Y(1:end-1,:)].';


Y_mean = mean(Y);
Y_mean_rep = Y_mean' * ones(1,nboot);

X_sim = zeros(k, T, nboot);

A_b=zeros(k,k,nboot);


% VAR(1)
% 1. obtain residuals
% use bootstrapped residuals
resid = zeros(k, T-1);
for t = 2:T
    resid(:,(t-1)) = Y(t,:)' - Y_mean' - (Gamma_hat_var * (Y(t-1,:) - Y_mean)' );
end
% 2. generate series
% randomly select initial values from the data Y
ind_start = randi(T,nboot,1);
X_sim(:,1,:) = Y( ind_start,:)';
for t = 2:T
    u_sim = resid(:, randi(T-1,nboot,1));
    X_sim(:,t,:) = Y_mean_rep + Gamma_hat_var * (squeeze(X_sim(:,t-1,:)) - Y_mean_rep) + u_sim;
end
    
for jj=1:nboot
    
    temp.Y_b=X_sim(:,:,jj)';
    temp.Y_b=temp.Y_b-ones(T,1)*mean(temp.Y_b); % deamn X so can focus on Gamma

    [temp.gamma_b,~,~, ~, ~]= olsvarc_noint(temp.Y_b,1);
    
    A_b(:,:,jj)=temp.gamma_b;
    clear temp
end

A_b_avg=mean(A_b,3);

inside=0;
delta=1;
while inside==0
    temp.A_boot=(1+delta).*Gamma_hat_var-delta.*A_b_avg;
    temp.eig=eig(temp.A_boot);
    if max(real(temp.eig))>=1
        delta=delta-0.01;
    else
        inside=1;
        Gamma_hat_boot=temp.A_boot;
        alpha_hat_boot=Y_mean'-Gamma_hat_boot*Y_mean';
        V_hat_boot=Y(2:end,:) - ([alpha_hat_boot Gamma_hat_boot]*Z).';
    end
    clear temp
end

    
    