      PROGRAM DOCFSAMP
C FORTRAN CODE FOR EPPS-SINGLETON TEST                            

c amy, this was the first change I made.
c my PC couldn't handle NMAX=150000 and my 
c data set was small--so I changed it to 
c NMAX=15000.
      PARAMETER(NMAX=15000)
      DIMENSION X(NMAX,2),XT(8*NMAX),X1X2(2*NMAX),
     +          WSQ(2),W(NMAX,2),N(2),NOBS(2),IQR(2)
      CHARACTER FILENAME*30
c
c     Reference:      R.V. Burkhauser, A.C. Cutts, M.C. Daly, and
c                     S.P. Jenkins
c                     "Testing the Significance of Income Distribution
c                     Changes Over the 1980s Business Cycle: A Cross-
c                     national Comparison"
c                     Journal of Applied Economtrics, forthcoming
c
c     Programmer:     Amy Crews Cutts            
c                     Housing Economics and Financial Research
c                     Freddie Mac
c                     8200 Jones Branch Drive, Mailstop 484
c                     McLean, VA 22102
c                     (703) 903-2321
c                     (703) 903-4045 (fax)
c                     amy_crews@freddiemac.com
c                     amycrews@erols.com
c                                                     
c     Version:        April, 1998
c
c     Original Reference:
c                     T.W. Epps and K.J. Singleton     
c                     J. Statist. Comp. and Simulation, 
c                     V.26 (1986), pp. 177-203.      
c                                                     
c     Original Programmer:
c                     T.W. Epps                        
c                     Economics Department, Rouss Hall  
c                     University of Virginia             
c                     Charlottesville, VA  22901          
c                                               
c     Original Version: October, 1992 
c
c     Purpose:        Test that two samples of sizes N1 and N2
c                     have the same distribution.  Variables  
c                     may be either continuous or count data, 
c                     but observations on each must be i.i.d. 
c                     The test would rarely be appropriate    
c                     for time series.                        
c                                                             
c                     Modified version takes account of sample
c                     weights.  Reference paper is 
c                     "Where in the World is the Middle Class:
c                     A CrossNational Comparison of the Vanishing
c                     Middle Class Using Kernel Density Estimates"
c                     by Burkhauser, Crews, Daly, and Jenkins,
c                     (1996) Center for Policy Research, Syracuse
c                     University.
c
c                     It is assumed that the sample weights will 
c                     sum to n in each sample
c
c     NOTES:          The data file and parameter file are input at
c                     prompts.  The datafile should contain the 
c                     random variable and sample weight for sample 1
c                     followed by the random variable and sample
c                     weight for sample 2.  The parameter file
c                     should contain, in order, N1 and N2 which are 
c                     the number of observations in the data file;
c                     NOBS1 and NOBS2 which are the number of 
c                     observations in the raw data files.  If the 
c                     data were not compressed by summing sample
c                     weights for multiple obs on a value of x
c                     then n=nobs;
c                     IQR1 and IQR2 which are the lower and upper
c                     obs. #s of the interquartile range of the 
c                     combined samples.  These values should use the
c                     sampling weights in the calculation.
c                     WSQ1 and WSQ2 which are the sum
c                     of sample weights squared for each sample;
c                     The parameter file must be in the format
c                     (1x,6(i5,1x),2(f8.0));
c                     
c                     
c                                     
c     Steps:       1. Data must be in 2(f8.0,1x) format, a record 
c                     comprising one observation on the x1, w1,
c                     x2, w2  whose distributions are to 
c                     be compared.  Entries with decimals may 
c                     appear  anywhere  inside the respective 
c                     8-col.  fields.  Integer  data must be 
c                     right  justified.  Leave one space between
c                     each column of data.  If sample sizes for 
c                     the two variables differ,  leave blanks 
c                     to  fill out the shorter column.   This 
c                     executable program allows  up to  5,000 
c                     observations per variable.              
c                                                             
c                  2. Enter the complete path to the parameter
c                     file
c                                                             
c                  3. Enter complete path of data file.       
c                                                             
c     Output:         Program returns value of test statistic,
c                     CF2, distributed as chi-squared(4) under
c                     Ho.  Reject at level à if CF2 ò Xý(4;à).
c
      PRINT*,'                                                        '
 1    PRINT*,'Enter the path of the parameter file:          '
      READ(5,'(A30)') FILENAME
      OPEN(1,FILE=FILENAME)
        READ(1,10) N(1),N(2),NOBS(1),NOBS(2),IQR(1),IQR(2),
     +     WSQ(1),WSQ(2)
   10 FORMAT(1x,4(i5,1x),2(i6,1x),2(f15.4,1x))
      PRINT*,'N(1) and N(2) are: '
      WRITE(6,'(3x,i5,2x,i5)')N(1),N(2)
      PRINT*,'NOBS(1) and NOBS(2) are: '
      WRITE(6,'(3x,i5,2x,i5)')NOBS(1),NOBS(2)
      PRINT*,'IQR(1) and IQR(2) are: '
      WRITE(6,'(3x,i5,2x,i5)')IQR(1),IQR(2)
      PRINT*,'WSQ(1) is :'
      WRITE(6,'(3x,f15.0)') WSQ(1)
      PRINT*,'WSQ(2) is :'
      WRITE(6,'(3x,f15.0)')WSQ(2)
      NNOBS=NOBS(1)+NOBS(2)
      NN=N(1)+N(2)
      MAXN=MAX0(N(1),N(2))
      PRINT*,'Enter complete path of file where sample 1 data are.'
      READ(5,'(A30)') FILENAME
      OPEN(2,FILE=FILENAME)
      PRINT*,'Reading the data...'
      XBAR=0.
      DO I=1,N(1)
         READ(2,*) X(I,1),W(I,1)
         XBAR=XBAR+X(I,1)*W(I,1)/FLOAT(NOBS(1))
      END DO
      PRINT*,'The mean value of data in sample 1 is:              '
      WRITE(6,'(3x,f15.8)')XBAR
      PRINT*,'                                                     '
      PRINT*,'Enter complete path of file where sample 2 data are.'
      READ(5,'(A30)') FILENAME
      OPEN(3,FILE=FILENAME)
      PRINT*,'Reading the data...'
      XBAR=0.
      DO I=1,N(2)
         READ(3,*) X(I,2),W(I,2)
