function [b_unb,ee]=estimate_VARboot(Y,X,L,REPS,SCALE,printinfo)
%VAR estimates with bootstrap bias correction

N=cols(Y);
T=rows(Y);
%OLS
bols=X\Y;
resid=Y-X*bols;
%starting values
history=cell(T,1);
for t=1:L
    history{t}=Y(1:L,:);
end
for t=L+1:T
    tmp=Y(t-L+1:t,:);
    history{t}=tmp;
end

%bootstrap
bsave=zeros(REPS,N*L+1,N);
for r=1:REPS
  if printinfo
       disp(strcat('REPS=',num2str(r)));
  end

    %re-sample resids
    id=randsample(T,T+L,true);
    residb=resid(id,:);
    %random start
    ids=randsample(T,1,true);
    start=history{ids};
    
    Yb=zeros(T+L,N);
    Yb(1:L,:)=start;
    for t=L+1:T+L
    Xb=ones(1,N*L+1);
    jj=1;
    for j=1:L
        Xb(:,jj:jj+N-1)=Yb(t-j,:);
        jj=jj+N;
    end
    Yb(t,:)=Xb*bols+residb(t,:);
    end
    
    %OLS estimation on bootstrap data
    xx=prepare(Yb,L);
    yy=Yb(L+1:end,:);
    xx=xx(L+1:end,:);
    bols_b=xx\yy;
    bsave(r,:,:)=bols_b;
end

%bias
bias=squeeze(mean(bsave))-bols;
b_unb=bols-bias;
FF=comp1(vec(b_unb),N,L,1);
ee=max(abs(eig(FF)));
%make sure VAR is stable
if ee>1
trys=0;
% (e.g. Kilian's stationarity adjustment)
  bias0 = bias;
  delta = 1;
  while (ee>1 && trys<1000000)
    delta = delta - SCALE;%.001;
    bias1 = delta*bias0;
    b_unb = bols - bias1;
    FF=comp1(vec(b_unb),N,L,1);
   ee=max(abs(eig(FF)));
   if printinfo
       disp(strcat('trys=',num2str([trys ee])));
  end
   trys=trys+1;
  end
end 
