PROGRAM TRIPLES

C*****************************************************************************
C    RANDAL J. VERBRUGGE
C    DEPT. OF ECONOMICS
C    VPI&SU, BLACKSBURG, VA  24061
C    address@hidden
C
C    NOV. 10, 1995
C    REVISED JUNE 17, 1996
C    SECOND REVISION JUNE 21, 1996
C    REVISIONS INCORPORATED INTO CURRENT PROGRAM: MARCH 31, 1997
C
C
C    THIS PROGRAM TAKES DATA FROM C:\DATA.DAT AND CONDUCTS THE TRIPLES
C    ASYMMETRY TEST ON IT
C    THIS TEST IS DESCRIBED IN RANDLES, FLIGNER, POLICELLO AND WOLFE
C    'AN ASYMPTOTICALLY DISTRIBUTION-FREE TEST FOR SYMMETRY VS. ASYMMETRY'
C    JASA 75 (MARCH 1980) 168-172
C
C    IMPORTANT NOTES: 
C    1) AN ATTEMPTED PORT OF THIS PROGRAM TO A WORKSTATION
C            FAILED FOR UNKNOWN REASONS. 
C    2) TO TEST WHETHER THE CODE IS WORKING PROPERLY, ENTER THE FOLLOWING
C            NUMBERS IN THE FILE C:\DATA.DAT:
C                       2.373
C                       3.339
C                       1.980
C                       3.102
C                       0.000
C                       3.335
C            OUTPUT OF THE TEST SHOULD BE:
C               ETA = -.23333
C               VARIANCE = .013333
C               TEST STATISTIC = -2.0207
C
C    THE NULL OF THE TEST IS A SYMMETRIC DISTRIBUTION
C    THIS TEST CANNOT DETECT ASYMMETRIC DISTRIBUTIONS WITH MEDIAN=MEAN
C    THE TEST REQUIRES A MINIMUM OF 5 DATA POINTS
C
C    ETA IS THE ESTIMATED VARIABLE
C    VARIANCE IS ITS ESTIMATED VARIANCE (I ALREADY DIVIDE BY N)
C    N = NUMBER OF OBSERVATIONS
C    TEST STATISTIC IS ETA/SQRT(VARIANCE), WHICH IS ASYMPTOTICALLY STD. N
C
C    FSTAR IS A FUNCTION THAT CALCULATES FSTAR(X(I),X(J),X(K))
C        ARGUMENTS I,J,K
C    FX1I IS A FUNCTION THAT CALCULATES FSTAR(1)(X(I))
C        ARGUMENTS: I
C    FX2IJ IS A FUNCTION THAT CALCULATES FSTAR(2)(X(I),X(J))
C        ARGUMENTS: I,J
C
C    IF YOUR DATA SET IS LARGER THAN 600, NEED TO CHANGE MAXN IN PARAMETER
C    STATEMENT BELOW AND IN ALL 3 FUNCTIONS
C*****************************************************************************

      PARAMETER (MAXN=600)

      INTEGER N,I,J,K,HALFN

      REAL X(MAXN),NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,KSI1,KSI2,
     +KSI3,TEMP,TEMP2,VARIANCE,STAT

      EXTERNAL FSTAR,FX1I,FX2IJ

      COMMON /A/ N,NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,X

      OPEN(UNIT=10,FILE='C:\DATA.DAT')

C     READ IN, AND COUNT NUMBER OF, DATA POINTS

      N=1
  10  READ(10,*,END=100) X(N)
      N=N+1
      GOTO 10

C     INITIALIZE, DO COMBINATORIAL CALCULATIONS

  100 N=N-1
      WRITE(*,1003) N
      HALFN=INT(REAL(N)/2.0)
      NN=REAL(N)
      NNMIN2=NN-2.0
      NNCHUS3=(NN*(NN-1.0)*(NN-2.0))/6.0
      WRITE(*,*)'N CHOOSE 3 IS ',NNCHUS3
      NNCHUS2=(NN*(NN-1.0))/2.0
      NNMIN1C2=((NN-1.0)*(NN-2.0))/2.0

C     CALCULATE ETA

      ETA=0.0
      DO 200 I=1,N-2
        DO 300 J=I+1,N-1
          DO 400 K=J+1,N
            ETA=ETA+FSTAR(I,J,K)
  400     CONTINUE
  300   CONTINUE
  200 CONTINUE
      ETA=ETA/NNCHUS3

      WRITE(*,1000) ETA

