function [outm,epsilon]=svolAR1(data,reps,burn,skip,varargin)

REPS=reps;
BURN=burn;
Sindex=BURN+1:skip:REPS;
fsize=length(Sindex);
T0=20;
Y=data;
X=[lagmatrix(Y,1) ones(rows(Y),1)];
X(1,:)=X(2,:);
T=rows(Y);
N=cols(Y);
%set priors for SVOL

sigma0=std(diff(Y(1:T0)))^2;

MU0=log(diag(sigma0)); %prior mean
SS0=10;

g0=0.1;  %prior scale parameter for inverse gamma
Tg0=1;     %prior degrees of freedom
g=g0;
alpha=0;
delta=1;



rho0=zeros(2,1);
SigmaR0=eye(2);
rho=X\Y;


hlast=(diff(Y).^2)+0.0001;   %rough guess for stochastic volatility 
hlast=[hlast(1:2,:);hlast];  %rough intial guess for svol

epsilon=Y-X*rho; %starting values for VAR residuals
%50 iterations to smooth initial estimates of stochastic volatility

for m=1:50
  hnew=zeros(T+1,N);
  for i=1:N
      htemp=getsvolxx(hlast(:,i),g0,MU0(i),SS0,epsilon(:,i),0,1);
      hnew(:,i)=htemp;
      
  end
hlast=hnew;

end
out=zeros(fsize,T,N);
jgibbs=1;
for igibbs=1:REPS
    
    igibbs
   
    
   
    %sample VOL
    htemp=getsvolxx(hlast,g,MU0,SS0,epsilon,alpha,delta); 
    hlast=htemp;
    
    gerrors=log(hlast)-alpha-lag0(log(hlast),1)*delta;
    g=IGx(Tg0,g0,gerrors(2:end));
    
  
    
    ystar=Y./sqrt(hlast(2:end));
    xstar=X./sqrt(hlast(2:end));
    rho=getARxx(ystar,xstar,rho0,SigmaR0,1);
    epsilon=Y-X*rho;
    
 
if igibbs>BURN
    if sum(Sindex==igibbs)>0
    out(jgibbs,:)=(hlast(2:end,:))./(1-rho(1)^2);
    jgibbs=jgibbs+1;
    end
end
end

outm=prctile(out,[50 16 84]);
