function [betamat,Emat]=gs_sur(X,Xcell,y,ycell,k1,k2,k3,k,M,n,r1,r2,mu0,...
    V0,betadraw,v0,S0,Edraw);

R = r1+r2;
betamat=zeros(k,r2); 
Emat=zeros(M*(M+1)/2,r2); %collects unique elements of E as a vector

thoucount=1000;
% start main loop
i=1;
for i=1:R
   
   %%%%%%%%%%%%%%%%%%%%%%
   % draw beta
   %%%%%%%%%%%%%%%%%%%%%
   
   int1=0; 
   int2=0;
   j=1;
   for j=1:n
       Xi=cell2mat(Xcell(j)); % note: round brackets needed for this
       yi=cell2mat(ycell(j));
       int1=int1+Xi'*inv(Edraw)*Xi;
       int2=int2+Xi'*inv(Edraw)*yi;
   end
   V1=inv(inv(V0)+int1);
   mu1=V1*(inv(V0)*mu0+int2);
   betadraw=mvnrnd(mu1,V1)';    
   if i>r1
       betamat(:,i-r1)= betadraw;
   end   
    
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % draw E
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 v1=v0+n;
 int=y-X*betadraw;%MxN by 1
 int=reshape(int,M,n)';%n x M
 S1=S0+(int'*int); 
 Edraw=iwishrnd(S1,v1);
  
 if i>r1
      Emat(:,i-r1)=var_vec(Edraw);%collect unique elements of E as a vector
  end

  
     
  if i== thoucount
      i
      thoucount=thoucount+1000;
  end
        
    
end


