Statalist The Stata Listserver


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: RE: durbin watson


From   "Steichen, Thomas J." <SteichT@rjrt.com>
To   <statalist@hsphsun2.harvard.edu>
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
        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

-----------------------------------------
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
recipient, please notify the sender by replying to this message and
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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index