SUBROUTINE B2MNB(TITLE,FI,FY,P,XA,XY,V,VCP,COVFIP) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER IDIM,IDIM2 PARAMETER (IDIM=50,IDIM2=IDIM*IDIM+2*IDIM) INTEGER I,J,K,L,S,MAXITR,YR(IDIM),IPRNT CHARACTER*128 TITLE DIMENSION P(IDIM),W(IDIM+2,IDIM),VCP(IDIM,IDIM),COVFIP(IDIM), , XA(IDIM),XY(IDIM),X(IDIM+2,IDIM),XM(IDIM),COVPM(IDIM,IDIM), , COV(IDIM,IDIM),V(IDIM+2,IDIM+2),COVFYP(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 EQUIVALENCE (COVPM,COVFYP) DATA AVGSE,XM,X/.0D0,IDIM*.0D0,IDIM2*.0D0/ C C THIS SUBROUTINE COMPUTES ESTIMATES OF M,N, AND B FOR Model B2. C C STATEMENT FUNCTIONS U(I)=(N(I)-M(I))/P(I) XN(I)=XM(I)+U(I) XB(I)=U(I+1)-FI*(.1D1-P(I))*U(I)+SY(I)*FY Q(I)=.1D1-P(I) C XM(S)=(M(S)+B(S))/P(S) AVG=XM(S) X(1,1)=-Q(2)*XA(2)*SY(1)/XY(1)*FY/FI ALFAJ=.1D1 DO 20 J=2,S-1 XM(J)=(M(J)+B(J))/(.1D1-Q(J)*XA(J)) AVG=AVG+XM(J) X(J,J-1)=-Q(J)*XA(J)/(.1D1-XA(J-1))*((XM(J-1)*Q(J-1)+NS(J-1))* * FI/XA(J-1)+SY(J-1)*FY/XY(J-1)) ALFAJ=ALFAJ*FI*Q(J+1) ALFA=ALFAJ DO 10 I=1,J X(1,J)=X(1,J)+(XM(I)/XA(I)-FY/FI*SY(I)/XY(I))*ALFA X(S+1,J)=X(S+1,J)+SY(I)*ALFA/XY(I) 10 ALFA=ALFA/FI/Q(J-I+1) X(1,J)=X(1,J)*Q(J+1)*XA(J+1) X(S+1,J)=X(S+1,J)*Q(J+1)*XA(J+1) ALFA=.1D1 DO 20 I=J,S-1 X(J,I)=Q(I+1)*XA(I+1)*ALFA/XA(J)*(XM(J-1)*Q(J-1)+NS(J-1)* * FI*(.1D1-FI)/XA(J-1)+SY(J-1)*FY*(.1D1-FY)/XY(J-1)) 20 ALFA=ALFA*FI*Q(I+1) C CALL MMULT(V,X,W,S+1,S+1,S-1,IDIM+2,IDIM+2,IDIM) C------------------------------------------------------ C COMPUTE VAR & COV OF M, AND OUTPUT IF(PRNT.GT.1)WRITE(7,30)TITLE 30 FORMAT(/1X,A128//' Model B2 Var-cov matrix of M (i=2,s;j=i,s)') DO 80 I=2,S QXAI=.1D1-Q(I)*XA(I) ALFA=.1D1 DO 60 J=I,S QXAJ=.1D1-Q(J)*XA(J) COV(I,J)=Q(J)*XA(J)*XM(I)*ALFA/QXAI+XM(I)*W(1,J-1)/FI+ + XM(I)*W(I,J-1)/(.1D1-XA(I-1))+XM(J)/FI*(XM(I)*V(1,1)/FI+ + W(1,I-1)+XM(I)*V(1,I)/(.1D1-XA(I-1)))+XM(J)/(.1D1-XA(J-1))* * (W(J,I-1)+XM(I)*V(1,J)/FI+XM(I)/(.1D1-XA(I-1))*V(I,J)) IF(J.LT.S)ALFA=ALFA*FI*Q(J+1) 60 AVGSE=AVGSE+.2D1*COV(I,J) IF(PRNT.GT.1)WRITE(7,70)I,(COV(I,J),J=I,S) 70 FORMAT(' i=',I4,10(1X,F11.2)/7X,10(1X,F11.2)) 80 AVGSE=AVGSE-COV(I,I) 90 FORMAT(I5,4X,F11.2,3X,F12.2,1X,F9.2,6X,F11.2,' - ',F11.2) WRITE(7,40) 40 FORMAT(/38X,'Standard'/' Period',11X,'M',8X,'Variance',6X,'error' ,,7X,'95% confidence interval'/1X,75('-')) DO 92 I=2,S SE=XSQRT(COV(I,I)) C1=XM(I)-.196D1*SE C2=XM(I)+.196D1*SE WRITE(7,90)YR(I),XM(I),COV(I,I),SE,C1,C2 92 IF(MOD(YR(I),5).EQ.0)WRITE(7,90) 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,100)AVG,AVGSE,SE,C1,C2 100 FORMAT(1X,76('-')/' MEAN',4X,F11.2,3X,F12.2,1X,F9.2,6X,F11.2, ,' - ',F11.2/1X,76('-')) C----------------------------------------------------------- C COMPUTE COVAR(P(I),M(J)) IF(PRNT.GT.1)WRITE(7,110) 110 FORMAT(/' Model B2 Covariance( P(i),M(j) ) (i=2,s;j=2,s)') DO 130 I=2,S-1 QXAI=.1D1-Q(I)*XA(I) DO 120 J=2,S 120 COVPM(I,J)=(Q(I)*W(I+1,J-1)-W(I,J-1)/FI-QXAI/FI*W(1,J-1)+ + XM(J)/FI*(Q(I)*V(1,I+1)-V(1,I)/FI-QXAI/FI*V(1,1))+ + XM(J)/(.1D1-XA(J-1))*(Q(I)*V(I+1,J)-V(I,J)/FI-QXAI/FI* * V(1,J)))/XA(I) 130 IF(PRNT.GT.1)WRITE(7,140)I,(COVPM(I,J),J=2,S) 140 FORMAT(I5,(10(1X,F12.6))) DO 150 J=2,S 150 COVPM(S,J)=-W(S,J-1)/FI-P(S)/FI*W(1,J-1)-XM(J)/FI* * (P(S)/FI*V(1,1)+V(1,S)/FI)-XM(J)/(.1D1-XA(J-1))* * (P(S)/FI*V(1,J)+V(S,J)/FI) IF(PRNT.GT.1)WRITE(7,140)S,(COVPM(S,J),J=2,S) C--------------------------------------------------------- C COMPUTE VAR & COV OF N AND OUTPUT IF(PRNT.GT.1)WRITE(7,160) 160 FORMAT(/' Model B2 N Variance-Covariance matrix (i=2,s;j=i,s)') AVG=.0D0 AVGSE=.0D0 DO 190 I=2,S DO 180 J=I,S COV(I,J)=COV(I,J)-U(I)/P(I)*COVPM(I,J)-U(J)/P(J)*COVPM(J,I)+ + U(I)/P(I)*U(J)/P(J)*VCP(I,J) 180 AVGSE=AVGSE+.2D1*COV(I,J) IF(I.LT.S)COVFYP(I)=(Q(I)*V(I+1,S+1)-V(I,S+1)/FI- - (.1D1-Q(I)*XA(I))/FI*V(1,S+1))/XA(I) AVGSE=AVGSE-.2D1*COV(I,I) COV(I,I)=COV(I,I)+U(I)/P(I)*Q(I) AVGSE=AVGSE+COV(I,I) AVG=AVG+XN(I) 190 IF(PRNT.GT.1)WRITE(7,70)I,(COV(I,J),J=I,S) WRITE(7,170) 170 FORMAT(/38X,'Standard'/' Period',11X,'N',8X,'Variance',6X,'error' ,,7X,'95% Confidence interval'/1X,75('-')) DO 192 I=2,S SE=XSQRT(COV(I,I)) C1=XN(I)-.196D1*SE C2=XN(I)+.196D1*SE WRITE(7,90)YR(I),XN(I),COV(I,I),SE,C1,C2 192 IF(MOD(YR(I),5).EQ.0)WRITE(7,90) COVFYP(S)=-(V(S,S+1)+P(S)*V(1,S+1))/FI 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,100)AVG,AVGSE,SE,C1,C2 C----------------------------------------------------------- C COMPUTE VAR & COV OF B AND OUTPUT IF(PRNT.GT.1)WRITE(7,200) 200 FORMAT(/' Model B2 B VAR-COVAR MATRIX (I=2,S-1;J=I,S-1)') AVG=.0D0 AVGSE=.0D0 DO 230 I=2,S-1 AVG=AVG+XB(I) VCP(I+1,I)=VCP(I,I+1) DO 220 J=I,S-1 C** , VCP(I,J),VCP(I,J+1),VCP(I+1,J),VCP(I+1,J+1) COV(I,J)=Q(I)*U(I)*(Q(J)*U(J)*V(1,1)-FI*U(J)/P(J)*COVFIP(J)+ + U(J+1)/P(J+1)*COVFIP(J+1))+FI*U(I)/P(I)*(-Q(J)*U(J)* * COVFIP(I)+FI*U(J)/P(J)*VCP(I,J)-U(J+1)/P(J+1)*VCP(I,J+1))+ + U(I+1)/P(I+1)*(Q(J)*U(J)*COVFIP(I+1)-FI*U(J)/P(J)*VCP(I+1,J)+ + U(J+1)/P(J+1)*VCP(I+1,J+1))+SY(I)*SY(J)*V(S+1,S+1)+ + SY(I)*(FI*U(J)/P(J)*COVFYP(J)-U(J)*Q(J)*V(1,S+1)- - U(J+1)/P(J+1)*COVFYP(J+1))+SY(J)*(FI*U(I)/P(I)*COVFYP(I)- - U(I)*Q(I)*V(1,S+1)-U(I+1)/P(I+1)*COVFYP(I+1)) IF(J.EQ.I+1)COV(I,J)=COV(I,J)-Q(J)*FI*U(J)/P(J) IF(J.EQ.I)COV(I,I)=COV(I,I)+U(I+1)*Q(I+1)/P(I+1)+ + Q(I)*FI*U(I)*(.1D1+Q(I)*FI/P(I)) 220 AVGSE=AVGSE+.2D1*COV(I,J) 230 IF(PRNT.GT.1)WRITE(7,70)I,(COV(I,J),J=I,S-1) WRITE(7,210) 210 FORMAT(/38X,'Standard'/' Period',11X,'B',8X,'Variance',6X,'error' ,,7X,'95% Confidence interval'/1X,75('-')) DO 232 I=2,S-1 SE=XSQRT(COV(I,I)) C1=XB(I)-.196D1*SE C2=XB(I)+.196D1*SE WRITE(7,90)YR(I),XB(I),COV(I,I),SE,C1,C2 232 IF(MOD(YR(I),5).EQ.0)WRITE(7,90) 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,100)AVG,AVGSE,SE,C1,C2 RETURN END