Hi Rick,

I'm doing a very basic model using SAS/IML for nonlinear optimization by maximizing the likelihood function. Below is my code. It turns out that the model does not converge. The model do converge if I run the it with proc mixed. I was trying to figure out why this happened but I couldn't find out the reason. It is appreciated if you could take a look of my code and let me know what's wrong with it. Thanks in advance!

proc iml;
use x;
read all into x;
use y;
read all into y;
use nsub;
read all into nsub;
use ntime;
read all into ntime;
use initial;
read all into initial;

nrun=0;
/***loglikelihood funtion***/
start loglik(parameters) global(x,y,nsub,ntime,nrun);
beta=parameters[1:2];
print beta;
rho=parameters[3];
print rho;
sigma2=parameters[4];
print sigma2;

ll=0;
ind=1;
do while (ind<=nsub);
xi=x[(ind-1)*ntime+1:ind*ntime,];
yi=y[(ind-1)*ntime+1:ind*ntime];
pi=constant('PI');

/*remove missing values (monotone missing)*/
if (yi[nrow(yi)]=.) then do; /*if yi is not complete*/
j=nrow(yi);
do while (yi[j-1]=.);
j=j-1;
end;
di=j;
yobs=yi[1:di-1];
end;
else do;
yobs=yi;
end;
ll1=0;
ni=nrow(yobs);
mui=xi[1:ni,]*beta;
help1=0:(ni-1);
help2=shape(help1,ni,ni);
h=help2-t(help2);
h=rho##abs(h); /*elementwise power*/
vi=sigma2*h;
ll1=-(ni/2)*log(2*pi)-0.5*log(det(vi))-0.5*t(yobs-mui)*inv(vi)*(yobs-mui);
/*calculate ll*/
ll=ll+ll1;
print ll;
ind=ind+1;
end;
print ll;
nrun=nrun+1;
file log;
put nrun;
return(ll);
finish loglik;

opt={1 11};
con={. . -1 0 ,
. . 1 . };
call nlpnrr(rc, xr, "loglik", initial, opt, con);
call nlpfdd(maxlik,grad,hessian,"loglik",est);
inf=-hessian;
covar=inv(inf);
var=vecdiag(covar);
stde=sqrt(var);
create result var {est stde};
append;
quit;