


Rudy Fichtenbaum <rudy.fichtenbaum@wright.edu> asked,

> Does Stata have a way of calculating the p value for a Durbin-Watson
> statistic? SAS does this and it is a lot easier for students because
> they don't have to rely on a Durbin-Watson table which can result in the
> test being inconclusive.

We at Stata are not fans of the original Durbin-Watson test because the
test's p value is known to be heavily dependent on the
normality-of-the-residuals assumption.  A far better test is Durbin's
alternative test, available in Stata by tying -estat durbina- after
estimation by -regress-.

As others on the list have noted, they have requested p values and we have
not responded.  Other than (1) doing homework or (2) reproducing the results
of others, we know of no good reason why anyone should want the p value that
assumes normality of the Durbin-Watson statistic.  If we are wrong on this,
we would like to hear about it.  We suspect some researchers think that
assuming normality is reasonable to obtain an exact statistic, but we would
prefer a nonexact statistic in such cases.

Our concerns about the Durbin-Watson not withstanding, we have written -dwe-
to provide the p value of the Durbin-Watson statistic after -regress-,
thanks largely to Tom Streichen <SteichT@rjrt.com>, who posted well known
FORTRAN code on the list to do a part of the problem.

-dwe- takes no arguments and no options, and is for use after -regress-:

	. regress y x

	. dwe

This is work in progress, but we feel reasonably confident that -dwe- will
produce numerically correct results if the sample is not too large, say _N<100
(we have tested that) and probably _N<500 (we are testing that now).

Putting aside numerical issues, you should only use the D-W test if (1) you
ave an interecept in the model, (2) all variables are strictly exogenous
which, among other things, means you must not include lagged values of the
dependent variable among the explanatory variables.

You can obtain -dwe- by typing 

	. net from http://www.stata.com
	. net cd users/ddrukker
	. net install dwe 

-dwe- replicates the results reported in Judge, et al., 1988, pages 399-400.


Background
----------

We admit our interest in this was more about Mata and less about the
Durbin-Watson statistic.

The FORTRAN code posted by Tom STreichen, the amended version of the Applied
Statistics Algorithm AS 153 (AS R52) by R. W.  Farebrowther, amended by Clint
Cummins of TSP fame, can be found in Attachment 1.

We translated the routine to Mata and tested it; see Attachment 2.
We suspect many Mata users will be interested in that step.

After having proved our translation was correct, we wrote an ado-file to
calculate the Durbin-Watson exact p value.  The FORTRAN routine requires the
eigenvalues of a matrix be supplied as input.  This was easy to do in Mata.

