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 ) 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 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 PTTRND( Y,XSS,VSS,N,MJ,M,TAU2,SIG2,FF,AIC ) C C ... Plot Original time series, trend ans residuals ... C C Inputs: C Y(I): Time Series C XSS: Smoothed state C VSS: Covariance matrices C N: Data length C MJ: Adjustable dimension of F and G C M: Model order C TAU2: Variance of the model C FF: Log-Likelihood of the mode C AIC: AIC of the model C IMPLICIT REAL*8(A-H,O-Z) CHARACTER VNAME*8 DIMENSION Y(N) DIMENSION XSS(MJ,N), VSS(MJ,MJ,N) DIMENSION DATA(1500), VNAME(20), VALUE(20) C VNAME(1) = 'M1 = ' VNAME(2) = 'TAU2 = ' VNAME(3) = 'SIG2 = ' VNAME(4) = 'FF = ' VNAME(5) = 'AIC = ' VALUE(1) = M VALUE(2) = TAU2 VALUE(3) = SIG2 VALUE(4) = FF VALUE(5) = AIC CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 10.2: TREND ESTIMATION ',38,5, * VNAME,VALUE ) WX = 15.0 WY = 4.5 DX =200.0 FN = N IY = 10 C C ... Plot original data ... C CALL MAXMIN( Y,N,YMIN1,YMAX1,DY ) CALL NEWPEN( 2 ) CALL PLOT( 3.0,16.0-SNGL(WY),-3 ) CALL SYMBOL( 0.25,SNGL(WY)+0.25,0.25,'TIME SERIES',0.0,11 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY,0.0D0,FN,YMIN1,YMAX1,DX,DY, * 0.2D0,1,IY,2 ) CALL PLOTY( Y,N,YMIN1,YMAX1,WX,WY,1,1 ) C C ... PLOT TREND ... C CALL PLOT( 0.0,-SNGL(WY+1.0),-3 ) CALL NEWPEN( 2 ) CALL SYMBOL( 0.25,SNGL(WY)+0.25,0.25,'TREND',0.0,5 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY,0.0D0,FN,YMIN1,YMAX1,DX,DY, * 0.2D0,1,IY,2) C DO 120 J=-1,1 DO 100 I=1,N 100 DATA(I) = XSS(1,I) + DSQRT( VSS(1,1,I)*SIG2 )*J CALL NEWPEN( 1 ) IF(J.EQ.0) CALL NEWPEN( 2 ) CALL PLOTY( DATA,N,YMIN1,YMAX1,WX,WY,1,1 ) 120 CONTINUE C C ... Plot residuals ... C DO 180 J=4,4,3 DO 160 I=1,N 160 DATA(I) = Y(I) - XSS(1,I) CALL MAXMIN( DATA,N,DMIN,DMAX,DDY ) YMIN3 = DY*INT(DMIN/DY) YMAX3 = DY*INT(DMAX/DY) IF( DMIN.LT.YMIN3 ) YMIN3 = YMIN3-DY IF( DMAX.GT.YMAX3 ) YMAX3 = YMAX3+DY WY3 = WY*(YMAX3-YMIN3)/(YMAX1-YMIN1) C CALL PLOT( 0.0,-SNGL(WY3+1.0),-3 ) CALL NEWPEN( 2 ) CALL SYMBOL( 0.25,SNGL(WY3)+0.25,0.25,'RESIDUAL',0.0,8 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY3,0.0D0,FN,YMIN3,YMAX3,DX,DY, * 0.2D0,1,IY,2) CALL NEWPEN( 1 ) IF(J.EQ.4) CALL NEWPEN( 2 ) CALL PLOTY2( DATA,N,YMIN3,YMAX3,WX,WY3,1,1,0.0D0 ) 180 CONTINUE C CALL PLOTE C call plot( 0.0,0.0,999 ) C RETURN 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 AXISXY(X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC,IBOX,IX,IY) C C ... This subroutine draws X and Y axes. C C Inputs: C X,Y: location of the left bottom of the figure C WX,WY: width and height of the figure C X0,X1: lower and upper bounds of the X axis C Y0,Y1: lower and upper bounds of the Y axis C DX,DY: step width in X and Y axes C DWC: size of characters C IBOX: = 0 draw X and Y axes only C = 1 draw window doundary C IX,IY: number of additional click in each step C Modified: 8/31/90, 11/21/90 C REAL*8 X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC C FX = X FY = Y FWX = WX FWY = WY FWC = DWC CALL NEWPEN( 2 ) CALL PLOT( FX,FY,-3 ) CALL PLOT( 0.0,FWY,3 ) CALL PLOT( 0.0,0.0,2 ) CALL PLOT( FWX,0.0,2 ) IF ( IBOX.EQ.1 ) THEN CALL PLOT( FWX,FWY,2 ) CALL PLOT( 0.0,FWY,2 ) END IF C C ... draw X axis ... C WC = DWC NX = (X1-X0)/DX + 0.001 NY = (Y1-Y0)/DY + 0.001 MX = -DLOG10( DX ) + 0.001 IF(DX.LT.1.0D0) MX = -DLOG10( DX ) + 0.999 IF( MX.LE.0 ) MX = -1 D = WX*DX/(X1-X0) DO 10 I=1,NX+1 XN = X0 + DX*(I-1) XX = D*(I-1) CALL PLOT( XX,-FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF LX = 0 IF(XN.GT.0.0D0) LX = ALOG10( XN ) + 0.001 IF( MX.GT.0 ) XXX = WC*(LX+MX)*0.5 + WC*4/7.0 IF( MX.LT.0 ) XXX = WC*LX*0.5 + WC*4/7.0 CALL NUMBER( XX-XXX,-(1.5*FWC+0.1),FWC,XN,0.0,MX ) 10 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 20 I=1,NX+1 DO 20 J=1,IX XX = D*(I-1) + (D*J)/IX IF( XX.LE.WX ) THEN CALL PLOT( XX,FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF END IF 20 CONTINUE CALL NEWPEN( 2 ) END IF C C ... draw Y axis ... C D = WY*DY/(Y1-Y0) MY = -DLOG10( DY ) + 0.99999 IF( MY.LE.0 ) MY = -1 DO 30 I=1,NY+1 YN = Y0 + DY*(I-1) YY = D*(I-1) CALL PLOT( -FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF LY = 1 IF(YN.LT. 0.0) LY = 2 IF(YN.GE. 1.0) LY = ALOG10( YN )+1.00001 IF(YN.LE.-1.0) LY = ALOG10(-YN )+2.00001 IF( MY.GT.0 ) YYY = WC*(LY+MY+1) + 0.20 IF( MY.LT.0 ) YYY = WC*LY + 0.20 CALL NUMBER( -YYY,YY-FWC/2,FWC,YN,0.0,MY ) 30 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 40 I=1,NY+1 DO 40 J=1,IY YY = D*(I-1) + (D*J)/IY IF( YY.LE.WY ) THEN CALL PLOT( FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF END IF 40 CONTINUE CALL NEWPEN( 2 ) END IF C RETURN E N D SUBROUTINE HEADER( PNAME,NP,M,VNAME,VALUE ) C C ... This subroutine draw title of data, program name, names C and values of parameters ... C C Inputs: C PNAME: program name C NP: number of characters of program name C M: number of parameters C VNAME(I): name of the I-th parameter C VALUE(I): value of the I-th parameter C Input via common: C TITLE: title of the data set C CHARACTER PNAME*80, TITLE*72 CHARACTER*8 VNAME(M) REAL*8 VALUE DIMENSION VALUE(M), DAY(2), TIME(3) COMMON /CMDATA/ TITLE C CALL NEWPEN( 2 ) CALL DATE( DAY ) CALL CLOCK( TIME,1 ) C CALL SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 ) CALL SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP ) CALL NEWPEN( 1 ) CALL SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 ) CALL SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 ) YY = 17.9 DO 10 I=1,M IF( MOD(I,4).EQ.1 ) THEN YY = YY - 0.40 X1 = 2.0 X2 = 3.5 ELSE X1 = X1 + 5.0 X2 = X2 + 5.0 END IF CALL SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 ) 10 CALL NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 ) C RETURN E N D SUBROUTINE MAXMIN( X,N,XMIN0,XMAX0,DXL ) C C ... This subroutine determines the minimum, the maximum and C the step width ... C C Inputs: C X(I): data C N: data length C Outputs: C XMIN0: the minimum bound for the figure C XMAX0: the maximum bound for the figure C DXL: step width C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N) C XMIN = 1.0D30 XMAX =-1.0D30 C DO 10 I=1,N IF( X(I) .LT. XMIN ) XMIN = X(I) 10 IF( X(I) .GT. XMAX ) XMAX = X(I) DX = XMAX-XMIN IF( DLOG10(DX) .GE. 0.0D0 ) DXL = INT( DLOG10(DX) ) IF( DLOG10(DX) .LT. 0.0D0 ) DXL = INT( DLOG10(DX) )-1.0 DXL = 10.0**DXL IF( DX/DXL.GT.6.0D0 ) DXL = DXL*2.0 DIF = INT( DX/DXL ) XMIN0 = INT( XMIN/DXL )*DXL XMAX0 = XMIN0 + DIF*DXL IF( XMIN0 .GT. XMIN ) XMIN0 = XMIN0 - DXL 30 IF( XMAX0 .GE. XMAX ) GO TO 40 XMAX0 = XMAX0 + DXL GO TO 30 40 CONTINUE C RETURN E N D SUBROUTINE PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C C ... This subroutine draws time vs data plot of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 connect data by straight line C = 2 connect data by dashed line C = 3 express data by character (character code=ISW) C Modified: 8/31/90 C REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IPOS) YY = (Y(I) - YMIN)*DY IF( I.EQ.1 .AND. ISW.LT.3 ) THEN CALL PLOT( XX,YY,3 ) ELSE IF(ISW.EQ.1) CALL PLOT( XX,YY,2 ) IF(ISW.EQ.2) CALL DASHP( XX,YY,0.2 ) IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) END IF 100 CONTINUE C RETURN E N D SUBROUTINE PLOTY2( Y,N,YMIN,YMAX,WX,WY,IOFF,ISW,YLEVEL ) C C ... This subroutine draws box car graph of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 draw box with thick line C = 2 draw box with thin line C YLEVEL: ground level for the box car C REAL*8 Y(N), YMIN, YMAX, WX, WY, YLEVEL C DX = WX/(N-1+IOFF) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 3 ) IF( N.GT.75 ) CALL NEWPEN( 2 ) IF( N.GT.150) CALL NEWPEN( 1 ) IF( ISW.EQ.2 ) CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IOFF) YY = (Y(I) - YMIN)*DY Y0 = (YLEVEL-YMIN)*DY CALL PLOT( XX,Y0,3 ) CALL PLOT( XX,YY,2 ) 100 CONTINUE CALL NEWPEN( 2 ) C RETURN E N D SUBROUTINE READTS( IDEV,Y,N ) REAL*8 Y(N) CHARACTER*72 TITLE COMMON /CMDATA/ TITLE C C OPEN( IDEV,FILE='temp.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N READ( IDEV,* ) (Y(I),I=1,N) C 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