Rudy Fichtenbaum <[email protected]> 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 <[email protected]>, 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 intercept 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
[email protected] [email protected]
References
----------
Judge G. G. et al (1988) "Introduction to the Theory and Practice of
Econometrics, 2d" New York: John Wiley and Sons.
----------------------------- 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>
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