The FORTRAN code suggested using the eigenvalues of M*A, where 
M=I-X(X'X)^(-1)X'.  In fact, convention is to use the eigenvalues of
A-D*X*(X'X)^(-1)*(D*X)', where D is the matrix of values that performs the
first-difference operation, and A=D*D', and that is what -dwe- uses.

Routine -dwe- can be found below in Attachment 3.

If you -net install dwe-, you will obtain files dwe.ado and dwe.hlp.
The help file is nothing more than this Statalist posting.

-- David                   -- Bill
   ddrukker@stata.com         ddrukker@stata.com


----------------------------- Attachment 1 --- file gradsol.f --- BEGIN ---
      DOUBLE PRECISION FUNCTION GRADSOL(A, M, C, N)
C
C     TRANSLATION OF AMENDED VERSION OF APPLIED STATISTICS ALGORITHM
C     AS 153 (AS R52), VOL. 33, 363-366, 1984.
C     BY R.W. FAREBROTHER  (ORIGINALLY NAMED GRADSOL OR PAN)
C
C     GRADSOL EVALUATES THE PROBABILITY THAT A WEIGHTED SUM OF
C     SQUARED STANDARD NORMAL VARIATES DIVIDED BY X TIMES THE UNWEIGHTED
C     SUM IS LESS THAN A GIVEN CONSTANT, I.E. THAT
C     A1.U1**2 + A2.U2**2 + ... + AM.UM**2 <
C     X*(U1**2 + U2**2 + ... + UM**2) + C
C     WHERE THE U'S ARE STANDARD NORMAL VARIABLES.
C     FOR THE DURBIN-WATSON STATISTIC, X = DW, C = 0, AND
C     A ARE THE NON-ZERO EIGENVALUES OF THE "M*A" MATRIX.
C
C     THE ELEMENTS A(I) MUST BE ORDERED.  A(0) = X
C     N = THE NUMBER OF TERMS IN THE SERIES.   THIS DETERMINES THE
C     ACCURACY AND ALSO THE SPEED.   NORMALLY N SHOULD BE ABOUT 10-15.
C     --------------
C     ORIGINALLY FROM STATLIB.  REVISED 5/3/1996 BY CLINT CUMMINS:
C     1. DIMENSION A STARTING FROM 0  (FORTRAN 77)
C        IF THE USER DOES NOT INITIALIZE A(0) = X,
C        THERE WOULD BE UNPREDICTABLE RESULTS, SINCE A(0) IS ACCESSED
C        WHEN J2=0 FOR THE FINAL DO 60 LOOP.
C     2. USE X VARIABLE TO AGREE WITH PUBLISHED CODE
C     3. FIX BUG 2 LINES BELOW  DO 60 L2 = J2, NU, D
C        PROD = A(J2)  -->  PROD = A(L2)
C        (PRIOR TO THIS FIX, ONLY THE TESTS WITH M=3 WORKED CORRECTLY)
C     4. TRANSLATE TO UPPERCASE AND REMOVE TABS
C     TESTED SUCCESSFULLY ON THE FOLLOWING BENCHMARKS:
C     1. FAREBROTHER 1984 TABLE (X=0):
C        A           C   PROBABILITY
C        1,3,6       1   .0542
C        1,3,6       7   .4936
C        1,3,6      20   .8760
C        1,3,5,7,9   5   .0544
C        1,3,5,7,9  20   .4853
C        1,3,5,7,9  50   .9069
C        3,4,5,6,7   5   .0405
C        3,4,5,6,7  20   .4603
C        3,4,5,6,7  50   .9200
C     2. DURBIN-WATSON 1951/71 SPIRITS DATASET, FOR X=.2,.3,...,3.8, C=0
C        COMPARED WITH BETA APPROXIMATION  (M=66), A SORTED IN REVERSE ORDER
C     3. JUDGE, ET AL 2ND ED. P.399 DATASET, FOR X=.2,.3,...,3.8, C=0
C        COMPARED WITH BETA APPROXIMATION  (M=8), A SORTED IN EITHER ORDER
C
      INTEGER M, N
      DOUBLE PRECISION A(0:M), C, X
C
C     LOCAL VARIABLES
C
      INTEGER D, H, I, J1, J2, J3, J4, K, L1, L2, NU, N2
      DOUBLE PRECISION NUM, PIN, PROD, SGN, SUM, SUM1, U, V, Y
      DOUBLE PRECISION ZERO, ONE, HALF, TWO
      DATA ZERO/0.D0/, ONE/1.D0/, HALF/0.5D0/, TWO/2.D0/
C
C     SET NU = INDEX OF 1ST A(I) >= X.
C     ALLOW FOR THE A'S BEING IN REVERSE ORDER.
C
      IF (A(1) .GT. A(M)) THEN
        H = M
        K = -1
        I = 1
      ELSE
        H = 1
        K = 1
        I = M
      ENDIF
      X = A(0)
      DO 10 NU = H, I, K
        IF (A(NU) .GE. X) GO TO 20
   10 CONTINUE
C
C     IF ALL A'S ARE -VE AND C >= 0, THEN PROBABILITY = 1.
C
      IF (C .GE. ZERO) THEN
        GRADSOL = ONE
        RETURN
      ENDIF
C
C     SIMILARLY IF ALL THE A'S ARE +VE AND C <= 0, THEN PROBABILITY = 0.
C
   20 IF (NU .EQ. H .AND. C .LE. ZERO) THEN
        GRADSOL = ZERO
        RETURN
      ENDIF
C
      IF (K .EQ. 1) NU = NU - 1
      H = M - NU
      IF (C .EQ. ZERO) THEN
        Y = H - NU
      ELSE
        Y = C * (A(1) - A(M))
      ENDIF
C
      IF (Y .GE. ZERO) THEN
        D = 2
        H = NU
        K = -K
        J1 = 0
        J2 = 2
        J3 = 3
        J4 = 1
      ELSE
        D = -2
        NU = NU + 1
        J1 = M - 2
        J2 = M - 1
        J3 = M + 1
        J4 = M
      ENDIF
      PIN = TWO * DATAN(ONE) / N
      SUM = HALF * (K + 1)
      SGN = K / DBLE(N)
      N2 = N + N - 1
C
C       FIRST INTEGRALS
C
      DO 70 L1 = H-2*(H/2), 0, -1
        DO 60 L2 = J2, NU, D
          SUM1 = A(J4)
C         FIX BY CLINT CUMMINS 5/3/96
C          PROD = A(J2)
          PROD = A(L2)
          U = HALF * (SUM1 + PROD)
          V = HALF * (SUM1 - PROD)
          SUM1 = ZERO
          DO 50 I = 1, N2, 2
            Y = U - V * DCOS(DBLE(I)*PIN)
            NUM = Y - X
            PROD = DEXP(-C/NUM)
            DO 30 K = 1, J1
              PROD = PROD * NUM / (Y - A(K))
   30       CONTINUE
            DO 40 K = J3, M
              PROD = PROD * NUM / (Y - A(K))
   40       CONTINUE
            SUM1 = SUM1 + DSQRT(DABS(PROD))
   50       CONTINUE
          SGN = -SGN
          SUM = SUM + SGN * SUM1
          J1 = J1 + D
          J3 = J3 + D
          J4 = J4 + D
   60   CONTINUE
C
C       SECOND INTEGRAL.
C
        IF (D .EQ. 2) THEN
          J3 = J3 - 1
        ELSE
          J1 = J1 + 1
        ENDIF
        J2 = 0
        NU = 0
   70 CONTINUE
C
      GRADSOL = SUM
      RETURN
      END
------------------------------- Attachment 1 --- file gradsol.f --- END ---


---------------------------- Attachment 2 --- file gradsol.do --- BEGIN ---
!! header Tom passed was truncated dmd has full header to put into gradsol
cscript

mata:

mata set matastrict on

/*
C
C     TRANSLATION OF AMENDED VERSION OF APPLIED STATISTICS ALGORITHM
C     AS 153 (AS R52), VOL. 33, 363-366, 1984.
C     BY R.W. FAREBROTHER  (ORIGINALLY NAMED GRADSOL OR PAN)
C
C     GRADSOL EVALUATES THE PROBABILITY THAT A WEIGHTED SUM OF
C     SQUARED STANDARD NORMAL VARIATES DIVIDED BY X TIMES THE UNWEIGHTED
C     SUM IS LESS THAN A GIVEN CONSTANT, I.E. THAT
C     A1.U1**2 + A2.U2**2 + ... + AM.UM**2 <
C     X*(U1**2 + U2**2 + ... + UM**2) + C
C     WHERE THE U'S ARE STANDARD NORMAL VARIABLES.
C     FOR THE DURBIN-WATSON STATISTIC, X = DW, C = 0, AND
C     A ARE THE NON-ZERO EIGENVALUES OF THE "M*A" MATRIX.
C
C     THE ELEMENTS A(I) MUST BE ORDERED.  A(0) = X
C     N = THE NUMBER OF TERMS IN THE SERIES.   THIS DETERMINES THE
C     ACCURACY AND ALSO THE SPEED.   NORMALLY N SHOULD BE ABOUT 10-15.
C     --------------
C     ORIGINALLY FROM STATLIB.  REVISED 5/3/1996 BY CLINT CUMMINS:
C     1. DIMENSION A STARTING FROM 0  (FORTRAN 77)
C        IF THE USER DOES NOT INITIALIZE A(0) = X,
C        THERE WOULD BE UNPREDICTABLE RESULTS, SINCE A(0) IS ACCESSED
C        WHEN J2=0 FOR THE FINAL DO 60 LOOP.
C     2. USE X VARIABLE TO AGREE WITH PUBLISHED CODE
C     3. FIX BUG 2 LINES BELOW  DO 60 L2 = J2, NU, D
C        PROD = A(J2)  -->  PROD = A(L2)
C        (PRIOR TO THIS FIX, ONLY THE TESTS WITH M=3 WORKED CORRECTLY)
C     4. TRANSLATE TO UPPERCASE AND REMOVE TABS
C     TESTED SUCCESSFULLY ON THE FOLLOWING BENCHMARKS:
C     1. FAREBROTHER 1984 TABLE (X=0):
C        A           C   PROBABILITY
C        1,3,6       1   .0542
C        1,3,6       7   .4936
C        1,3,6      20   .8760
C        1,3,5,7,9   5   .0544
C        1,3,5,7,9  20   .4853
C        1,3,5,7,9  50   .9069
C        3,4,5,6,7   5   .0405
C        3,4,5,6,7  20   .4603
C        3,4,5,6,7  50   .9200
C     2. DURBIN-WATSON 1951/71 SPIRITS DATASET, FOR X=.2,.3,...,3.8, C=0
C        COMPARED WITH BETA APPROXIMATION  (M=66), A SORTED IN REVERSE ORDER
C     3. JUDGE, ET AL 2ND ED. P.399 DATASET, FOR X=.2,.3,...,3.8, C=0
C        COMPARED WITH BETA APPROXIMATION  (M=8), A SORTED IN EITHER ORDER
C
*/


real scalar gradsol(real vector A, 	/* A[0:m] -> A[0+1 : m+1] */
		    real scalar m, 	/* integer */
		    real scalar c,
		    real scalar n	/* integer */
		   )
{
	real scalar /*integer */ d, h, i, j1, j2, j3, j4, k, l1, l2, nu, n2
	real scalar /* double */ num, pin, prod, sgn, sum, sum1, u, v, y
	real scalar /* double */ x

/*
C     SET NU = INDEX OF 1ST A(I) >= X.
C     ALLOW FOR THE A'S BEING IN REVERSE ORDER.
*/
	if (A[1+1]>A[m+1]) {
		h = m 
		k = -1
		i = 1 
	}
	else {
		h = k = 1
		i = m 
	}

	x = A[0+1]
	if (k>0) {
		for(nu=h; nu<=i; nu++) {
			if (A[nu+1]>=x) goto L20
		}
	}
	else {
		for(nu=h; nu>=i; nu--) {
			if (A[nu+1]>=x) goto L20
		}
	}

/*
C     IF ALL A'S ARE -VE AND C >= 0, THEN PROBABILITY = 1.
*/
	if (c>=0) return(1)

/*
C     SIMILARLY IF ALL THE A'S ARE +VE AND C <= 0, THEN PROBABILITY = 0.
*/
L20: 	
	if (nu==h & c<=0) return(0) 

	if (k==1) nu-- 
	h = m - nu
	if (c==0) y = h-nu ; 
	else 	  y = c * (A[1+1] - A[m+1])

	if (y>=0) {
		d = 2 
		h = nu 
		k = -k 
		j1 = 0 
		j2 = 2 
		j3 = 3 
		j4 = 1
	}
	else {
		d = -2 
		nu++ 
		j1 = m-2 
		j2 = m-1 
		j3 = m+1
		j4 = m
	}

	pin = 2*atan(1) / n
	sum = (k+1)/2
	sgn = k / n			// DBLE(N)
	n2  = n + n - 1 

/*
C       FIRST INTEGRALS
*/
	for (l1=h-2*floor(h/2); l1>=0; l1--) {		/* L70 */
		for (l2=j2; (d>0) ? (l2<=nu) : (l2>=nu); l2=l2+d) {
			sum1 = A[j4+1]
			prod = A[l2+1]
			u = (sum1+prod)/2
			v = (sum1-prod)/2
			sum1 = 0 
			for (i=1; i<=n2; i=i+2) {	/* L50 */
				y = u - v * cos(i*pin)
				num = y - x 
				prod = exp(-c/num)
				for (k=1; k<=j1; k++) {
					prod = prod*num/(y-A[k+1])
				}
				for (k=j3; k<=m; k++) {
					prod = prod*num/(y-A[k+1])
				}
				sum1 = sum1 + sqrt(abs(prod))
			}
			sgn = -sgn 
			sum = sum + sgn*sum1
		
			j1 = j1 + d 
			j3 = j3 + d 
			j4 = j4 + d
		} 
	
/*
C       SECOND INTEGRAL.
*/
		if (d==2)  j3 = j3 - 1 
		else	   j1 = j1 + 1
		j2 = 0 
		nu = 0 
	}
	return(sum)
}


real scalar testgradsol(real vector A, real scalar c)
{
	return(gradsol(A, length(A)-1, c, 10))
}

assert(reldif(testgradsol((0,1,3,6),  1), .0542)<.0001)
assert(reldif(testgradsol((0,1,3,6),  7), .4936)<.0001)
assert(reldif(testgradsol((0,1,3,6), 20), .8760)<.0001)

assert(reldif(testgradsol((0,1,3,5,7,9),  5), .0544)<.0001)
assert(reldif(testgradsol((0,1,3,5,7,9), 20), .4853)<.0001)
assert(reldif(testgradsol((0,1,3,5,7,9), 50), .9069)<.0001)

assert(reldif(testgradsol((0,3,4,5,6,7),  5), .0405)<.0001)
assert(reldif(testgradsol((0,3,4,5,6,7), 20), .4603)<.0001)
assert(reldif(testgradsol((0,3,4,5,6,7), 50), .9200)<.0001)

end
---------------------------- Attachment 2 --- file gradsol.do --- END -----


---------------------------- Attachment 3 --- file dwe.ado ---- BEGIN -----
*! version 1.0.0  08may2006
program define dwe, rclass
	version 9.2 

	if "`e(cmd)'" != "regress" {
		di as err "{cmd:dwe} only works after {cmd:regress}"
		exit 498
	}

	local N = e(N)

	qui tsset

	tempvar e top bot touse cons
	tempname xpxi tops bots dw

	mat `xpxi' = e(V)
	mat `xpxi' = (1/(e(rmse)^2))*`xpxi'

	qui predict double `e' , res

	qui gen double `top' = (D.`e' )^2 
	qui gen `touse'      = `top' < .
	qui gen double `bot' = `e'^2

	qui sum `top' 
	scalar `tops'        = r(sum)

	qui sum `bot' 
	scalar `bots'        = r(sum) 

	scalar `dw'          = `tops'/`bots'

// now get p-value 	

	tempvar  dcons 
	tempname b p
	
	mat `b' = e(b)
	local bnames : colnames  `b'

	local bnames : subinstr local bnames "_cons" "", 		///
		all word count(local ncons)

	if `ncons' == 0 {
		di as err "{cmd:dwe} requires a constant in the model"
		exit 498
	}

	if `ncons' > 1 {
		di as err "error parsing variable names"
		exit 498
	}

	qui tsrevar `bnames'
	local bnames = r(varlist)
	foreach var of local bnames {
		tempvar d`var'
		local dvars `dvars' `d`var''
		qui gen double `d`var'' = D.`var' if `touse'
	}

	qui gen double `cons'  = 1
	tempvar d`cons'
	local dvars `dvars' `d`cons''
	qui gen double `d`cons'' = 0 if `touse'

	local bnames `bnames' `cons'

	mata: _dwe_getpvalue("`bnames'", "`dvars'", "`touse'", 	///
		`N', "`dw'", "`xpxi'", "`p'" )

	local len 52
	di
	di as txt "Durbin-Watson test with normal p value"
	di as txt "{hline `len'}
	di as txt "{col 8}dw{col 22}Prob < dw{col 41}Prob > dw"
	di as txt "{hline `len'}
	di as res _col(5) %7.6g `dw' _col(22) %7.4f `p' _col(41) %7.4f 1-`p'
	di as txt "{hline `len'}

	ret clear
	ret scalar dw   = `dw'
	ret scalar p_l  = `p'
	ret scalar p_u  = 1-`p'
	ret scalar  N    = `N'
end


mata:

mata set matastrict on


void _dwe_getpvalue(	string scalar vars,  string scalar dvars, 
			string scalar touse, real   scalar n,
			string scalar dw,    string scalar xpxi_s,
			string scalar p) 
{			
	
	string vector	vnames, dvnames
	real matrix 	A, xpxi
	real vector	v, pv
	real scalar	tol 

	tol = 1e-14
	
	vnames  = tokens(vars)
	dvnames = tokens(dvars)
	
	st_view(xmat=.,  . , vnames,  touse)
	st_view(dxmat=., . , dvnames, touse)

	A = _dwe_mkddp(n)

	xpxi = st_matrix(xpxi_s)

	MA = (A - dxmat*xpxi*dxmat')

	v =symeigenvalues(MA)
	
	tol = tol * v[1]
	for(i=1; i<=cols(v); i++) {
		if (v[i]<=tol)	{
			break 
		}
	}
	i--
	v = v[1..i]
	v = (st_numscalar(dw) , v)

	pv = gradsol(v, i, 0, 50)

	st_numscalar(p, pv)
}	

real matrix _dwe_mkddp(real scalar n) 
{
	real matrix A
	real scalar i, n1, n2
	
	n1 = n-1
	n2 = n-2
	A = J(n1, n1, 0)
	
	A[1,1]   =  2
	A[1,2]   = -1
	A[n1,n1] =  2
	A[n1,n2] = -1

	for(i=2; i<=n2; i++) {
		A[i,i-1] = -1
		A[i,i]   =  2
		A[i,i+1] = -1
	}
	return(A) 
}
	
/*
C
C     TRANSLATION OF AMENDED VERSION OF APPLIED STATISTICS ALGORITHM
C     AS 153 (AS R52), VOL. 33, 363-366, 1984.
C     BY R.W. FAREBROTHER  (ORIGINALLY NAMED GRADSOL OR PAN)
C
C     GRADSOL EVALUATES THE PROBABILITY THAT A WEIGHTED SUM OF
C     SQUARED STANDARD NORMAL VARIATES DIVIDED BY X TIMES THE UNWEIGHTED
C     SUM IS LESS THAN A GIVEN CONSTANT, I.E. THAT
C     A1.U1**2 + A2.U2**2 + ... + AM.UM**2 <
C     X*(U1**2 + U2**2 + ... + UM**2) + C
C     WHERE THE U'S ARE STANDARD NORMAL VARIABLES.
C     FOR THE DURBIN-WATSON STATISTIC, X = DW, C = 0, AND
C     A ARE THE NON-ZERO EIGENVALUES OF THE "M*A" MATRIX.
C
C     THE ELEMENTS A(I) MUST BE ORDERED.  A(0) = X
C     N = THE NUMBER OF TERMS IN THE SERIES.   THIS DETERMINES THE
C     ACCURACY AND ALSO THE SPEED.   NORMALLY N SHOULD BE ABOUT 10-15.
C     --------------
C     ORIGINALLY FROM STATLIB.  REVISED 5/3/1996 BY CLINT CUMMINS:
C     1. DIMENSION A STARTING FROM 0  (FORTRAN 77)
C        IF THE USER DOES NOT INITIALIZE A(0) = X,
C        THERE WOULD BE UNPREDICTABLE RESULTS, SINCE A(0) IS ACCESSED
C        WHEN J2=0 FOR THE FINAL DO 60 LOOP.
C     2. USE X VARIABLE TO AGREE WITH PUBLISHED CODE
C     3. FIX BUG 2 LINES BELOW  DO 60 L2 = J2, NU, D
C        PROD = A(J2)  -->  PROD = A(L2)
C        (PRIOR TO THIS FIX, ONLY THE TESTS WITH M=3 WORKED CORRECTLY)
C     4. TRANSLATE TO UPPERCASE AND REMOVE TABS
C     TESTED SUCCESSFULLY ON THE FOLLOWING BENCHMARKS:
C     1. FAREBROTHER 1984 TABLE (X=0):
C        A           C   PROBABILITY
C        1,3,6       1   .0542
C        1,3,6       7   .4936
C        1,3,6      20   .8760
C        1,3,5,7,9   5   .0544
C        1,3,5,7,9  20   .4853
C        1,3,5,7,9  50   .9069
C        3,4,5,6,7   5   .0405
C        3,4,5,6,7  20   .4603
C        3,4,5,6,7  50   .9200
C     2. DURBIN-WATSON 1951/71 SPIRITS DATASET, FOR X=.2,.3,...,3.8, C=0
C        COMPARED WITH BETA APPROXIMATION  (M=66), A SORTED IN REVERSE ORDER
C     3. JUDGE, ET AL 2ND ED. P.399 DATASET, FOR X=.2,.3,...,3.8, C=0
C        COMPARED WITH BETA APPROXIMATION  (M=8), A SORTED IN EITHER ORDER
C
*/



real scalar gradsol(real vector A, 	/* A[0:m] -> A[0+1 : m+1] */
		    real scalar m, 	/* integer */
		    real scalar c,
		    real scalar n	/* integer */
		   )
{
	real scalar /*integer */ d, h, i, j1, j2, j3, j4, k, l1, l2, nu, n2
	real scalar /* double */ num, pin, prod, sgn, sum, sum1, u, v, y
	real scalar /* double */ x


/*
C     SET NU = INDEX OF 1ST A(I) >= X.
C     ALLOW FOR THE A'S BEING IN REVERSE ORDER.
*/
	if (A[1+1]>A[m+1]) {
		h = m 
		k = -1
		i = 1 
	}
	else {
		h = k = 1
		i = m 
	}

	x = A[0+1]
	if (k>0) {
		for(nu=h; nu<=i; nu++) {
			if (A[nu+1]>=x) goto L20
		}
	}
	else {
		for(nu=h; nu>=i; nu--) {
			if (A[nu+1]>=x) goto L20
		}
	}

/*
C     IF ALL A'S ARE -VE AND C >= 0, THEN PROBABILITY = 1.
*/
	if (c>=0) return(1)

/*
C     SIMILARLY IF ALL THE A'S ARE +VE AND C <= 0, THEN PROBABILITY = 0.
*/
L20: 	
	if (nu==h & c<=0) return(0) 

	if (k==1) nu-- 
	h = m - nu
	if (c==0) y = h-nu ; 
	else 	  y = c * (A[1+1] - A[m+1])

	if (y>=0) {
		d = 2 
		h = nu 
		k = -k 
		j1 = 0 
		j2 = 2 
		j3 = 3 
		j4 = 1
	}
	else {
		d = -2 
		nu++ 
		j1 = m-2 
		j2 = m-1 
		j3 = m+1
		j4 = m
	}

	pin = 2*atan(1) / n
	sum = (k+1)/2
	sgn = k / n			// DBLE(N)
	n2  = n + n - 1 

/*
C       FIRST INTEGRALS
*/
	for (l1=h-2*floor(h/2); l1>=0; l1--) {		/* L70 */
		for (l2=j2; (d>0) ? (l2<=nu) : (l2>=nu); l2=l2+d) {
			sum1 = A[j4+1]
			prod = A[l2+1]
			u = (sum1+prod)/2
			v = (sum1-prod)/2
			sum1 = 0 
			for (i=1; i<=n2; i=i+2) {	/* L50 */
				y = u - v * cos(i*pin)
				num = y - x 
				prod = exp(-c/num)
				for (k=1; k<=j1; k++) {
					prod = prod*num/(y-A[k+1])
				}
				for (k=j3; k<=m; k++) {
					prod = prod*num/(y-A[k+1])
				}
				sum1 = sum1 + sqrt(abs(prod))
			}
			sgn = -sgn 
			sum = sum + sgn*sum1
		
			j1 = j1 + d 
			j3 = j3 + d 
			j4 = j4 + d
		} 
	
/*
C       SECOND INTEGRAL.
*/
		if (d==2)  j3 = j3 - 1 
		else	   j1 = j1 + 1
		j2 = 0 
		nu = 0 
	}
	return(sum)
}

end

---------------------------- Attachment 3 --- file dwe.ado ------ END -----

<end>
