APPENDIX C  THE MULTIPLE DYNAMIC WEIGHTED LEAST SQUARES PROGRAM


       PROGRAM MDWLS
C       MULTIPLE WEIGHTED D.L.S.
       COMMON DATA(7,20),X(19),Y(19),IT,B1(5,1),SMSQER,ID,NEXOG
       SMSQER=0.0
       ID=0
       NEXOG=3       	       ! INCLUDING THE CONSTANT
       DO 1 I=1,19
1       READ(11,*)(DATA(J,I),J=1,3)
       DO 2 I=1,19
       Y(I)=DATA(3,I)
       DATA(3,I)=ALOG(DATA(2,I))
       DATA(2,I)=ALOG(DATA(1,I))
       DATA(1,I)=1.0
2       CONTINUE
       WRITE(55,1000)
1000       FORMAT('    T        B1        B2        B3         Y        FY',
       1'     SQERR   SMSQERR      RMSE')
       DO 3 IT=NEXOG,19
       CALL REGRES
       CALL FORECAST
3       CONTINUE
       STOP
       END


       SUBROUTINE REGRES
       COMMON DATA(7,20),X(19),Y(19),IT,B1(5,1),SMSQER,ID,NEXOG
       DIMENSION WT(19,19),XX(19,5),XT(5,19),T(5,19),XWX(5,5),YY(19,1)
       1,S(5,19),B2(5,1),DXWX(5,5)
       DOUBLE PRECISION DXWX,DET
       DO 2 I=1,IT
       DO 3 J=1,IT
3       WT(I,J)=0.0
       WT(I,I)=1.00**(IT-I)
       YY(I,1)=Y(I)
       DO 2 J=1,NEXOG
       XX(I,J)=DATA(J,I)
2       XT(J,I)=DATA(J,I)
       CALL MXMULT(XT,NEXOG,IT,WT,IT,T,5,19,19)        !X'W
       CALL MXMULT(T,NEXOG,IT,XX,NEXOG,XWX,5,19,5)         !X'WX
       DO 4 I=1,NEXOG
       DO 4 J=1,NEXOG
4       DXWX(I,J)=XWX(I,J)
       CALL MATINV(DXWX,NEXOG,DET,5)
       IF(DET.EQ.0.0) STOP 'SINGULARITY'
       DO 5 I=1,NEXOG
       DO 5 J=1,NEXOG
5       XWX(I,J)=DXWX(I,J)
       CALL MXMULT(T,NEXOG,IT,YY,1,B2,5,19,1)           !X'WY
       CALL MXMULT(XWX,NEXOG,NEXOG,B2,1,B1,5,5,1)  !X'WX**-1 X'WY
       RETURN
       END
       SUBROUTINE MXMULT(A1,I1,J1,A2,K1,A3,I2,J2,K2)
       DIMENSION A1(I2,J2),A2(J2,K2),A3(I2,K2)
       DO 1 I=1,I1
       DO 1 K=1,K1
       A3(I,K)=0.0
       DO 1 J=1,J1
1       A3(I,K)=A3(I,K)+A1(I,J)*A2(J,K)
       RETURN
       END


      SUBROUTINE MATINV(A,M,DET,SIZE)   
C       
C     GAUSS REDUCTION INVERSION WITHOUT ROW AND COLUMN INTERCHANGES.    
C     M IS THE ORDER OF THE SQUARE MATRIX , A.  
C     A-INVERSE IS RETURNED IN A.       
C     DETERMINANT IS RETURNED IN DET.   
C       
      INTEGER SIZE      
      DOUBLE PRECISION A,DET,PVT,W      
      DIMENSION A(SIZE,SIZE)    
      DET=1.0   
      DO 1 J=1,M
      PVT=A(J,J)
      DET=DET*PVT       
      A(J,J)=1.0
      DO 2 K=1,M
C     DIVIDE THE PIVOT ROW BY PIVOT ELEMENT.    
    2 A(J,K)=A(J,K)/PVT 
      DO 1 K=1,M
C     REDUCE THE NON-PIVOT ROWS.
      IF(K-J) 3,1,3     
    3 W=A(K,J)  
      A(K,J)=0.0
      DO 4 L=1,M
    4 A(K,L)=A(K,L)-A(J,L)*W    
    1 CONTINUE  
      RETURN    
      END       


       SUBROUTINE FORECAST
       COMMON DATA(7,20),X(19),Y(19),IT,B1(5,1),SMSQER,ID,NEXOG
       IF(IT.EQ.19)WRITE(55,1000)J,(B1(I,1),I=1,NEXOG)
       IF(IT.EQ.19) STOP
       J=IT+1
       ID=ID+1
       FY=0.0
       DO 1 I=1,NEXOG
