C PROGRAM 8.1 LSAR1 C C ... Decomposition of time interval to stationary subintevals ... C C Inputs: C LAG: Highest order od AR model C NS: Basic local span C The following inputs are required in the subroutine READTS. C TITLE: Caption of the data set C N: Data length C Y(I): Time series, (i=1,...,N) C Parameters: C IDEV: Input device for time series C NMAX: Adjustable dimension of Y (NMAX.GE.N) C KMAX,MJ2: Adjustable dimensions C NF: Number of frequencies for computing spectrum C PARAMETER( NMAX=3000,KMAX=20,MJ2=KMAX+1,MJ1=100,NF=100,IDEV=1) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), X(MJ1,MJ2), D(MJ1), U(MJ2,MJ2) DIMENSION AS(KMAX), A(KMAX) EXTERNAL SETXAR C C ISW = 0 READ( 5,* ) LAG, NS C C ... Read time series ... C CALL READTS( IDEV,Y,N ) C IF = 0 C NBLOCK = N/NS DO 100 II=1,NBLOCK C L = NS*(II-1) WRITE(6,600) L IF( II.EQ.NBLOCK ) NS = N - NS*(II-1) - LAG LK = L + LAG C C ... Locally stationary time series ... C CALL LOCAL( SETXAR,Y,X,U,D,LAG,L,NS,LAG,IF,MJ1,MJ2,A,MF,SIG2) C IF( IF .EQ. 2 ) LK0 = LK + 1 IF( IF.EQ.2 .AND. II.GT.1 ) THEN C CALL ARMASP( AS,MFS,B,0,SIG2S,NF,SP ) C CALL PTLSP( Y,N,SP,NF,LK0S,LK2S ) END IF MFS = MF SIG2S = SIG2 LK0S = LK0 LK2S = LK + NS DO 20 I=1,MF 20 AS(I) = A(I) C 100 CONTINUE C CALL ARMASP( AS,MFS,B,0,SIG2S,NF,SP ) C CALL PTLSP( Y,N,SP,NF,LK0S,LK2S ) C CALL PLOTE C call plot( 0.0,0.0,999 ) C STOP 600 FORMAT( 1H ,'L =',I5 ) E N D SUBROUTINE LOCAL( SETX,Z,X,U,D,LAG,N0,NS,K,IF,MJ1,MJ2, * A,MF,SDF ) C C ... Locally stationary AR model ... C C Inputs: C SETX: Name of the subroutine for making X(I,J) C Z(I): Data vector C D,U: Working area C LAG: Highest order of the model C N0: Time point of the previous set ofobservations C NS: Number of new observations C MJ1,MJ2: Adjustable dimension C Output: C A(I): AR coefficients of the current model C MF: Order of the current model C SDF: Innovation variance of the current model C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Z(1) DIMENSION X(MJ1,1), U(MJ2,1), D(1), A(1) DIMENSION B(20), AA(20,20), AIC(0:20), SIG2(0:20) EXTERNAL SETX C K1 = K + 1 K2 = K1*2 NN0 = N0 + LAG + 1 NN1 = N0 + LAG + NS WRITE(6,600) NN0, NN1 C CALL REDUCT( SETX,Z,D,NS,N0,K,MJ1,X ) CALL REGRES( X,K,NS,MJ1,20,AA,SIG2,AIC,MS ) C SDS = SIG2(MS) DO 10 I=1,MS 10 B(I) = AA(I,MS) IF( IF.EQ.0 ) THEN CALL COPY( X,K1,0,0,MJ1,MJ2,U ) AICS = AIC(MS) AIC0 = AIC(MS) WRITE(6,610) NS, MS, SDS, AICS ELSE C AICS = AIC(MS) + AICF AIC0 = AIC(MS) WRITE(6,620) NF, NS, MS, SDS, AICS CALL COPY( X,K1,0,K2,MJ1,MJ1,X ) CALL COPY( U,K1,0,K1,MJ2,MJ1,X ) CALL HUSHLD( X,D,MJ1,K2,K1 ) NP = NF + NS CALL REGRES( X,K,NP,MJ1,20,AA,SIG2,AIC,MP ) C AICP = AIC(MP) SDP = SIG2(MP) DO 20 I=1,MP 20 A(I) = AA(I,MP) WRITE(6,630) NP, MP, SDP, AICP IF( AICS.GE.AICP ) GO TO 40 WRITE(6,640) CALL COPY( X,K1,K2,0,MJ1,MJ2,U ) END IF C IF = 2 NF = NS MF = MS AICF = AIC0 DO 30 I=1,MF 30 A(I) = B(I) SDF = SDS C GO TO 50 40 IF = 1 CALL COPY( X,K1,0,0,MJ1,MJ2,U ) WRITE(6,650) SDF = SDP MF = MP AICF = AICP NF = NF + NS 50 CONTINUE C RETURN 600 FORMAT( 1H0,'<<< NEW DATA (N =',I4,' ---',I4,') >>>' ) 610 FORMAT( 1H ,' INITIAL MODEL: NS =',I4,/,10X,'MS =',I2,3X, 1 'SDS =',D13.6,3X,'AICS =',F12.3 ) 620 FORMAT( 1H ,' SWITCHED MODEL: (NF =',I4,', NS =',I4,1H), 2 /,10X,'MS =',I2,3X,'SDS =',D13.6,3X,'AICS =',F12.3 ) 630 FORMAT( 1H ,' POOLED MODEL: (NP =',I4,1H), 3 /,10X,'MP =',I2,3X,'SDP =',D13.6,3X,'AICP =',F12.3 ) 640 FORMAT( 1H ,30X,'*** SWITCHED MODEL ACCEPTED ***' ) 650 FORMAT( 1H ,30X,'*** POOLED MODEL ACCEPTED ***' ) 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 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 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 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 COPY( X,K,II,JJ,MJ1,MJ2,Y ) C C ... Make a copy of X on Y C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION X(MJ1,1) , Y(MJ2,1) C DO 10 I=1,K I1 = I + II I2 = I + JJ DO 10 J=1,K 10 Y(I2,J) = X(I1,J) C RETURN E N D