C PROGRAM 8.2 LSAR2 C C ... Estimation of the change point ... C C Inputs: C K: Highest order od AR model C [N0,NE]: Time interval used for model fitting C [N1,N2]: Candidate for change point 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 IDEV: Input device for time series C NMAX: Adjustable dimension of Y (NMAX.GE.N) C MJ,K,NSPAN: Adjustable dimensions C MODIFIED 2/15/93 C PARAMETER(NMAX=3000,MJ=100,NSPAN=1000,KMAX=20,K1=KMAX+1,IDEV=1) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX) DIMENSION X(MJ,K1), AIC1(NSPAN), AIC2(NSPAN), AICS(NSPAN) C DIMENSION DATA(200) C READ( 5,* ) K, N0, N1, N2, NE NS = 1 M = N2 - N1 C C ... Read time series ... C CALL READTS( IDEV,Y,N ) C C ... First models and AIC1's ... C CALL UPDATE( X,Y,N0,N1,M,NS,K,MJ,AIC1 ) C C ... Second models and AIC2's ... C CALL BUPDAT( X,Y,N2,NE,M,NS,K,MJ,AIC2 ) C C ... AICs of the locally stationary AR models ... C DO 20 I=1,M 20 AICS(I) = AIC1(I) + AIC2(I) C AICMIN = 1.0D30 DO 60 I=1,M IF( AICS(I) .GT. AICMIN ) GO TO 60 AICMIN = AICS(I) MMIN = I 60 CONTINUE C C ... Print out and plot results ... C CALL PRLSAR( AICS,AICMIN,MMIN,M,N0,N1,N2,NE,K ) C CALL PTLSAR( Y,AICS,AICMIN,MMIN,N,M,N0,N1,N2,NE,K ) C STOP E N D SUBROUTINE UPDATE( X,Z,N0,N1,M,NS,K,MJ,AIC ) C C ... Fit AR models to first part of data ... C C Inputs: C X: Working area C Z: Time series C N0,N1,M,NS: Fit AR models on [N0,N1],[N0,N1+NS],..., C [N0,N1+(M-1)*NS] C K: Highest order of AR models C MJ: Adjustable dimension C Output: C AIC(I): AIC of the AR model fitted on [N0,N1+(I-1)*NS] C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ,1), AIC(1), D(200), A(20,20), AICS(0:20) DIMENSION Z(N1), SIG2(0:20) EXTERNAL SETXAR C NMK = N1 - K - N0 C CALL REDUCT( SETXAR,Z,D,NMK,N0,K,MJ,X ) C DO 100 I=1,M II = N1 + (I-1)*NS CALL REGRES( X,K,II-K-N0,MJ,20,A,SIG2,AICS,IMIN ) AIC(I) = AICS(IMIN) C CALL SETXAR( Z,II-K,NS,K,MJ,1,X ) CALL HUSHL2( X,D,MJ,K+1+NS,K+1 ) C 100 CONTINUE RETURN E N D SUBROUTINE BUPDAT( X,Z,N2,N,M,NS,K,MJ,AIC ) C C ... Fit AR models to the second part of data ... C C Inputs: C X: Working area C Z: Time series C N0,N2,N,M,NS: Fit AR models on [N2-(M-1)*NS,N],..., C [N2-NS,S],[N2,N] C K: Highest order of AR models C MJ: Adjustable dimension C Output: C AIC(I): AIC of the AR model fitted on [N0,N1+(I-1)*NS] C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ,1), AIC(1), D(200), A(20,20), AICS(0:20) DIMENSION Z(1), SIG2(0:20) EXTERNAL SETXAR C NMK = N - N2 J = M C CALL REDUCT( SETXAR,Z,D,NMK,N2-K-NS,K,MJ,X ) C DO 100 I=1,M II = N2 - (I-2)*NS CALL REGRES( X,K,N-II,MJ,20,A,SIG2,AICS,IMIN ) AIC(J) = AICS(IMIN) J = J - 1 C CALL SETXAR( Z,II-K-NS,NS,K,MJ,1,X ) CALL HUSHL2( X,D,MJ,K+1+NS,K+1 ) C 100 CONTINUE RETURN E N D SUBROUTINE HUSHL2( X,D,MJ1,N,K ) C C ... Householder transformation for adding new observations ... C C Inputs: C X: Data matrix (Householder reduced form augmented with C new observations) C D: Working area C MJ1: Adjustable dimension C N: Number of rows of X C K: Number of columns of X C Output: C X: Householder reduced form C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(MJ1,K) , D(MJ1) C TOL = 1.0D-30 C DO 100 II=1,K D1 = X(II,II) H = D1**2 DO 10 I=K+1,N D(I) = X(I,II) 10 H = H + D(I)**2 IF( H .GT. TOL ) GO TO 20 G = 0.0D00 GO TO 100 20 G = DSQRT( H ) IF( D1 .GE. 0.0D00 ) G = -G H = H - D1*G D1 = D1 - G C IF( II .EQ. K ) GO TO 100 II1 = II+1 DO 60 J=II1,K S = D1*X(II,J) DO 40 I=K+1,N 40 S = S + D(I)*X(I,J) S = S/H X(II,J) = X(II,J) - D1*S DO 50 I=K+1,N 50 X(I,J) = X(I,J) - D(I)*S 60 CONTINUE 100 X(II,II) = G C RETURN E N D SUBROUTINE PRLSAR( AICS,AICMIN,MMIN,M,N0,N1,N2,NE,K ) C C ... Print out AIC's and change point ... C C Inputs: C AICS(I): AIC's of the LSAR models C AICMIN: Minimum of AIC's C MMIN: Estimated change point C M: Number of models C N0,N1,N2,NE: C K: Highest order of AR models C IMPLICIT REAL*8(A-H,O-Z) DIMENSION AICS(M) C NMIN = N1 + MMIN WRITE(6,600) WRITE(6,610) K, N0, N1, N2, NE WRITE(6,620) NMIN, AICMIN WRITE(6,630) WRITE(6,640) (AICS(I),I=1,M) C RETURN 600 FORMAT( 1H ,'PROGRAM 8.2' ) 610 FORMAT( 1H ,' K =',I3,3X,'N0 =',I4,3X,'N1 =',I4,3X, * 'N2 =',I4,3X,'NE =',I4 ) 620 FORMAT( 1H ,'NMIN =',I5,5X,'AICMIN =',F12.3 ) 630 FORMAT( 1H ,'AIC(1),...,AIC(M)' ) 640 FORMAT( 1H ,10F8.2 ) 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 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