/* ........... PROCEDURES FOR PDF ............................

                Anurag N Banerjee
                Center For Economic Research
                Tilburg University
                Netherlands
                                        Date 4/11/97
..............................................................
This program is in the public domain.  While the author disclaims
any responsibility for the performance of this software, he
would appreciate receiving any comments.

This written by Anurag Banerjee and may be distributed as freeware
for public non-commercial use. Please provide appropriate
acknowledgment if this supports supports published work.
..................................................................
                  IMHOF1(A,B,x)

  This  proccedure for computing the distribution
  function

          P(u'Au/u'Bu < x)

  where u dist as N(0,I) and A sym matrix and B is a
  psd matrix.

                   IMHOFINV(A,B,X1,PROB);

  This  proccedure for computing the inverse of the
  distribution function

          P(u'Au/u'Bu < x) = prob

  where u dist as N(0,I) and A sym matrix and B is a
  psd matrix.
             X1 is a 2X1 vector, having the search
  ranges of x.
...............................................................
References:

1) Abrahamse, A.P.J and Koerts,J. (1969)  On the Theory and
application of the General Linear Model , Rotterdam University Press.

2) Imhof, J.P (1961), '' Computing the Distribution of Quadratic
forms of Normal Variables'', Biometrika 48,419-426.
.................................................................
Comment: These proccedues has been used in various journal articles and
Anurag N Banerjee's PhD thesis as well

Useful Applications include calculation of significance levels of DW
type tests and their power.
....................................................................*/

DECLARE MATRIX _LAMDA;
DECLARE MATRIX _TOL;


PROC IMHOFINV(A,B,X1,PROB);
LOCAL ICOUNT,EPS,B1,FB1,B2,FB2,WT,Z,PROB1;
B1 = MINC(X1);
B2 = MAXC(X1);
PROB1 = 0;
FB1 = IMHOF1(A,B,B1);
FB2 = IMHOF1(A,B,B2);
EPS = 1E-8;
ICOUNT = 1;
DO WHILE (ABS(PROB1-PROB) GT EPS);
   WT = (FB2-PROB)/(FB2-FB1);
   Z = B1*WT+B2*(1-WT);
   PROB1 = IMHOF1(A,B,Z);
   IF (PROB1 > PROB);
   B2 = Z;
   FB2 = PROB1;
   ELSE;
   B1 = Z;
   FB1 = PROB1;
   ENDIF;
ICOUNT = ICOUNT +1;
IF (ICOUNT > 1000 );
PRINT "FAILED TO CONVERGE ITERATION MORE THAN 1000 ";
RETP(Z);
ENDIF;

ENDO;
PRINT ICOUNT;

RETP(Z);
ENDP;

PROC IMHOF1(A,B,X);
LOCAL TOL,Z,X1,EV;

EV = EIGH(A-X*B);

TOL = 1E-20;
_lamda = EV;
_TOL   = TOL;
_INTORD = 40;
X1 = {1,
      0};
Z = INTQUAD1(&INTEGRAD0,X1);
Z = 0.5 - Z/PI;
RETP(Z);
ENDP;


PROC INTEGRAD0(U);
RETP(INTEGRAD1(U,_LAMDA,_TOL));
ENDP;

PROC INTEGRAD1(U,LAMDA,TOL);
LOCAL IU,UBAR,TETA,RHO,FF0,FF1,FF2,FF,EVALUE,N,I,SUM1,SUM2,EPS;

EPS = TOL/100;
EVALUE = LAMDA;
N = ROWS(EVALUE);
SUM1  = 0;
SUM2  = 0 ;
FF0  = 0.5*SUMC(EVALUE);

IU = (( U .> TOL) .AND ( U .< (1- TOL))) ;
UBAR = U.*IU + TOL.*(1-IU);

I =1;
DO WHILE I <=N;
 SUM1 = SUM1 + ATAN(EVALUE[I,1]*(1-UBAR)./UBAR);
 SUM2 = SUM2 + LN(UBAR^2 + (EVALUE[I,1].*(1-UBAR))^2);
I = I+1;
ENDO;
  TETA = 0.5 .* SUM1;
  RHO  = (N/2-1).*LN(UBAR)-LN((1-UBAR)) - (0.25 .* SUM2);

  FF1 = SIN(TETA).*EXP(RHO) ;
  FF2 =TOL;
  FF = FF0.*( U .> 1- tol) + IU.*FF1 + FF2 .* (U .< TOL) ;

RETP(FF);
ENDP;