C     PROGRAM  4.4   BOXCOX
C
C  ...  Box-Cox transformation
C
C     The following inputs are required in READTS
C        TITLE:   title of the data set
C        N:       data length
C        Y(I):    data
C     12/27/90 Y.I. 8/7/91 G.K.
C
      PARAMETER(MJ=1000)
      IMPLICIT REAL*8(A-H,O-Z )
      CHARACTER  TITLE*72
      DIMENSION  Y(MJ), Z(MJ)
      COMMON  /CMDATA/  TITLE
C
      CALL  READTS( 1,Y,N )
      CALL  GAUSSM( Y,N,YMEAN,YVAR,FFY,AICY )
      WRITE(6,600)
      WRITE(6,610)  TITLE
      WRITE(6,620)
      DO 200 II=10,-10,-1
      A = II/10.0D0
      CALL  BOXCOX( Y,N,A,Z,ZJACOB )
C
      CALL  GAUSSM( Z,N,ZMEAN,ZVAR,FFZ,AICZ )
C
      FFZT = FFZ + ZJACOB
      AICZT = AICZ-2*ZJACOB
      WRITE(6,630)  A, AICZT, FFZT, AICZ, FFZ, ZMEAN, ZVAR
  200 CONTINUE
      STOP
  600 FORMAT( 1H ,'PROGRAM 4.4:   BOX-COX TRANSFORMATION' )
  610 FORMAT( 1H ,A72 )
  620 FORMAT( 1H ,'LAMBDA',5X,'AIC''',8X,'LL''',7X,'AIC',9X,'LL',
     *            9X,'MEAN',9X,'VARIANCE' )
  630 FORMAT( 1H ,F5.2,4F11.2,2D15.6 )
      E N D
      SUBROUTINE BOXCOX( Y,N,A,Z,ZJACOB )
C
C  ...  Box-Cox transformation:  Z = (Y**A - 1)/A  ...
C
C     Inputs:
C        Y(I):   data
C        N:      data length
C        A:      lambda
C     Outputs:
C        Z(I):   transformed data
C        ZJACOB: log(Jacobian), log |dz/dy|
C
      IMPLICIT REAL*8(A-H,O-Z )
      DIMENSION  Y(N), Z(N)
C
      SUM = 0.0D0
      DO 10 I=1,N
      IF( A .NE. 0.0D0 ) THEN
         Z(I) = (Y(I)**A - 1)/A
         SUM = SUM + (A-1)*DLOG( DABS( Y(I) ) )
      ELSE
         Z(I) = DLOG( Y(I) )
         SUM = SUM - DLOG( DABS( Y(I) ) )
      END IF
   10 CONTINUE
      ZJACOB = SUM
C
      RETURN
      E N D
      SUBROUTINE GAUSSM( Y,N,YMEAN,YVAR,FF,AIC )
C
C  ...  This subroutine fits Gaussian distribution model  ...
C
C     Inputs:
C        Y(I):   data
C        N:      data length
C     Outputs:
C        YMEAN:  mean
C        YVAR:   variance
C        FF:     log-likelihood
C        AIC:    AIC = -2FF + 4
C
      IMPLICIT REAL*8( A-H,O-Z )
      DIMENSION Y(N)
      DATA  PI/3.1415926535D0/
C
      SUM = 0.0D0
      DO 10 I=1,N
   10 SUM = SUM + Y(I)
      YMEAN = SUM/N
C
      SUM = 0.0D0
      DO 20 I=1,N
   20 SUM = SUM + (Y(I)-YMEAN)**2
      YVAR = SUM/N
      FF  = -0.5D0*N*(DLOG( 2*PI*YVAR ) + 1)
      AIC = -2*FF + 4
C
      RETURN
      E N D
      SUBROUTINE  READTS( IDEV,Y,N )
      REAL*8        Y(N)
      CHARACTER*72  TITLE
      COMMON  /CMDATA/  TITLE
C
C     OPEN( IDEV,FILE='temp.dat' )
      READ( IDEV,1 )  TITLE
      READ( IDEV,* )  N
      READ( IDEV,* )  (Y(I),I=1,N)
C     CLOSE( IDEV )
      RETURN
    1 FORMAT( A72 )
      E N D