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