C PROGRAM 13.1 TVVAR C C ... TIME-VARYING VARIANCE ... C C Inputs: C M: Trend Order C IOPT: Search Method C TAU20: Initial Estimate of TAU2 C DELTA: Search Width C Parameters: C NMAX: Adjustable dimension of Y C MJ: Adjustable dimension of XF, VF, etc. C MAXM: Adjustable dimension of A, B, C C NC: Number of components C @TEST.FILTER2: SEP.08,1990 C MODIFIED 2/15/93 C PARAMETER( MJ1=3000,NMAX=MJ1/2,MJ=2,K=1,NDIM=NMAX ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(MJ1), Y2(NMAX) DIMENSION F(MJ,MJ), G(MJ), H(MJ), Q(K,K) DIMENSION XPS(MJ,NMAX), XFS(MJ,NMAX), XSS(MJ,NMAX) DIMENSION VPS(MJ,MJ,NMAX), VFS(MJ,MJ,NMAX), VSS(MJ,MJ,NMAX) DIMENSION XF(MJ), VF(MJ,MJ) DATA OUTMIN, OUTMAX /-1.0D30, 1.0D30/ C IDEV = 1 SIG2 =1.0D0 READ( 5,* ) M, IOPT IF( IOPT.EQ.1 ) READ( 5,* ) TAU20, DELTA C C ... Read Time Series ... C CALL READTS( IDEV,Y,N0 ) sum = 0.0d0 do 5 i=1,n0 5 sum = sum + y(i) ymean = sum/n0 do 6 i=1,n0 6 y(i) = y(i) - ymean N = N0/2 YMIN = 1.0D30 DO 10 I=1,N J = 2*I Y2(I) = (Y(J-1)**2 + Y(J)**2)/2 10 IF( Y2(I).GT.0.0D0 .AND. Y2(I).LT.YMIN ) YMIN = Y2(I) DO 20 I=1,N 20 Y2(I) = DLOG( DMAX1(Y2(I),YMIN/2) ) C NS = 1 NFE = N NPE = N CALL MOMENT( Y2,N/10,YMEAN,VAR ) C FFMAX = -1.0D30 C DO 100 II=1,19 IF( IOPT.NE.0 ) TAU2 = TAU20 + DELTA*(II-9) IF( IOPT.EQ.0 .AND. M.EQ.1 ) TAU2 = 2.0D0**(-II) IF( IOPT.EQ.0 .AND. M.GT.1 ) TAU2 = 2.0D0**(-II-5) C C ... Set Up State Space Model for Trend Estimation ... C CALL SETTRN( M,MJ,F,G,H,R ) R = (3.14159265D0)**2/6 C C ... Set initial values ... C CALL ISTATE( M,MJ,YMEAN,VAR,XF,VF ) C C ... Kalman Filter ... C Q(1,1) = TAU2 CALL FILTER( Y2,XF,VF,F,G,H,Q,R,M,K,1,NS,NFE,NPE,MJ,NMAX, * NDIM,OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 ) IF( FF.GT.FFMAX ) THEN FFMAX = FF TAUMAX = TAU2 SIG2M = SIG2 END IF 100 WRITE(6,600) TAU2, SIG2, FF AIC = -2*FFMAX + 2*(M+2) WRITE(6,610) TAUMAX, SIG2M, FFMAX, AIC C C ... Filtering with the best parameter ... C CALL ISTATE( M,MJ,YMEAN,VAR,XF,VF ) Q(1,1) = TAUMAX CALL FILTER( Y2,XF,VF,F,G,H,Q,R,M,K,1,NS,NFE,NPE,MJ,NMAX, * NDIM,OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 ) C C ... Fixed Interval Smoother ... C CALL SMOOTH( F,M,MJ,NDIM,NS,NFE,NPE,VFS,VPS,XFS,XPS, * VSS,XSS ) C C ... Plot and print data, trend and noise ... C C CALL PTTRND( Y2,XSS,VSS,N,MJ,M,TAUMAX,SIG2,FF,AIC ) CALL PRVAR( Y,M,TAUMAX,SIG2,FF,AIC,XSS,MJ,N,N0 ) C STOP 600 FORMAT( 1H ,5X,F15.10,F12.6,F13.5 ) 610 FORMAT( 1H ,'TAUMAX =',F15.10,3X,'SIG2 =',F15.10,3X, * 'FF =',F13.5,3X,'AIC =',F13.4 ) E N D SUBROUTINE PRVAR( Y,M,TAU2,SIG2,FF,AIC,XSS,MJ,N,N0 ) C C ... Print out changing variance and normalized data ... C C Inputs: C Y: Time Series C M: Trend order C TAU2: Variance of the system noise C SIG2: Variance of the observational noise C FF: Log-Likelihood of the mode C AIC: AIC of the model C XSS: Smoothed state C MJ: Adjustable dimension of F and G C N: Data LENGTH C IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION XSS(MJ,*), Y(*), DATA(1500) COMMON /CMDATA/ TITLE C WRITE(6,600) WRITE(6,610) TITLE WRITE(6,620) M, TAU2, SIG2, FF, AIC DO 10 I=1,N 10 DATA(I) = DEXP(XSS(1,I)+0.57721D0) WRITE(6,630) WRITE(6,640) (DATA(I),I=1,N) DO 20 I=1,N0 J = (I+1)/2 20 Y(I) = Y(I)/DSQRT( DATA(J) ) WRITE(6,650) WRITE(6,640) (Y(I),I=1,N0) RETURN C 600 FORMAT( 1H1,'PROGRAM 12.3: CHANGING VARIANCE ESTIMATION' ) 610 FORMAT( 1H ,'DATA: ',A72 ) 620 FORMAT( 1H ,'M =',I3,3X,'TAU2 =',D12.5,3X,'SIG2 =',D12.5, * 3X,'LOG-L =',F12.4,3X,'AIC =',F12.4 ) 630 FORMAT( 1H ,'*** TIME VARYING VARIANCE ***' ) 640 FORMAT( 1H ,5D13.6 ) 650 FORMAT( 1H ,'*** NORMALIZED DATA ***' ) E N D SUBROUTINE SETTRN( M,MJ,F,G,H,R ) C C ... State space model for trend estimation ... C C Input: C M: Order of the trend C MJ: Adjustable Dimension of F C Outputs: C F: M*M Transition Matrix C G,H: M Vector C R: Observational Noise (=1) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(MJ,MJ), G(MJ), H(MJ) C DO 10 J=1,MJ G(J) = 0.0D0 H(J) = 0.0D0 DO 10 I=1,MJ 10 F(I,J) = 0.0D0 C IF( M.EQ.1 ) THEN F(1,1) = 1.0D0 END IF IF( M.EQ.2 ) THEN F(1,1) = 2.0D0 F(1,2) =-1.0D0 F(2,1) = 1.0D0 END IF IF( M.EQ.3 ) THEN F(1,1) = 3.0D0 F(1,2) =-3.0D0 F(1,3) = 1.0D0 F(2,1) = 1.0D0 F(3,2) = 1.0D0 END IF G(1) = 1.0D0 H(1) = 1.0D0 R = 1.0D0 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 FILTER( Y,XF,VF,F,G,H,Q,R,M,K,ISW,NS,NFE,NPE,MJ,NMAX, * NDIM,OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 ) C C ... Kalman Filter (General Form, L=1) ... C C Inputs: C Y: time series C NS: Start position of filtering C NFE: End position of filtering C NPE: End position of prediction C XF: Initial state vector C VF: Initial covariance matrix C M: Dimension of the state vector C K: Dimension of the system noise C F: M*M matrix C G: M*K matrix C H: M vector C Q: K*K matrix, system noise covariance C R: observation variance C ISW: = 0; R will be estimated C = 1; R is given C MJ: Adjustable dimension of XF, VF C NDIM: Adjustable dimension of XFS, XPS, VFS, VPS C = 0 XF, XP, VF, VP are not stored C > 0 They are stored for smoothing C NMAX Adjustable dimension of Y C OUTMIN: Lower limit for detecting outliers C OUTMAX: Upper limit for detecting outliers C Outputs: C VFS: Covariance matrices of the filter C VPS: Covariance matrices of the predictor C XFS: Mean vectors of the filter C XPS: Mean vectors of the predictor C FF: Log likelihood C SIG2: Estimated observational noise variance C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX) DIMENSION F(MJ,MJ), G(MJ,K), H(MJ), Q(K,K) DIMENSION XF(MJ), VF(MJ,MJ), XP(40), VP(40,40) DIMENSION XFS(MJ,NDIM), XPS(MJ,NDIM) DIMENSION VFS(MJ,MJ,NDIM), VPS(MJ,MJ,NDIM) DIMENSION WRK(40,40), VH(40), GAIN(40) DATA PI /3.1415926535D0/ C SIG2 = 0.0D0 SDET = 0.0D0 NSUM = 0 C DO 500 II=NS,NPE C C ... ONE STEP AHEAD PREDICTION ... C DO 20 I=1,M SUM = 0.0D0 DO 10 J=1,M 10 SUM = SUM + F(I,J)*XF(J) 20 XP(I) = SUM C DO 40 I=1,M DO 40 J=1,M SUM = 0.0D0 DO 30 JJ=1,M 30 SUM = SUM + F(I,JJ)*VF(JJ,J) 40 WRK(I,J) = SUM C DO 60 I=1,M DO 60 J=1,M SUM = 0.0D0 DO 50 JJ=1,M 50 SUM = SUM + WRK(I,JJ)*F(J,JJ) 60 VP(I,J) = SUM C DO 80 I=1,M DO 80 J=1,K SUM = 0.0D0 DO 70 JJ=1,K 70 SUM = SUM + G(I,JJ)*Q(JJ,J) 80 WRK(I,J) = SUM C DO 100 I=1,M DO 100 J=1,M SUM = VP(I,J) DO 90 JJ=1,K 90 SUM = SUM + WRK(I,JJ)*G(J,JJ) 100 VP(I,J) = SUM C C ... FILTERING ... C IF( Y(II).GT.OUTMIN.AND.Y(II).LT.OUTMAX.AND. II.LE.NFE ) THEN C DO 210 I=1,M SUM = 0.0D0 DO 200 J=1,M 200 SUM = SUM + VP(I,J)*H(J) 210 VH(I) = SUM C PERR = Y(II) PVAR = R DO 220 I=1,M PERR = PERR - H(I)*XP(I) 220 PVAR = PVAR + H(I)*VH(I) C DO 250 I=1,M 250 GAIN(I) = VH(I)/PVAR C DO 290 I=1,M 290 XF(I) = XP(I) + GAIN(I)*PERR C DO 310 I=1,M DO 310 J=1,M 310 VF(I,J) = VP(I,J) - GAIN(I)*VH(J) C SIG2 = SIG2 + PERR**2/PVAR SDET = SDET + DLOG(PVAR) NSUM = NSUM + 1 C C ... MISSING OBSERVATION ... C ELSE DO 350 I=1,M XF(I) = XP(I) DO 350 J=1,M 350 VF(I,J) = VP(I,J) END IF C C ... SAVE MEAN AND COVARIANCE ... C IF( NDIM.GT.1 ) THEN DO 360 I=1,M XPS(I,II) = XP(I) XFS(I,II) = XF(I) DO 360 J=1,M VPS(I,J,II) = VP(I,J) 360 VFS(I,J,II) = VF(I,J) END IF C 500 CONTINUE SIG2 = SIG2/NSUM IF(ISW.EQ.0) FF = -0.5D0*(NSUM*(DLOG(PI*2*SIG2) + 1) + SDET) IF(ISW.EQ.1) FF = -0.5D0*(NSUM*(DLOG(PI*2) + SIG2) + SDET) C RETURN E N D SUBROUTINE SMOOTH( F,M,MJ,NDIM,NS,NFE,NPE,VFS,VPS,XFS,XPS, * VSS,XSS ) C C ... Fixed Interval Smoother (General Form) ... C C Inputs: C NS: Start position of filtering C NFE: End position of filtering C NPE: End position of prediction C M: Dimension of the state vector C F: M*M matrix C MJ: Adjustable dimension of XF, VF C NDIM: Adjustable dimension of XFS, XPS, VFS, VPS C NMAX Adjustable dimension of Y C VFS: Covariance matrices of the filter C VPS: Covariance matrices of the predictor C XFS: Mean vectors of the filter C XPS: Mean vectors of the predictor C Outputs: C VSS: Covariance matrices of the smoother C XSS: Mean vectors of the smoother C IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(MJ,MJ) DIMENSION XS(40), VS(40,40), VP(40,40) DIMENSION XFS(MJ,NDIM), XPS(MJ,NDIM), XSS(MJ,NDIM) DIMENSION VFS(MJ,MJ,NDIM), VPS(MJ,MJ,NDIM), VSS(MJ,MJ,NDIM) DIMENSION WRK(40,40), SGAIN(40,40) C C C ... SMOOTHING ... C DO 10 II=NFE,NPE DO 10 I=1,M XSS(I,II) = XFS(I,II) DO 10 J=1,M 10 VSS(I,J,II) = VFS(I,J,II) DO 20 I=1,M XS(I) = XFS(I,NFE) DO 20 J=1,M 20 VS(I,J) = VFS(I,J,NFE) C DO 500 II=NFE-1,NS,-1 C NZERO = 0 DO 100 I=1,M 100 IF( VFS(I,I,II).GT.1.0D-12 ) NZERO = NZERO + 1 C IF( NZERO.EQ.0 ) THEN DO 110 I=1,M XS(I) = XFS(I,II) XSS(I,II) = XFS(I,II) DO 110 J=1,M VS(I,J) = VFS(I,J,II) 110 VSS(I,J,II) = VFS(I,J,II) C ELSE DO 410 I=1,M DO 410 J=1,M 410 VP(I,J) = VPS(I,J,II+1) C CALL GINVRS( VP,VDET,M,40 ) C DO 425 I=1,M DO 425 J=1,M SUM = 0.0D0 DO 420 IJ=1,M 420 SUM = SUM + VFS(I,IJ,II)*F(J,IJ) 425 WRK(I,J) = SUM C DO 440 I=1,M DO 440 J=1,M SUM = 0.0D0 DO 430 IJ=1,M 430 SUM = SUM + WRK(I,IJ)*VP(IJ,J) 440 SGAIN(I,J) = SUM C DO 450 I=1,M XS(I) = XFS(I,II) DO 450 J=1,M WRK(I,J) = 0.0D0 450 VS (I,J) = VFS(I,J,II) C DO 460 J=1,M DO 460 I=1,M 460 XS(I) = XS(I) + SGAIN(I,J)*(XSS(J,II+1) - XPS(J,II+1)) C DO 470 J=1,M DO 470 IJ=1,M DO 470 I=1,M 470 WRK(I,J)=WRK(I,J) + SGAIN(I,IJ)*(VSS(IJ,J,II+1)-VPS(IJ,J,II+1)) C DO 480 J=1,M DO 480 IJ=1,M DO 480 I=1,M 480 VS(I,J) = VS(I,J) + WRK(I,IJ)*SGAIN(J,IJ) DO 485 I=1,M 485 IF( VS(I,I).LT.0.0D0 ) VS(I,I) = 0.0D0 C DO 490 I=1,M XSS(I,II) = XS(I) DO 490 J=1,M 490 VSS(I,J,II) = VS(I,J) END IF C 500 CONTINUE C RETURN E N D SUBROUTINE GINVRS( A,DET,M,MJ ) C C ... Generalized inverse of a square matrix A ... C C Inputs: C A: M*M matrix C M: Dimension of A C MJ: Adjustable dimension of A C Outputs: C A: Generalize inverse of A C DET: Determinant of A C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(MJ,MJ), IND(50) EPS = 1.0D-10 C DO 10 I=1,M 10 IND(I) = I DO 60 L=1,M AMAX = 0.0D0 DO 20 I=L,M IF( A(IND(I),IND(I)).GT.AMAX ) THEN AMAX = A(IND(I),IND(I)) I0 = I END IF 20 CONTINUE IF( AMAX.GT.EPS*A(IND(1),IND(1)) ) THEN IMAX = IND(I0) DO 30 I=I0,L+1,-1 30 IND(I) = IND(I-1) IND(L) = IMAX LMAX = L DO 40 I=L+1,M A(IND(I),IMAX) = -A(IND(I),IMAX)/A(IMAX,IMAX) DO 40 J=L+1,M 40 A(IND(I),IND(J)) = A(IND(I),IND(J)) * + A(IND(I),IMAX)*A(IMAX,IND(J)) ELSE DO 50 I=L,M DO 50 J=L,M 50 A(IND(I),IND(J)) = 0.0D0 GO TO 70 END IF 60 CONTINUE C 70 DET = 1.0D0 DO 80 I=1,M C DET = DET*A(IND(I),IND(I)) C IF( A(IND(I),IND(I)).GT.EPS*AMAX ) THEN IF( A(IND(I),IND(I)).GT.0.0D0 ) THEN A(IND(I),IND(I)) = 1.0D0/A(IND(I),IND(I)) ELSE A(IND(I),IND(I)) = 0.0D0 END IF 80 CONTINUE C MS = MIN0( M-1,LMAX ) DO 200 L=MS,1,-1 DO 100 J=L+1,M SUM = 0.0D0 DO 90 I=L+1,M 90 SUM = SUM + A(IND(I),IND(L))*A(IND(I),IND(J)) 100 A(IND(L),IND(J)) = SUM SUM = A(IND(L),IND(L)) DO 110 I=L+1,M 110 SUM = SUM + A(IND(L),IND(I))*A(IND(I),IND(L)) A(IND(L),IND(L)) = SUM DO 120 I=L+1,M 120 A(IND(I),IND(L)) = A(IND(L),IND(I)) IMAX = IND(L) DO 130 I=L+1,M IF( IMAX.GT.IND(I) ) THEN IND(I-1) = IND(I) IND(I) = IMAX END IF 130 CONTINUE 200 CONTINUE C RETURN E N D SUBROUTINE ISTATE( M,MJ,XMEAN,XVAR,XF,VF ) C C ... Initial state ... C C Inputs: C M: Dimension of the state vector C MJ: Adjustable dimension of F C XMEAN: Mean value C XVAR: Varaince C Outputs: C XF: State vector, X(0|0) C VF: State covarance matrix, V(0|0) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION XF(MJ), VF(MJ,MJ) C DO 10 I=1,MJ DO 10 J=1,MJ 10 VF(I,J) = 0.0D0 C DO 20 I=1,M 20 XF(I) = XMEAN C DO 30 I=1,M 30 VF(I,I) = XVAR C RETURN E N D SUBROUTINE MOMENT( Y,N,YMEAN,VAR ) C C ... Mean and variance of the data ... C C Inputs: C Y(I): data C N: data length C Outputs: C YMEAN: mean C VAR: variance C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N) C SUM = 0.0D0 DO 10 I=1,N 10 SUM = SUM + Y(I) YMEAN = SUM/N SUM = 0.0D0 DO 20 I=1,N 20 SUM = SUM + (Y(I)-YMEAN)**2 VAR = SUM/N C RETURN E N D