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 PTMCOR( R,L,LAG,MJ1,N ) 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 PTMCOR( R,ID,K,MJ1,N ) C C ... Plot cross-correlation function ... C C Inputs: 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 N: data length C IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 VNAME REAL*4 SWX DIMENSION R(0:K,MJ1,MJ1), VNAME(3), VALUE(3) C WX = 16.0/ID -0.5 SWX = WX VNAME(1) = 'N =' VNAME(2) = 'ID =' VNAME(3) = 'LAG =' VALUE(1) = N VALUE(2) = ID VALUE(3) = K C CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 2.3: CROSS-CORRELATION FUNCTION' , * 40,3,VNAME,VALUE ) CALL PLOT( 2.0,16.5-SWX,-3 ) C DO 10 I=1,ID IF( I.GE.2 ) CALL PLOT( -(SWX+0.5)*(ID-1),-SWX-0.5,-3 ) DO 10 J=1,ID IF( J.GE.2 ) CALL PLOT( SWX+0.5,0.0,-3 ) IPEN = 1 IF( I.EQ.J ) IPEN = 2 CALL PLTBOX( WX,WX,0.0D0,DBLE(K),10.0D0,-1.D0,1.D0,1.D0,IPEN ) CALL PLOT( 0.0,SWX/2,3 ) CALL PLOT( SWX,SWX/2,2 ) CALL PLOTY( R(0,I,J),K+1,-1.0D0,1.0D0,WX,WX,0,1 ) 10 CONTINUE CALL PLOTE C call plot( 0.0,0.0,999 ) C RETURN 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 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) C CLOSE( IDEV ) C RETURN 1 FORMAT( A72 ) E N D SUBROUTINE HEADER( PNAME,NP,M,VNAME,VALUE ) C C ... THIS SUBROUTINE DRAW TITLE OF DATA, PROGRAM NAME, NAMES C AND VALUES OF PARAMETERS ... C C INPUTS: C PNAME: PROGRAM NAME C NP: NUMBER OF CHARACTERS OF PROGRAM NAME C M: NUMBER OF PARAMETERS C VNAME(I): NAME OF THE I-TH PARAMETER C VALUE(I): VALUE OF THE I-TH PARAMETER C INPUT VIA COMMON: C TITLE: TITLE OF THE DATA SET C CHARACTER PNAME*80, TITLE*72 CHARACTER*8 VNAME(M) REAL*8 VALUE DIMENSION VALUE(M), DAY(2), TIME(3) COMMON /CMDATA/ TITLE C CALL NEWPEN( 2 ) CALL DATE( DAY ) CALL CLOCK( TIME,1 ) C CALL SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 ) CALL SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP ) CALL NEWPEN( 1 ) CALL SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 ) CALL SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 ) YY = 17.9 DO 10 I=1,M IF( MOD(I,4).EQ.1 ) THEN YY = YY - 0.40 X1 = 2.0 X2 = 4.0 ELSE X1 = X1 + 5.5 X2 = X2 + 5.5 END IF CALL SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 ) 10 CALL NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 ) C RETURN E N D SUBROUTINE PLTBOX ( WX,WY,XMIN,XMAX,XSTEP,YMIN,YMAX,YSTEP,IPEN ) C C ... THIS SUBROUTINE DRAWS BOX WINDOW ... C C INPUTS: C WX,WY: WIDTH AND HEIGHT OF THE BOX WINDOW C XMIN,XMAX: C XSTEP: C YMIN,YMAX: C YSTEP: C IPEN: THICKNESS OR COLOUR OF THE WINDOW C REAL*8 WX, WY, XMIN, XMAX, XSTEP, YMIN, YMAX, YSTEP C CALL NEWPEN( IPEN ) CALL PLOT( 0.0,0.0,3 ) CALL PLOT( SNGL(WX),0.0,2 ) CALL PLOT( SNGL(WX),SNGL(WY),2 ) CALL PLOT( 0.0,SNGL(WY),2 ) CALL PLOT( 0.0,0.0,2 ) C CALL NEWPEN(1) DX = XSTEP/(XMAX-XMIN) NX = 1.0D0/DX DX = DX*WX NX = ( XMAX-XMIN )/XSTEP DO 10 I=1,NX CALL PLOT( DX*I,0.0,3 ) 10 CALL PLOT( DX*I,0.1,2 ) DY = YSTEP/( YMAX-YMIN ) NY = 1.0D0/DY DY = DY*WY DO 20 I=1,NY CALL PLOT( 0.0,DY*I,3 ) 20 CALL PLOT( 0.1,DY*I,2 ) RETURN E N D SUBROUTINE PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C C ... THIS SUBROUTINE DRAWS TIME VS DATA PLOT OF THE DATA ... C C INPUTS: C Y(I): DATA C N: DATA LENGTH C YMIN: MINIMUM BOUND FOR THE Y AXIS C YMAX: MAXIMUM BOUND FOR THE Y AXIS C WX,WY: WIDTH AND HEIGHT OF THE FIGURE C IPOS: STARTING POSITION OF THE DATA C ISW: = 1 CONNECT DATA BY STRAIGHT LINE C = 2 CONNECT DATA BY DASHED LINE C = 3 EXPRESS DATA BY CHARACTER (CHARACTER CODE=ISW) C MODIFIED: 8/31/90 C REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IPOS) YY = (Y(I) - YMIN)*DY IF( I.EQ.1 .AND. ISW.LT.3 ) THEN CALL PLOT( XX,YY,3 ) ELSE IF(ISW.EQ.1) CALL PLOT( XX,YY,2 ) IF(ISW.EQ.2) CALL DASHP( XX,YY,0.2 ) IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) END IF 100 CONTINUE C RETURN 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