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

%%% Written by Jonathan Hambur

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

[t,n]=size(Y);
A_b=zeros(n,n+1,nboot);


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

A=[alpha_hat_var, Gamma_hat_var];
Z = [ones(t-1,1), Y(1:end-1,:)].';

residuals = Y(2:end,:) - (A*Z).';

residuals = detrend(residuals,'constant'); %recentre residuals

for jj=1:nboot
    
    r=randi(t-1,t-1,1);
    start_r=randi(t,1);
    
    temp.Y_b=zeros(t,n);
    temp.Y_b(1,:)=Y(start_r,:);
    
    for ii=2:t
        temp.Y_b(ii,:)=(A*[1,temp.Y_b(ii-1,:)].').'+residuals(r(ii-1),:);
    end
    
    [temp.gamma_b,~,~, temp.alpha_b, ~]= olsvarc(temp.Y_b,1);
    
    A_b(:,:,jj)=[temp.alpha_b,temp.gamma_b];
    clear temp
end

A_b_med=median(A_b,3);

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

    
    