C PROGRAM 12.1 SEASON C C ... SEASONAL ADJUSTMENT BY STATE SPACE MODELING ... C C Inputs: C M1: Trend Order (M1=0,1,2,3) C M2: Seasonal Order (M1=0,1,2) C M3: AR Order C M4: Trading day effect (M4=0 or 6) C PERIOD: Number of seasons in one period C (=12 for monthly data, =4 for quarterly data) C IPARAM: =0 Use defalt initail values C =1 Read intial values C OUTMIN: Lower limit of the observations C OUTMAX: Upper limit of the observations C JYEAR: Starting year of the data (e.g. 1967) C MONTH: Starting month of the data (e.g. 1) C TAU2(I): System noise variances (I=1,4) C AR(I): AR coefficients (I=1,M3) C NS: Start position of filtering C NFE: End point of filtering C NPE: End point of prediction C IOPT: = 1 To get MLE 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.FILTER2O NOV.29,1990, SEP.02,1992 C PARAMETER( NMAX=160,MJ=20,MAXM=22,NC=10) IMPLICIT REAL*8(A-H,O-Z) INTEGER*4 PERIOD DIMENSION M(10), TAU2(10), TAU20(10) DIMENSION XMEAN(10), XVAR(MAXM,10), AA(20), AR(10), PAR(10) DIMENSION A(MAXM,10), B(MAXM,10), C(MAXM,10), Q(10,10) 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), MTYPE(10) COMMON /C92827/ M1, M2, M3, M4, PERIOD COMMON /C92825/ OUTMIN, OUTMAX COMMON /C92826/ Y(NMAX), REG(NMAX,7) COMMON /C92907/ ALIMIT COMMON /CMMODL/ M, XMEAN, XVAR, N, NS, NFE, NPE COMMON / CCC / ISW, IPR, ISMT, IDIF EXTERNAL FFSEAS C C ... Read Time Series ... C CALL READTS( 1,Y,N ) C C ... Read Model Orders ... C READ( 5,* ) M1, M2, M3, M4, PERIOD C C ... Set Defalt Parameters ... C IPR = 7 SIG2 = 1.0D0 ALIMIT = 0.9D0 CALL SPARAM( M1,M2,M3,M4,N,TAU2,AR,OUTMIN,OUTMAX,NS,NFE,NPE, * LOGT,IOPT ) C READ( 5,* ) IPARAM IF( IPARAM.EQ.1 ) THEN READ( 5,* ) (TAU2(I),I=1,3) READ( 5,* ) NS, NFE, NPE READ( 5,* ) LOGT, IOPT READ( 5,* ) OUTMIN, OUTMAX IF(M3.GT.0) READ( 5,* ) (AR(I),I=1,M3) END IF WRITE(6,*) M1,M2,M3,M4,PERIOD WRITE(6,*) (TAU2(I),I=1,3) WRITE(6,*) (AR(I),I=1,M3) C IF( M4.GT.0 ) READ( 5,* ) JYEAR, MONTH C C ... Log Transformation ... C IF( LOGT.EQ.1 ) THEN DO 5 I=1,N 5 Y(I) = DLOG10( Y(I) ) END IF C C ... Trading day factor ... C IF( M4.GT.0 ) CALL TRADE( JYEAR,MONTH,N,NMAX,REG ) DO 10 I=1,N DO 10 J=1,6 10 REG(I,J) = REG(I,J) - REG(I,7) C CALL PARCOR( AR,M3,PAR ) NP = ID(M1) + ID(M2) + ID(M3) DO 20 I=1,NP 20 AA(I) = DLOG( TAU2(I)/(1.0D0-TAU2(I)) ) DO 30 I=1,M3 30 AA(NP+I) = DLOG( (ALIMIT+PAR(I))/(ALIMIT-PAR(I)) ) DO 40 I=1,NC 40 TAU20(I) = TAU2(I) C C ... Maximum Likelihood Method ... C IF( IOPT.EQ.1 ) CALL DAVIDN( FFSEAS,AA,NP+M3,2 ) C DO 50 I=1,NP 50 TAU2(I) = DEXP( AA(I) )/(1.0D0 + DEXP( AA(I) ) ) DO 60 I=1,M3 60 PAR(I) = ALIMIT*(DEXP(AA(NP+I))-1.0D0)/(DEXP(AA(NP+I))+1.0D0) CALL ARCOEF( PAR,M3,AR ) C C ... The Maxmum Likelihood Model ... C CALL SETABC( M,L,NC,MTYPE,TAU2,MAXM,MM,AR,Y,N,A,B,C,Q,XMEAN, * XVAR ) CALL ISTAT1( L,M,MJ,MAXM,A,XMEAN,XVAR,XF,VF ) C C ... The Kalman Filter ... C CALL FILTR1( Y,XF,VF,A,B,C,Q,SIG2,REG,MTYPE,NMAX,M,MAXM,NC,L, * NS,NFE,NPE,MJ,N, OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,OVAR ) C C ... Fixed Interval Smoothing ... C CALL SMOTH1( A,M,MAXM,L,NS,NFE,NPE,MJ, * VFS,VPS,VSS,XFS,XPS,XSS ) AIC = -2*FF + 2*(NP+M3+MM) C C ... Print out and draw estimates ... C CALL PRSEAS( M1,M2,M3,M4,PERIOD,TAU2,TAU20,OVAR,FF,AIC,AR, * XSS,N,REG,NMAX,MJ ) CALL PTSEAS( Y,XSS,VSS,N,NPE,MJ,M1,M2,M3,M4, * PERIOD,TAU2,OVAR,FF,AIC,REG,NMAX ) C STOP E N D SUBROUTINE SPARAM( M1,M2,M3,M4,N,TAU2,AR,OUTMIN,OUTMAX,NS,NFE, * NPE,LOGT,IOPT ) C C ... Set Default Parameters ... C C Inputs: C M1,M2,M3,M4: Model orders C N: Data length C Outputs: C TAU2: System noise variance C AR: AR coefficients C OUTMIN: Lower bound for outliers C OUTMAX: Upper bound for outliers C NS: Start Point C NFE: End point of filtering C NPE: End point of prediction C LOGT: (=1 Log transformation) C IOPT: (=1 MLE by numerical optimization) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION TAU2(4), AR(10), PAR(10) C L = 0 IF( M1.GT.0 ) THEN L = 1 TAU2(1) = 0.0500D0 END IF IF( M2.GT.0 ) THEN L = L+1 TAU2(L) = 0.11D-7 END IF IF( M3.GT.0 ) THEN L = L+1 TAU2(L) = 0.9D0 DO 10 I=1,M3 10 PAR(I) = -(-0.8D0)**I CALL ARCOEF( PAR,M3,AR ) END IF IF( M4.GT.0 ) THEN L = L+1 TAU2(L) = 0.0D0 END IF C OUTMIN = -1.0D30 OUTMAX = 1.0D30 NS = 1 NFE = N NPE = N LOGT = 0 IOPT = 1 C RETURN E N D SUBROUTINE PRSEAS( M1,M2,M3,M4,LPER,TAU2,TAU20,OVAR,FF,AIC, * AR,XSS,N,REG,NMAX,MJ ) C C ... Print out Estimates ... C C Inputs: C M1-M4: Model orders C LPER: Period C TAU2: Variance of the system noise C TAU20: Variance of the system noise (initial estimat) C OVAR: Variance of the observational noise C FF: Log-Likelihood of the mode C AIC: AIC of the model C AR: AR coefficients C XSS: Smoothed state vector C N: Data length C REG: Number of trading days C NMAX: Adjustable dimension of REG C MJ: Adjustable dimension of XSS C IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION TAU2(*), TAU20(*), AR(*) DIMENSION XSS(MJ,NMAX), REG(NMAX,7) COMMON /CMDATA/ TITLE C NP = ID(M1) + ID(M2) + ID(M3) M12 = M1 + M2*(LPER-1) M123 = M12 + M3 WRITE(6,600) WRITE(6,610) TITLE WRITE(6,620) M1, M2, M3, M4, LPER WRITE(6,630) WRITE(6,640) WRITE(6,650) (TAU20(I),I=1,NP) WRITE(6,660) WRITE(6,640) WRITE(6,650) (TAU2(I),I=1,NP) IF( M3.GT.0 ) THEN WRITE(6,690) WRITE(6,700) (AR(I),I=1,M3) END IF WRITE(6,680) OVAR, FF, AIC IF( M1.GT.0 ) THEN WRITE(6,720) WRITE(6,710) (XSS(1,I),I=1,N) END IF IF( M2.GT.0 ) THEN WRITE(6,730) WRITE(6,710) (XSS(M1+1,I),I=1,N) END IF IF( M3.GT.0 ) THEN WRITE(6,740) WRITE(6,710) (XSS(M12+1,I),I=1,N) END IF IF( M4.GT.0 ) THEN DO 210 I=1,N SUM = 0.0D0 DO 200 J=1,6 200 SUM = SUM + XSS(M123+J,I)*REG(I,J) 210 REG(I,7) = SUM WRITE(6,750) WRITE(6,710) (REG(I,7),I=1,N) END IF RETURN 600 FORMAT( 1H1,'PROGRAM 10.1: SEASONAL ADJUSTMENT' ) 610 FORMAT( 1H ,A72 ) 620 FORMAT( 1H ,'M1 =',I2,3X,'M2 =',I2,3X,'M3 =',I2,3X,'M4 =',I2, * 'PERIOD =',I3 ) 630 FORMAT( 1H ,5X,'*** INITIAL ESTIMATE ***' ) 640 FORMAT( 1H ,'TAU2(I),I=1,NC:' ) 650 FORMAT( 1H ,5D13.6 ) 660 FORMAT( 1H ,5X,'*** FINAL ESTIMATE ***' ) 680 FORMAT( 1H ,'SIG2 =',D13.6,3X,'LOG-LIKELIHOOD =',F13.4,3X, * 'AIC =',F13.4 ) 690 FORMAT( 1H ,'AR COEFFICIENTS' ) 700 FORMAT( 1H ,5F10.6 ) 710 FORMAT( 1H ,6F12.5 ) 720 FORMAT( 1H ,'TREND COMPONENT' ) 730 FORMAT( 1H ,'SEASONAL COMPONENT' ) 740 FORMAT( 1H ,'AR COMPONENT' ) 750 FORMAT( 1H ,'TRADING DAY EFFECT' ) E N D SUBROUTINE FILTR1( Y,XF,VF,A,B,C,Q,SIG2,REG,MTYPE,NMAX,M,MMAX, * NCM,NC,NS,N,NPE,MJ,NDIM,OUTMIN,OUTMAX, * VFS,VPS,XFS,XPS,FF,OVAR ) C C ... Kalman filter ... C C Inputs: C Y: time series C N: data length C NS: Start position of filtering C NE: End position of prediction C XF: Initial state vector C VF: Initial covariance matrix C A: Parameters of the matrix F C B: Parameters of the matrix G C C: Parameters of the matrix H C Q: Covariance matrix of the system noise C M: Dimension of the state vector C MMAX: Upper limit of the dimension of each component C NC: Number of components 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 OUTMIN: Lower limit for detecting outliers C OUTMAX: Upper limit for detecting outliers C SIG2: Variance of the observational noise 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 OVAR: Estimated variance C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), REG(NMAX,*) DIMENSION A(MMAX,NCM), B(MMAX,NCM), C(MMAX,NCM), M(NCM) DIMENSION XF(MJ), VF(MJ,MJ), XP(40), VP(40,40) DIMENSION XFS(MJ,N), XPS(MJ,N), Q(NCM,NCM) DIMENSION VFS(MJ,MJ,NDIM), VPS(MJ,MJ,NDIM) DIMENSION WRK(40,40), VH(40), GAIN(40) DIMENSION I0(10), MTYPE(NCM) DATA PI /3.1415926535D0/ C OVAR = 0.0D0 SDET = 0.0D0 NSUM = 0 I0(1) = 0 DO 10 I=2,NC 10 I0(I) = I0(I-1) + M(I-1) MM = I0(NC) + M(NC) C DO 300 II=NS,NPE C C ... ONE STEP AHEAD PREDICTION ... C DO 100 L=1,NC XP(I0(L)+M(L)) = A(M(L),L)*XF(I0(L)+1) DO 100 I=1,M(L)-1 100 XP(I0(L)+I) = A(I,L)*XF(I0(L)+1) + XF(I0(L)+I+1) C DO 110 J=1,MM DO 110 L=1,NC WRK(I0(L)+M(L),J) = A(M(L),L)*VF(I0(L)+1,J) DO 110 I=1,M(L)-1 110 WRK(I0(L)+I,J) = A(I,L)*VF(I0(L)+1,J) + VF(I0(L)+I+1,J) C DO 120 I=1,MM DO 120 L=1,NC VP(I,I0(L)+M(L)) = WRK(I,I0(L)+1)*A(M(L),L) DO 120 J=1,M(L)-1 120 VP(I,I0(L)+J) = WRK(I,I0(L)+1)*A(J,L) + WRK(I,I0(L)+J+1) C DO 130 J=1,NC DO 130 L=1,NC DO 130 I=1,M(L) 130 WRK(I0(L)+I,J) = B(I,L)*Q(L,J) C DO 140 I=1,MM DO 140 L=1,NC DO 140 J=1,M(L) 140 VP(I,I0(L)+J) = VP(I,I0(L)+J) + WRK(I,L)*B(J,L) C C ... FILTERING ... C IF( Y(II).GT.OUTMIN .AND. Y(II).LT.OUTMAX .AND. II.LE.N ) THEN C DO 210 I=1,MM SUM = 0.0D0 DO 205 L=1,NC IF( MTYPE(L).EQ.0 ) THEN DO 200 J=1,M(L) 200 SUM = SUM + VP(I,I0(L)+J)*C(J,L) ELSE SUM = SUM + VP(I,I0(L)+1)*REG(II,MTYPE(L)) END IF 205 CONTINUE 210 VH(I) = SUM C PVAR = SIG2 DO 225 L=1,NC IF( MTYPE(L).EQ.0 ) THEN DO 220 J=1,M(L) 220 PVAR = PVAR + VH(I0(L)+J)*C(J,L) ELSE PVAR = PVAR + VH(I0(L)+1)*REG(II,MTYPE(L)) END IF 225 CONTINUE C DO 230 I=1,MM 230 GAIN(I) = VH(I)/PVAR C PERR = Y(II) DO 245 L=1,NC IF( MTYPE(L).EQ.0 ) THEN DO 240 J=1,M(L) 240 PERR = PERR - C(J,L)*XP(I0(L)+J) ELSE PERR = PERR - REG(II,MTYPE(L))*XP(I0(L)+1) END IF 245 CONTINUE C DO 250 I=1,MM 250 XF(I) = XP(I) + GAIN(I)*PERR C DO 260 J=1,MM DO 260 I=1,MM 260 VF(I,J) = VP(I,J) - GAIN(I)*VH(J) C OVAR = OVAR + PERR**2/PVAR SDET = SDET + DLOG(PVAR) NSUM = NSUM + 1 C C ... MISSING OBSERVATION ... C ELSE DO 270 I=1,MM XF(I) = XP(I) DO 270 J=1,MM 270 VF(I,J) = VP(I,J) END IF C C ... SAVE MEAN AND COVARIANCE ... C IF( NDIM.GT.1 ) THEN DO 280 I=1,MM XPS(I,II) = XP(I) XFS(I,II) = XF(I) DO 280 J=1,MM VPS(I,J,II) = VP(I,J) 280 VFS(I,J,II) = VF(I,J) END IF C 300 CONTINUE OVAR = OVAR/NSUM FF = -0.5D0*(NSUM*DLOG(PI*2*OVAR) + SDET + NSUM) C WRITE(6,*) 'FF =',FF, OVAR, SDET, NSUM C RETURN E N D SUBROUTINE SMOTH1( A,M,MMAX,NC,NS,N,NE,MJ, * VFS,VPS,VSS,XFS,XPS,XSS ) C C ... Fixed Interval Smoother (General Form) ... C C Inputs: C A: Parameter of matrix F C M: Dimension of the state vector C MMAX: Adjustable dimension of A C NC: Number of components C N: Data length C NS: Start position of filtering C NE: End position of filtering C MJ: Adjustable dimension of XF, VF 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 A(MMAX,NC), M(NC) DIMENSION XS(40), VS(40,40), VP(40,40) DIMENSION XFS(MJ,N), XPS(MJ,N), XSS(MJ,N) DIMENSION VFS(MJ,MJ,N), VPS(MJ,MJ,N), VSS(MJ,MJ,N) DIMENSION WRK(40,40), SGAIN(40,40) DIMENSION I0(10) C I0(1) = 0 DO 10 I=2,NC 10 I0(I) = I0(I-1) + M(I-1) MM = I0(NC) + M(NC) C C ... SMOOTHING ... C NSS = MIN0( N,NE ) DO 20 I=1,MM XS(I) = XFS(I,NSS) XSS(I,NSS) = XFS(I,NSS) DO 20 J=1,MM VS(I,J) = VFS(I,J,NSS) 20 VSS(I,J,NSS) = VFS(I,J,NSS) DO 30 II=NSS+1,NE DO 30 I=1,MM XSS(I,II) = XFS(I,II) DO 30 J=1,MM 30 VSS(I,J,II) = VFS(I,J,II) C DO 500 II=NSS-1,NS,-1 C NZERO = 0 DO 100 I=1,MM 100 IF( VFS(I,I,II).GT.1.0D-12 ) NZERO = NZERO + 1 C IF( NZERO.EQ.0 ) THEN DO 110 I=1,MM XS(I) = XFS(I,II) XSS(I,II) = XFS(I,II) DO 110 J=1,MM VS(I,J) = VFS(I,J,II) 110 VSS(I,J,II) = VFS(I,J,II) C ELSE DO 410 I=1,MM DO 410 J=1,MM 410 VP(I,J) = VPS(I,J,II+1) C CALL GINVRS( VP,VDET,MM,40 ) C DO 420 I=1,MM DO 420 L=1,NC WRK(I,I0(L)+M(L)) = VFS(I,I0(L)+1,II)*A(M(L),L) DO 420 J=1,M(L)-1 420 WRK(I,I0(L)+J) = VFS(I,I0(L)+1,II)*A(J,L) + VFS(I,I0(L)+J+1,II) C DO 440 I=1,MM DO 440 J=1,MM SUM = 0.0D0 DO 430 IJ=1,MM 430 SUM = SUM + WRK(I,IJ)*VP(IJ,J) 440 SGAIN(I,J) = SUM C DO 450 I=1,MM XS(I) = XFS(I,II) DO 450 J=1,MM WRK(I,J) = 0.0D0 450 VS (I,J) = VFS(I,J,II) C DO 460 J=1,MM DO 460 I=1,MM 460 XS(I) = XS(I) + SGAIN(I,J)*(XSS(J,II+1) - XPS(J,II+1)) C DO 470 J=1,MM DO 470 IJ=1,MM DO 470 I=1,MM 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,MM DO 480 IJ=1,MM DO 480 I=1,MM 480 VS(I,J) = VS(I,J) + WRK(I,IJ)*SGAIN(J,IJ) DO 485 I=1,MM 485 IF( VS(I,I).LT.0.0D0 ) VS(I,I) = 0.0D0 C DO 490 I=1,MM XSS(I,II) = XS(I) DO 490 J=1,MM 490 VSS(I,J,II) = VS(I,J) END IF C 500 CONTINUE C RETURN E N D SUBROUTINE ISTAT1( NC,M,MJ,MAXM,A,XMEAN,XVAR,XF,VF ) C C ... Initial state ... C C Inputs: C NC: Number of components C M: Dimension of the state vector C MJ: Adjustable dimension of F C MAXM: Adjustable dimension of A C A: Parameter of matrix 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 M(NC), A(MAXM,NC), XMEAN(NC), XVAR(MAXM,NC) DIMENSION XF(MJ), VF(MJ,MJ), I0(10) C I0(1) = 0 DO 5 I=2,NC 5 I0(I) = I0(I-1) + M(I-1) DO 10 I=1,MJ DO 10 J=1,MJ 10 VF(I,J) = 0.0D0 C DO 70 L=1,NC C XF(I0(L)+1) = XMEAN(L) SUM = 0.0D0 DO 20 I=M(L),2,-1 SUM = SUM + A(I,L) 20 XF(I0(L)+I) = SUM*XMEAN(L) C VF(I0(L)+1,I0(L)+1) = XVAR(1,L) DO 40 I=2,M(L) SUM = 0.0D0 DO 30 J=I,M(L) 30 SUM = SUM + A(J,L)*XVAR(J,L) VF(I0(L)+1,I0(L)+I) = SUM 40 VF(I0(L)+I,I0(L)+1) = SUM DO 60 I=2,M(L) DO 60 J=I,M(L) SUM = 0.0D0 DO 50 I1=I,M(L) DO 50 J1=J,M(L) 50 SUM = SUM + A(I1,L)*A(J1,L)*XVAR(IABS(J1-J-I1+I)+1,L) VF(I0(L)+I,I0(L)+J) = SUM 60 VF(I0(L)+J,I0(L)+I) = SUM 70 CONTINUE C RETURN E N D SUBROUTINE PTSEAS( Y,XSS,VSS,N,NE,MJ,M1,M2,M3,M4,MPER, * TAU2,SIG2,FF,AIC,REG,NMAX ) 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 NE: Final point of smoothing C MJ: Adjustable dimension of F and G C M1-M4: Model orders C MPER: Period C TAU2: System noise variance C SIG2: Observational noise variance C FF: Log-Likelihood of the mode C AIC: AIC of the model C REG: Exogenious variable (number of trading days) C NMAX: Adjustable dimension of REG C IMPLICIT REAL*8(A-H,O-Z) CHARACTER VNAME*8 DIMENSION Y(N), TAU2(10), REG(NMAX,*) DIMENSION XSS(MJ,N), VSS(MJ,MJ,N) DIMENSION DATA(1000), VNAME(20), VALUE(20) C VNAME(1) = 'M1 =' VNAME(2) = 'M2 =' VNAME(3) = 'M3 =' VNAME(4) = 'M4 =' VNAME(5) = 'PERIOD =' VNAME(6) = 'TAU1 =' VNAME(7) = 'TAU2 =' VNAME(8) = 'TAU3 =' VNAME(9) = 'SIG2 =' VNAME(10)= 'LOG-LH =' VNAME(11)= 'AIC =' VALUE(1) = M1 VALUE(2) = M2 VALUE(3) = M3 VALUE(4) = M4 VALUE(5) = MPER VALUE(6) = TAU2(1) VALUE(7) = TAU2(2) VALUE(8) = TAU2(3) VALUE(9) = SIG2 VALUE(10)= FF VALUE(11)= AIC CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 11.1: SEASONAL ADJUSTMENT ',36,11, * VNAME,VALUE ) WY = 5.0 WX = 12.0 DX = 12.0 FN = N IY = 12 CALL MAXMIN( Y,N,YMIN1,YMAX1,DY ) YMIN3 = -DY YMAX3 = DY CALL NEWPEN( 2 ) CALL PLOT( 1.5,17.0,-3 ) C C ... Trend Component ... C IF( M1.GT.0 ) THEN CALL PLOT( 0.0,-SNGL(WY)-1.0,-3 ) CALL SYMBOL( 0.2,SNGL(WY)+0.1,0.2,'TREND',0.0,5 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY,0.0D0,FN,YMIN1,YMAX1,DX,DY, * 0.2D0,1,IY,2 ) DO 120 J=-1,1 DO 100 I=1,NE 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 END IF C C ... Seasonal Component ... C WY2 = -1.0 IF( M2.GT.0 ) THEN DO 130 I=1,NE 130 DATA(I) = XSS(M1+1,I) CALL MAXMIN( DATA,N,DMIN,DMAX,DDY ) YMAX2 = DY*INT( DMAX/DY ) IF( YMAX2.LT.DMAX ) YMAX2 = YMAX2 + DY YMIN2 = DY*INT( DMIN/DY ) IF( YMIN2.GT.DMIN ) YMIN2 = YMIN2 - DY WY2 = WY*(YMAX2-YMIN2)/(YMAX1-YMIN1) C CALL PLOT( 0.0,-SNGL(WY2+1.0),-3 ) CALL NEWPEN( 2 ) CALL SYMBOL( 0.2,SNGL(WY2)+0.1,0.2,'SEASONAL',0.0,8 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY2,0.0D0,FN,YMIN2,YMAX2,DX,DY, * 0.2D0,1,IY,2) C CALL PLOTY( DATA,N,YMIN2,YMAX2,WX,WY2,1,1 ) CALL PLOTB( DATA,N,YMIN2,YMAX2,WX,WY2,1 ) END IF C C ... Noise Component ... C WY3 = WY*(YMAX3-YMIN3)/(YMAX1-YMIN1) CALL PLOT( 0.0,-SNGL(WY3+1.0),-3 ) CALL NEWPEN( 2 ) CALL SYMBOL( 0.2,SNGL(WY3)+0.1,0.2,'IRREGULAR',0.0,9 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY3,0.0D0,FN,YMIN3,YMAX3,DX,DY, * 0.2D0,1,IY,2) C M12 = M1+ M2*(MPER-1) DO 160 I=1,NE DATA(I) = Y(I) - XSS(1,I) IF(M2.GT.0) DATA(I) = DATA(I) - XSS(M1+1,I) IF(M3.GT.0) DATA(I) = DATA(I) - XSS(M12+1,I) IF(M4.GT.0) THEN DO 150 J=1,M4 150 DATA(I) = DATA(I) - REG(I,J)*XSS(M12+M3+J,I) END IF 160 CONTINUE C CALL NEWPEN( 1 ) IF(J.EQ.4) CALL NEWPEN( 2 ) C CALL PLOTY( DATA,N,YMIN3,YMAX3,WX,WY3,1,1 ) CALL PLOTB( DATA,N,YMIN3,YMAX3,WX,WY3,1 ) 180 CONTINUE C C ... PLOT ORIGINAL DATA ... C CALL PLOT( SNGL(WX)+2.0,SNGL(WY2+WY3+2.0),-3 ) CALL NEWPEN( 2 ) CALL SYMBOL( 0.2,SNGL(WY)+0.1,0.2,'OBSERVED',0.0,8 ) 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 CALL PLOTY5( Y,N,YMIN1,YMAX1,WX,WY,1,3 ) C C ... AR Component ... C IF( M3.GT.0 ) THEN DO 210 I=1,NE 210 DATA(I) = XSS(M12+1,I) C CALL MAXMIN( DATA,N,DMIN,DMAX,DDY ) YMAX4 = DY*INT( DMAX/DY ) IF( YMAX4.LT.DMAX ) YMAX4 = YMAX4 + DY YMIN4 = DY*INT( DMIN/DY ) IF( YMIN4.GT.DMIN ) YMIN4 = YMIN4 - DY WY4 = WY*(YMAX4-YMIN4)/(YMAX1-YMIN1) CALL PLOT( 0.0,-SNGL(WY4+1.0),-3 ) CALL NEWPEN( 2 ) CALL SYMBOL( 0.2,SNGL(WY4)+0.1,0.2,'STATIONARY',0.0,10 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY4,0.0D0,FN,YMIN4,YMAX4,DX,DY, * 0.2D0,1,IY,2) C CALL PLOTY( DATA,N,YMIN4,YMAX4,WX,WY4,1,1 ) CALL PLOTB( DATA,N,YMIN4,YMAX4,WX,WY4,1 ) 230 CONTINUE END IF C C ... Trading Day Effect ... C IF( M4.GT.0 ) THEN CALL PLOT( 0.0,-SNGL(WY3+1.0),-3 ) CALL NEWPEN( 2 ) CALL SYMBOL( 0.2,SNGL(WY3)+0.1,0.2,'TRADING DAY EFFECT',0.0,18) CALL AXISXY( 0.0D0,0.0D0,WX,WY3,0.0D0,FN,YMIN3,YMAX3,DX,DY, * 0.2D0,1,IY,2) C DO 250 I=1,NE SUM = 0.0D0 DO 240 J=1,M4 240 SUM = SUM + REG(I,J)*XSS(M12+M3+J,I) 250 DATA(I) = SUM C C CALL PLOTY( DATA,N,YMIN3,YMAX3,WX,WY3,1,1 ) CALL PLOTB( DATA,N,YMIN3,YMAX3,WX,WY3,1 ) END IF C CALL PLOTE C call plot( 0.0,0.0,999 ) RETURN E N D SUBROUTINE PLOTB( Y,N,YMIN,YMAX,WX,WY,IOFF ) REAL*8 Y, YMIN, YMAX, WX, WY DIMENSION Y(N) C DX = WX/(N + IOFF-1) DY = WY/(YMAX - YMIN) C DO 100 I=1,N XX = DX*(I + IOFF -1) YY = (Y(I) - YMIN)*DY Y0 = (0.0 - YMIN)*DY IF( YY.GT.WY ) THEN CALL SYMBOL( XX,SNGL(WY),0.2,A,0.0,-1 ) YY = WY END IF IF( YY.LT.0.0 ) THEN CALL SYMBOL( XX,0.0,0.2,A,0.0,-1 ) YY = 0.0 END IF CALL PLOT( XX,Y0,3 ) CALL PLOT( XX,YY,2 ) 100 CONTINUE C RETURN E N D SUBROUTINE PLOTY5( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C MODIFIED 8/31/90 REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) CALL NEWPEN( 1 ) C 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 ) C IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) C IF(ISW.GE.3) CALL CIRC1( XX,YY,0.05,1 ) END IF IF(I.GT.40.AND.I.LE.70) CALL CIRC1( XX,YY,0.05,1 ) IF(I.GT.100.AND.I.LE.120) CALL CIRC1( XX,YY,0.05,1 ) 100 CONTINUE C RETURN E N D SUBROUTINE SETABC( M,L,NC,MTYPE,TAU2,MAXM,MM,AR,Y,N,A,B,C,Q, * XMEAN,XVAR ) C C ... State space model for Seasonal Adjustment ... C C Input: C M: Order of the trend C L: Number of components C NC: Adjustable dimension C MTYPE: Model type C TAU2: System noise variance C MAXM: Adjustable dimension of A, B, C C MM: State dimension C AR: AR coefficient C Y(I): Time series C N: Data length C Outputs: C A,B,C: Parameters of F, G, H C Q: System noise covarance matrix C XMEAN: Initial value of each component C XVAR: Initial variance of each component C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(MAXM,NC), B(MAXM,NC), C(MAXM,NC), Q(NC,NC) DIMENSION M(NC), TAU2(NC), AR(NC), MTYPE(NC) DIMENSION Y(N), XMEAN(NC), XVAR(MAXM,NC), DUM(1), COV(0:10) COMMON /C92827/ M1, M2, M3, M4, MPER C L = 0 CALL MOMENT( Y,N/4,YMEAN,YVAR ) YVAR = 1.0D4 C DO 10 J=1,NC DO 10 I=1,MAXM XVAR(I,J) = 0.0D0 C(I,J) = 0.0D0 A(I,J) = 0.0D0 10 B(I,J) = 0.0D0 DO 20 I=1,NC DO 20 J=1,NC 20 Q(I,J) = 0.0D0 C C ... TREND MODEL ... C IF( M1.GT.0 ) THEN L = 1 M(L) = M1 MTYPE(L) = 0 IF( M1.EQ.1 ) THEN A(1,1) = 1.0D0 END IF IF( M1.EQ.2 ) THEN A(1,1) = 2.0D0 A(2,1) = -1.0D0 END IF IF( M1.EQ.3 ) THEN A(1,1) = 3.0D0 A(2,1) = -3.0D0 A(3,1) = 1.0D0 END IF B(1,1) = 1.0D0 C(1,1) = 1.0D0 Q(1,1) = TAU2(1) XMEAN(1) = YMEAN XVAR(1,1) = YVAR END IF C C ... SEASONAL MODEL ... C IF( M2.GT.0 ) THEN L = L+1 M(L) = M2*(MPER-1) MTYPE(L) = 0 IF( M2.NE.2 ) THEN DO 30 I=1,M(L) 30 A(I,L) = -1.0D0 END IF IF( M2.EQ.2 ) THEN DO 40 I=1,MPER-1 A(I,L) = -(I+1) 40 A(M(L)-I+1,L) = -I END IF B(1,L) = 1.0D0 C(1,L) = 1.0D0 Q(L,L) = TAU2(L) XMEAN(L) = 0.0D0 XVAR(1,L) = YVAR END IF C C ... AR MODEL ... C IF( M3.GT.0 ) THEN L = L+1 M(L) = M3 MTYPE(L) = 0 DO 50 I=1,M3 50 A(I,L) = AR(I) B(1,L) = 1.0D0 C(1,L) = 1.0D0 Q(L,L) = TAU2(L) CALL ARMCOV( M3,0,AR,DUM,1.0D0,M3,COV ) XMEAN(L) = 0.0D0 DO 60 I=1,M3 60 XVAR(I,L) = COV(I-1) END IF C C ... TRADING DAY EFFECT MODEL ... C IF( M4.GT.0 ) THEN DO 70 I=L+1,L+6 M(I) = 1 MTYPE(I) = I-L XMEAN(I) = 0.0D0 XVAR(1,I) = YVAR A(1,I) = 1.0D0 B(1,I) = 1.0D0 C(1,I) = 1.0D0 70 Q(I,I) = 0.0D0 L = L + 6 END IF C C ... STATE DIMENSION ... C MSUM = 0 DO 80 I=1,L 80 MSUM = MSUM + M(I) MM = MSUM C RETURN E N D SUBROUTINE FFSEAS( K,AA,FF,IFG ) C C ... Function for Maximum likelihood estimation (seasonal model) ... C C Inputs: C K: Number of parameters C AA(I): Parameter vector C Outputs: C FF: -(Log likelihood) C IFG: =1 if some conditions are violated C IMPLICIT REAL*8(A-H,O-Z) DIMENSION M(10), TAU2(10), PAR(10) DIMENSION XMEAN(10), XVAR(22,10), AA(20), AR(10) DIMENSION A(22,9), B(22,9), C(22,9), Q(9,9) DIMENSION XPS(25,300), XFS(25,300) DIMENSION VPS(25,25,300), VFS(25,25,300) DIMENSION XF(40), VF(40,40), MTYPE(10) COMMON /C92827/ M1, M2, M3, M4, NPER COMMON /C92826/ Y(160), REG(160,7) COMMON /C92825/ OUTMIN, OUTMAX COMMON /C92907/ ALIMIT COMMON /CMMODL/ M, XMEAN, XVAR, N, NS, NFE, NPE NMAX= 300 MJ = 25 MAXM= 22 NC = 9 NP = ID(M1) + ID(M2) + ID(M3) C DO 10 I=1,K 10 IF( DABS(AA(I)).GT.20.0D0 ) GO TO 100 DO 20 I=1,NP 20 TAU2(I) = DEXP( AA(I) )/(1.0D0 + DEXP( AA(I) )) DO 30 I=1,M3 30 PAR(I) = ALIMIT*(DEXP(AA(NP+I))-1.0D0)/(DEXP(AA(NP+I))+1.0D0) CALL ARCOEF( PAR,M3,AR ) SIG2 = 1.0D0 IFG = 0 C CALL SETABC( M,L,NC,MTYPE,TAU2,MAXM,MM,AR,Y,N,A,B,C,Q, * XMEAN,XVAR ) CALL ISTAT1( L,M,MJ,MAXM,A,XMEAN,XVAR,XF,VF ) CALL FILTR1( Y,XF,VF,A,B,C,Q,SIG2,REG,MTYPE,NMAX,M,MAXM,NC,L, * NS,NFE,NPE,MJ,N,OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,OVAR ) WRITE(6,*) FF FF = -FF RETURN C 100 IFG = 1 FF = 1.0D20 RETURN E N D SUBROUTINE TRADE( JYEAR,MONTH,N,MJ,TDAY ) C C ... This subroutine computes the number of days of the week C in each month, C C Inputs: C JYEAR: Starting year of the data set C MONTH: Starting month of the data set C N: Data length C MJ: Adjustable dimension of TDAY C Output: C TDAY(I,J): Number of J-th days of the week in I-th month C Written by F.K. Nov.1981 C REAL*8 TDAY DIMENSION TDAY(MJ,7), IX(12) DATA IX /3,0,3,2,3,2,3,3,2,3,2,3/ C JS = JYEAR - 1900 I0 = MOD( JS+(JS-1)/4,7 ) + 1 JJ = 2 - MONTH II = 0 5 II = II + 1 I1 = II + JS - 1 IX(2) = 0 IF( MOD(I1,4) .EQ. 0 ) IX(2) = 1 IF( MOD(I1,100).EQ.0 ) IX(2) = 0 IF( MOD(I1,400).EQ.0 ) IX(2) = 1 DO 30 J=1,12 DO 10 I=1,7 10 IF( JJ.GT.0 ) TDAY(JJ,I) = 4 C IE = IX(J) IF( IE .EQ. 0 ) GO TO 30 DO 20 I=1,IE I2 = I0 + I IF( I2 .GT. 7 ) I2 = I2 - 7 20 IF( JJ.GT.0 ) TDAY(JJ,I2) = 5 I0 = I2 30 JJ = JJ + 1 IF( JJ .GT.N ) RETURN GO TO 5 C E N D FUNCTION ID( K ) ID = 0 IF( K.GT.0 ) ID = 1 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 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 IMPULS( M,L,A,B,K,G ) C C ... Impulse Response Function ... C C Inputs: C M: AR order C L: MA order C A(I): AR coefficient C B(I): MA coefficient C K: Required maximum lag of impulse respose C Output: C G(I): Impulse response function C Y.I. IMPLICIT REAL*8( A-H,O-Z ) DIMENSION A(*), B(*), G(0:K) C G(0) = 1.0 DO 20 I=1,K SUM = 0.0D0 IF(I.LE.L) SUM = -B(I) DO 10 J=1,I 10 IF(J.LE.M) SUM = SUM + A(J)*G(I-J) 20 G(I) = SUM C RETURN E N D SUBROUTINE ARMCOV( M,L,A,B,SIG2,K,COV ) C C ... Autocovariance Function of ARMA model ... C C Inputs: C M: AR order C L: MA order C A(I): AR coefficient C B(I): MA coefficient C SIG2: innovation variance C K: Required maximum lag of autocovariance C Output: C COV(I): Autocovariance function C Y.I. IMPLICIT REAL*8( A-H,O-Z ) DIMENSION A(*), B(*), COV(0:K), G(0:100), X(30,30) DIMENSION Z(100), UL(30,30), IPS(100) C KMAX = MAX(M,L,K) CALL IMPULS( M,L,A,B,KMAX,G) C DO 10 I=1,M+1 DO 10 J=1,M+1 10 X(I,J) = 0.0D0 DO 20 I=1,M+1 20 X(I,I) = 1.0D0 DO 30 I=1,M DO 30 J=2,M-I+2 30 X(I,J) = X(I,J) - A(I+J-2) DO 40 I=2,M+1 DO 40 J=1,I-1 40 X(I,J) = X(I,J) - A(I-J) C CALL DECOM( M+1,X,30,UL,IPS ) C SUM = 1.0D0 DO 50 J=1,L 50 SUM = SUM - B(J)*G(J) Z(1)= SIG2*SUM DO 70 I=2,M+1 SUM = 0.0D0 DO 60 J=I-1,L 60 SUM = SUM - B(J)*G(J-I+1) 70 Z(I) = SIG2*SUM C CALL SOLVE( M+1,UL,30,Z,COV,IPS) C DO 100 J=M+1,K SUM = 0.0D0 DO 80 I=1,M 80 SUM = SUM + A(I)*COV(J-I) DO 90 I=J,L 90 SUM = SUM - B(I)*G(I-J)*SIG2 100 COV(J) = SUM C RETURN E N D SUBROUTINE DECOM( N,A,MJ,UL,IPS ) C C ... UL decomposition: A = L*U ... C C Inputs: C N: Dimension of the matrix A C A(I,J): N*N positive definite matrix C MJ: Adjustable dimension of A and UL C Outputs: C UL(I,J): L-I and U C IPS: Index vector C Y.I. IMPLICIT REAL*8(A-H,O-Z ) DIMENSION A(MJ,*),UL(MJ,*),SCALES(100),IPS(100) C DO 20 I=1,N IPS(I) = I RNORM = 0.0D0 DO 10 J=1,N UL(I,J) = A(I,J) IF(RNORM.LT.ABS(UL(I,J)) ) RNORM = ABS( UL(I,J) ) 10 CONTINUE IF( RNORM .NE. 0.0D0 ) THEN SCALES(I) = 1/RNORM ELSE SCALES(I) = 0.0D0 CALL SING(0) END IF 20 CONTINUE C DO 60 K=1,N-1 BIG = 0.0D0 DO 30 I=K,N SIZE = ABS( UL(IPS(I),K) )*SCALES( IPS(I) ) IF( BIG.LT.SIZE ) THEN BIG = SIZE INDEX = I END IF 30 CONTINUE IF( BIG.EQ. 0.0D0 ) THEN CALL SING(1) GO TO 60 END IF IF( INDEX.NE.K ) THEN J = IPS(K) IPS(K) = IPS(INDEX) IPS(INDEX) = J END IF C PIVOT = UL(IPS(K),K) DO 50 I=K+1,N TM = UL( IPS(I),K)/PIVOT UL( IPS(I),K) = TM IF( TM.NE. 0.0D0 ) THEN DO 40 J = K+1,N 40 UL( IPS(I),J ) = UL( IPS(I),J)-TM*UL( IPS(K),J) C WRITE(6,*) (UL(IPS(I),J),J=1,N) END IF 50 CONTINUE 60 CONTINUE C IF( UL(IPS(N),N) .EQ. 0.0D0 ) CALL SING(2) RETURN E N D SUBROUTINE SOLVE( N,UL,MJ,B,X,IPS ) C C ... Solve Ax=b using UL obtained by DECOM ... C C Inputs: C N: Dimension of UL and B C UL: LU decomposition of A C MJ: Adjustable dimension of A C B: C IPS: index vector C Output: C X: Solution C Y.I. IMPLICIT REAL*8( A-H,O-Z ) DIMENSION UL(MJ,*),B(*),X(*),IPS(100) C DO 20 I=1,N SUM = 0.0D0 DO 10 J=1,I-1 10 SUM = SUM + UL(IPS(I),J)*X(J) 20 X(I) = B(IPS(I)) - SUM C DO 40 I=N,1,-1 SUM = 0.0D0 DO 30 J=I+1,N 30 SUM = SUM + UL(IPS(I),J)*X(J) 40 X(I) = ( X(I)-SUM )/UL(IPS(I),I) RETURN 600 FORMAT(1H ,'N=',I10,/,(5X,'IPS=',I10 ) ) E N D SUBROUTINE SING( IFG ) C C ... ERROR MESSAGE ... C IF( IFG.EQ.0) THEN WRITE(6,100) ELSE IF( IFG.EQ.1 ) THEN WRITE(6,200) ELSE WRITE(6,300) END IF END IF RETURN 100 FORMAT(' MATRIX WITH ZERO ROW IN DECOMPOSE.') 200 FORMAT(' SINGULAR MATRIX IN DECOMPOSE.ZERO DIDIVIDE IN SOLVE.') 300 FORMAT(' CONVERGENCE IN IMPRUV.MATRIX IS NEARLY SINGULAR.') E N D SUBROUTINE ARCOEF( PAR,K,A ) C C ... This subroutine computes AR coefficients from PARCOR ... C C Inputs: C PAR(I): PARCOR C K: Order of AR model C Output: C A(I): AR coefficient C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PAR(K), A(K), AA(50) C DO 30 II=1,K A(II) = PAR(II) AA(II) = PAR(II) IF( II-1.LE.0 ) GO TO 30 DO 10 J=1,II-1 10 A(J) = AA(J) - PAR(II)*AA(II-J) IF( II.LT.K ) THEN DO 20 J=1,II-1 20 AA(J) = A(J) END IF 30 CONTINUE RETURN E N D SUBROUTINE PARCOR( A,K,PAR ) C C ... This subriutine computes PARCOR for AR coefficients ... C C Inputs: C A: Vector of AR coefficients C K: Order of the model C C Output: C PAR: Vector of partial autocorrelations (PARCOR) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(K), PAR(K), G(50) C DO 10 I=1,K 10 PAR(I) = A(I) C DO 40 II=K-1,1,-1 S = 1.D0 - PAR(II+1)**2 DO 20 I=1,II 20 G(I) = (PAR(I) + PAR(II+1)*PAR(II-I+1))/S I2 = (II+1)/2 IF( MOD( II,2 ).EQ.1 ) G(I2) = PAR(I2)/(1.D0 - PAR(II+1)) DO 30 I=1,II 30 PAR(I) = G(I) 40 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 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 DAVIDN( FUNCT,X,N,NDIF ) C C ... 6/20/83, 12/19/92 C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION X(40), DX(40), G(40), G0(40), Y(40) DIMENSION H(40,40), WRK(40), S(40) COMMON / CCC / ISW, IPR, ISMT, IDIF COMMON / DDD / XM , AIC , SD EXTERNAL FUNCT DATA TAU2 / 1.0D-6 / DATA EPS1, EPS2 / 1.0D-6 , 1.0D-6 / RAMDA = 0.5D0 CONST1 = 1.0D-30 IPR = 0 C C INITIAL ESTIMATE OF INVERSE OF HESSIAN C ICOUNT = 0 1000 DO 20 I=1,N DO 10 J=1,N 10 H(I,J) = 0.0D0 S(I) = 0.0D0 DX(I) = 0.0D0 20 H(I,I) = 1.0D0 ISW = 0 C IF( NDIF.EQ.0 ) CALL FUNCT( N,X,XM,G,IG ) IF( NDIF.GE.1 ) CALL FUNCND( FUNCT,N,X,XM,G,IG ) C IF( IPR .GE. 2 ) WRITE( 6,640 ) XM, SD, AIC WRITE(6,650) (X(I),I=1,N), XM, (G(I),I=1,N) C ICC = 0 C ITERATION 2000 CONTINUE ICC = ICC + 1 IF( ICC .EQ. 1 ) GO TO 120 C DO 40 I=1,N 40 Y(I) = G(I) - G0(I) DO 60 I=1,N SUM = 0.0D0 DO 50 J=1,N 50 SUM = SUM + Y(J)*H(I,J) 60 WRK(I) = SUM S1 = 0.0D0 S2 = 0.0D0 DO 70 I=1,N S1 = S1 + WRK(I)*Y(I) 70 S2 = S2 + DX(I) *Y(I) IF( S1.LE.CONST1 .OR. S2.LE.CONST1 ) GO TO 900 C IF( S1 .LE. S2 ) GO TO 100 C C UPDATE THE INVERSE OF HESSIAN MATRIX C C --- BROYDEN-FLETCHER-GOLDFARB-SHANNO TYPE CORRECTION - C 100 CONTINUE STEM = S1 / S2 + 1.0D0 DO 110 I=1,N DO 110 J=I,N H(I,J) = H(I,J)-(DX(I)*WRK(J)+WRK(I)*DX(J)-DX(I)*DX(J)*STEM)/S2 110 H(J,I) = H(I,J) C 120 CONTINUE SS = 0.0D0 DO 140 I=1,N SUM = 0.0D0 DO 130 J=1,N 130 SUM = SUM + H(I,J)*G(J) SS = SS + SUM**2 140 S(I) = -SUM C S1 = 0.0D0 S2 = 0.0D0 DO 150 I=1,N S1 = S1 + S(I)*G(I) 150 S2 = S2 + G(I)*G(I) C DS2 = DSQRT(S2) C GTEM = DABS(S1)/DS2 C IF( GTEM.LE.TAU1 .AND. DS2.LE.TAU2 ) GO TO 900 IF( S1.LT.0.0D0 ) GO TO 200 DO 170 I=1,N DO 160 J=1,N 160 H(I,J) = 0.0D0 H(I,I) = 1.0D0 170 S(I) = -S(I) 200 CONTINUE C ED = XM C C LINEAR SEARCH C CALL LINEAR( FUNCT,X,S,RAMDA,ED,N,IG ) C IF( IPR .GE. 2 ) WRITE( 6,630 ) RAMDA, ED, SD, AIC C S1 = 0.0D0 DO 210 I=1,N DX(I) = S(I)*RAMDA S1 = S1 + DX(I)*DX(I) G0(I) = G(I) 210 X(I) = X(I) + DX(I) XMB = XM ISW = 0 C IF( NDIF.EQ.0 ) CALL FUNCT( N,X,XM,G,IG ) IF( NDIF.GE.1 ) CALL FUNCND( FUNCT,N,X,XM,G,IG ) WRITE(6,650) (X(I),I=1,N), XM, (G(I),I=1,N) C S2 = 0.D0 DO 220 I=1,N 220 S2 = S2 + G(I)*G(I) IF( DSQRT(S2) .LT. TAU2 ) GO TO 900 IF( XMB/XM-1.D0.LT.EPS1 .AND. DSQRT(S1).LT.EPS2 ) GO TO 900 C IF( ICC .GE. 5 ) GO TO 900 GO TO 2000 900 CONTINUE IF( IPR .LE. 0 ) RETURN WRITE( 6,600 ) WRITE( 6,610 ) (X(I),I=1,N) WRITE( 6,620 ) WRITE( 6,610 ) (G(I),I=1,N) C ICOUNT = ICOUNT + 1 S2 = 0.D0 DO 910 I=1,N 910 S2 = S2 + G(I)*G(I) IF( S2.GT.1.0D0.AND.ICOUNT.LE.1 ) GO TO 1000 RETURN 600 FORMAT( 1H0,'----- X -----' ) 610 FORMAT( 1H ,10D13.5 ) 620 FORMAT( 1H0,'*** GRADIENT ***' ) 630 FORMAT( 1H ,'LAMBDA =',D15.7,3X,'(-1)LOG LIKELIHOOD =',D23.15, * 3X,'SD =',D22.15,5X,'AIC =',D23.15 ) 640 FORMAT( 1H ,26X,'(-1)LOG-LIKELIHOOD =',D23.15,3X,'SD =',D22.15, * 5X,'AIC =',D23.15 ) 650 FORMAT( 10X,5F12.7 ) E N D SUBROUTINE FUNCND( FUNCT,M,A,F,G,IFG ) C C ... FUNCTION EVALUATION AND NUMERICAL DIFFERENCING ... C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION A(M) , G(M) , B(20),GD(5) COMMON / CCC / ISW , IPR, ISMT, IDIF COMMON /CMFUNC/ DJACOB,FC,SIG2,AIC,FI,SIG2I,AICI,GI(20),GC(20) C DATA ICNT /0/ CONST = 0.00001D0 C CALL FUNCT( M,A,F,GD,IFG ) FB = F IF( ISW .GE. 1 ) RETURN C C WRITE( 6,600 ) (A(I),I=1,M) DO 10 I=1,M 10 B(I) = A(I) C DO 30 II=1,M B(II) = A(II) + CONST CALL FUNCT( M,B,FF,GD,IFG ) IF( IDIF .EQ. 1 ) GO TO 20 B(II) = A(II) - CONST CALL FUNCT( M,B,FB,GD,IFG ) 20 G(II) = (FF-FB)/(CONST*IDIF) IF( G(II) .GT. 1.0D20 ) G(II) = (F-FB)/CONST IF( G(II) .LT.-1.0D20 ) G(II) = (FF-F)/CONST IF( FB.GT.F .AND. FF.GT.F ) G(II) = 0.0D0 30 B(II) = A(II) C RETURN C 600 FORMAT( 3X,'--- PARAMETER ---',(/,3X,5D13.5) ) 610 FORMAT( 3X,'--- GRADIENT ---',(/,3X,5D13.5) ) E N D SUBROUTINE LINEAR( FUNCT,X,H,RAM,EE,K,IG ) C C ... LINEAR SEARCH ... C IMPLICIT REAL*8( A-H,O-Z ) INTEGER RETURN,SUB DIMENSION X(1), H(1), X1(40) DIMENSION G(40) COMMON / CCC / ISW , IPR, ISMT, IDIF C IPR = 7 C C1 = 1.0D-7 C C2 = 1.0D-30 C ISW = 1 IG = 0 RAM = 0.1D0 CONST2 = 1.0D-60 14 DO 15 I=1,K 15 IF( DABS(H(I)) .GT. 1.0D10 ) GO TO 16 GO TO 18 16 DO 17 I=1,K 17 H(I) = H(I)*1.0D-10 GO TO 14 18 CONTINUE HNORM = 0.D0 DO 10 I=1,K 10 HNORM = HNORM + H(I)**2 HNORM = DSQRT( HNORM ) C RAM2 = RAM E1 =EE RAM1 = 0.D0 C DO 20 I=1,K 20 X1(I) = X(I) + RAM2*H(I) CALL FUNCT( K,X1,E2,G,IG ) IF(IPR.GE.7) WRITE(6,2) RAM2,E2 C IF( IG .EQ. 1 ) GO TO 50 IF( E2 .GT. E1 ) GO TO 50 30 RAM3 = RAM2*4.D0 DO 40 I=1,K 40 X1(I) = X(I) + RAM3*H(I) CALL FUNCT( K,X1,E3,G,IG ) IF( IG.EQ.1 ) GO TO 500 IF( IPR.GE.7 ) WRITE(6,3) RAM3,E3 IF( E3 .GT. E2 ) GO TO 70 IF(RAM3.GT.1.0D10 .AND. E3.LT.E1) GO TO 45 IF(RAM3.GT.1.0D10 .AND. E3.GE.E1) GO TO 46 RAM1 = RAM2 RAM2 = RAM3 E1 = E2 E2 = E3 GO TO 30 C 45 RAM = RAM3 EE = E3 RETURN C 46 RAM = 0.0D0 RETURN C 50 RAM3 = RAM2 E3 = E2 RAM2 = RAM3*0.1D0 IF( RAM2*HNORM .LT. CONST2 ) GO TO 400 DO 60 I=1,K 60 X1(I) = X(I) + RAM2*H(I) CALL FUNCT( K,X1,E2,G,IG ) IF(IPR.GE.7) WRITE(6,4) RAM2,E2 IF( E2.GT.E1 ) GO TO 50 C 70 ASSIGN 80 TO RETURN GO TO 200 C 80 DO 90 I=1,K 90 X1(I) = X(I) + RAM*H(I) CALL FUNCT( K,X1,EE,G,IG ) IF(IPR.GE.7) WRITE(6,5) RAM,EE C IFG = 0 ASSIGN 300 TO SUB ASSIGN 200 TO SUB 95 ASSIGN 130 TO RETURN IF( RAM .GT. RAM2 ) GO TO 110 IF( EE .GE. E2 ) GO TO 100 RAM3 = RAM2 RAM2 = RAM E3 =E2 E2 =EE GO TO SUB,( 200,300 ) C 100 RAM1 = RAM E1 = EE GO TO SUB,( 200,300 ) C 110 IF( EE .LE. E2 ) GO TO 120 RAM3 = RAM E3 = EE GO TO SUB,( 200,300 ) C 120 RAM1 = RAM2 RAM2 = RAM E1 = E2 E2 = EE GO TO SUB,( 200,300 ) C 130 DO 140 I=1,K 140 X1(I) = X(I) + RAM*H(I) CALL FUNCT( K,X1,EE,G,IG ) IF( IPR.GE.7 ) WRITE(6,6) RAM,EE ASSIGN 200 TO SUB IFG = IFG+1 600 FORMAT( 1H ,6D20.13 ) C IF( HNORM*(RAM3-RAM1) .GT. 1.0D-12 ) GO TO 95 IF( RAM3-RAM1 .GT. RAM2*1.0D-1 ) GO TO 95 C IF( E1-EE .GT. 1.0D-15 ) GO TO 95 C IF( E3-EE .GT. 1.0D-15 ) GO TO 95 C IF(DX.LE.C1*XN+C2 .AND. DF.LE.C1*DABS(EE)+C2) RETURN C IF( IFG .LE. 2 ) GO TO 95 C IF( E2 .LT. EE ) RAM = RAM2 RETURN C C ------- INTERNAL SUBROUTINE SUB1 ------- 200 IF( RAM3-RAM2 .GT. 5.0D0*(RAM2-RAM1) ) GO TO 202 IF( RAM2-RAM1 .GT. 5.0D0*(RAM3-RAM2) ) GO TO 204 A1 = (RAM3-RAM2)*E1 A2 = (RAM1-RAM3)*E2 A3 = (RAM2-RAM1)*E3 B2 = (A1+A2+A3)*2.D0 B1 = A1*(RAM3+RAM2) + A2*(RAM1+RAM3) + A3*(RAM2+RAM1) IF( B2 .EQ. 0.D0 ) GO TO 210 RAM = B1 /B2 C IF( RAM .LE. RAM1 ) RAM = (RAM1 + RAM2)/2.0D0 IF( RAM .GE. RAM3 ) RAM = (RAM2 + RAM3)/2.0D0 IF( DABS(RAM-RAM2) .LE. 1.0D-15 ) RAM = (RAM2*4.D0 + RAM3)/5.0D0 GO TO RETURN ,( 80,130 ) 202 RAM = (4.0D0*RAM2 + RAM3)/5.0D0 GO TO RETURN, (80,130) 204 RAM = (RAM1 + 4.0D0*RAM2)/5.0D0 GO TO RETURN, (80,130) C 210 IG = 1 RAM = RAM2 RETURN C C ------- INTERNAL SUBROUTINE SUB2 ------- C 300 IF( RAM3-RAM2 .GT. RAM2-RAM1 ) GO TO 310 RAM = (RAM1+RAM2)*0.5D0 GO TO RETURN ,( 80,130 ) C 310 RAM = (RAM2+RAM3)*0.5D0 GO TO RETURN ,( 80,130 ) C 400 RAM = 0.D0 RETURN C 500 RAM = (RAM2+RAM3)*0.5D0 510 DO 520 I=1,K 520 X1(I) = X(I) + RAM*H(I) CALL FUNCT( K,X1,E3,G,IG ) IF( IPR.GE.7 ) WRITE(6,7) RAM,E3 IF( IG.EQ.1 ) GO TO 540 IF( E3.GT.E2 ) GO TO 530 RAM1 = RAM2 RAM2 = RAM E1 = E2 E2 = E3 GO TO 500 C 530 RAM3 = RAM GO TO 70 C 540 RAM = (RAM2+RAM)*0.5D0 GO TO 510 C C ------------------------------------------------------------ 1 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E1 =',D25.17 ) 2 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E2 =',D25.17 ) 3 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E3 =',D25.17 ) 4 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E4 =',D25.17 ) 5 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E5 =',D25.17 ) 6 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E6 =',D25.17 ) 7 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E7 =',D25.17 ) E N D