*! 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