1       FY=FY+DATA(I,J)*B1(I,1)
       ERR=FY-Y(J)
       SQERR=ERR*ERR
       SMSQER=SMSQER+SQERR
       RMSE=SQRT(SMSQER/ID)
       WRITE(55,1000)J,(B1(I,1),I=1,NEXOG),Y(J),FY,SQERR,SMSQER,RMSE
1000       FORMAT(I5,10F10.5)
       RETURN
       END



'Two added to one - if that could but be done,'
It said, 'with one's fingers and thumbs!'
The Hunting of the Snark, Lewis Carroll

APPENDIX D  THE KALMAN FILTER PROGRAM FOR CHAPTER 9


       PROGRAM BAYESCH9
C  FORECASTING PROGRAM USING RECURSIVE BAYESIAN ESTIMATION      
C  SEE 'BAYESIAN FORECASTING', HARRISON AND STEVENS,3000
C  J.R.STATIST.SOC,B,38,PP205-247       
C  NOTATION USED IS THEIRS. LIMITS (FAIRLY EASY TO CHANGE)      
C  ARE M.LE.20 N.LE.24  
C  ITIME IS TIME, NTIME IS THE NUMBER OF DATA POINTS    
C  M IS THE NUMBER OF DEPENDENT VARIABLES       
C  N IS THE NUMBER OF PARAMETERS
C  DATA WILL NOT BE STORED, BUT READ IN AS NEEDED       
       IMPLICIT REAL*8 (A-H,O-Z)
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
C  CURRENT TIME, TIME MAX, M, N, NO OF FORECASTS, OUTPUT OPTIONS
C  THETAM IS THEIR M(T), THETAC IS THEIR C(T)   
C  Y IS THE OBSERVED VECTOR, F IS THE MATRIX OF INDEPENDENT VARIABLES   
C  G IS THE MATRIX DESCRIBING THE EVOLUTION OF THE PARAMETERS   
C  G=I FOR CONSTANT PARAMETERS.  YHAT IS THE FORECAST OF Y, BIGYHAT     
C  IS THE VARIANCE-COVARIANCE MATRIX OF YHAT.   
C  V IS THE ERROR IN OBSERVING Y, MEAN ZERO     
C  W IS THE ERROR IN UPDATING THE PARAMETERS, MEAN ZERO 
C  BIGV IS THE VAR-COV MATRIX OF V; BIGW IS THE SAME FOR W      
C  IF BIGW=0, THE CERTAINTY WITH WHICH YOU KNOW PARAMETERS NEVER
C  REDUCES. IF BIGW.NE.0 THEN THE CERTAINTY REDUCES AS TIME PASSES.     
C  THE DYNAMIC LINEAR MODEL IS: 


C       
C  Y=F.THETA+v  
C  THETA=G.THETA(-1)+w
      MMAX=20   
      NMAX=24   
      ITIME=0   
      CALL   GETDATA    
C      DOES INITIALIZATION FOR A GIVEN MODEL    
C       GETDATA IS  PROBLEM SPECIFIC    
C      PICKS UP NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS AND G    
C      INITIALIZES THETA,THETAM AND THETAC      
      CALL INITIALIZE   
100   ITIME=ITIME+1     
      CALL GETDATA      
C      PICKS UP CURRENT Y,F,BIGV,BIGW   
      CALL KALMAN       
C      UPDATES THETAM AND THETAC BY THE KALMAN FILTER   
       CALL UPERR
      IF(IFMATS.EQ.1)CALL OUTMATS       
      IF(ITIME.LT.NTIME) GOTO 100       
      IF(IFMATS.EQ.0) CALL OUTMATS      
101   CONTINUE  
      IF(NFORC.EQ.0) STOP 'BAYES PROGRAM TERMINATES NORMALLY'   
      NFORC=NFORC-1     
      ITIME=ITIME+1     
      CALL GETDATA      
C      SAME AS GETDATA, EXCEPT THERE ISNT ANY Y, AND THERE      
C      CAN BE UNCERTAINTY IN F  
      CALL FORECAST     
      IF(IFFMATS.EQ.1) CALL OUTMATS     
      IF(NFORC.EQ.0.AND.IFFMATS.EQ.0) CALL OUTMATS      
      GOTO 101  
