function b_unb=estimate_VARSA_C(data,L,REPS,BURN,REPSBOOT,SCALE)
%VAR estimates with SA bias correction
N=cols(data);

Y=data(L+1:end,:);
X=prepare(data,L);
X=X(L+1:end,1:end);
reps=REPS+BURN;
gamma_i = .5;


%OLS
bols=X\Y;
bolsvec=vec(bols);
NN=rows(bolsvec);
theta=zeros(reps,NN);
theta(1,:)=bolsvec';
%SA algorithm
for r=1:reps
    r
    thetanew=estimate_VAR_mean_c(theta(r,:)',Y,L,REPSBOOT,Y,X);
    d=bolsvec-thetanew;
    theta(r+1,:)=theta(r,:)+gamma_i*d';
end
b_unb=reshape(mean(theta(BURN+1:end,:)),N*L+1,N);

%bias
bias=bols-b_unb;
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<10000)
    delta = delta - SCALE;%.001;
    bias1 = delta*bias0;
    b_unb = bols - bias1;
    FF=comp1(vec(b_unb),N,L,1);
   ee=max(abs(eig(FF)));
   trys=trys+1;
  end
end 


