SUBROUTINE D2MNB(TITLE,FI,FY,P,XA,XY,DLT,V) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER IDIM PARAMETER (IDIM=50) INTEGER I,J,K,L,S,MAXITR,YR(IDIM),IPRNT CHARACTER*128 TITLE DIMENSION THETA(3),THETA1(3),H(3,3),W(3,IDIM),V(3,3), , G(3),XA(IDIM),XY(IDIM),DLT(IDIM),X(3,IDIM),XM(IDIM),COVPM(IDIM), , COV(IDIM,IDIM) COMMON/BLK1/M(IDIM),N(IDIM),NY(IDIM),B(IDIM),R(IDIM),RY(IDIM), , NS(IDIM),SY(IDIM),YR COMMON /BLK3/S,FIMA,FIMD,FYMA,FYMD,EPSLON,MAXITR,IPRNT C C THIS SUBROUTINE COMPUTES ESTIMATES OF M,N, AND B FOR Model D2. C C STATEMENT FUNCTIONS U(I)=(N(I)-M(I))/P XN(I)=XM(I)+U(I) XB(I)=U(I+1)-FI*Q*U(I)+SY(I)*FY C Q=.1D1-P XM(S)=(M(S)+B(S))/P AVG=XM(S) DO 30 J=1,S-1 DO 10 I=1,3 10 X(I,J)=.0D0 XM(J)=(M(J)+B(J))/(.1D1-Q*XA(J)) AVG=AVG+XM(J) DO 20 I=1,J X(1,J)=X(1,J)+(Q*FI)**(J-I) /(.1D1-XA(I))*(DLT(I)*NS(I)/XA(I)+ + FY*Q*DLT(I+1)*SY(I)/XY(I)+P*Q*DLT(I)*XM(I)/(.1D1-Q*XA(I))) X(2,J)=X(2,J)+(Q*FI)**(J-I) /(.1D1-Q*XA(I+1))*((DLT(I)- - (.1D1-XA(I))/P)*(NS(I)/XA(I)+FY/FI*SY(I)/XY(I))/Q+ + DLT(I)*XM(I)*P/(.1D1-Q*XA(I))) 20 X(3,J)=X(3,J)+(Q*FI)**(J-I) *SY(I)/XY(I) X(1,J)=X(1,J)*Q*XA(J+1) X(2,J)=-X(2,J)*Q*XA(J+1) 30 X(3,J)=X(3,J)*Q*XA(J+1) C CALL MMULT(V,X,W,3,3,S-1,3,3,IDIM) C---------------------------------------------------------- C COMPUTE VAR & COV OF M IF(PRNT.GT.1)WRITE(7,50)TITLE 50 FORMAT(/1X,A128//' Model D2 M Variance-Covariance matrix ', ,'(i=2,s;j=i,s)') AVGSE=.0D0 DO 90 I=2,S QXAI=.1D1-Q*XA(I) PPI=DLT(I)-QXAI/P COVPM(I)=W(2,I-1)-Q*DLT(I)*XM(I)*V(1,2)/(.1D1-XA(I-1))+ + PPI*XM(I)*V(2,2)/QXAI DO 70 J=I,S QXAJ=.1D1-Q*XA(J) PPJ=DLT(J)-QXAJ/P COV(I,J)=Q*XA(J)*XM(I)/QXAI*(Q*FI)**(J-I)- - XM(J)/QXAJ*(Q*DLT(J)/FI*W(1,I-1)-PPJ*W(2,I-1))- - Q*DLT(I)*XM(I)/FI/QXAI*(W(1,J-1)-Q*DLT(J)*XM(J)/FI/QXAJ* * V(1,1)+PPJ*XM(J)*V(1,2)/QXAJ)+ + PPI*XM(I)/QXAI*(W(2,J-1)-Q*DLT(J)*XM(J)/FI/QXAJ*V(1,2)+ + PPJ*XM(J)*V(2,2)/QXAJ) 70 AVGSE=AVGSE+.2D1*COV(I,J) IF(PRNT.GT.1)WRITE(7,80)I,(COV(I,J),J=I,S) 80 FORMAT(' i=',I4,10(1X,F11.2)/7X,10(1X,F11.2)) 90 AVGSE=AVGSE-COV(I,I) WRITE(7,60) 60 FORMAT(/38X,'Standard'/' Period',11X,'M',8X,'Variance',6X,'error' ,,7X,'95% confidence interval Covariance(p,M(i))'/1X,96('-')) DO 92 I=2,S SE=XSQRT(COV(I,I)) C1=XM(I)-.196D1*SE C2=XM(I)+.196D1*SE WRITE(7,100)YR(I),XM(I),COV(I,I),SE,C1,C2,COVPM(I) 92 IF(MOD(YR(I),5).EQ.0)WRITE(7,100) 100 FORMAT(I5,4X,F11.2,3X,F12.2,2X,F8.2,6X,F11.2,' - ',F11.2,7X,F15.2) AVG=AVG/DBLE(FLOAT(S-1)) AVGSE=AVGSE/DBLE(FLOAT(S-1)) SE=XSQRT(AVGSE) C1=AVG-.196D1*SE C2=AVG+.196D1*SE WRITE(7,110)AVG,AVGSE,SE,C1,C2 110 FORMAT(1X,76('-')/' MEAN',4X,F11.2,3X,F12.2,2X,F8.2,6X,F11.2, ,' - ',F11.2/1X,76('-')) C---------------------------------------------------------- C COMPUTE VAR & COV OF N IF(PRNT.GT.1)WRITE(7,120) 120 FORMAT(/' Model D2 N Variance-Covariance matrix (i=2,s;j=i,s)') AVG=.0D0 AVGSE=.0D0 DO 150 I=2,S DO 140 J=I,S COV(I,J)=COV(I,J)-U(I)/P*COVPM(J)-U(J)/P*COVPM(I)+ + U(I)*U(J)/P/P*V(2,2) 140 AVGSE=AVGSE+.2D1*COV(I,J) AVGSE=AVGSE-COV(I,I) AVG=AVG+XN(I) COV(I,I)=COV(I,I)+U(I)*Q/P 150 IF(PRNT.GT.1)WRITE(7,80)I,(COV(I,J),J=I,S) WRITE(7,130) 130 FORMAT(/38X,'Standard'/' Period',11X,'N',8X,'Variance',6X,'error' ,,7X,'95% Confidence interval'/1X,75('-')) DO 152 I=2,S SE=XSQRT(COV(I,I)) C1=XN(I)-.196D1*SE C2=XN(I)+.196D1*SE WRITE(7,100)YR(I),XN(I),COV(I,I),SE,C1,C2 152 IF(MOD(YR(I),5).EQ.0)WRITE(7,100) AVG=AVG/DBLE(FLOAT(S-1)) AVGSE=AVGSE/DBLE(FLOAT(S-1)) SE=XSQRT(AVGSE) C1=AVG-.196D1*SE C2=AVG+.196D1*SE WRITE(7,110)AVG,AVGSE,SE,C1,C2 C---------------------------------------------------------- C COMPUTE VAR & COV OF B WRITE(7,160) 160 FORMAT(/38X,'Standard'/' Period',11X,'B',8X,'Variance',6X,'error' ,,7X,'95% Confidence interval Covariance(B(i-1),B(i))'/ /1X,100('-')) AVG=.0D0 AVGSE=.0D0 DO 180 I=2,S-1 AVG=AVG+XB(I) TMP=(U(I+1)-FI*U(I))/P VARB=Q*U(I+1)/P+Q*FI*U(I)*(.1D1+Q*FI/P)+ + TMP*TMP*V(2,2)+Q*U(I)*Q*U(I)*V(1,1)+ + .2D1*TMP*Q*U(I)*V(1,2)+SY(I)*SY(J)*V(3,3)+ + SY(I)*((FI*U(J)-U(J+1))/P*V(2,3)-U(J)*Q*V(1,3))+ + SY(J)*((FI*U(I)-U(I+1))/P*V(2,3)-U(I)*Q*V(1,3)) AVGSE=AVGSE+VARB COVB=.1D30 IF(I.LE.2)GO TO 170 TMP1=(U(I)-FI*U(I-1))/P COVB=-Q*FI*U(I)/P+TMP1*TMP*V(2,2)+Q*Q*U(I-1)*U(I)*V(1,1)+ + (Q*U(I-1)*TMP+Q*U(I)*TMP1)*V(1,2) AVGSE=AVGSE+.2D1*COVB SE=XSQRT(VARB) C1=XB(I)-.196D1*SE C2=XB(I)+.196D1*SE 170 WRITE(7,100)YR(I),XB(I),VARB,SE,C1,C2,COVB 180 IF(MOD(YR(I),5).EQ.0)WRITE(7,100) AVG=AVG/DBLE(FLOAT(S-2)) AVGSE=AVGSE/DBLE(FLOAT(S-2)) SE=XSQRT(AVGSE) C1=AVG-.196D1*SE C2=AVG+.196D1*SE WRITE(7,110)AVG,AVGSE,SE,C1,C2 RETURN END