st: RE: durbin watson

 From "Steichen, Thomas J." To Subject st: RE: durbin watson Date Fri, 5 May 2006 10:49:42 -0400

```Rudy Fichtenbaum writes:
>
> 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.

FORTRAN algorithm AS 153 is available for the p-value (copied below). Are you (or anyone else) brave enough to translate it to Mata
code?  It is known to be slow, even in compiled code!

Tom

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
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
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
RETURN
END

-----------------------------------------
CONFIDENTIALITY NOTE: This e-mail message, including any
attachment(s), contains information that may be confidential,
protected by the attorney-client or other legal privileges, and/or
proprietary non-public information. If you are not an intended
recipient of this message or an authorized assistant to an intended
then delete it from your system. Use, dissemination, distribution,
or reproduction of this message and/or any of its attachments (if
any) by unintended recipients is not authorized and may be
unlawful.

*
*   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/
```