/* bgl.g
Procedure that implements the Bera Jarque Lee test on normality in a probit
model, see 'Testing the Normality Assumption in Limited Dependent Variable
Models', A.K. Bera, C.M. Jarque and L-F. Lee, International Economic Review
(1984), vol 23, pp. 563--578. This procedure calculates the test statistic
for a test of normality of the disturbances against any other (than normal)
member of the Pearson family of distributions.

usage:  t=bjl(i,yhat,x);
input:  i:      n-vector with observed response variables (0's and 1's)
        yhat:   n-vector with predicted latent responses (b'x)
        x:      n x k matrix with regressors
output: t:      scalar, value of test statistic (equation below 4.15 on p.571)

*/

proc bjl(i,yhat,x);
 local k,n,g,imatrix,omega,j,yhat2,xx,pdfni,cdfni;

 /* input checking */

 if ((cols(i) ne 1) or (cols(yhat) ne 1));
  errorlog("input error bjl: i and yhat must have 1 column only");
  retp(-1);
 endif;
 if ((maxc(i) ne 1) or (minc(i) ne 0));
  errorlog("data error bjl: i must consist of 0's and 1's only");
  retp(-1);
 endif;
 if (det(x'x)==0);
  errorlog("data error bjl: x matrix must be of full column rank");
  retp(-1);
 endif;
 if ((rows(i) ne rows(yhat)) or (rows(i) ne rows(x)));
  errorlog("data error bjl: i, yhat, and x must have same number of rows");
  retp(-1);
 endif;



 /* start calculations */

 k=cols(x);
 n=rows(x);
 pdfni=pdfn(yhat);
 cdfni=cdfn(yhat);
 yhat2=yhat.^2;
 xx=x~(yhat2-1)~(yhat.*(3+yhat2));
 omega=pdfni./(cdfni.*(1-cdfni));
 g=(omega.*(yhat2~(yhat.*(3+yhat2))))'(i-cdfni)/sqrt(n);
 imatrix=(pdfni.*omega.*xx)'xx/n;
 if (det(imatrix)<=0);
  errorlog("error bjl: estimate of information matrix is not pd");
  retp(-1);
 endif;
 j=zeros(k+2,2);j[k+1,1]=1;j[k+2,2]=1;
 retp( g'*j'*inv(imatrix)*j*g );
endp;