C PROGRAM 7.1 ARFIT C C ... AR model fitting ... C C Inputs: C IDEV: Input device for time series C LAG: Highest order od AR model C ISW: Estimation procedure C = 1: Yule-Walker method C = 2: Least squares (Householder) method C = 3: Partial autoregression method C = 4: PARCOR method C = 5: Burg's algorithm (MEM) 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 MJ,MJ2: Adjustable dimensions C NF: Number of frequencies for computing spectrum C @TEST.PN71: 12/12/90,1/7/91 Y.I. and G.K. 9/4/91 C PARAMETER( NMAX=200,MJ=20,MJ2=100,NF=200,IDEV=1 ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), COV(0:MJ,4), SIG2(0:MJ), SP(0:NF) DIMENSION D(MJ2), X(MJ2,MJ+1) DIMENSION AIC(0:MJ), A(MJ,MJ), B(MJ), PAR(MJ) DIMENSION FE(NMAX), BE(NMAX) EXTERNAL SETXAR DATA OUTMIN/-1.0D30/, OUTMAX/1.0D30/ C READ(5,*) LAG, ISW CALL READTS( IDEV,Y,N ) CALL MEAN( Y,N,-1.0D30,1.0D30,NSUM,YMEAN ) C C ... Yule-Walker method ... C IF( ISW.EQ.1 ) THEN CALL UNICOR( Y,N,LAG,OUTMIN,OUTMAX,COV ) CALL ARYULE( COV,N,LAG,SIG2,AIC,PAR,A,MAR) END IF C C ... Householder method ... C IF( ISW.EQ.2 ) THEN CALL REDUCT( SETXAR,Y,D,N-LAG,0,LAG,MJ2,X ) CALL REGRES( X,LAG,N-LAG,MJ2,MJ,A,SIG2,AIC,MAR ) CALL PARCOR( A(1,MAR),MAR,PAR ) END IF C C ... PARCOR method ... C IF( ISW.GE.3 ) THEN CALL ARPCOR( Y,FE,BE,SIG2,AIC,LAG,N,PAR,ISW-2,MAR ) DO 10 I=1,LAG 10 CALL ARCOEF( PAR,I,A(1,I) ) END IF C C ... Power spectrum ... C CALL ARMASP( A(1,MAR),MAR,B,0,SIG2(MAR),NF,SP ) C CALL PRARSP( N,LAG,A,MAR,PAR,SIG2,AIC,SP,MJ,NF,ISW ) C CALL PTAR( PAR,AIC,SP,LAG,NF,MAR,SIG2,ISW ) STOP E N D SUBROUTINE ARYULE( C,N,MAXM,SIG2,AIC,PARCOR,A,MAR ) C C ... Yule-Walker method ... C C Inputs: C C(I): Autocovariance function C N: Data length C MAXM: Highest AR order C Outputs: C SIG2(I): Innovation variance C AIC(I): AIC C PARCOR(I): PARCOR C AMIN: AR coefficients of the best model C MAR: Selected order of the model C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION C(0:MAXM), SIG2(0:MAXM), AIC(0:MAXM) DIMENSION PARCOR(MAXM), A(MAXM,MAXM) CONST = N*(DLOG(2*3.1415926535D0) + 1) C SIG2(0) = C(0) AIC(0) = CONST + N*DLOG(SIG2(0)) + 2 AICMIN = AIC(0) MAR = 0 C DO 50 M=1,MAXM SUM = C(M) DO 10 I=1,M-1 10 SUM = SUM - A(I,M-1)*C(M-I) A(M,M) = SUM /SIG2(M-1) DO 20 J=1,M-1 20 A(J,M) = A(J,M-1)-A(M,M)*A(M-J,M-1) SIG2(M) = SIG2(M-1)*(1.0D0-A(M,M)**2) AIC(M) = CONST + N*DLOG(SIG2(M)) + 2*(M+1) PARCOR(M) = A(M,M) IF( AIC(M).LT.AICMIN ) THEN AICMIN = AIC(M) MAR = M END IF 50 CONTINUE RETURN E N D SUBROUTINE SETXAR( Z,N0,L,K,MJ1,JSW,X ) C C ... Data matrix for AR model ... C C Inputs: C Z(I): Time series C N0: Origin of the current observations C L: Number of current observations C K: Number of regressors C MJ1: Adjustable dimension of X C JSW=0: Make initial data matrix C =1: Apend L*(K+1) data matrix below the triangular one C Output: C X(I,J): Data matrix C REAL*8 X(MJ1,1), Z(1) C I0 = 0 IF( JSW .EQ. 1 ) I0 = K+1 DO 10 I=1,L II = I + I0 JJ = N0 + K + I X(II,K+1) = Z(JJ) DO 10 J=1,K JJ = JJ - 1 10 X(II,J) = Z(JJ) C RETURN E N D SUBROUTINE ARPCOR( Y,FE,BE,SIG2,AIC,K,N,PARCOR,ISW,MAR ) C C ... PARCOR method ... C C Inputs: C Y(I): Time series C N: Data length C K: Highest AR order C ISW: Estimation procedure C = 1: Partial autoregression C = 2: PARCOR C = 3: Burg's algorithm (MEM) C FE,BE: Working area C Outputs: C SIG2(I): Innovation variance C AIC(I): AIC C PARCOR(I): PARCOR C MAR: MAICE order C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N), FE(N), BE(N) DIMENSION SIG2(0:K), AIC(0:K), PARCOR(K) DIMENSION A(50), B(50), FA(50), BA(50) PI = 3.1415926535D0 C SUM = 0.0D0 DO 5 I=K+1,N 5 SUM = SUM + Y(I)**2 SIG2(0) = SUM/(N-K) AIC(0) = (N-K)*( DLOG(2*PI) + 1+ DLOG(SIG2(0))) + 2 MAR = 0 AICM = AIC(0) C DO 10 I=1,N FE(I) = Y(I) 10 BE(I) = Y(I) C DO 100 M=1,K FB = 0.0D0 FF = 0.0D0 BB = 0.0D0 DO 20 I=M+1,N FB = FB + FE(I)*BE(I-M) FF = FF + BE(I-M)**2 20 BB = BB + FE(I)**2 C IF( ISW.EQ.1) THEN A(M) = FB/FF B(M) = FB/BB END IF IF( ISW.EQ.2 ) THEN A(M) = FB/(DSQRT(FF*BB)) B(M) = FB/(DSQRT(FF*BB)) END IF IF( ISW.EQ.3 ) THEN A(M) = FB/((FF+BB)/2) B(M) = FB/((FF+BB)/2) END IF C DO 50 L=1,M-1 A(L) = FA(L)-A(M)*BA(M-L) 50 B(L) = BA(L)-B(M)*FA(M-L) DO 60 L=1,M FA(L) = A(L) 60 BA(L) = B(L) DO 70 I=M+1,N FE0 = FE(I) FE(I) = FE(I)-A(M)*BE(I-M) 70 BE(I-M) = BE(I-M)-B(M)*FE0 C SUM = 0.0D0 DO 80 I=K+1,N 80 SUM = SUM + FE(I)**2 X = SUM PARCOR(M) = A(M) SIG2(M) = X/(N-K) AIC(M) = (N-K)*( DLOG(2*PI) + 1+ DLOG(SIG2(M)))+ 2*(M+1) IF( AIC(M).LT.AICM ) THEN MAR = M AICM = AIC(M) END IF C 100 CONTINUE C RETURN E N D SUBROUTINE PRARSP( N,LAG,A,MAR,PARCOR,SIG2,AIC,SP,MJ,MJ2,ISW ) C C ... Print OUT SIG2, PARCOR, AIC and Power spectrum ... C C Inputs: C N: Data length C LAG: Highest AR order C A: AR coefficients C MAR: Selected AR order C PARCOR(I): PARCOR C SIG2(I): Innovation variance C AIC(I): AIC C SP(I): Power spectrum C MJ: Adjustable dimension C MJ2: Number of frequencies C ISW: Specification of estimation procedure C IMPLICIT REAL*8(A-H,O-Z) CHARACTER*72 TITLE DIMENSION A(MJ,MJ), PARCOR(MJ), SIG2(0:MJ), AIC(0:MJ) DIMENSION SP(0:MJ2) COMMON /CMDATA/ TITLE C WRITE(6,600) WRITE(6,610) TITLE WRITE(6,620) N, LAG IF( ISW.EQ.1 ) WRITE(6,630) IF( ISW.EQ.2 ) WRITE(6,640) IF( ISW.GE.3 ) WRITE(6,650) ISW WRITE(6,660) DO 10 I=0,LAG 10 WRITE(6,670) I, SIG2(I), AIC(I) WRITE(6,680) MAR, SIG2(MAR), AIC(MAR) WRITE(6,690) (A(I,MAR),I=1,MAR) WRITE(6,700) (PARCOR(I),I=1,MAR) IF(MJ2.GT.0) WRITE(6,710) (SP(I),I=0,MJ2) C 600 FORMAT( 1H ,'PROGRAM 7.1: AR MODEL FITTING' ) 610 FORMAT( 1H ,' DATA: ',A72 ) 620 FORMAT( 1H ,' N =',I5,3X,'MAX. ORDER =',I3 ) 630 FORMAT( 1H ,'**** YULE WALKER METHOD ****' ) 640 FORMAT( 1H ,'**** HOUSEHOLDER METHOD ****' ) 650 FORMAT( 1H ,'**** PARCOR METHOD (ISW =',I1,') ****' ) 660 FORMAT( 1H ,' M',10X,'SIG2',10X,'AIC' ) 670 FORMAT( 1H ,I3,D17.6,F12.3 ) 680 FORMAT( 1H ,' MAICE ORDER =',I2,3X,'SIG2 =',D13.5,3X,'AIC =', * F12.4 ) 690 FORMAT( 1H0,'** AR COEFFICIENTS: A(1),...,A(MAR) **',/ * (1H ,5F10.6) ) 700 FORMAT( 1H0,'** PARCOR: PAR(1),...,PAR(LAG) **',/ * (1H ,5F10.6) ) 710 FORMAT( 1H0,'** POWER SPECTRUM (IN LOG SCALE: SP(0),...,SP(MJ2) ***',/ (1H ,10F8.4) ) 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 SUBROUTINE REGRES( X,K,N,MJ1,MJ2,A,SIG2,AIC,IMIN ) C C ... Regression model fitting ... C ... Order of the model is selected by AIC ... C C Inputs: C X: Householder reduced form (upper triangular matrix) C K: Maximum number of regressors C N: Number of data C MJ1: Adjustable dimension of X C MJ2: Adjustable dimension of A, ISG2 and AIC C Outputs: C A(I,M): Regression coefficients of the model with order M C SIG2: Residual variances C AIC: AIC's C IMIN: MAICE order C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,1), A(MJ2,MJ2), SIG2(0:MJ2), AIC(0:MJ2) C CALL COMAIC( X,N,K,MJ1,SIG2,AIC ) C IMIN = 0 AICM = AIC(0) DO 10 M=1,K IF( AIC(M).LT.AICM ) THEN IMIN = M AICM = AIC(M) END IF CALL RECOEF( X,M,K,MJ1,A(1,M) ) 10 CONTINUE C RETURN E N D SUBROUTINE REDUCT( SETX,Z,D,NMK,N0,K,MJ1,X ) C C ... Successive Householder reduction ... C C Inputs: C SETX: Name of the subroutine for making X(I,J) C Z(I): Data vector C D(I): Working area C NMK: Number of actually used observations C N0: Time point of the previous set ofobservations C K: Heighest order of the model C MJ1: Adjustable dimension of X C Output: C X(I,J): data matrix C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION X(MJ1,1) , D(1), Z(1) C L = MIN0( NMK,MJ1 ) K1 = K + 1 N1 = L C CALL SETX( Z,N0,L,K,MJ1,0,X ) CALL HUSHLD( X,D,MJ1,L,K1 ) IF( N1 .GE. NMK ) RETURN C 10 L = MIN0( NMK-N1,MJ1-K1 ) C LK = L + K1 N2 = N0 + N1 CALL SETX( Z,N2,L,K,MJ1,1,X ) CALL HUSHLD( X,D,MJ1,LK,K1 ) N1 = N1 + L IF( N1.LT.NMK ) GO TO 10 C RETURN C E N D SUBROUTINE COMAIC( X,N,K,MJ1,SIG2,AIC ) C C ... This subroutine computes residual variances and AIC's ... C C Inputs: C X(I,J): Householder reduced form C N: Data length C K: Heighest order C MJ1: Adjustable dimension of X C Outputs: C SIG2(I): residual variance of the model with order I C AIC(I): AIC of the model with order I C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,1), AIC(0:K), SIG2(0:K) DATA PI2/6.28318531D0/ C PVAR = 0.0D0 DO 10 I=K,0,-1 PVAR = PVAR + X(I+1,K+1)**2 SIG2(I) = PVAR / N 10 AIC(I) = N*DLOG( PI2*SIG2(I) ) + N + 2*(I+1) C RETURN E N D SUBROUTINE HUSHLD( X,D,MJ1,N,K ) C C ... Householder transformation ... C C Inputs: C X(I,J): Original matrix C D(I): Working area C MJ1: Adjustable dimension of X C N: Number of rows of X C K: Number of columns of X C Output: C X(I,J): Householder reduced form (upper triangular form) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,1), D(MJ1) TOL = 1.0D-60 C DO 100 II=1,K H = 0.0D0 DO 10 I=II,N D(I) = X(I,II) 10 H = H + D(I)**2 IF( H .GT. TOL ) GO TO 20 G = 0.0D0 GO TO 100 20 G = DSQRT( H ) F = X(II,II) IF( F .GE. 0.0D0 ) G = -G D(II) = F - G H = H - F*G C DO 30 I=II+1,N 30 X(I,II) = 0.0D0 DO 60 J=II+1,K S = 0.0D0 DO 40 I=II,N 40 S = S + D(I)*X(I,J) S = S/H DO 50 I=II,N 50 X(I,J) = X(I,J) - D(I)*S 60 CONTINUE 100 X(II,II) = G C RETURN E N D SUBROUTINE RECOEF( X,M,K,MJ,A ) C C ... Regression coefficients ... C C Inputs: C X(I,J): Householder reduced form C M: Number of actually used regressors C K: Heighest order C MJ: Adjustable dimension of X C Output: C A(I): Vector of regression coefficients C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(MJ,1), A(1) C A(M) = X(M,K+1)/X(M,M) DO 20 I=M-1,1,-1 SUM = X(I,K+1) DO 10 J=I+1,M 10 SUM = SUM - A(J)*X(I,J) 20 A(I) = SUM/X(I,I) C RETURN E N D SUBROUTINE ARCOEF( PAR,K,A ) C C ... This subroutine computes AR coefficients from PARCOR ... C C Inputs: C PAR(I): PARCOR C K: Order of AR model C Output: C A(I): AR coefficient C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PAR(K), A(K), AA(50) C DO 30 II=1,K A(II) = PAR(II) AA(II) = PAR(II) IF( II-1.LE.0 ) GO TO 30 DO 10 J=1,II-1 10 A(J) = AA(J) - PAR(II)*AA(II-J) IF( II.LT.K ) THEN DO 20 J=1,II-1 20 AA(J) = A(J) END IF 30 CONTINUE RETURN E N D SUBROUTINE PARCOR( A,K,PAR ) C C ... This subriutine computes PARCOR for AR coefficients ... C C Inputs: C A: Vector of AR coefficients C K: Order of the model C C Output: C PAR: Vector of partial autocorrelations (PARCOR) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(K), PAR(K), G(50) C DO 10 I=1,K 10 PAR(I) = A(I) C DO 40 II=K-1,1,-1 S = 1.D0 - PAR(II+1)**2 DO 20 I=1,II 20 G(I) = (PAR(I) + PAR(II+1)*PAR(II-I+1))/S I2 = (II+1)/2 IF( MOD( II,2 ).EQ.1 ) G(I2) = PAR(I2)/(1.D0 - PAR(II+1)) DO 30 I=1,II 30 PAR(I) = G(I) 40 CONTINUE C RETURN E N D SUBROUTINE ARMASP( A,M,B,L,SIG2,NF,SP ) C C ... Logarithm of the power spectrum of the ARMA model ... C C Inputs: C M: AR order C L: MA order C A(I): AR coefficient C B(I): MA coefficient C SIG2: Innovation variance C NF: Number of frequencies C Output: C SP(I): Power spectrum (in log scale) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M), B(L) DIMENSION SP(0:NF), H(0:500), FR(0:500), FI(0:500) C H(0) = 1.0D0 DO 10 I=1,M 10 H(I) = -A(I) C CALL FOURIE( H,M+1,NF+1,FR,FI ) C DO 20 I=0,NF 20 SP(I) = SIG2/( FR(I)**2 + FI(I)**2 ) C H(0) = 1.0D0 DO 30 I=1,L 30 H(I) = -B(I) CALL FOURIE( H,L+1,NF+1,FR,FI ) DO 40 I=0,NF 40 SP(I) = SP(I)*( FR(I)**2 + FI(I)**2 ) C DO 50 I=0,NF 50 SP(I) = DLOG10( SP(I) ) C RETURN E N D SUBROUTINE FOURIE( X,N,M,FC,FS ) C C ... Discrete Fourier transformation by Goertzel method ... C C Inputs: C X(I): data (I=1,N) C N: data length C M: number of Fourier components C FC(J): Fourier cosine transform (J=1,M) C FS(J): Fourier sine transform (J=1,M) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N), FC(M), FS(M) DATA PI/3.14159265358979D0/ C W = PI/(M-1) DO 20 I=1,M CI = DCOS(W*(I-1)) SI = DSIN(W*(I-1)) T1 = 0.0 T2 = 0.0 DO 10 J=N,2,-1 T0 = 2*CI*T1 - T2 + X(J) T2 = T1 10 T1 = T0 FC(I) = CI*T1 - T2 + X(1) 20 FS(I) = SI*T1 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 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 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