function[beta2,error,roots,problem,M,V]=...
    GLSdrawnew(Y,X,sigma,sig1,beta0,P00,L,CHECK,maxdraws,EX,mid)
T=rows(Y);
N=cols(Y);
%%Step 2a Set up matrices for the Kalman Filter

ns=cols(beta0);

terma=0;
termb=0;
iR1=invpd(sigma);
iR2=invpd(sig1);
for i=1:T
   if mid(i)==1 %instrument missing
    Rinv=iR1;
   else %instrument available so there is covariance
    Rinv=iR2;
   end
   terma=terma+X(i,:)'*(Y(i,:))*Rinv;
   termb=termb+kron(Rinv,X(i,:)'*X(i,:));
end
iP00=(P00);
V=invpd(iP00+termb);
M=V*(vec(terma)+iP00*beta0');
cV=cholx(V);
chck=-1;
problem=0;
trys=1;
while chck<0 && trys<=maxdraws
    
wa=randn(1,ns);


beta2=M'+(wa*cV);   %draw for beta in period t from N(beta_tt,ptt)
error=Y-X*reshape(beta2,N*L+EX,N);  %var residuals
roots=stability(beta2',N,L,EX);

if CHECK
if sum(roots)==0
    chck=1;
else
    trys=trys+1;
end
else
 chck=1;
end
end
if CHECK
    if chck<0
        problem=1;
    end
end