200   STOP 'BAYES PROGRAM TERMINATES IMPOSSIBLY'
      END       


      SUBROUTINE INITIALIZE     
       IMPLICIT REAL*8 (A-H,O-Z)
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
       IF(MVARS.GT.MMAX.OR.NPARS.GT.NMAX)STOP'TOO MANY VARIABLES
     1OR PARAMETERS'    
       DO 2 I=1,NPARS
       DO 2 J=1,NPARS
       IF(THETAC(I,J).NE.THETAC(J,I))STOP ' ASYMMETRIC THETAC 0'
2       CONTINUE
C   READ THE INITIAL VALUES OF THETA IN GETDATA
C   AND THETAC
1002  FORMAT(16F5.0)    
C  DONT NEED MUCH ACCURACY; ZEROS OR ONES WILL OFTEN BE APPROPRIATE     
C  AND NOW THE VAR-COV MATRIX. I SHALL ASSUME THAT THETAC IS    
C  INITIALLY DIAGONAL (I.E. PARAMETERS ARE INDEPENDENT). IT IS  
       DO 4 I=1,MVARS    
      SUMSQ(I)=0.0      
      SUMNAI(I)=0.0     
4     YLAST(I,1)=0.0    
      WRITE(7,1004)     
      DO 5 I=1,NPARS    
5     WRITE(7,1003)(G(I,J),J=1,NPARS)   
1003  FORMAT(1X,8E10.3) 
1004  FORMAT(' G-MATRIX')       
1005  FORMAT(' THETAM AND THETAC STARTING VALUES')      
      WRITE(7,1005)     
      WRITE(7,1003)(THETAM(I,1),I=1,NPARS)      
      DO 7 I=1,NPARS    
7     WRITE(7,1003)(THETAC(I,J),J=1,NPARS)      
      RETURN    
      END       


      SUBROUTINE KALMAN 
       IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION T7,DDD   
      DIMENSION T1(24,24),T2(20,24),T3(24,20),T4(20,20),
     1T5(24,24),T6(24,20),T7(20,20),T8(24,1),T9(24,1),  
     2T10(20,24)
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
      M=MMAX    
      N=NMAX    
C YHAT AND E    
      CALL MXMULT(F,M,N,G,N,T2) 
      CALL MXMULT(T2,M,N,THETAM,1,YHAT) 
      CALL MXADD(Y,M,1,YHAT,E,-1.0)     
C      R
      CALL MXTRANS(G,N,N,T1)    
      CALL MXMULT(THETAC,N,N,T1,N,T5)   
      CALL MXMULT(G,N,N,T5,N,T1)
      CALL MXADD(T1,N,N,BIGW,R,1.0)     
C      BIGYHAT  
      CALL MXTRANS(F,M,N,T3)    
      CALL MXMULT(R,N,N,T3,M,T6)
      CALL MXMULT(F,M,N,T6,M,T4)
      CALL MXADD(T4,M,M,BIGV,BIGYHAT,1.0)       
C  INVERT BIGYHAT. AS A PRECAUTION, DO THIS IN DOUBLE PRECISION 
      DO 1 I=1,MVARS    
      DO 1 J=1,MVARS    
1     T7(I,J)=BIGYHAT(I,J)      
      CALL MATINV(T7,MVARS,DDD,M)       
      DET=DDD   
      DO 2 I=1,MVARS    
      DO 2 J=1,MVARS    


2     YINV(I,J)=T7(I,J) 
C       A       
      CALL MXTRANS(F,M,N,T3)    
      CALL MXMULT(T3,N,M,YINV,M,T6)     
      CALL MXMULT(R,N,N,T6,M,A) 
C       THETAM  
      CALL MXMULT(G,N,N,THETAM,1,T8)    
      CALL MXMULT(A,N,M,E,1,T9) 
      CALL MXADD(T8,N,1,T9,THETAM,1.0)  
C       THETAC  
      CALL MXTRANS(A,N,M,T2)    
      CALL MXMULT(BIGYHAT,M,M,T2,N,T10) 
      CALL MXMULT(A,N,M,T10,N,T5)       
      CALL MXADD(R,N,N,T5,THETAC,-1.0)  
CTHATS ALL FOLKS
      RETURN    
      END       


      SUBROUTINE MATINV(A,M,DET,SIZE)   
       IMPLICIT REAL*8 (A-H,O-Z)
