C PROGRAM 9.1 SMOOTH C C ... Prediction and Interpolation of Time Series ... C C Inputs: C NFE: End point of filtering C NPE: End point of prediction C M: Order of the AR model C TAU2: Innovation variance of the AR model C AR(I): AR coefficients (II=1,M) C OUTMIN,OUTMAX: Lower and upper limits of observations C NMISS: Number of missed intervals C N0(I): Start position of missed intervals (I=1,NMISS) C NN(I): Number of missed observations (I=1,NMISS) C Inputs required in subroutine READTS C TITLE: Title of the data set C N: Data length C Y(I): Tiem series (I=1,N) C Parameters: C IDEV: Input device specification C NMAX: Adjustable dimension of Y, YMISS (NMAX >= N) C MJ: Highest dimension of the state (MJ >= M) C K: Dimension of the system noise C ISW: =1 (R is specified) C @TEST.FILTER2: SEP.08,1990, SEP.02,1992 C MODIFIED 2/15/93 C PARAMETER( NMAX=1000,MJ=20,K=1,ISW=1,IDEV=1 ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), YMISS(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) DIMENSION N0(10), NN(10), AR(MJ) C READ( 5,* ) NFE, NPE READ( 5,* ) M READ( 5,* ) TAU2 READ( 5,* ) (AR(I),I=1,M) READ( 5,* ) OUTMIN, OUTMAX READ( 5,* ) NMISS DO 10 I=1,NMISS 10 READ( 5,* ) N0(I), NN(I) NS = 1 C C ... Read Time Series ... C CALL READTS( IDEV,Y,N ) CALL MOMENT( Y,N,YMEAN,YVAR ) C C ... SET MISSING OBSERVATIONS ... CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 20 I=1,N 20 YMISS(I) = Y(I) - YMEAN DO 30 J=1,NMISS DO 30 I=1,NN(J) 30 YMISS(N0(J)+I-1) = OUTMIN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ... Set State Space Model ... C CALL SETFGH( M,MJ,K,AR,TAU2,F,G,H,Q,R ) C C ... Set Initial State ... C CALL ISTATE( M,MJ,0.0D0,YVAR,XF,VF ) C C ... Kalman Filtering ... C CALL FILTER( YMISS,XF,VF,F,G,H,Q,R,M,K,ISW,NS,NFE,NPE,MJ,NMAX, * NMAX,OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 ) C C ... Fixed Interval Smoothing ... C CALL SMOOTH( F,M,MJ,NDIM,NS,NFE,NPE,VFS,VPS,XFS,XPS, * VSS,XSS ) FLK = -FF AIC = -2*FF + 2*(M+1) C C ... Plot and Print ... C C CALL PTSTAT( Y,M,TAU2,FLK,AIC,XSS,VSS,N,NFE,NPE,MJ, C * YMEAN,NMISS,N0,NN ) CALL PRSTAT( M,TAU2,FLK,AIC,XSS,VSS,YMEAN,MJ,NFE,NPE, * NMISS,N0,NN ) C STOP 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 PRSTAT( M,TAU2,FF,AIC,XSS,VSS,YMEAN,MJ,NFE,NPE, * NMISS,N0,NN ) C C ... Print out Estimates ... C C Inputs: C M: Model order C TAU2: Variance of the model C FF: Log-Likelihood of the mode C AIC: AIC of the model C XSS: Smoothed state C VSS: Covariance matrices C YMEAN: Mean of the original time series C MJ: Adjustable dimension of F and G C NFE: End position of filtering C NPE: End position of prediction C NMISS: Number of missed intervals C N0(I): Start position of I-th missed interval C NN(I): Numbers of missed observations in I-th interval C IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION XSS(MJ,*), VSS(MJ,MJ,*), N0(10), NN(10) COMMON /CMDATA/ TITLE C WRITE(6,600) WRITE(6,610) TITLE WRITE(6,620) M, TAU2, FF, AIC WRITE(6,630) DO 10 J=1,NMISS DO 10 I=1,NN(J) II = N0(J)+I-1 XM = XSS(1,II) + YMEAN 10 WRITE(6,640) II, XM, VSS(1,1,II) DO 20 I=NFE+1,NPE XM = XSS(1,I) + YMEAN 20 WRITE(6,640) I, XM, VSS(1,1,I) RETURN C 600 FORMAT( 1H ,'PROGRAM 9.1: PREDICTION AND INTERPOLATION' ) 610 FORMAT( 1H ,'DATA: ',A72 ) 620 FORMAT( 1H ,'M =',I3,3X,'TAU2 =',D12.5,3X,'LOG-L =',F12.4, * 3X,'AIC =',F12.4 ) 630 FORMAT( 1H ,' I',10X,'MEAN',12X,'VARANCE' ) 640 FORMAT( 1H ,I4,3X,F15.6,5X,D13.6 ) E N D SUBROUTINE SETFGH( M,MJ,K,AR,TAU2,F,G,H,Q,R ) C C ... Specify the state space representation of AR model ... C C Inputs: C M: Order of AR model C MJ: Adjustable dimension of F C AR(I): AR coefficents C TAU2: Innovation variance C Outputs: C F: State transition matrix C G,H: M-vector C Q: System noise varaince C R: Observational noise varaince (=0) C MODIFIED 2/15/93 C IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(MJ,MJ), G(MJ), H(MJ), AR(MJ), Q(K,K) 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 DO 20 I=1,M 20 F(1,I) = AR(I) DO 30 I=2,M 30 F(I,I-1) = 1.0D0 G(1) = 1.0D0 H(1) = 1.0D0 Q(1,1)=TAU2 R = 0.0D0 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 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