C Amy, This was the second change, I changed it to 
C Float(NOBS(2)) since I had a differnt sample size in 
C the second group.
         XBAR=XBAR+X(I,2)*W(I,2)/FLOAT(NOBS(2))
      END DO
      PRINT*,'The mean value of data in sample 2 is:              '
      WRITE(6,'(3x,f15.8)')XBAR
      PRINT*,'                                                     '
      PRINT*,'  Computing test statistic...'
 2    WRITE(6,'(2(/))')
      CALL CF2SAMP(N,NN,NOBS,NNOBS,NMAX,X,XT,X1X2,W,WSQ,IQR,CF2)
      STOP
      END
C
      SUBROUTINE CF2SAMP(N,NN,NOBS,NNOBS,NMAX,X,XT,X1X2,W,WSQ,IQR,CF2)
      INTEGER N(2),NOBS(2),IWK(4),IQR(2)
      REAL X(NMAX,2),W(NMAX,2),T(2),XT(2,4,NMAX),S(4,4),SS(2,4,4),
     +     G(2,4),X1X2(2*NMAX),WSQ(2),OMINV(4,4)
      DATA T/.4,.8/
      DO I=1,N(1)
         X1X2(I)=X(I,1)
      END DO
      DO I=1,N(2)
         X1X2(N(1)+I)=X(I,2)
      END DO
      CALL SORT(NN,X1X2)
      XSCALE=.25*(X1X2(IQR(1)-1)+X1X2(IQR(1))-X1X2(IQR(2)+1)
     +       -X1X2(IQR(2)))
      IF(XSCALE.LE.0.) XSCALE=1.
      DO J=1,2
         T(J)=T(J)/XSCALE
      END DO
      DO NS=1,2
         DO J=1,4
            G(NS,J)=0.
            DO JJ=1,4
               SS(NS,J,JJ)=0.
            END DO
         END DO
      END DO
      DO J=1,2
         DO NS=1,2
            DO I=1,N(NS)
               XT(NS,2*J-1,I) = COS(T(J)*X(I,NS))
               XT(NS,2*J  ,I) = SIN(T(J)*X(I,NS))
            END DO
         END DO
      END DO
      DO NS=1,2
         DO I=1,N(NS)
            DO J=1,4
               G(NS,J) = G(NS,J) + XT(NS,J,I)*W(I,NS)/FLOAT(NOBS(NS))
               DO JJ=1,J
                  SS(NS,J,JJ) = SS(NS,J,JJ) + XT(NS,J,I)*XT(NS,JJ,I)
     +                          *W(I,NS)/FLOAT(NOBS(NS))
                  SS(NS,JJ,J) = SS(NS,J,JJ)
               END DO
            END DO
         END DO
      END DO
      DO J=1,4
         DO JJ=1,J
            S(J,JJ) = (WSQ(1)*(SS(1,JJ,J)-G(1,J)*G(1,JJ))
     +                +WSQ(2)*(SS(2,JJ,J)-G(2,J)*G(2,JJ)))/FLOAT(NNOBS)
     +               *(1./FLOAT(NOBS(1))+1./FLOAT(NOBS(2)))
            S(JJ,J)=S(J,JJ)
         END DO
      END DO
      CALL MATINV(S,4,4,OMINV,IWK)
      CF2=0.
      DO J=1,4
         DO JJ=1,4
            CF2=CF2+(G(1,J)-G(2,J))*OMINV(J,JJ)*(G(1,JJ)-G(2,JJ))
         END DO
      END DO