C       
C     GAUSS REDUCTION INVERSION WITHOUT ROW AND COLUMN INTERCHANGES.    
C     M IS THE ORDER OF THE SQUARE MATRIX , A.  
C     A-INVERSE IS RETURNED IN A.       
C     DETERMINANT IS RETURNED IN DET.   
C       
      INTEGER SIZE      
      DOUBLE PRECISION A,DET,PVT,W      
      DIMENSION A(SIZE,SIZE)    
      DET=1.0   
      DO 1 J=1,M
      PVT=A(J,J)
      DET=DET*PVT       
      A(J,J)=1.0
      DO 2 K=1,M
C     DIVIDE THE PIVOT ROW BY PIVOT ELEMENT.    
    2 A(J,K)=A(J,K)/PVT 
      DO 1 K=1,M
C     REDUCE THE NON-PIVOT ROWS.
      IF(K-J) 3,1,3     
    3 W=A(K,J)  
      A(K,J)=0.0
      DO 4 L=1,M
    4 A(K,L)=A(K,L)-A(J,L)*W    
    1 CONTINUE  
      RETURN    
      END       


      SUBROUTINE MXADD(ARR1,II,JJ,ARR2,ARROUT,SIGN)     
       IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ARR1(II,JJ),ARR2(II,JJ),ARROUT(II,JJ)   
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
      IL=NPARS  
      JL=NPARS  
      IF(II.EQ.MMAX)IL=MVARS    
      IF(JJ.EQ.MMAX)JL=MVARS    
      IF(II.EQ.1)IL=1   
       IF(JJ.EQ.1)JL=1  
C  IF SIGN = 1.0, ADDS. IF -1.0, SUBTRACTS      
      DO 1 I=1,IL       
      DO 1 J=1,JL       
1     ARROUT(I,J)=ARR1(I,J)+SIGN*ARR2(I,J)      
      RETURN    
      END       
      SUBROUTINE MXTRANS(ARR,II,JJ,ARROUT)      
       IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ARR(II,JJ),ARROUT(JJ,II)
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
      IL=NPARS  
      JL=NPARS  
      IF(II.EQ.MMAX)IL=MVARS    
      IF(JJ.EQ.MMAX)JL=MVARS    
      IF(II.EQ.1)IL=1   
       IF(JJ.EQ.1)JL=1  
      DO 1 I=1,IL       
      DO 1 J=1,JL       
1     ARROUT(J,I)=ARR(I,J)      
      RETURN    
      END       


      SUBROUTINE MXMULT(ARR1,II,JJ,ARR2,KK,ARROUT)      
       IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ARR1(II,JJ),ARR2(JJ,KK),ARROUT(II,KK)   
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
      IL=NPARS  
      JL=NPARS  
      IF(II.EQ.MMAX)IL=MVARS    
      IF(JJ.EQ.MMAX)JL=MVARS    
      IF(II.EQ.1)IL=1   
       IF(JJ.EQ.1)JL=1  
      KL=NPARS  
      IF(KK.EQ.MMAX)KL=MVARS    
      IF(KK.EQ.1)KL=1   
      DO 1 I=1,IL       
      DO 1 K=1,KL       
      ARROUT(I,K)=0.0   
      DO 1 J=1,JL       
1     ARROUT(I,K)=ARROUT(I,K)+ARR1(I,J)*ARR2(J,K)       
      RETURN    
      END       
      SUBROUTINE TOVAR(R)       
       IMPLICIT REAL*8 (A-H,O-Z)
      IF(R/100. + 1.0.LE.0.0) GOTO 2    
      R=DLOG(R/100.0+1.0)       
      R=R*R     
C  TRANSFORMS AVERAGE % ERROR TO VARIANCE       
      RETURN    
2     R=99999.0 
      RETURN    
      END       


      SUBROUTINE FROMVAR(R)     
       IMPLICIT REAL*8 (A-H,O-Z)
      IF(R.LT.0.0) GOTO 2       
      IF(R.GT.100.) GOTO 1
      R=100.0*(EXP(SQRT(R))-1.0)
C  TRANSFORMS VARIANCE TO AVERAGE % ERROR       
      RETURN    
1     R=99999.  
      RETURN    
2     CONTINUE  
      RETURN    
      END       
      SUBROUTINE CHECKDATA      
       IMPLICIT REAL*8 (A-H,O-Z)
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
C   CHECKS THAT VAR-COV MATRIX IS SYMMETRIC     
      DO 1 I=1,MVARS    
      DO 1 J=1,MVARS    
      IF(BIGV(I,J).NE. BIGV(J,I)) STOP 'ASSYMETRIC BIGV'
