SUBROUTINE TIME(X,N) C C PURPOSE--THIS SUBROUTINE PERFORMS A TIME SERIES ANALYSIS C ON THE DATA IN THE INPUT VECTOR X. C THE ANALYSIS CONSISTS OF THE FOLLOWING-- C 1) A PLOT OF AUTOCORRELATION VERSUS LAG NUMBER; C 2) A TEST FOR WHITE NOISE (ASSUMING NORMALITY); C 3) A 'PILOT' SPECTRUM; AND C 4) 4 OTHER ESTIMATED SPECTRA--EACH BASED C ON A DIFFERING BANDWIDTH. C C IN ORDER THAT THE RESULTS OF THE TIME SERIES ANALYSIS C BE VALID AND PROPERLY INTERPRETED, THE INPUT DATA C IN X SHOULD BE EQUI-SPACED IN TIME C (OR WHATEVER VARIABLE CORRESPONDS TO TIME). C C THE HORIZONTAL AXIS OF THE SPECTRA PRODUCED C BY THIS SUBROUTINE IS FREQUENCY. C THIS FREQUENCY IS MEASURED IN UNITS OF C CYCLES PER 'DATA POINT' OR, MORE PRECISELY, IN C CYCLES PER UNIT TIME WHERE C 'UNIT TIME' IS DEFINED AS THE C ELAPSED TIME BETWEEN ADJACENT OBSERVATIONS. C THE RANGE OF THE FREQUENCY AXIS IS 0.0 TO 0.5. C C INPUT ARGUMENTS--X = THE SINGLE PRECISION VECTOR OF C (UNSORTED) OBSERVATIONS. C N = THE INTEGER NUMBER OF OBSERVATIONS C IN THE VECTOR X. C OUTPUT--7 TO 11 PAGES (DEPENDING ON C THE INPUT SAMPLE SIZE) OF C AUTOMATIC PRINTOUT-- C 1) A PLOT OF AUTOCORRELATION VERSUS LAG NUMBER; C THIS PLOT MAY TAKE AS LITTLE AS 1 C OR AS MANY AS 5 PAGES C (THE EXACT NUMBER DEPENDING ON C THE INPUT SAMPLE SIZE N); C 2) A TEST FOR WHITE NOISE (ASSUMING NORMALITY); C 3) A 'PILOT' SPECTRUM; AND C 4) AN ESTIMATED SPECTRUM BASED ON A C BANDWIDTH DERIVED FROM THE DATA SET; C 5) AN ESTIMATED SPECTRUM BASED ON A C BANDWIDTH ONLY 1/2 AS WIDE AS IN 4; C 6) AN ESTIMATED SPECTRUM BASED ON A C BANDWIDTH ONLY 1/4 AS WIDE AS IN 4; C 7) AN ESTIMATED SPECTRUM BASED ON A C BANDWIDTH ONLY 1/8 AS WIDE AS IN 4; C PRINTING--YES. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C --THE SAMPLE SIZE N MUST BE GREATER C THAN OR EQUAL TO 3. C OTHER DATAPAC SUBROUTINES NEEDED--PLOTC0, PLOTSP, AND CHSPPF. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--THE 'FAST FOURIER TRANSFORM' IS NOT USED C IN THIS VERSION OF TIME, BUT WILL BE C IMPLEMENTED IN A FUTURE VERSION. C --THE USUAL MAXIMUM NUMBER OF LAGS C FOR WHICH THE AUTOCORRELATION IS C COMPUTED IS N/4 WHERE N IS C THE SAMPLE SIZE (LENGTH OF THE C DATA RECORD IN THE VECTOR X). C THIS RULE IS OVERRIDDEN IN C LARGE DATA SETS AND IS REPLACED C BY THE RULE THAT THE MAXIMUM C NUMBER OF LAGS = 500 C (WHICH CORRESPONDS TO AN C AUTOCORRELATION PLOT COVERING C 5 COMPUTER PAGES). C IF MORE LAGS ARE DESIRED, C CHANGE THE VALUE OF THE C VARIABLE MAXLAG C WITHIN THIS SUBROUTINE C FROM 500 TO WHATEVER DESIRED, C AND ALSO CHANGE THE DIMENSION OF C THE VECTOR R FROM ITS PRESENT 500 TO HOWEVER C MANY LAGS ARE DESIRED. C --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED C TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME, C THEN THE FREQUENCY AXIS OF THE RESULTING C SPECTRA WOULD BE IN UNITS OF HERTZ C (= CYCLES PER SECOND). C --THE FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE C IN THE DATA OF INFINITE (= 1/(0.0)) C LENGTH OR PERIOD. C THE FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE C IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS. C --ANY EQUI-SPACED TIME SERIES ANALYSIS IS C INTRINSICALLY LIMITED TO DETECTING FREQUENCIES C NO LARGER THAN 0.5 CYCLES PER DATA POINT; C THIS CORRESPONDS TO THE FACT THAT THE C SMALLEST DETECTABLE CYCLE IN THE DATA C IS 2 DATA POINTS PER CYCLE. C REFERENCES--JENKINS AND WATTS, ESPECIALLY PAGE 290. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-2315 C ORIGINAL VERSION--JUNE 1972. C UPDATED --NOVEMBER 1975. C UPDATED --FEBRUARY 1977. C C--------------------------------------------------------------------- C DIMENSION X(1) DIMENSION R(500) DIMENSION S(125) DIMENSION PSSQ(6),PMSQ(6),PS(6),P(5),L(4) DATA PI/3.14159265358979/ C IPR=6 ILOWER=3 MAXLAG=500 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.ILOWER)GOTO50 HOLD=X(1) DO65I=2,N IF(X(I).NE.HOLD)GOTO90 65 CONTINUE WRITE(IPR, 9)HOLD RETURN 50 WRITE(IPR,17)ILOWER WRITE(IPR,47)N RETURN 90 CONTINUE 9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT (A VECTOR) TO THE TIME SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 17 FORMAT(1H , 96H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 TIME SUBROUTINE IS OUTSIDE THE ALLOWABLE (,I6,11H,INFINITY) , 1 14HINTERVAL *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C-----START POINT----------------------------------------------------- C AN=N C C COMPUTE THE SAMPLE MEAN C SUM=0.0 DO100I=1,N SUM=SUM+X(I) 100 CONTINUE XBAR=SUM/AN C C COMPUTE THE SAMPLE VARIANCE AND THE SUM OF SQUARED DEVIATIONS C SUM=0.0 DO200I=1,N SUM=SUM+(X(I)-XBAR)*(X(I)-XBAR) 200 CONTINUE SSQ=SUM VARB=SSQ/AN VAR=SSQ/(AN-1.0) SD=SQRT(VAR) C C COMPUTE THE SAMPLE AUTOCORRELATIONS C REFERENCE--JENKINS AND WATTS, PAGES 290 AND 259 (7.1.6) C KMAX=N/4 IF(N.LE.32)KMAX=N/2 IF(N.LE.16)KMAX=N IF(KMAX.GT.MAXLAG)KMAX=MAXLAG DO300K=1,KMAX SUM=0.0 NMK=N-K DO400I=1,NMK J=I+K SUM=SUM+(X(I)-XBAR)*(X(J)-XBAR) 400 CONTINUE R(K)=SUM/SSQ 300 CONTINUE C C PLOT THE SAMPLE AUTOCORRELATIONS C CALL PLOTCO(R,KMAX) WRITE(IPR,1705) WRITE(IPR,1710)N,KMAX C C DO A WHITE NOISE ANALYSIS C SDR=1.0/SQRT(AN) R975=1.96*SDR IF(R975.GT.1.0)R975=1.0 R025=-R975 NUMOUT=0 DO410K=1,KMAX ABSR=R(K) IF(ABSR.LT.0.0)ABSR=-ABSR IF(ABSR.GT.R975)NUMOUT=NUMOUT+1 410 CONTINUE PEROUT=FLOAT(NUMOUT)/FLOAT(KMAX) PEROUT=100.0*PEROUT WRITE(IPR,999) WRITE(IPR,420)R025,R975 WRITE(IPR,430) WRITE(IPR,440)NUMOUT,KMAX,PEROUT DO415I=1,5 WRITE(IPR,999) 415 CONTINUE WRITE(IPR,1715)N WRITE(IPR,1720)XBAR WRITE(IPR,1725)VAR WRITE(IPR,1730)SD WRITE(IPR,1735)VARB C C COMPUTE THE PILOT SPECTRUM FOR THE REDUCED (2**J) SAMPLE C REFERENCE--JENKINS AND WATTS, PAGE 288 C DO450I=1,20 NDIV=N/(2**I) IF(NDIV.EQ.0)I2=I-1 IF(NDIV.EQ.0)GOTO460 450 CONTINUE 460 IF(7.LT.I2)I2=7 N2=2**I2 AN2=N2 DO500K=1,I2 SUM=0.0 IMIN=2**K JMAX=IMIN/2 DO600I=IMIN,N2,IMIN SUM1=0.0 SUM2=0.0 DO700J=1,JMAX IARG1=I+J-JMAX IARG2=IARG1-JMAX SUM1=SUM1+X(IARG1) SUM2=SUM2+X(IARG2) 700 CONTINUE SUM=SUM+(SUM1-SUM2)*(SUM1-SUM2) 600 CONTINUE PSSQ(K)=SUM/FLOAT(IMIN) PMSQ(K)=PSSQ(K)/AN2 PS(K)=FLOAT(2*IMIN)*PMSQ(K) PS(K)=PS(K)/VARB 500 CONTINUE C C FORM THE PILOT SPECTRUM PLOT C DO900I=1,I2 IREV=I2-I+1 JMIN=(120/(2**I))+1 IF(I.EQ.I2)JMIN=1 JMAX=120/(2**(I-1)) DO1000J=JMIN,JMAX S(J)=PS(I) 1000 CONTINUE 900 CONTINUE CALL PLOTSP(S,120,0) WRITE(IPR,999) WRITE(IPR,1750) C C DEFINE 4 LAG WINDOW TRUNCATION POINTS C REFERENCE--JENKINS AND WATTS, PAGES 290 AND 260 C P(1)=.2 P(2)=.1 P(3)=.05 P(4)=.025 P(5)=.01 LMAX=0 DO1100I=1,5 DO1200K=1,KMAX KREV=KMAX-K+1 RK=R(KREV) IF(RK.LT.0.0)RK=-RK IF(RK.GE.P(I))LMAX=KREV IF(RK.GE.P(I))GOTO1350 1200 CONTINUE 1100 CONTINUE IF(LMAX.NE.0)GOTO1350 RMAX=ABS(R(1)) DO1300K=1,KMAX RK=R(K) IF(RK.LT.0.0)RK=-RK IF(RK.GE.RMAX)LMAX=K IF(RK.GE.RMAX)RMAX=RK 1300 CONTINUE 1350 ALMAX=LMAX L(1)=(3.0/2.0)*ALMAX IF(L(1).LE.32)LMAX=32 IF(L(1).LE.32)ALMAX=32.0 IF(L(1).LE.32)L(1)=32 IF(L(1).GE.KMAX)LMAX=KMAX IF(L(1).GE.KMAX)ALMAX=KMAX IF(L(1).GE.KMAX)L(1)=KMAX L(2)=(ALMAX/2.0)+0.1 L(3)=(ALMAX/4.0)+0.1 L(4)=(ALMAX/8.0)+0.1 IF(L(4).GE.3)NUMSP=4 IF(L(4).GE.3)GOTO1380 IF(L(3).GE.3)NUMSP=3 IF(L(3).GE.3)GOTO1380 IF(L(2).GE.3)NUMSP=2 IF(L(2).GE.3)GOTO1380 IF(L(1).GE.3)NUMSP=1 IF(L(1).GE.3)GOTO1380 WRITE(IPR,1390)N RETURN 1380 CONTINUE C C COMPUTE THE 4 SPECTRUM ESTIMATES C REFERENCE--JENKINS AND WATTS, PAGES 260 AND 244 C C COMPUTE BANDWIDTHS C REFERENCE--JENKINS AND WATTS, PAGES 257 AND 252 C C COMPUTE DEGREES OF FREEDOM FOR THE SPECTAL DENSITY ESTIMATE AT INDIVIDUAL C FREQUENCIES C REFERENCE--JENKINS AND WATTS, PAGES 254 AND 252 C C COMPUTE 95 PERCENT CONFIDENCE INTERVAL LENGTHS FOR THE LOG SPECTRAL C DENSITY ESTIMATES C REFERENCE--JENKINS AND WATTS, PAGES 255 AND 252 C C WRITE OUT THE 4 SPECTRUM PLOTS C DO1400I=1,NUMSP AL=L(I) LM1=L(I)-1 DO1500LLP1=1,121 LL=LLP1-1 ALL=LL SUM=0.0 DO1600K=1,LM1 AK=K ARG1=PI*AK/AL ARG2=PI*AK*ALL/120.0 WK=0.0 IF(K.LE.L(I))WK=0.5*(1.0+COS(ARG1)) SUM=SUM+R(K)*WK*COS(ARG2) 1600 CONTINUE SUM=2.0+4.0*SUM S(LLP1)=SUM 1500 CONTINUE BW=(4.0/3.0)/FLOAT(L(I)) DF=(8.0/3.0)*AN/FLOAT(L(I)) IDF=DF+0.5 CALL PLOTSP(S,121,IDF) DFROUN=IDF WRITE(IPR,999) WRITE(IPR,1760)L(I),BW,DFROUN 1400 CONTINUE C 420 FORMAT(1H ,103HUNDER THE NULL HYPOTHESIS OF WHITE NOISE (AND NORMA 1LITY), A 2-SIDED 95 PERCENT ACCEPTANCE INTERVAL IS (,F6.4,1H,,F6.4 1,1H)) 430 FORMAT(1H ,125HUNDER THE NULL HYPOTHESIS, ONLY 5 PERCENT (ON THE A 1VERAGE) OF THE OBSERVED AUTOCORRELATIONS SHOULD FALL OUTSIDE THIS 1INTERVAL) 440 FORMAT(1H ,20HIT IS OBSERVED THAT ,I5,12H OUT OF THE ,I5,11H (THAT 1 IS, ,F5.1,72H PERCENT) OF THE COMPUTED AUTOCORRELATIONS FALL OUTS 1IDE OF THIS INTERVAL) 1390 FORMAT(1H ,45HDUE TO THE SMALL NUMBER OF OBSERVATIONS (N = ,I6,59H 1), THERE ARE NOT ENOUGH LAGS TO PRODUCE A RELIABLE SPECTRUM) 1705 FORMAT(1H ,30X, 63HAUTOCORRELATION PLOT--PLOT OF SAMPLE AUTOCORREL 1ATION VERSUS LAG) 1710 FORMAT(1H ,10X,29HTHE NUMBER OF OBSERVATIONS = ,I6,10X,56HTHE NUMB 1ER OF COMPUTED (AND PLOTTED) AUTOCORRELATIONS = ,I6) 1715 FORMAT(1H ,18HTHE SAMPLE SIZE = ,I6) 1720 FORMAT(1H ,18HTHE SAMPLE MEAN = ,E15.8) 1725 FORMAT(1H ,22HTHE SAMPLE VARIANCE = ,E15.8) 1730 FORMAT(1H ,32HTHE SAMPLE STANDARD DEVIATION = ,E15.8) 1735 FORMAT(1H ,29HTHE BIASED SAMPLE VARIANCE = ,E15.8) 1750 FORMAT(1H ,50X,14HPILOT SPECTRUM) 1760 FORMAT(1H ,17HNUMBER OF LAGS = ,I5,10X,11HBANDWIDTH =,F10.3,10X,21 1HDEGREES OF FREEDOM = ,F10.3) 999 FORMAT(1H ) C RETURN END