function Ysave=getbootsamples(bols,Y,X,L,REPS)
%VAR estimates with bootstrap bias correction

N=cols(Y);
T=rows(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
Ysave=zeros(T+L,N,REPS);
for r=1:REPS
    
    %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
    
    
    Ysave(:,:,r)=Yb;
end