1     CONTINUE  
      DO 2 I=1,NPARS    
      DO 2 J=1,NPARS    
      IF(BIGW(I,J).NE.BIGW(J,I)) STOP 'ASSYMETRIC BIGW' 
2     CONTINUE  
      RETURN    
      END       


      SUBROUTINE OUTMATS
       IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION T(24),S(20)     
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
C  OUTPUTS THETAM, THETAC, Y, YHAT, E, R, BIGYHAT, YINV, A      
      WRITE(7,999) ITIME
999   FORMAT('1TIME PERIOD NUMBER',40X,I5)      
      WRITE(7,2000)     
2000  FORMAT(' INDEPENDENT VARIABLES (F)')      
       DO 8 I=1,MVARS
8     WRITE(7,1001)(F(I,J),J=1,NPARS)   
      WRITE(7,2001)     
2001  FORMAT(/' VAR-COV OF DEPENDENT, INDEPENDENT (BIGV, BIGW)')
      DO 11 I=1,MVARS   
11    WRITE(7,1001)(BIGV(I,J),J=1,MVARS)
      DO 12 I=1,NPARS   
12    WRITE(7,1001)(BIGW(I,J),J=1,NPARS)
      WRITE(7,1005)     
1005  FORMAT(/' DEPENDENT VARIABLE; LAST YEAR, ACTUAL, FORECAST,'/      
     1' ERROR, %ERROR, %ERROR NAIVE MODEL, SUMSQ ERRORS, SUMSQ NAIVE')  
      IF(ITIME.LE.NTIME)WRITE(7,1021)(YLAST(I,1),I=1,MVARS)     
      IF(ITIME.LE.NTIME) WRITE(7,1031)(Y(I,1),I=1,MVARS)
       WRITE(7,1041)(YHAT(I,1),I=1,MVARS)       
      IF(ITIME.LE.NTIME)WRITE(7,1051)(E(I,1),I=1,MVARS) 
      DO 3 I=1,MVARS    
      IF(Y(I,1).EQ.0.0)S(I)=99999.      
      IF(Y(I,1).NE.0.0)S(I)=100.0*E(I,1)/Y(I,1) 
      T(I)=BIGYHAT(I,I) 
3     T(I)=SQRT(ABS(T(I)))      
      IF(ITIME.GT.NTIME) GOTO 3000      
      WRITE(7,1063)(S(I),I=1,MVARS)     


      DO 21 I=1,MVARS   
      IF(Y(I,1).NE.0.0)S(I)=100.0*(Y(I,1)-YLAST(I,1))/Y(I,1)    
      IF(Y(I,1).EQ.0.0)S(I)=99999.      
21    CONTINUE  
      WRITE(7,1073)(S(I),I=1,MVARS)     
      WRITE(7,1081)(SUMSQ(I),I=1,MVARS) 
      WRITE(7,1091)(SUMNAI(I),I=1,MVARS)
1021  FORMAT(' PREVIOUS YEAR  ',8E10.3) 
1031  FORMAT(' THIS YEAR      ',8E10.3) 
1041  FORMAT(' FORECAST       ',8E10.3) 
1051  FORMAT(' ERROR          ',8E10.3) 
1063  FORMAT(' %ERROR         ',8F10.2) 
1073  FORMAT(' %ERROR NAIVE   ',8F10.2) 
1081  FORMAT(' SUM OF SQ.     ',8E10.3) 
1091  FORMAT(' SUMSQ NAIVE    ',8E10.3) 
3000  CONTINUE  
      WRITE(7,1006)     
1006  FORMAT(/' STANDARD DEVIATION ON FORECAST')
      WRITE(7,1003)(T(I),I=1,MVARS)     
       WRITE(7,1010)    
1010  FORMAT(/' FORECAST VARIANCE-COVARIANCE MATRIX FOR THETA (R)')     
      DO 6 I=1,NPARS    
6     WRITE(7,1001)(R(I,J),J=1,NPARS)   
      WRITE(7,1007)     
1007  FORMAT(/' VARIANCE-COVARIANCE MATRIX FOR FORECAST (BIGYHAT)')     
      DO 4 I=1,MVARS    
4     WRITE(7,1001)(BIGYHAT(I,J),J=1,MVARS)     
      IF(ITIME.GT.NTIME)RETURN  
      WRITE(7,1008)     
1008  FORMAT(/' INVERSE OF VARIANCE-COVARIANCE MATRIX (BIGYHAT-1)')     
      DO 5 I=1,MVARS    
5     WRITE(7,1001)(YINV(I,J),J=1,MVARS)
      WRITE(7,1009) DET 
