PROGRAM MAIN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC THIS PROGRAM PROVIDES BOOTSTRAPPED P-VALUES FOR SB, MT AND THE CC CC CHI-SQUARE TEST. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT DOUBLE PRECISION (A-H,Q-Z) REAL*4 RAN3 INTEGER Y(400),ITER,N,NOB DIMENSION C(3) DIMENSION XC(400),XH(400),BLOGL(20),Z3Z3(20,20),BLOG2(2) DIMENSION X(10000),OUT(3,10000) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC EXPECT, XOBSER AND XCUT ARE USED TO CONSTUCT THE CHI-SQUARE TEST CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DIMENSION EXPECT(15),XOBSER(15),XCUT(15) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC THE FILE OUTPUT2 IS GENERATED BY THE PROGRAM DISCRETE.FOR . CC CC OUTPUT2 REPORTS THE ECONOMETRIC SAMPLE SIZE, THE EMPIRICAL VALUES CC CC FOR SB AND MT, AND THE GEOMETRIC PROBABILITY. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC OPEN(13,FILE='OUTPUT2') OPEN(15,FILE='OUTPUT3') CC THE MONTE CARLO REALIZATIONS CAN BE RECORDED, IF DESIRED. OPEN(3,FILE='OUT1') OPEN(4,FILE='OUT2') OPEN(5,FILE='OUT3') IDUM = (-1.0D0)*8677206.0D0 4 READ(13,*,END=7000)N,C(1),C(2),XP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC ASSIGN C(3) THE EMPIRICAL VALUE OF THE CHI-SQUARE TEST. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C(3) = 3.0000D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC THE NUBER OF ITERATIONS IS SET TO 10000. THE USER CAN CHANGE THIS. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ITER = 10000 XN = N DO 5000 I = 1,ITER DO 5 JJ = 1,15 XOBSER(JJ) = 0.0D0 EXPECT(JJ) = 0.0D0 XCUT(JJ) = 0.0D0 5 CONTINUE DO 10 JJ = 1,N XC(JJ) = 0.0 XH(JJ) = 0.0 Y(JJ) = 0 10 CONTINUE DO 90 JJ = 1,N 15 RNN1 = RAN3(IDUM) XX = 0.0D0 DO 20 III = 1,1000 J = III-1 XX = XP*(1.0D0-XP)**J + XX IF(RNN1.LE.XX)GOTO 80 20 CONTINUE WRITE(*,110) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC A FAILURE USUALLY MEANS THAT THE GEOMETRIC PROBABILITY IS SO CC CC SMALL THAT IT IS DIFFICULT TO GENERATE DATA FROM THE GEOMETRIC CC CC DISTRIBUTION. THE PROGRAM MAY TAKE A LONG TIME TO COMPLETE, OR CC CC OR IT MIGHT NOT BE ABLE TO COMPLETE THE TASK. LOWER THE CC CC NUMBER OF ITERATIONS TO 100 TO SEE IF THE PROGRAM IS ABLE TO CC CC COMPLETE THE COMPUTATIONS IN A REASONABLE AMOUNT OF TIME. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 110 FORMAT(1H ,'FAILURE TO OBTAIN A Y VALUE') GOTO 15 80 Y(JJ) = J 90 CONTINUE YBAR = 0.0 DO 100 J = 1,N XY = Y(J) YBAR = YBAR + XY/XN 100 CONTINUE IF(YBAR.EQ.0.0D0)YBAR = 0.10D0 XPHAT = 1.0 / (YBAR + 1.0) SIG2 = 0.0 XNUM = 0.0 XDEN = 0.0 SUM1 = 0.0 SUM2 = 0.0 DO 120 JJ = 1,N XY = Y(JJ) XH(JJ) = (YBAR - XY)*(1.0 + (1.0/YBAR)) XC(JJ) = (XY - (1.0-XPHAT)/XPHAT)**2 - (1.0-XPHAT)/XPHAT**2 XNUM = XNUM + XH(JJ)*XC(JJ) XDEN = XDEN + XH(JJ)*XH(JJ) SIG2 = SIG2 + (1/XN)*(Y(JJ) - YBAR)**2 SUM1 = SUM1 + (1/XN)*(Y(JJ)**2) SUM2 = SUM2 + (1/XN)*(Y(JJ)**0.5) 120 CONTINUE IF(XDEN.EQ.0)XDEN = 0.10D0 IF(SIG2.EQ.0)SIG2 = 0.10D0 BHAT = XNUM/XDEN ROOTN = XN**0.5 XGMM = SIG2 - YBAR**2 - YBAR XSIG2 = 0.0 DO 200 JJ = 1,N XSIG2 = ((XC(JJ)-XGMM-BHAT*XH(JJ))**2)/(XN-2.0) + XSIG2 200 CONTINUE XSIG = (XSIG2**0.5)/ROOTN IF(XSIG.EQ.0)XSIG = 0.10D0 Z3 = XGMM/XSIG XMT = Z3 IF(Z3.LT.-999.0D0)Z3 = -999.0D0 IF(Z3.GT.999.0D0)Z3 = 999.0D0 WRITE(4,9020)XMT OUT(2,I) = XMT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC THE DISCRETE-TIME TESTS DO NOT REQUIRED THE DATA TO BE ORDERED, BUT SOME CC CC CONTINUOUS-TIME TESTS DO. THIS PROGRAM CAN BE MODIFIED FOR THESE TESTS. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 325 DO 300 II = 1,N YMIN = Y(II) DO 250 JJ = II,N IF(Y(JJ).GE.YMIN)GOTO 250 YMIN = Y(JJ) Y(JJ) = Y(II) Y(II) = YMIN 250 CONTINUE 300 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC COMPUTE THE CHI-SQUARE TEST STATISTIC CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUM = 0.0D0 SUM2 = 0.0D0 IYEND = Y(N) IBI = 1 DO 101 J=1,IYEND XPROB = XPHAT*(1.0D0-XPHAT)**(J-1) XPRSA = XPROB*XN SUM = XPRSA + SUM SUM2 = XPRSA + SUM2 CC WE HAVE AN EXECTED FREQUENCY OF A LITTLE MORE THAN 6 IF(SUM.LT.6.0D0)GOTO 101 EXPECT(IBI) = SUM SUM = 0.0D0 XCUT(IBI) = J-1 IBI = IBI + 1 101 CONTINUE IBI = IBI - 1 XTERM = SUM + (XN-SUM2) CC THE LAST BIN MUST HAVE AN EXPECTED FREQUENCY OF 5 OR MORE, CC OR ELSE IT WILL BE COMBINED WITH THE PENULTIMATE BIN IF(XTERM.GE.5.0D0)IBI = IBI + 1 XCUT(IBI) = 100000.0D0 EXPECT(IBI) = EXPECT(IBI) + SUM + (XN-SUM2) DO 202 J=1,N IXC = 1 DO 201 JBJ = 1,IBI IF(Y(J).GT.XCUT(JBJ))IXC = IXC + 1 201 CONTINUE XOBSER(IXC) = XOBSER(IXC) + 1 202 CONTINUE XCHI = 0.0D0 DO 203 JBJ = 1,IBI XCHI = XCHI + ( EXPECT(JBJ) - XOBSER(JBJ) )**2 / EXPECT(JBJ) 203 CONTINUE WRITE(5,9020)XCHI OUT(3,I) = XCHI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC BEGIN THE SB REGRESSION TEST CC CC WE REGRESS ST ON A CONSTANT ON DT, WHERE ST IS 0 IF THE CYCLE CC CC ENDS, AND 1 OTHERWISE. DT IS THE LAGGED DURATION OF THE CYCLE. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC XX1 = 1.0D0 DO 510 JJ = 1,2 BLOGL(JJ) = 0.0D0 DO 505 III = 1,2 Z3Z3(JJ,III) = 0.0D0 505 CONTINUE 510 CONTINUE XNOB = 0.0D0 DO 47 JJ = 1,N NOB = Y(JJ) + 1 XX2 = 0.0D0 DO 45 IT = 1,NOB XNOB = XNOB + 1.0D0 YS = 1.0D0 XX2 = XX2 + 1.0D0 IF(IT.EQ.NOB)YS = 0.0D0 Z3Z3(1,1) = Z3Z3(1,1) + XX1*XX1 Z3Z3(1,2) = Z3Z3(1,2) + XX1*XX2 Z3Z3(2,1) = Z3Z3(1,2) Z3Z3(2,2) = Z3Z3(2,2) + XX2*XX2 BLOGL(1) = BLOGL(1) + XX1*YS BLOGL(2) = BLOGL(2) + XX2*YS 45 CONTINUE 47 CONTINUE NNN = 2 CALL DMATEQ(Z3Z3,NNN,BLOG2,DET) XAY2 = Z3Z3(1,1)*BLOGL(1) + Z3Z3(1,2)*BLOGL(2) XBY2 = Z3Z3(2,1)*BLOGL(1) + Z3Z3(2,2)*BLOGL(2) 48 SIG11 = 0.0D0 DO 67 JJ = 1,N NOB = Y(JJ) + 1 XX2 = 0.0D0 DO 65 IT = 1,NOB YS = 1.0D0 XX2 = XX2 + 1.0D0 IF(IT.EQ.NOB)YS = 0.0D0 E1 = YS - XAY2*XX1 - XBY2*XX2 SIG11 = SIG11 + (E1**2)/(XNOB-2.0D0) 9030 FORMAT(1H ,4(F15.5,1X)) 65 CONTINUE 67 CONTINUE CC CONSTRUCT TEST STATISTIC IF(SIG11.EQ.0)SIG11 = 0.10D0 IF(Z3Z3(2,2).EQ.0)Z3Z3(2,2) = 0.10D0 SB = XBY2 / ( (SIG11*Z3Z3(2,2))**0.5 ) WRITE(3,9020)SB OUT(1,I) = SB 9020 FORMAT(F15.8) 5000 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC THIS ROUTINE SORTS A UNIVARIATE SERIES OF 10000 OBSERVATIONS AND CC CC RETURNS A P-VALUE. IT IS ASSSUMED THAT ITER = 10000 FROM ABOVE. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 2301 JJX = 1,3 XCRIT = C(JJX) DO 6200 I=1,10000 X(I) = OUT(JJX,I) 6200 CONTINUE DO 6300 I = 1,10000 XMIN = X(I) DO 6250 J = I,10000 IF(X(J).GE.XMIN)GOTO 6250 XMIN = X(J) X(J) = X(I) X(I) = XMIN 6250 CONTINUE 6300 CONTINUE DO 6400 I=1,10000 IF(XCRIT.GT.X(I))GOTO 6400 XI = I IF(JJX.EQ.1)WRITE(15,2001) 2001 FORMAT(1H ,'THE SB P-VALUE IS '/) IF(JJX.EQ.2)WRITE(15,2002) 2002 FORMAT(1H ,'THE MT P-VALUE IS '/) IF(JJX.EQ.3)GOTO 2302 IF(I.GT.5000)GOTO 6350 XI = I XP1 = 2.0D0 * XI / 10000.0D0 WRITE(15,2000)XP1,N,XCRIT,XP WRITE(15,2100) 2000 FORMAT(1H ,'P-VALUE = ',F10.5,' N =',I4,' XCRIT =',F10.5, & ' GEOMETRIC PROB =',F10.7) 2100 FORMAT(1H ,'LOWER TAIL'/) GOTO 2301 6350 XI = I XP1 = 2.0D0*(10000.0D0 - XI)/10000.0D0 WRITE(15,2000)XP1,N,XCRIT,XP WRITE(15,2200) 2200 FORMAT(1H ,'UPPER TAIL'/) GOTO 2301 6400 CONTINUE IF(JJX.EQ.3)WRITE(15,2003) XP1 = 0.0D0 WRITE(15,2000)XP1,N,XCRIT,XP WRITE(15,2200) 1000 FORMAT(1H ,F20.8) 2301 CONTINUE GOTO 4 2302 WRITE(15,2003) 2003 FORMAT(1H ,'THE CHI-SQUARE P-VALUE IS ',/) XP1 = (10000.0D0 - XI)/10000.0D0 WRITE(15,2000)XP1,N,XCRIT,XP WRITE(15,2200) GOTO 4 7000 CONTINUE STOP END CC UNIFORM RANDOM NUMBER GENERATOR FUNCTION RAN3(IDUM) C IMPLICIT REAL*4(M) C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=2.5E-7) PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.E-9) DIMENSION MA(55) DATA IFF /0/ IF(IDUM.LT.0.OR.IFF.EQ.0)THEN IFF=1 MJ=MSEED-IABS(IDUM) MJ=MOD(MJ,MBIG) MA(55)=MJ MK=1 DO 11 I=1,54 II=MOD(21*I,55) MA(II)=MK MK=MJ-MK IF(MK.LT.MZ)MK=MK+MBIG MJ=MA(II) 11 CONTINUE DO 13 K=1,4 DO 12 I=1,55 MA(I)=MA(I)-MA(1+MOD(I+30,55)) IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG 12 CONTINUE 13 CONTINUE INEXT=0 INEXTP=31 IDUM=1 ENDIF INEXT=INEXT+1 IF(INEXT.EQ.56)INEXT=1 INEXTP=INEXTP+1 IF(INEXTP.EQ.56)INEXTP=1 MJ=MA(INEXT)-MA(INEXTP) IF(MJ.LT.MZ)MJ=MJ+MBIG MA(INEXT)=MJ RAN3=MJ*FAC RETURN END SUBROUTINE DMATEQ (A,N,Y,DET) TOB11970 IMPLICIT DOUBLE PRECISION(A-H,O-Z) TOB11980 DIMENSION A(20,20),Y(20) TOB11990 DOUBLE PRECISION DABS TOB12000 C DMATEQ SUBROUTINE (/360 FILE NO. 310.4.501) TOB12010 C /360 REAL*8 VERSION OF MATEQ (7074 FILE NO. 10.4.502) TOB12020 C A = REAL*8 MATRIX TO BE INVERTED TOB12030 C N = INTEGER*4 VARIALBLE FOR ROW SIZE OF MATRIX A TOB12040 C Y = REAL*8 VECTOR TO SOLVE AX=Y TOB12050 C DET = REAL*8 VARIABLE FOR DETERMINANT OF A TOB12060 C USER SHOULD CHECK DET IMMEDIATELY FOR SINGULARITY TOB12070 C A BECOMES AINVERSE, YOFAX=Y BECOMES X, DET IS DET A TOB12080 C DMATEQ SOLVES AX=Y FOR X, COMPUTES A INVERSE, AND CALCULATES THE TOB12090 C DETERMINANT USING A VARIATION OF GAUSSIAN ELIMINATION. TOB12100 C P J EBERLEIN, REVISED BY E C REINEMAN FOR /360 TOB12110 C UNIVERSITY OF ROCHESTER COMPUTING CENTER TOB12120 C AUGUST 3, 1967 TOB12130 C KEYPUNCH IS 029 TOB12140 DIMENSION ICHG(20) TOB12150 DET=1.0D0 TOB12160 DO 118 K=1,N TOB12170 AMX=DABS(A(K,K)) TOB12180 IMX=K TOB12190 DO 100 I=K,N TOB12200 IF (DABS(A(I,K)).LE.AMX) GO TO 100 TOB12210 AMX=DABS(A(I,K)) TOB12220 IMX=I TOB12230 100 CONTINUE TOB12240 IF (DABS(AMX).GT.0.1D-70) GO TO 102 TOB12250 DET=0.0D0 TOB12260 GO TO 124 TOB12270 102 IF (IMX.EQ.K) GO TO 106 TOB12280 DO 104 J=1,N TOB12290 TEMP=A(K,J) TOB12300 A(K,J)=A(IMX,J) TOB12310 104 A(IMX,J)=TEMP TOB12320 ICHG(K)=IMX TOB12330 TEMP=Y(K) TOB12340 Y(K)= Y(IMX) TOB12350 Y(IMX)= TEMP TOB12360 DET=-DET TOB12370 GO TO 108 TOB12380 106 ICHG(K)=K TOB12390 108 DET=DET*A(K,K) TOB12400 A(K,K)=1.0D0/A(K,K) TOB12410 DO 110 J=1,N TOB12420 IF (J.NE.K) A(K,J)=A(K,J)*A(K,K) TOB12430 110 CONTINUE TOB12440 Y(K) = Y(K)*A(K,K) TOB12450 DO 114 I=1,N TOB12460 DO 112 J=1,N TOB12470 IF (I.EQ.K) GO TO 114 TOB12480 IF (K.NE.J) A(I,J)=A(I,J)-A(I,K)*A(K,J) TOB12490 112 CONTINUE TOB12500 Y(I) = Y(I)-A(I,K)*Y(K) TOB12510 114 CONTINUE TOB12520 DO 116 I=1,N TOB12530 IF (I.NE.K) A(I,K)=-A(I,K)*A(K,K) TOB12540 116 CONTINUE TOB12550 118 CONTINUE TOB12560 DO 122 K=1,N TOB12570 L=N+1-K TOB12580 KI=ICHG(L) TOB12590 IF (L.EQ.KI) GO TO 122 TOB12600 DO 120 I=1,N TOB12610 TEMP = A(I,L) TOB12620 A(I,L) = A(I,KI) TOB12630 120 A(I,KI) = TEMP TOB12640 122 CONTINUE TOB12650 124 RETURN TOB12660 END TOB12670