function [outm,epsilon]=svol0x(data,reps,burn,skip,varargin)
if isempty(varargin)
    tight=0;
else
    tight=varargin{1};
end
REPS=reps;
BURN=burn;
Sindex=BURN+1:skip:REPS;
fsize=length(Sindex);
T0=20;
Y=data;
T=rows(Y);
N=cols(Y);
%set priors for SVOL

sigma0=std(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;
if tight
B0=[0;1];
Sigma0=eye(2)*0.000001;
else
B0=[0.95;0];
Sigma0=0.1;
Sigma0(1,1)=0.1;
end

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

epsilon=Y; %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));
    yy=log(hlast);
    xx=[lag0(yy,1) ones(rows(yy),1)];
    yy=yy(2:end,:);
    xx=xx(2:end,:);
    bdraw=getARxx(yy,xx,B0,Sigma0,g);
    delta=bdraw(1);
    alpha=bdraw(2);
    [bdraw' g]
 
if igibbs>BURN
    if sum(Sindex==igibbs)>0
    out(jgibbs,:)=log(hlast(2:end,:));
    out1(jgibbs,:)=[alpha/(1-delta) delta sqrt(g)]; 
    jgibbs=jgibbs+1;
    end
end
end

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