1009  FORMAT(' DETERMINANT= ',E10.3)    
      WRITE(7,1011)     
1011  FORMAT(/' SMOOTHING FACTOR (A)')  


      DO 7 I=1,NPARS    
7     WRITE(7,1001)(A(I,J),J=1,MVARS)   
      WRITE(7,1000)     
1000  FORMAT(/' MAXIMUM LIKELIHOOD PARAMETER ESTIMATES (THETAM)')       
      WRITE(7,1001)(THETAM(I,1),I=1,NPARS)      
1001  FORMAT(8(1X,E10.3))       
      WRITE(7,1002)     
1002  FORMAT(/' STANDARD DEVIATION ON PARAMETER ESTIMATES')     
      DO 1 I=1,NPARS    
      T(I)=THETAC(I,I)  
1     T(I)=SQRT(ABS(T(I)))      
      WRITE(7,1003)(T(I),I=1,NPARS)     
1003  FORMAT(1X,8F10.2) 
      WRITE(7,1004)     
1004  FORMAT(/' VARIANCE-COVARIANCE MATRIX FOR THETA (THETAC)') 
      DO 2 I=1,NPARS    
2     WRITE(7,1001)(THETAC(I,J),J=1,NPARS)      
      RETURN    
      END       


      SUBROUTINE FORECAST       
       IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION T1(24,1),T2(24,24),T3(24,24),T4(24,20),T5(24,20),       
     1T6(20,20) 
C  THIS ASSUMES THAT F IS KNOWN PRECISELY       
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
C  FIRST UPDATE THETAM AND THETAC USING EQUATION 2.3.3 AND 2.3.4
      N=NMAX    
      M=MMAX    
      CALL MXMULT(G,N,N,THETAM,1,T1)    
      DO 1 I=1,N
1     THETAM(I,1)=T1(I,1)       
      CALL MXTRANS(G,N,N,T2)    
      CALL MXMULT(THETAC,N,N,T2,N,T3)   
      CALL MXMULT (G,N,N,T3,N,T2)       
      CALL MXADD(T2,N,N,BIGW,THETAC,1.0)
C  NOW THETAM AND THETAC ARE UPDATED    
C  NOW FORECAST USING EQUATIONS 2.3.5 AND 2.3.6 
      CALL MXMULT(F,M,N,THETAM,1,YHAT)  
      CALL MXTRANS(F,M,N,T4)    
      CALL MXMULT(THETAC,N,N,T4,M,T5)   
      CALL MXMULT(F,M,N,T5,M,T6)
      CALL MXADD(T6,M,M,BIGV,BIGYHAT,1.0)       
C  NOW YHAT AND BIGYHAT ARE UPDATED     
      RETURN    
      END       


      SUBROUTINE UPERR
       IMPLICIT REAL*8 (A-H,O-Z)
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
      IF (ITIME.LE.NPARS+1.OR. ITIME .GT. NTIME) GOTO 4000    
      DO 8 I=1,MVARS    
      SUMSQ(I)=SUMSQ(I)+E(I,1)*E(I,1)   
      Q=Y(I,1)-YLAST(I,1)       
8      SUMNAI(I)=SUMNAI(I)+Q*Q   
4000  CONTINUE  
      DO 1 I=1,MVARS    
1     YLAST(I,1)=Y(I,1) 
       IF(IFMATS.EQ.1)RETURN
       WRITE(7,1000)(THETAM(I,1),I=1,NPARS),(SUMSQ(I),I=1,MVARS)
1000       FORMAT(' PARAMETER ESTIMATES AND ERRORS',T35,9F10.3)
       RETURN
       END


      SUBROUTINE GETDATA
       IMPLICIT REAL*8 (A-H,O-Z)
       DOUBLE PRECISION DDATA(23)
       CHARACTER FILENAME*50
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
       COMMON /SPEC/ DW(24),DATA(6,22),YP(22),DV(20)
C  IF ITIME=0, THEN WE INITIALIZE AND EXIT      
C  IF ITIME.GT.1 AND .LE. NTIME, NORMAL KALMAN  
C  IF ITIME .GT. NTIME, FORECASTING TIME
C  THIS SUBROUTINE IS PROBLEM SPECIFIC  
C       
C  WE PICK UP AND/OR CALCULATE Y,F,BIGV,BIGW.   
C       THIS VERSION IS FOR CHAPTER NINE
      IF(ITIME.GT.0) GOTO 2000  
