C PROGRAM 2.1 UNICOR C C ... This program computes sample autocovariance and sample C autocorrelation function ... C C The following inputs are required in the subroutine READTS. C TITLE: caption of the data set C N: data length C FORMAT: reading format C Y(I): time series, (i=1,...,N) C Parameters: C NMAX: adjustable dimension of Y (NMAX.GE.N) C LAG: maximum lag of autocovariance function C IDEV: input file specification C JDEV: output file specification C 5/16,1989 modified 12/05/90, 12/06/90, 7/18/92 C PARAMETER( NMAX=1000,LAG=30,IDEV=1,JDEV=6 ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), COV(0:LAG,4) DATA OUTMIN/-1.0D30/, OUTMAX/ 1.0D30/ C C ... read in data ... C CALL READTS( IDEV,Y,N ) C C ... compute autocovariance, autocorrelation and their error bound . C CALL UNICOR( Y,N,LAG,OUTMIN,OUTMAX,COV ) C C ... print out results ... C CALL PRACOR( JDEV,COV,N,LAG ) STOP E N D SUBROUTINE AUTCOR( C,MAXLAG,R ) C C ... This subroutine compute sample autocorrelation ... C C Inputs: C C(I): sample autocovariance C MAXLAG: maximum lag of autocovariance C Output: C R(I): sample autocorrelation C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION C(0:MAXLAG), R(0:MAXLAG) C DO 10 LAG=0,MAXLAG 10 R(LAG) = C(LAG)/C(0) RETURN E N D SUBROUTINE AUTCOV( Y,N,MAXLAG,OUTMIN,OUTMAX,C ) C C ... This subroutine computes sample autocovariance ... C C Inputs: C Y(I): time series C N: data length C MAXLAG: maximum lag of autocovariance C OUTMIN: bound for outliers in low side C OUTMAX: bound for outliers in high side C OUTPUT: C C(I): autocovariance function C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION Y(N), C(0:MAXLAG ) C C ... sample mean ... C CALL MEAN( Y,N,OUTMIN,OUTMAX,NSUM,YMEAN ) C C ... sample autocovariance ... C DO 20 LAG = 0,MAXLAG SUM = 0.0D0 DO 10 I=LAG+1,N IF( Y(I).GT.OUTMIN .AND. Y(I).LT.OUTMAX ) THEN IF( Y(I-LAG).GT.OUTMIN .AND. Y(I-LAG).LT.OUTMAX ) THEN SUM = SUM + ( Y(I)-YMEAN)*( Y(I-LAG) - YMEAN ) END IF END IF 10 CONTINUE 20 C(LAG) = SUM/NSUM C RETURN E N D SUBROUTINE ERRACF( C,N,MAXLAG,CERR,RERR ) C C ... This subroutine computes error bounds for sample autocovariance C and autocorrelation function ... C C Inputs: C C(I): sample autocovariance C N: data length C MAXLAG: maximum lag of autocovariance C Outputs: C CERR(I): error bound for I-th autocovariance C RERR(I): error bound for I-th autocorrelation C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION C(0:MAXLAG), CERR(0:MAXLAG), RERR(0:MAXLAG) C SUM = C(0)**2 CERR(0) = DSQRT( 2*SUM/N ) DO 10 LAG=1,MAXLAG IF(LAG.GE.2) SUM = SUM + 2*C(LAG-1)**2 10 CERR(LAG) = DSQRT( SUM/N ) C RERR(0) = 0.0D0 DO 20 LAG=1,MAXLAG 20 RERR(LAG) = CERR(LAG)/C(0) 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 SUBROUTINE PRACOR( JDEV,COV,N,LAG ) C C ... This subroutine print out sample autocovariances, C sample atocorrelations and their error bound ... C C Inputs: C JDEV: output file specification C C(I,J): J=1 sample autocovariance C =2 sample autocorrelation C =3 error bound for sample autocovariance C =4 error bound for sample autocorrelation C N: length of original time series C LAG: maximum lag of autocovariance function C IMPLICIT REAL*8(A-H,O-Z) CHARACTER*72 TITLE DIMENSION COV(0:LAG,4) COMMON /CMDATA/ TITLE C WRITE(JDEV,600) WRITE(JDEV,610) TITLE WRITE(JDEV,620) N, LAG WRITE(JDEV,630) (COV(I,1),I=0,LAG) WRITE(JDEV,640) (COV(I,2),I=1,LAG) C 600 FORMAT( 1H ,'PROGRAM 2.1: ' ) 610 FORMAT( 1H ,' DATA: ',A72 ) 620 FORMAT( 1H ,' N =',I5,3X,'LAG =',I3 ) 630 FORMAT( 1H0,'** AUTOCOVARIANCE FUNCTION: C(0),...,C(LAG) **',/ * (1H ,5D14.5) ) 640 FORMAT( 1H0,'** AUTOCORRELATION FUNCTION: R(1),...,R(LAG) **',/ * (1H ,10F8.4) ) RETURN E N D SUBROUTINE UNICOR( Y,N,MAXLAG,OUTMIN,OUTMAX,COV ) C C ... This subroutine computes sample autocovariance function, C sample autocorrelationfunction and their error bounds ... C C Inputs: C Y(I): time series C N: data length C MAXLAG: maximum lag of autocovariance C OUTMIN: bound for outliers in low side C OUTMAX: bound for outliers in high side C Output: C COV(I,J): I=0,...,MAXLAG C J = 1 sample autocovariance C = 2 sample autocorrelation C = 3 error bound for sample autocovariance C = 4 error bound for sample autocorrelation C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N), COV(0:MAXLAG,4) C C ... sample autocovariance function ... C CALL AUTCOV( Y,N,MAXLAG,OUTMIN,OUTMAX,COV ) C C ... sample autocorrelation function ... C CALL AUTCOR( COV,MAXLAG,COV(0,2) ) C C ... error bounds ... C CALL ERRACF( COV,N,MAXLAG,COV(0,3),COV(0,4) ) C RETURN E N D SUBROUTINE READTS( IDEV,Y,N ) REAL*8 Y(N) CHARACTER*72 TITLE COMMON /CMDATA/ TITLE C OPEN( IDEV,FILE='temp.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N READ( IDEV,* ) (Y(I),I=1,N) CLOSE( IDEV ) RETURN 1 FORMAT( A72 ) E N D