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)<logacc;
            rhodraw = rhocan;
            rhocount = rhocount +1;
        end
       

        if i>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];