C  INITIALIZE   
       OPEN(7,FILE='BAYESCH9.LIS',STATUS='NEW')
       NTIME=22
       NFORC=0
       IFFMATS=0
      MVARS=1   
      NPARS=5   
       READ(5,2222)FILENAME
2222       FORMAT(A)
       OPEN(11,FILE=FILENAME,STATUS='OLD',FORM='UNFORMATTED')
       READ(11)FILENAME(1:1)
       DO 4 I=1,20
       READ(11)(DDATA(J),J=1,23)
C       1.0, PAGG, INCOME, DEGDAYS OR B/S, Q(T-1)
       DO 5 J=2,23
       IF(I.EQ.1)DATA(1,J-1)=1.0
       IF(I.EQ.5)DATA(2,J-1)=DDATA(J)


       IF(I.EQ.15)DATA(3,J-1)=DDATA(J)
       IF(I.EQ.16.AND.FILENAME(10:12).EQ.'TDF')DATA(4,J-1)=DDATA(J)/10000.
       IF(I.EQ.14)DATA(5,J-1)=DDATA(J-1)
       IF(I.EQ.15.AND.FILENAME(10:12).NE.'TDF')
       1DATA(4,J-1)=DDATA(J)-DDATA(J-1)-.022
       IF(I.EQ.14)YP(J-1)=DDATA(J)
5       CONTINUE
4       CONTINUE
       CLOSE(11)
      DO 12 I=1,NPARS    
      DO 12 J=1,NPARS    
      G(I,J)=0.0
      IF(I.EQ.J)G(I,J)=1.0      
12     CONTINUE  
C  G WILL OFTEN BE THE IDENTITY MATRIX. WHERE IT ISNT, CHANGE IT
       READ(5,*)IFMATS
       READ(5,*)(THETAM(I,1),I=1,NPARS)
       DO 20 I=1,NPARS
20       READ(5,*)(THETAC(I,J),J=1,NPARS)
       READ(5,*) (DW(I),I=1,NPARS)
       READ(5,*) (DV(I),I=1,MVARS)
      RETURN    
2000  CONTINUE  
      IF(ITIME.GT.NTIME) GOTO 3000      
C  GET DATA     
      Y(1,1)=YP(ITIME)
       DO 55 J=1,NPARS
55       F(1,J)=DATA(J,ITIME)
      DO 1 I=1,NPARS
      DO 1 J=1,NPARS
      BIGW(I,J)=0.0     
       BIGW(I,I)=DW(I)
1     CONTINUE  
C     .01 MEANS SIGMA=0.1 = LIKELY SIZE OF CHANGE IN PARAMETER  
C   BIGW CONTROLS THE DEGRADATION IN PARAMETER ESTIMATES
      BIGV(1,1)=DV(1)


      RETURN    
3000  CONTINUE  
C  FORECAST TIME
      Y(1,1)=1.0
      BIGV(1,1)=0.0     
      F(1,1)=1.0
      F(1,2)=F(1,2)+.1  
      F(1,3)=F(1,3)+0.00
      IF(ITIME.GT.NTIME+1) GOTO 4000    
      F(1,2)=DLOG(PW)   
      F(1,3)=DLOG(PS)   
4000  CONTINUE  
      DO 11 I=1,3       
      DO 11 J=1,3       
      BIGW(I,J)=0.0     
      IF(I.EQ.J)BIGW(I,I)=0.0   
11    CONTINUE  
C  ASSUMES NO DEGRADATION OF PARAMETER ESTIMATES AND    
C  AND NO UNCERTAINTY ON F      
      RETURN    
      END       


APPENDIX D (CONTINUED) THE GETDATA SUBROUTINE FOR CHAPTER 10


      SUBROUTINE GETDATA
C       THIS IS THE GETDATA FOR CHAPTER 10
       IMPLICIT REAL*8 (A-H,O-Z)
       DOUBLE PRECISION DDATA(23)
       CHARACTER FILENAME*50,DATAF*50,ERRORS*50
      COMMON Y(20,1),F(20,24),THETA(24,1),G(24,24),BIGV(20,20), 
     1BIGW(24,24),THETAM(24,1),THETAC(24,24),YHAT(20,1),E(20,1),
     2R(24,24),BIGYHAT(20,20),YINV(20,20),A(24,20),DET  
     3,YLAST(20,1),SUMSQ(20),SUMNAI(20) 
      COMMON /BITS/ ITIME,NTIME,MVARS,NPARS,NFORC,IFMATS,IFFMATS
     1,MMAX,NMAX
       COMMON /SPEC/ DW(24),DATA(6,22),YP(22),DV(20)
