C COMPUTATION OF TEST STATISTICS FOR LINEAR REGRESSION C WITH AR(1) ERRORS USING REAL DATA C ESTIMATION UNDER H1 USING THE BEACH-MCKINNON METHOD C DATA FROM WESOLOWSKY (1976, P. 142) C Y: QUANTITY OF BEEF CONSUMED C X1: PRICE OF BEEF C X2: PER CAPITA DISPOSABLE PERSONAL INCOME C C INPUT IS IN BEEF18.INP C FOR N=16 CASE, INPUT IS IN BEEF16.INP implicit real*8 (a-h,o-z) parameter (n=18,np=3) C N: SAMPLE SIZE C NP: COLUMN DIMENSION OF REGRESSOR MATRIX X, INCLUDING THE FIRST C COLUMN OF ONES real*8 x(n,np),y(n,1) real*8 beta0(np,1),e0(n,1),amat(n,n) real*8 beta1(np,1) data zero/0.0d0/,one/1.0d0/ C INPUT DATA FROM FILE do i=1,n read(5,*)y(i,1),(x(i,j),j=1,np) enddo C INPUT ESTIMATES UNDER HO: beta0(j),j=1 to np and sigma0 do j=1,np read(5,*)beta0(j,1) enddo read(5,*)sigma0 C INPUT ESTIMATES UNDER H1: beta1(j),j=1 to np, sigma1 and rho1 do j=1,np read(5,*)beta1(j,1) enddo read(5,*)sigma1 read(5,*)rho1 C EVALUATION OF TEST STATISTICS C ASYMPTOTIC TEST, AST ast=float(n)*(rho1**2) C Likelihood Ratio Test, LRT alrt=float(n)*dlog(sigma0/sigma1)+dlog(1.0-rho1**2) C SCORE TEST, ST do i=1,n sum=zero do j=1,np sum=sum+beta0(j,1)*x(i,j) enddo e0(i,1)=y(i,1)-sum enddo sumst=zero do i=2,n sumst=sumst+e0(i,1)*e0(i-1,1) enddo sct=(sumst**2)/(float(n-1)*sigma0**2) C DURBIN WATSON STATISTIC DW rho=0.0 call durbinw(rho,y,x,beta0,e0,amat,dw) C MODIFIED DW TEST (JUDGE ET AL, P. 399) 50 call dwext(x,amat,dstar,ddstar) C MODIFIED SCORE STATISTIC call modscor(x,y,sigma0,beta0,sumst,sct,pscor) C MODIFIED PROFILE LIKELIHOOD RATIO TEST call modprof(x,sigma0,sigma1,rho1,alrt,plrt) write(6,*)'Test Statistics' write(6,*)'DW,MDW,AST,LRT,MLRT,ST,MST' write(6,30)dw,dstar,ast,alrt,plrt,sct,pscor 30 format(1x,7f10.4) stop end C C C subroutine durbinw(rho,y,x,beta0,e0,amat,dw) parameter(n=18,np=3) implicit real*8 (a-h,o-z) real*8 y(n,1),x(n,np),beta0(np,1),amat(n,n),e0(n,1) real*8 e0t(1,n),den(1,1),anum(1,1),t1(n,1) data zero/0.0d0/,one/1.0d0/,two/2.0d0/,onem/-1.0d0/ data four/4.0d0/ C FORM THE AMAT MATRIX IN D=E0'AE0/E0'E0 do 4 i=1,n do 4 j=1,n amat(i,j)=zero 4 continue amat(1,1)=one do 2 i=2,n-1 amat(i,i)=two 2 continue amat(n,n)=one do 3 i=1,n-1 j=i+1 amat(i,j)=onem amat(j,i)=onem 3 continue call trans(e0,e0t,n,1) call matmul(e0t,e0,den,1,n,1) call matmul(amat,e0,t1,n,n,1) call matmul(e0t,t1,anum,1,n,1) dw=anum(1,1)/den(1,1) return end C C C subroutine dwext(x,amat,dstar,ddstar) parameter(n=18,np=3) implicit real*8 (a-h,o-z) real*8 x(n,np),xt(np,n),temp(np,np),tempinv(np,np),amat(n,n) real*8 z1(n,np),z2(np,np),z3(np,np),z4(n,n),z5(n,np),z6(np,np) real*8 z7(np,np),z3sq(np,np) data zero/0.0d0/,one/1.0d0/,two/2.0d0/,three/3.0d0/ data four/4.0d0/ data dwl18/1.046/,dwu18/1.535/,edu18/2.257/,vdu18/0.18510/ data dwl16/0.982/,dwu16/1.539/,edu16/2.293/,vdu16/0.20125/ call trans(x,xt,n,np) call matmul(xt,x,temp,np,n,np) call dlinds(np,temp,np,tempinv,np) call matmul(amat,x,z1,n,n,np) call matmul(xt,z1,z2,np,n,np) call matmul(z2,tempinv,z3,np,np,np) trace=zero do 75 i=1,np trace=trace+z3(i,i) 75 continue qty1=two*(float(n)-one)-trace call matmul(amat,amat,z4,n,n,n) call matmul(z4,x,z5,n,n,np) call matmul(xt,z5,z6,np,n,np) call matmul(z6,tempinv,z7,np,np,np) trace2=zero do 76 i=1,np trace2=trace2+z7(i,i) 76 continue call matmul(z3,z3,z3sq,np,np,np) trace3=zero do 77 i=1,np trace3=trace3+z3sq(i,i) 77 continue qty2=two*(three*float(n)-four)-two*trace2+trace3 if (n .eq. 18) edu=edu18 if (n .eq. 16) edu=edu16 if (n .eq. 18) vdu=vdu18 if (n .eq. 16) vdu=vdu16 if (n .eq. 18) dwu=dwu18 if (n .eq. 16) dwu=dwu16 if (n .eq. 18) dwl=dwl18 if (n .eq. 16) dwl=dwl16 expd=qty1/(float(n-np)) vard=two*(qty2-qty1*expd)/(float(n-np)*float(n-np+2)) b=dsqrt(vard/vdu) a=expd-b*edu dstar=a+b*dwu ddstar=a+b*(4.0-dwl) return end C C C subroutine matmul(x1,y1,z1,n1,n2,n3) C THIS SUBROUTINE FORMS THE MATRIX PRODUCT OF THE TWO MATRICES X & Y C RESULT IN THE MATRIX Z. implicit real*8 (a-h,o-z) real*8 x1(n1,n2),y1(n2,n3),z1(n1,n3) data zero/0.0d0/ do 30 i=1,n1 do 20 j=1,n3 z1(i,j)=0.0 do 10 k=1,n2 z1(i,j)=z1(i,j)+x1(i,k)*y1(k,j) 10 continue 20 continue 30 continue return end C C C subroutine modprof(x,sigma0,sigma1,rho1,alrt,plrt) implicit real*8 (a-h,o-z) parameter (n=18,np=3) real*8 x(n,np),psi(n,n),psiinv(n,n) real*8 xt(np,n),xtx(np,np),temp1(np,n),temp2(np,np) real*8 eval1(np),evec1(np,np),eval2(np),evec2(np,np) external dlinds, devcsf data zero/0.0d0/,one/1.0d0/,two/2.0d0/ C FORM PSI MATRIX t1=one-rho1**2 do i=1,n do j=1,n dij=dabs(dfloat(i-j)) psi(i,j)=(rho1**dij)/t1 enddo enddo C INVERT PSI call dlinds(n,psi,n,psiinv,n) call trans(x,xt,n,np) call matmul(xt,x,xtx,np,n,np) call devcsf(np,xtx,np,eval1,evec1,np) detp1=one do j=1,np detp1=detp1*eval1(j) enddo detp1l=dlog(detp1) call matmul(xt,psiinv,temp1,np,n,n) call matmul(temp1,x,temp2,np,n,np) call devcsf(np,temp2,np,eval2,evec2,np) detp2=one do j=1,np detp2=detp2*eval2(j) enddo detp2l=dlog(detp2) sigratio=sigma0/sigma1 plrt=alrt-float(np)*dlog(sigratio)+detp1l-detp2l return end C C subroutine modscor(x,y,sigma0,beta0,sumst,sct,pscor) implicit real*8 (a-h,o-z) parameter (n=18,np=3) real*8 y(n,1),x(n,np),xt(np,n),xtx(np,np),xtxinv(np,np) real*8 xvec1(np,1),xvec2(np,1),xvec1t(1,np) real*8 xvec2t(1,np) real*8 temp3(np,np),temp1(np,np),temp2(np,np),temp4(np,np) real*8 beta0(np,1) external dlinds data zero/0.0d0/ call trans(x,xt,n,np) call matmul(xt,x,xtx,np,n,np) call dlinds(np,xtx,np,xtxinv,np) sum1=sumst/(float(n-1)*sigma0) do j1=1,np do k1=1,np temp3(j1,k1)=zero enddo enddo do i=2,n do j=1,np xvec1(j,1)=x(i-1,j) xvec2(j,1)=x(i,j) enddo call trans(xvec1,xvec1t,np,1) call trans(xvec2,xvec2t,np,1) call matmul(xvec1,xvec2t,temp1,np,1,np) call matmul(xvec2,xvec1t,temp2,np,1,np) do j=1,np do k=1,np temp3(j,k)=temp3(j,k)+temp1(j,k)+temp2(j,k) enddo enddo enddo call matmul(xtxinv,temp3,temp4,np,np,np) trc=zero do j=1,np trc=trc+temp4(j,j) enddo pscor=sct+trc*sum1 return end C C C subroutine trans(a,b,i,j) implicit real*8 (a-h,o-z) real*8 a(i,j),b(j,i) do 1 i1=1,i do 1 j1=1,j b(j1,i1)=a(i1,j1) 1 continue return end