C     CALCULATE KSI1

      KSI1=0.0
      DO 500 I=1,N
         TEMP=FX1I(I)
         TEMP2=(TEMP-ETA)*(TEMP-ETA)
         KSI1=KSI1+TEMP2
  500 CONTINUE
      KSI1=KSI1/NN

C     CALCULATE KSI2

      KSI2=0.0
      DO 600 I=1,N-1
        DO 700 J=I+1,N
         TEMP=FX2IJ(I,J)
         TEMP2=(TEMP-ETA)*(TEMP-ETA)
         KSI2=KSI2+TEMP2
  700   CONTINUE
  600 CONTINUE
      KSI2=KSI2/NNCHUS2

C     CALCULATE KSI3

      KSI3=(1.0/9.0)-(ETA*ETA)

C     CALCULATE VARIANCE

      VARIANCE=KSI1*3.0*(NN-3.0)*(NN-4.0)/2.0
      VARIANCE=VARIANCE+(KSI2*3.0*(NN-3.0))+KSI3
      VARIANCE=VARIANCE/NNCHUS3
      WRITE(*,1001) VARIANCE

C     CALCULATE TEST STATISTIC (WHICH IS STD. N)

      IF (VARIANCE .NE. 0.0) THEN
          STAT=ETA/SQRT(VARIANCE)
      ELSE
          STAT = 999999
          WRITE(*,*) 'TEST IS CLAIMING 0 VARIANCE'
          STOP
      ENDIF
      WRITE(*,1002) STAT

 1000 FORMAT(1X,'ETA = ',F9.6)
 1001 FORMAT(1X,'VARIANCE = ',F13.10)
 1002 FORMAT(1X,'TEST STATISTIC = ',F8.4)
 1003 FORMAT(1X,'NUMBER OF DATA POINTS = ',I4)

      END

C     **************************************************************
C
      REAL FUNCTION FSTAR(I,J,K)
      PARAMETER (MAXN=600)

      REAL NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,X(MAXN),T1,T2,T3,T4,
     +TMP
      INTEGER N,I,J,K
      COMMON /A/ N,NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,X

      T1=X(I)+X(J)-2.0*X(K)
      T2=X(I)+X(K)-2.0*X(J)
      T3=X(K)+X(J)-2.0*X(I)
      T4=T1*T2*T3

      IF (T4 .EQ. 0.0) THEN
        TMP = 0.0
      ELSE
        T1=SIGN(1.0,T1)
        T2=SIGN(1.0,T2)
        T3=SIGN(1.0,T3)
        TMP=(T1+T2+T3)/3.0
      ENDIF
      FSTAR=TMP
      RETURN
      END


C     **************************************************************
C
      REAL FUNCTION FX1I(I)
      PARAMETER (MAXN=600)

      REAL NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,X(MAXN),T1
      INTEGER N,I,J,K
      COMMON /A/ N,NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,X

      T1=0.0
      DO 5000 J=1,I-2
        DO 5010 K=J+1,I-1
            T1=T1+FSTAR(J,K,I)
 5010   CONTINUE           
 5000 CONTINUE
      DO 5015 J=1,I-1
         DO 5020 K=I+1,N
            T1=T1+FSTAR(J,I,K)
 5020   CONTINUE           
 5015 CONTINUE
      DO 5025 J=I+1,N-1
        DO 5030 K=J+1,N
            T1=T1+FSTAR(I,J,K)
 5030   CONTINUE           
 5025 CONTINUE

      FX1I=T1/NNMIN1C2
      RETURN
      END


C     **************************************************************
C
      REAL FUNCTION FX2IJ(I,J)
      PARAMETER (MAXN=600)

      REAL NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,X(MAXN),T1
      INTEGER N,I,J,K,II,JJ
      COMMON /A/ N,NN,NNMIN2,NNCHUS3,NNCHUS2,NNMIN1C2,X

      II=I
      JJ=J
      IF(I .GT. J) THEN
         II=J
         JJ=I
      ENDIF

      T1=0.0
      DO 6000 K=1,II-1
         T1=T1+FSTAR(K,II,JJ)
 6000 CONTINUE
      DO 6001 K=II+1,JJ-1
         T1=T1+FSTAR(II,K,JJ)
 6001 CONTINUE
      DO 6002 K=JJ+1,N
         T1=T1+FSTAR(II,JJ,K)
 6002 CONTINUE

      FX2IJ=T1/NNMIN2
      RETURN
      END