C  IF ITIME=0, THEN WE INITIALIZE M,C,G,BITS, OPEN FILES AND EXIT      
C  IF ITIME.GT.1 AND .LE. NTIME, NORMAL KALMAN  
C  IF ITIME .GT. NTIME, FORECASTING TIME
C  THIS SUBROUTINE IS PROBLEM SPECIFIC  
C       
C  WE PICK UP AND/OR CALCULATE Y,F,BIGV,BIGW.   
C  THIS VERSION IS FOR NDC AND ASSUMES THE DATAFILE IS ON FTN11 
C  MODEL IS THAT ESTIMATED IN 'MULTICOLLINEARITY AND A WOOL     
C  CONSUMPTION MODEL'   
C  (NDC/POP)/(PCE/POP)**C=A*PW**D*PS**E 
      IF(ITIME.GT.0) GOTO 2000  
C  INITIALIZE   
       OPEN(7,FILE='BAYESCH10.LIS',STATUS='NEW')
       NTIME=19
       NFORC=0
       IFMATS=0
       IFFMATS=1
      MVARS=2   
      NPARS=5   
       DO 1 I=1,NPARS
       DO 1 J=1,NPARS


       IF(I.EQ.J)G(I,J)=1.0
       IF(I.NE.J)G(I,J)=0.0
1       CONTINUE
1000  FORMAT(4I5,F10.0) 
       READ(5,4000)DATAF
       READ(5,4000)ERRORS
4000       FORMAT(A)
       READ(5,*)(DW(I),I=1,NPARS)
       READ(5,*)(THETAM(I,1),I=1,NPARS)
       READ(5,*)(THETAC(I,I),I=1,NPARS)
       WRITE(7,*)' THE DATA COMES FROM ',DATAF(1:20)
       1,' AND THE ERRORS FROM',ERRORS(1:20)
       OPEN(11,FILE=DATAF,STATUS='OLD')
       OPEN(12,FILE=ERRORS,STATUS='OLD')
      RETURN    
2000  CONTINUE  
      IF(ITIME.GT.NTIME) GOTO 3000      
C  GET DATA     
      IF(ITIME.EQ.1)READ(11,1001)PW,PS  
      READ(11,1100) POP,PCE,PRICE,PWN,PSN,WOOL,TOTFIB   
       PCE=PCE/1000.0
      READ(12,*)ERRNDC       
      WRITE(7,1200)POP,PCE,PRICE,PWN,PSN,WOOL,PW,PS,TOTFIB,ERRNDC      
1200  FORMAT(/' NDC  ',10F8.3)   
1100  FORMAT(7F8.0)     
1001  FORMAT(24X,2F8.0) 
C  USES THE FAMOUS FISK FILES   
      Y(1,1)=DLOG(WOOL/POP)     
      Y(2,1)=DLOG(TOTFIB/POP)   
      F(1,1)=1.0
      F(2,1)=0.0
      F(1,2)=0.0
      F(2,2)=1.0
      Z1=DLOG(PCE/(PRICE*POP))      
       IF(ITIME.GT.1)DW(3)=0.05*(Z1-F(1,3))**2
       F(1,3)=Z1


      F(2,3)=F(1,3)     
      Z2=DLOG(PW)   
       IF(ITIME.GT.1)DW(4)=0.05*(Z2-F(1,4))**2
       F(1,4)=Z2
      F(2,4)=0.0
      Z3=DLOG(PS)   
       IF(ITIME.GT.1)DW(5)=0.05*(Z3-F(1,5))**2
       F(1,5)=Z3
      F(2,5)=0.0
9000  CONTINUE  
      CPW=DLOG(PW/PWN)  
      CPS=DLOG(PS/PSN)  
      PW=PWN    
      PS=PSN    
C  ERRNDC IS READ AS AVERAGE % ERROR    
      CALL TOVAR(ERRNDC)
      BIGV(1,1)=ERRNDC  
      BIGV(2,2)=ERRNDC  
      DO 11 I=1,NPARS
      DO 11 J=1,NPARS
      BIGW(I,J)=0.0     
      IF(I.EQ.J)BIGW(I,I)=DW(I)
11     CONTINUE  
C     .01 MEANS SIGMA=0.1 = LIKELY SIZE OF CHANGE IN PARAMETER  
C   BIGW CONTROLS THE DEGRADATION IN PARAMETER ESTIMATES
      RETURN    
3000  CONTINUE  
C  FORECAST TIME
      STOP 'NO FORECASTS'       
      END