C     PROGRAM 2.3  CRSCOR
C
C  ...  This program computes sample cross-covariance and
C       sample cross-correlation functions  ...
C
C     The following inputs are required in the subroutine READMD.
C        TITLE:   title of the data
C        N:       data length
C        L:       dimension of the observation
C        IFM:     = 1   ((Y(I,J),J=1,L),I=1,N)
C                 = 2   ((Y(I,J),I=1,N),J=1,L)
C        Y(I,J):  multi-variate time series
C     Parameters:
C        MJ:      adjustable dimension of Y; (MJ.GE.N)
C        MJ1:     adjustable dimension of Y, C, R; (MJ1.GE.L)
C        LAG:     maximum lag of the cross-covariance function
C        IDEV:    input device specification
C     MAR.24,1989:  modified  12/20/90, 7/18/92
C
      PARAMETER( MJ=1000,MJ1=7,LAG=50,IDEV=1,JDEV=6 )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Y(MJ,MJ1), OUTMIN(10), OUTMAX(10)
      DIMENSION  C(0:LAG,MJ1,MJ1), R(0:LAG,MJ1,MJ1), YMEAN(10)
      DATA  OUTMIN/10*-1.0D30/, OUTMAX/10*1.0D30/
C
C  ...  read in multivariate time series  ...
C
      CALL  READMD( IDEV,MJ,Y,N,L )
C
C  ...  cross-covariance function  ...
C
      CALL  CRSCOR( Y,N,L,LAG,MJ,MJ1,OUTMIN,OUTMAX,C,R,YMEAN )
C
C  ...  print out and plot cross-correlation function  ...
C
      CALL  PRMCOR( JDEV,C,R,L,LAG,MJ1 )
C
      STOP
      E N D
      SUBROUTINE  CRSCOR( Y,N,ID,LAG,MJ,MJ1,OUTMIN,OUTMAX,C,R,YMEAN )
C
C ... cross correlation function computation ...
C
C     Inputs:
C        Y(I,J):   multi-variate time series
C        N:        data length
C        ID:       dimension of the observation
C        LAG:      maximum lag of cross-covariance
C        MJ:       adjustable dimension (MJ.GE.N)
C        MJ1:      adjustable dimension (MJ1.GE.ID)
C        OUTMIN:   bound for outliers in low side
C        OUTMAX:   bound for outliers in high side
C     Outputs:
C        C:        sample cross-covariance function
C        R:        sample cross-correlation function
C        YMEAN:    sample mean vector
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Y(MJ,MJ1), OUTMIN(MJ1), OUTMAX(MJ1)
      DIMENSION  C(0:LAG,MJ1,MJ1), R(0:LAG,MJ1,MJ1)
      DIMENSION  YMEAN(MJ1), NSUM(10)
C
      DO 10 J=1,ID
   10 CALL  MEAN( Y(1,J),N,OUTMIN(J),OUTMAX(J),NSUM(J),YMEAN(J) )
C
      DO 30 I=1,ID
      DO 30 J=1,ID
      WRITE(6,*)  I,J
      DO 30 L=0,LAG
      SUM = 0.0D0
      NNN = 0
      DO 20 II=L+1,N
      IF( Y(II,I).GT.OUTMIN(I).AND.Y(II,I).LT.OUTMAX(I) )  THEN
        IF( Y(II-L,J).GT.OUTMIN(J).AND.Y(II-L,J).LT.OUTMAX(J) )  THEN
          SUM = SUM + (Y(II,I)-YMEAN(I))*(Y(II-L,J)-YMEAN(J))
          NNN = NNN + 1
        END IF
      END IF
   20 CONTINUE
   30 C(L,I,J)=SUM/DSQRT( DBLE( NSUM(I)*NSUM(J) ) )
C
      DO 40 I=1,ID
      DO 40 J=1,ID
      DO 40 L=0,LAG
   40 R(L,I,J) = C(L,I,J)/DSQRT(C(0,I,I)*C(0,J,J))
      RETURN
C
      E N D
      SUBROUTINE  PRMCOR( JDEV,C,R,ID,K,MJ1 )
C
C ... print out cross-covariance and cross-correlation function ...
C
C     Inputs:
C        C(L,I,J):  cross-covariance function
C        R(L,I,J):  cross-correlation function
C        ID:        dimension of the time series
C        K:         maximum lag of cross-correlation
C        MJ1:       adjustable dimension (MJ1.GE.ID)
C        JDEV:      output device specification
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION   C(0:K,MJ1,MJ1), R(0:K,MJ1,MJ1)
C
      WRITE(JDEV,600)
      DO 10 I=1,ID
      DO 10 J=1,ID
      WRITE(JDEV,610)  I, J
   10 WRITE(JDEV,620)  (C(L,I,J),L=0,K)
C
      WRITE(JDEV,630)
      DO 20 I=1,ID
      DO 20 J=1,ID
      WRITE(JDEV,610)  I, J
   20 WRITE(JDEV,640)  (R(L,I,J),L=0,K)
C
      RETURN
  600 FORMAT( 1H ,'**  AUTO AND CROSS COVARIANCE FUNCTION  **' )
  610 FORMAT( 1H ,'I =',I2,3X,'J =',I2 )
  620 FORMAT( 1H ,5D15.7 )
  630 FORMAT( 1H ,'**  AUTO AND CROSS CORRELATION FUNCTION  **' )
  640 FORMAT( 1H ,10F8.4 )
      E N D
      SUBROUTINE  READMD( IDEV,NMAX,Y,N,L )
C
C   ...  THIS SUBROUTINE READS IN MULTI-VARIATE TIME SERIES  ...
C
C     INPUTS:
C         IDEV:    INPUT FILE SPECIFICATION
C         NMAX:    ADJUSTABLE DIMENSION OF Y(I,J)
C     OUTPUTS:
C         Y(I,J):  TIME SERIES (I=1,N;J=1,L)
C         N:       DATA LENGTH
C         L:       DIMENSION OF THE TIME SERIES
C
      CHARACTER*72  TITLE
      REAL*8        Y(NMAX,*)
      COMMON  /CMDATA/  TITLE
C
      OPEN( IDEV,FILE='temp2.dat' )
      READ( IDEV,1 )  TITLE
      READ( IDEV,* )  N, L, IFM
      IF( IFM.EQ.0 )  READ( IDEV,* )  ((Y(I,J),J=1,L),I=1,N)
      IF( IFM.GE.1 )  READ( IDEV,* )  ((Y(I,J),I=1,N),J=1,L)
      CLOSE( IDEV )
C
      RETURN
    1 FORMAT( A72 )
      E N D
      SUBROUTINE  MEAN( Y,N,OUTMIN,OUTMAX,NSUM,YMEAN )
C
C  ...  THIS SUBROUTINE COMPUTES SAMPLE MEAN  ...
C
C     INPUTS:
C        Y(I):    TIME SERIES
C        N:       DATA LENGTH
C        OUTMIN:  BOUND FOR OUTLIERS IN LOW SIDE
C        OUTMAX:  BOUND FOR OUTLIERS IN HIGH SIDE
C     OUTPUTS:
C        NSUM:    NUMBER OF NON-OUTLIER OBSERVAIONS
C        YMEAN:   SAMPLE MEAN
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(N)
C
      NSUM = 0
      YSUM = 0.0D0
      DO 10 I=1,N
      IF( Y(I) .GT.OUTMIN .AND. Y(I).LT.OUTMAX ) THEN
         NSUM = NSUM + 1
         YSUM = YSUM + Y(I)
      END IF
   10 CONTINUE
      YMEAN = YSUM/NSUM
C
      RETURN
      E N D