c           CF2=CF2*FLOAT(NNOBS)
      PRINT*,'(G prime S inverse G) is equal to:             '
      WRITE(6,'(3x,f15.8)') CF2
      PRINT*,'                                               '
c     COR=1.+ NNOBS**(-.45) + 10.1*(1./(FLOAT(NOBS(1)))**1.7 + 
c    + 1./(FLOAT(NOBS(2))**1.7)
c     CF2=CF2/COR
      WTILDE=(WSQ(1)+WSQ(2))/(FLOAT(NNOBS))
      PRINT*,'WTILDE is equal to:                            '
      WRITE(6,'(3x,f15.8)') WTILDE
      PRINT*,'                                               '
      CF2=CF2*FLOAT(NOBS(1))*FLOAT(NOBS(2))/WTILDE
      WRITE(6,100) CF2
 100  FORMAT('CF statistic, distributed as Xý(4) under Ho: ', F15.4)
      WRITE(6,110)NNOBS
 110  FORMAT('NNOBS equals ', I9)
      RETURN
      END
C Subroutines based on routines in Numerical Recipes (Cambridge Press, 1989)
      SUBROUTINE SORT(N,RA)
      DIMENSION RA(N)
      L=N/2+1
      IR=N
 10   CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          RRA=RA(L)
        ELSE
          RRA=RA(IR)
          RA(IR)=RA(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            RA(1)=RRA
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
 20     IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(RA(J).LT.RA(J+1))J=J+1
          ENDIF
          IF(RRA.LT.RA(J))THEN
            RA(I)=RA(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        RA(I)=RRA
      IIII=1
      IF(IIII.EQ.1) GO TO 10
      END
      SUBROUTINE MATINV(A,N,NP,AINV,IWK)
      DIMENSION A(NP,NP),AINV(NP,NP),IWK(NP)
      DO 2 I=1,N
         DO 1 J=1,N
         AINV(I,J)=0.
 1       CONTINUE
      AINV(I,I)=1.
 2    CONTINUE
      CALL LUDCMP(A,N,NP,IWK,D)
      DO 3 J=1,N
      CALL LUBKSB(A,N,NP,IWK,AINV(1,J))
 3    CONTINUE
      RETURN
      END
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      PARAMETER (NMAX=100,TINY=1.0E-20)
      DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
 11     CONTINUE
        IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
        VV(I)=1./AAMAX
 12   CONTINUE
      DO 19 J=1,N
        IF (J.GT.1) THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF (I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
 13           CONTINUE
              A(I,J)=SUM
            ENDIF
 14       CONTINUE
        ENDIF
        AAMAX=0.
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
 15         CONTINUE
            A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
 16     CONTINUE
        IF (J.NE.IMAX)THEN
          DO 17 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
 17       CONTINUE
          D=-D
          VV(IMAX)=VV(J)
        ENDIF
        INDX(J)=IMAX
        IF(J.NE.N)THEN
          IF(A(J,J).EQ.0.)A(J,J)=TINY
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
 18       CONTINUE
        ENDIF
 19   CONTINUE
      IF(A(N,N).EQ.0.)A(N,N)=TINY
      RETURN
      END
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
      DIMENSION A(NP,NP),INDX(N),B(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0)THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
 11       CONTINUE
        ELSE IF (SUM.NE.0.) THEN
          II=I
        ENDIF
        B(I)=SUM
 12   CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF(I.LT.N)THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
 13       CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
 14   CONTINUE
      RETURN
      END


