function bmean=estimate_VAR_mean(bolsvec,Y,L,REPS,y,x)
%VAR estimates with bootstrap bias correction
N=cols(Y);
T=rows(Y);
bols=reshape(bolsvec,N*L,N);
Y_mean=mean(Y);
resid=y-x*reshape(bolsvec,N*L,N);
%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,N);
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=zeros(1,N*L);
    jj=1;
    for j=1:L
        Xb(:,jj:jj+N-1)=Yb(t-j,:)-Y_mean;
        jj=jj+N;
    end
    Yb(t,:)=Y_mean+Xb*bols+residb(t,:);
    end
  
    %OLS estimation on bootstrap data
    data=detrend(Yb,'constant');
    xx=prepare(data,L);
    yy=data(L+1:end,:);
    xx=xx(L+1:end,1:end-1);
    bols_b=xx\yy;
    bsave(r,:,:)=bols_b;
end
%bias
bmean=vec(squeeze(mean(bsave)));
