*! version 1.0.0  09may2006
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