function [betamat,sig2mat,rhomat,acceptmat]=...
gs_AR1_ic(X,y,k,T,r1,r2,...
mu0,V0,betadraw,...
nu0,tau0,sig2draw,...
murho,Vrho,tc,tv,rhodraw);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data preparation for initial stuff that
% (i) will not be updated in GS, or
% (ii) will eventually be updated, but is needed before that point
% in the sampler
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R = r1+r2;
betamat=zeros(k,r2);
sig2mat=zeros(1,r2);
rhomat=zeros(1,r2);
rhocount=0;
% start main loop
thoucount=100;
i=1;
for i=1:R
%update full Omega
hat=abs(repmat([-(T-1):T-1],T,1));
int1=rhodraw.^hat;
Om=full(spdiags(int1,-(T-1):T-1,T,T));
%%%%%%%%%%%%%%%%%
% draw beta
%%%%%%%%%%%%%%%
int1=X'*((sig2draw*Om)\X);
int2=X'*((sig2draw*Om)\y);
V1=inv(inv(V0)+int1);
mu1=V1*(V0\mu0+int2);
betadraw=mvnrnd(mu1,V1)';%k by 1
if i>r1
betamat(:,i-r1)=betadraw;
end
%%%%%%%%%%%%%%%%%%
% draw sig2
%%%%%%%%%%%%%%%%%%
int3=(y-X*betadraw)'*((Om)\(y-X*betadraw)); %careful-no sig2draw in here!
nu1=(T+2*nu0)/2;
tau1=(1/2)*(2*tau0+int3);
sig2draw=1/gamrnd(nu1,1/tau1);
if i>r1
sig2mat(:,i-r1)=sig2draw;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw rho - MH algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x0=atanh(rhodraw); %this will keep rho between [-1,1]
options=optimset('LargeScale','off','display','off','TolFun',0.0001,'TolX',0.0001,...
'GradObj','off', 'Hessian','off','DerivativeCheck','off');
[x,fval,exitflag,output,gradient,H]=fminunc('rho_mode',x0,options,y,X,...
sig2draw,betadraw,murho,Vrho,T);
tstd = tc*sqrt(inv(H));
tmean=tanh(x); %convert back to actual rho -tanh is the inverse function fo atanh
rhocan=trnd_klaus(tmean,tstd,tv,1,1);
rhoold=rhodraw;
if abs(rhocan)<1 %- we only consider valid draws of rho
% CANDIDATE stuff:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%build Omega
hat=abs(repmat([-(T-1):T-1],T,1));
int1=rhocan.^hat;
Om=full(spdiags(int1,-(T-1):T-1,T,T));
int=-(1/2)* log(det(Om))-(1/2)*((y-X*betadraw)'*((sig2draw*Om)\(y-X*betadraw)));
priorchunk=tnormpdf(murho,sqrt(Vrho),-1,1,rhocan);
logcanrho=log(priorchunk)+ int - log(tpdf_klaus(rhocan,tmean,tstd,tv));
% OLD stuff:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%build Omega
hat=abs(repmat([-(T-1):T-1],T,1));
int1=rhoold.^hat;
Om=full(spdiags(int1,-(T-1):T-1,T,T));
int=-(1/2)* log(det(Om))-(1/2)*((y-X*betadraw)'*((sig2draw*Om)\(y-X*betadraw)));
priorchunk=tnormpdf(murho,sqrt(Vrho),-1,1,rhoold);
logoldrho=log(priorchunk)+ int - log(tpdf_klaus(rhoold,tmean,tstd,tv));
logacc = logcanrho-logoldrho;
% acceptance rule
lr=log(rand);
if log(rand)r1
rhomat(:,i-r1)=rhodraw;
end
end %end of abs(rho)<1 condition
if i== thoucount
i
thoucount=thoucount+100;
end
end
acceptmat=[rhocount/R];