time(3f) - [M_datapac:ANALYSIS] perform a time series analysis
(autocorrelation plot, a test for white noise, a "pilot" spectrum,
and 4 other estimated spectra based on differing bandwidth)
SUBROUTINE TIME(X,N)
time(3f) performs a time series analysis on the data in the input
vector x.
the analysis consists of the following--
1. a plot of autocorrelation versus lag number;
2. a test for white noise (assuming normality);
3. a 'pilot' spectrum; and
4. 4 other estimated spectra--each based
on a differing bandwidth.
in order that the results of the time series analysis be valid and
properly interpreted, the input data in x should be equi-spaced in time
(or whatever variable corresponds to time).
the horizontal axis of the spectra produced by time(3f) is frequency.
this frequency is measured in units of cycles per 'data point' or,
more precisely, in cycles per unit time where 'unit time' is defined
as the elapsed time between adjacent observations.
the range of the frequency axis is 0.0 to 0.5.
X description of parameter
Y description of parameter
Sample program:
program demo_time
use M_datapac, only : time
implicit none
! call time(x,y)
end program demo_time
Results:
The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.
John Urban, 2022.05.31
CC0-1.0
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | dimension(:) | :: | X | |||
| integer | :: | N |
SUBROUTINE TIME(X,N) REAL(kind=wp) :: absr, ak, al, all, almax, an, an2, arg1, arg2, bw, df, dfroun, hold, p, perout, pi, pmsq, ps, pssq, r REAL(kind=wp) :: r025, r975, rk, rmax, s, sd, sdr, ssq, sum, sum1, sum2, var, varb, wk, X, xbar INTEGER :: i, i2, iarg1, iarg2, idf, ilower, imin, irev, j, jmax, jmin, k, kmax, krev, l, ll, llp1, lm1, lmax INTEGER :: maxlag, N, n2, ndiv, nmk, numout, numsp ! ! INPUT ARGUMENTS--X = THE VECTOR OF ! (UNSORTED) OBSERVATIONS. ! N = THE INTEGER NUMBER OF OBSERVATIONS ! IN THE VECTOR X. ! OUTPUT--7 TO 11 pages (DEPENDING ON ! THE INPUT SAMPLE SIZE) OF ! AUTOMATIC PRINTOUT-- ! 1) A PLOT OF AUTOCORRELATION VERSUS LAG NUMBER; ! THIS PLOT MAY TAKE AS LITTLE AS 1 ! OR AS MANY AS 5 pages ! (THE EXACT NUMBER DEPENDING ON ! THE INPUT SAMPLE SIZE N); ! 2) A TEST FOR WHITE NOISE (ASSUMING NORMALITY); ! 3) A 'PILOT' SPECTRUM; AND ! 4) AN ESTIMATED SPECTRUM BASED ON A ! BANDWIDTH DERIVED FROM THE DATA SET; ! 5) AN ESTIMATED SPECTRUM BASED ON A ! BANDWIDTH ONLY 1/2 AS WIDE AS IN 4; ! 6) AN ESTIMATED SPECTRUM BASED ON A ! BANDWIDTH ONLY 1/4 AS WIDE AS IN 4; ! 7) AN ESTIMATED SPECTRUM BASED ON A ! BANDWIDTH ONLY 1/8 AS WIDE AS IN 4; ! PRINTING--YES. ! RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE ! OF N FOR THIS SUBROUTINE. ! --THE SAMPLE SIZE N MUST BE GREATER ! THAN OR EQUAL TO 3. ! COMMENT--THE 'FAST FOURIER TRANSFORM' IS NOT USED ! IN THIS VERSION OF TIME, BUT WILL BE ! IMPLEMENTED IN A FUTURE VERSION. ! --THE USUAL MAXIMUM NUMBER OF LAGS ! FOR WHICH THE AUTOCORRELATION IS ! COMPUTED IS N/4 WHERE N IS ! THE SAMPLE SIZE (LENGTH OF THE ! DATA RECORD IN THE VECTOR X). ! THIS RULE IS OVERRIDDEN IN ! LARGE DATA SETS AND IS REPLACED ! BY THE RULE THAT THE MAXIMUM ! NUMBER OF LAGS = 500 ! (WHICH CORRESPONDS TO AN ! AUTOCORRELATION PLOT COVERING ! 5 COMPUTER pages). ! IF MORE LAGS ARE DESIRED, ! CHANGE THE VALUE OF THE ! VARIABLE MAXLAG ! WITHIN THIS SUBROUTINE ! FROM 500 TO WHATEVER DESIRED, ! AND ALSO CHANGE THE DIMENSION OF ! THE VECTOR R FROM ITS PRESENT 500 TO HOWEVER ! MANY LAGS ARE DESIRED. ! --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED ! TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME, ! THEN THE FREQUENCY AXIS OF THE RESULTING ! SPECTRA WOULD BE IN UNITS OF HERTZ ! (= CYCLES PER SECOND). ! --THE FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE ! IN THE DATA OF INFINITE (= 1/(0.0)) ! LENGTH OR PERIOD. ! THE FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE ! IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS. ! --ANY EQUI-SPACED TIME SERIES ANALYSIS IS ! INTRINSICALLY LIMITED TO DETECTING FREQUENCIES ! NO LARGER THAN 0.5 CYCLES PER DATA POINT; ! THIS CORRESPONDS TO THE FACT THAT THE ! SMALLEST DETECTABLE CYCLE IN THE DATA ! IS 2 DATA POINTS PER CYCLE. ! !--------------------------------------------------------------------- ! DIMENSION X(:) DIMENSION r(500) DIMENSION s(125) DIMENSION pssq(6) , pmsq(6) , ps(6) , p(5) , l(4) DATA pi/3.14159265358979_wp/ ! ilower = 3 maxlag = 500 ! ! CHECK THE INPUT ARGUMENTS FOR ERRORS ! IF ( N<ilower ) THEN WRITE (G_IO,99001) ilower 99001 FORMAT (' ', & &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE TIME SUBROU& &TINE IS OUTSIDE THE ALLOWABLE (',I0,',INFINITY) ', & &'INTERVAL *****') WRITE (G_IO,99002) N 99002 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****') RETURN ELSE hold = X(1) DO i = 2 , N IF ( X(i)/=hold ) GOTO 100 ENDDO WRITE (G_IO,99003) hold 99003 FORMAT (' ', & &'***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUMENT (A VECTOR) & &TO THE TIME SUBROUTINE HAS ALL ELEMENTS = ',E15.8,' *****') RETURN ENDIF ! !-----START POINT----------------------------------------------------- ! 100 an = N ! ! COMPUTE THE SAMPLE MEAN ! sum = 0.0_wp DO i = 1 , N sum = sum + X(i) ENDDO xbar = sum/an ! ! COMPUTE THE SAMPLE VARIANCE AND THE SUM OF SQUARED DEVIATIONS ! sum = 0.0_wp DO i = 1 , N sum = sum + (X(i)-xbar)*(X(i)-xbar) ENDDO ssq = sum varb = ssq/an var = ssq/(an-1.0_wp) sd = SQRT(var) ! ! COMPUTE THE SAMPLE AUTOCORRELATIONS ! REFERENCE--JENKINS AND WATTS, pages 290 AND 259 (7.1.6) ! kmax = N/4 IF ( N<=32 ) kmax = N/2 IF ( N<=16 ) kmax = N IF ( kmax>maxlag ) kmax = maxlag DO k = 1 , kmax sum = 0.0_wp nmk = N - k DO i = 1 , nmk j = i + k sum = sum + (X(i)-xbar)*(X(j)-xbar) ENDDO r(k) = sum/ssq ENDDO ! ! PLOT THE SAMPLE AUTOCORRELATIONS ! CALL PLOTCO(r,kmax) WRITE (G_IO,99004) 99004 FORMAT (' ',30X, & & 'AUTOCORRELATION PLOT--PLOT OF SAMPLE AUTOCORRELATION VERSUS LAG'& & ) WRITE (G_IO,99005) N , kmax 99005 FORMAT (' ',10X,'THE NUMBER OF OBSERVATIONS = ',I0,10X, & & 'THE NUMBER OF COMPUTED (AND PLOTTED) AUTOCORRELATIONS = '& & ,I0) ! ! DO A WHITE NOISE ANALYSIS ! sdr = 1.0_wp/SQRT(an) r975 = 1.96_wp*sdr IF ( r975>1.0_wp ) r975 = 1.0_wp r025 = -r975 numout = 0 DO k = 1 , kmax absr = r(k) IF ( absr<0.0_wp ) absr = -absr IF ( absr>r975 ) numout = numout + 1 ENDDO perout = FLOAT(numout)/FLOAT(kmax) perout = 100.0_wp*perout WRITE (G_IO,99017) WRITE (G_IO,99006) r025 , r975 ! 99006 FORMAT (' ', & &'UNDER THE NULL HYPOTHESIS OF WHITE NOISE (AND NORMALITY), A 2-SID& &ED 95 PERCENT ACCEPTANCE INTERVAL IS (',F6.4,',',F6.4,')') WRITE (G_IO,99007) 99007 FORMAT (' ', & &'UNDER THE NULL HYPOTHESIS, ONLY 5 PERCENT (ON THE AVERAGE) OF THE& & OBSERVED AUTOCORRELATIONS SHOULD FALL OUTSIDE THIS INTERVAL') WRITE (G_IO,99008) numout , kmax , perout 99008 FORMAT (' ','IT IS OBSERVED THAT ',I5,' OUT OF THE ',I5, & & ' (THAT IS, ',F5.1, & &' PERCENT) OF THE COMPUTED AUTOCORRELATIONS FALL OUTSIDE OF THIS I& &NTERVAL') DO i = 1 , 5 WRITE (G_IO,99017) ENDDO WRITE (G_IO,99009) N 99009 FORMAT (' ','THE SAMPLE SIZE = ',I0) WRITE (G_IO,99010) xbar 99010 FORMAT (' ','THE SAMPLE MEAN = ',E15.8) WRITE (G_IO,99011) var 99011 FORMAT (' ','THE SAMPLE VARIANCE = ',E15.8) WRITE (G_IO,99012) sd 99012 FORMAT (' ','THE SAMPLE STANDARD DEVIATION = ',E15.8) WRITE (G_IO,99013) varb 99013 FORMAT (' ','THE BIASED SAMPLE VARIANCE = ',E15.8) ! ! COMPUTE THE PILOT SPECTRUM FOR THE REDUCED (2**J) SAMPLE ! REFERENCE--JENKINS AND WATTS, page 288 ! DO i = 1 , 20 ndiv = N/(2**i) IF ( ndiv==0 ) i2 = i - 1 IF ( ndiv==0 ) EXIT ENDDO IF ( 7<i2 ) i2 = 7 n2 = 2**i2 an2 = n2 DO k = 1 , i2 sum = 0.0_wp imin = 2**k jmax = imin/2 DO i = imin , n2 , imin sum1 = 0.0_wp sum2 = 0.0_wp DO j = 1 , jmax iarg1 = i + j - jmax iarg2 = iarg1 - jmax sum1 = sum1 + X(iarg1) sum2 = sum2 + X(iarg2) ENDDO sum = sum + (sum1-sum2)*(sum1-sum2) ENDDO pssq(k) = sum/FLOAT(imin) pmsq(k) = pssq(k)/an2 ps(k) = FLOAT(2*imin)*pmsq(k) ps(k) = ps(k)/varb ENDDO ! ! FORM THE PILOT SPECTRUM PLOT ! DO i = 1 , i2 irev = i2 - i + 1 jmin = (120/(2**i)) + 1 IF ( i==i2 ) jmin = 1 jmax = 120/(2**(i-1)) DO j = jmin , jmax s(j) = ps(i) ENDDO ENDDO CALL PLOTSP(s,120,0) WRITE (G_IO,99017) WRITE (G_IO,99014) 99014 FORMAT (' ',50X,'PILOT SPECTRUM') ! ! DEFINE 4 LAG WINDOW TRUNCATION POINTS ! REFERENCE--JENKINS AND WATTS, pages 290 AND 260 ! p(1) = 0.2_wp p(2) = 0.1_wp p(3) = 0.05_wp p(4) = 0.025_wp p(5) = 0.01_wp lmax = 0 DO i = 1 , 5 DO k = 1 , kmax krev = kmax - k + 1 rk = r(krev) IF ( rk<0.0_wp ) rk = -rk IF ( rk>=p(i) ) lmax = krev IF ( rk>=p(i) ) GOTO 200 ENDDO ENDDO IF ( lmax==0 ) THEN rmax = ABS(r(1)) DO k = 1 , kmax rk = r(k) IF ( rk<0.0_wp ) rk = -rk IF ( rk>=rmax ) lmax = k IF ( rk>=rmax ) rmax = rk ENDDO ENDIF 200 almax = lmax l(1) = (3.0_wp/2.0_wp)*almax IF ( l(1)<=32 ) lmax = 32 IF ( l(1)<=32 ) almax = 32.0_wp IF ( l(1)<=32 ) l(1) = 32 IF ( l(1)>=kmax ) lmax = kmax IF ( l(1)>=kmax ) almax = kmax IF ( l(1)>=kmax ) l(1) = kmax l(2) = (almax/2.0_wp) + 0.1_wp l(3) = (almax/4.0_wp) + 0.1_wp l(4) = (almax/8.0_wp) + 0.1_wp IF ( l(4)>=3 ) numsp = 4 IF ( l(4)<3 ) THEN IF ( l(3)>=3 ) numsp = 3 IF ( l(3)<3 ) THEN IF ( l(2)>=3 ) numsp = 2 IF ( l(2)<3 ) THEN IF ( l(1)>=3 ) numsp = 1 IF ( l(1)<3 ) THEN WRITE (G_IO,99015) N 99015 FORMAT (' ', & & 'DUE TO THE SMALL NUMBER OF OBSERVATIONS (N = '& & ,I0, & & '), THERE ARE NOT ENOUGH LAGS TO PRODUCE A RELIABLE SPECTRUM'& & ) RETURN ENDIF ENDIF ENDIF ENDIF ! ! COMPUTE THE 4 SPECTRUM ESTIMATES ! REFERENCE--JENKINS AND WATTS, pages 260 AND 244 ! ! COMPUTE BANDWIDTHS ! REFERENCE--JENKINS AND WATTS, pages 257 AND 252 ! ! COMPUTE DEGREES OF FREEDOM FOR THE SPECTAL DENSITY ESTIMATE AT INDIVIDUAL ! FREQUENCIES ! REFERENCE--JENKINS AND WATTS, pages 254 AND 252 ! ! COMPUTE 95 PERCENT CONFIDENCE INTERVAL LENGTHS FOR THE LOG SPECTRAL ! DENSITY ESTIMATES ! REFERENCE--JENKINS AND WATTS, pages 255 AND 252 ! ! WRITE OUT THE 4 SPECTRUM PLOTS ! DO i = 1 , numsp al = l(i) lm1 = l(i) - 1 DO llp1 = 1 , 121 ll = llp1 - 1 all = ll sum = 0.0_wp DO k = 1 , lm1 ak = k arg1 = pi*ak/al arg2 = pi*ak*all/120.0_wp wk = 0.0_wp IF ( k<=l(i) ) wk = 0.5_wp*(1.0_wp+COS(arg1)) sum = sum + r(k)*wk*COS(arg2) ENDDO sum = 2.0_wp + 4.0_wp*sum s(llp1) = sum ENDDO bw = (4.0_wp/3.0_wp)/FLOAT(l(i)) df = (8.0_wp/3.0_wp)*an/FLOAT(l(i)) idf = df + 0.5 CALL PLOTSP(s,121,idf) dfroun = idf WRITE (G_IO,99017) WRITE (G_IO,99016) l(i) , bw , dfroun 99016 FORMAT (' NUMBER OF LAGS = ',I5,10X,'BANDWIDTH =',F10.3,10X,'DEGREES OF FREEDOM = ',F10.3) ENDDO 99017 FORMAT (' ') ! END SUBROUTINE TIME