      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL DFP,FUNC1
      CHARACTER*8 ALABEL(114),ALAB2(8)
      DIMENSION X(114)
      PARAMETER (M=177,NX=8,NXP=8,NZ=8,NZP=8,NC=3,NCP=3,NPX=114,
     1  NPSQ=NPX*NPX,DLT=0.001)      
C
      PARAMETER(LG1=6,LG2=12,LG3=18,LG4=24,N11=LG1+2,N12=LG2+2,
     1  N13=LG3+2,N14=LG4+2,N21=LG1,N22=LG2,N23=LG3,N24=LG4,
     2  N41=N21,N42=N22,N43=N23,N44=N24,N51=LG1+1,N52=LG2,N53=LG3,
     3  N54=LG4,N61=LG1,N62=LG2,N63=LG3,N64=LG4,N71=LG1,
     4  N72=LG2,N73=LG3,N74=LG4,N81=LG1+8,N82=LG2+8,N83=LG3+8,
     5  N84=LG4+8,NS1=LG1,NS2=LG2,NS3=LG3,NS4=LG4)
C
      DIMENSION HES(NPSQ),HES2(NPX,NPX),STD(NPX),TRAT(NPX),NDF(8),
     1  TEST(8),PVAL(8),RH11(N11),RH12(N12),RH13(N13),RH14(N14),
     2  RH21(N21),RH22(N22),RH23(N23),RH24(N24),RH31(N21),RH32(N22),
     3  RH33(N23),RH34(N24),RH41(N41),RH42(N42),RH43(N43),RH44(N44),
     4  RH51(N51),RH52(N52),RH53(N53),RH54(N54),RH61(N61),RH62(N62),
     5  RH63(N63),RH64(N64),RH71(N71),RH72(N72),RH73(N73),RH74(N74),
     6  RH81(N81),RH82(N82),RH83(N83),RH84(N84),RHS11(NS1),
     7  RHS12(NS2),RHS13(NS3),RHS14(NS4),RHS21(NS1),RHS22(NS2),
     8  RHS23(NS3),RHS24(NS4),RHS31(NS1),RHS32(NS2),RHS33(NS3),
     9  RHS34(NS4)
C
      DIMENSION DWIBS(M),DCATF(M),DHOGF(M),DCF(M),DCATS(M),DHOGS(M),
     1  DCS(M),AT(M),ES1(M),ES2(M),ES3(M),EV1(M),EV2(M),EV3(M)
C
      DIMENSION GRD1(NPX),GRD2(NPX),GRD3(NPX),GRD4(NPX),GRD5(NPX),
     1  GRD6(NPX),GRD7(NPX),GRD8(NPX),DRV1(32),DRV2(32),DRV3(32),
     2  DRV4(32),DRV5(32),DRV6(32),DRV7(32),DRV8(32),WGHTS(NZ),
     3  WGHTT(NZ),WGHTW(NZ),STORWT(NZ),XX(32),GRDMT(NZ,NPX),
     4  GRDM(NPX,NZ),TEMP1(NZ,NPX),TEMP2(NZ,NZ),TEMP3(NZ,NZ),XK(32),
     5  DIFWT(NZ),CHOLO(NZ,NZ),TEM1(NX,NX),GRDMT2(4,NPX),GRDM2(NPX,4),
     6  TEMP4(4,NPX),TEMP5(4,4),TEMP6(4,4),GRDMT3(3,NPX),GRDM3(NPX,3),
     7  TEMP7(3,NPX),TEMP8(3,3),TEMP9(3,3),GRDMT4(3,NPX),GRDM4(NPX,3),
     8  TEMP10(3,NPX),TEMP11(3,3),TEMP12(3,3),GRDMT5(2,NPX),
     9  GRDM5(NPX,2),TEMP13(2,NPX),TEMP14(2,2),TEMP15(2,2),CHIW(NZ),
     1  PVW(NZ)
C
      DIMENSION CHOLE(NC,NC),CHOLET(NC,NC),CHOLS(NC,NC),CHOLST(NC,NC),
     1  PSI(NC,NC),PSIT(NC,NC),BLK1(NC,NC),BLK2(NC,NC),QA(NC),
     2  QQ(NC,NC),STORE(NX,NX),INDX(NX),OUT(NX,NX),OMEGA2(NZ,NZ),
     3  OMEGA2I(NZ,NZ),INDX2(NZ),SIGZY(NZ),STOR2(NZ,NZ),OUT2(NZ,NZ),
     4  INDX3(4),WGHT3(4),INDX4(2),WGHT4(2),WGHT5(2),WGHT6(2),INDX6(2)
C
      DIMENSION YEAR(M),CATF(M),CATS(M),HOGF(M),HOGS(M),CF(M),CS(M),
     1  SMF(M),SMS(M),WIBF(M),WIBS(M),WIBGD(M),SEAS1(M),SEAS2(M),
     2  SEAS3(M),SEAS4(M),SEAS5(M),HA(M),BP(M),TS(M),TRND(M),FC(M),
     3  EP1(M),EP2(M),EP3(M),EP4(M),EP5(M),EP6(M),EP7(M),EP8(M),
     4  OMEGA(NX,NX),OMEGAI(NX,NX),WIBSE(M),CATFE(M),HOGFE(M),
     5  CFE(M),SMFE(M),CATSE(M),HOGSE(M),CSE(M),SMSE(M),BPE(M),
     6  WIBSX(M),WT1(M),WT2(M),WT3(M),WT4(M),WGHT(NZ),WT(NC),
     7  ETA(NC,NC),SIGMA(NC,NC),ALPHA(M),EEP8(M),UEP8(M)
C
      COMMON/USER1/YEAR,CATF,CATS,HOGF,HOGS,CF,CS,SMF,SMS,WIBF,WIBS,
     1  WIBGD,SEAS1,SEAS2,SEAS3,SEAS4,SEAS5,HA,BP,TS,TRND,FC,EP1,EP2,
     2  EP3,EP4,EP5,EP6,EP7,EP8,OMEGA,OMEGAI,WIBSE,CATFE,HOGFE,CFE,
     3  SMFE,CATSE,HOGSE,CSE,SMSE,BPE,WIBSX,WT1,WT2,WT3,WT4,WGHT,WT,
     4  WIBSV,ETA,SIGMA,ALPHA,EEP8,UEP8
C
      DATA ALABEL/'A0','A1','A2','A3','A4','A5','A6','A7','A8','A9',
     1  'A10','A11','A12','TAU1','B0','B1','B2','B3','B4','B5','B6',
     2  'B7','B8','B9','B10','B11','B12','TAU2','C0','C1','C2','C3',
     3  'C4','C5','C6','C7','C8','C9','C10','C11','C12','TAU3','D0',
     4  'D1','D2','D3','D4','D5','D6','D7','D8','D9','D10','D11','D12',
     5  'TAU4','F0','F1','TAU5','G0','G1','TAU6','H0','H1','TAU7',
     8  'T0','T1','T2','T3','T4','T5','T6','T7','T8','T9','T10','T11',
     9  'T12','T13','T14','T15','T16','T17','SIGA','SIGB','QA1','QA2',
     1  'QA3','CE1','CE21','CE2','CE31','CE32','CE3','CS1','CS21',
     2  'CS2','CS31','CS32','CS3','P1','P12','P13','P21','P2','P23',
     3  'P31','P32','P3','QN','CQS1','CQS2','CQS3','SIGN'/
C
      DATA ALAB2/'WT1','WT2','WT3','WT4','WT5','WT6','WT7','WT8'/
C
      COMMON/BOPT/IVER,LT,IFP,ISP,NLOOP,IST,ILOOP
      COMMON/BFIDIF/FDFRAC,FDMIN
      COMMON/BLNSR/STEP1,STPACC,NLNSR
      COMMON/BDFP/STPMIN,FOPT
      COMMON/BSTACK/AINT(265000)
      COMMON/BSTAK/NQ,NTOP
      COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
      COMMON/BSTOP/NVAR1,ISTOP(3)
      COMMON/BTRAT/ITRFLG
      COMMON/BINPUT/INFLG
      COMMON/BREAD/NREAD
      COMMON/BNOBS/NOBS
C
C     HMLST40A.PRN  - CONTAINS THE STARTING VALUES FOR THE PARAMETERS
C                     TO BE USED TO INITIALIZE THE FIML ESTIMATION
C
C     BROILD3.PRN  - CONTAINS THE BI-MONTHLY BROILER PRICE AND PROD
C                    DATA, ALONG WITH LIVESTOCK FUTURES AND SPOT 
C                    PRICES TO BE USED IN THE ECONOMETRIC (FIML) 
C                    ESTIMATION OF HAMILTON-TYPE MODEL OF BROILER 
C                    PRICE EXPECTATIONS, 1966-95.
C
      OPEN(6,FILE='C:\DATA\HMLST40A.PRN',STATUS='OLD')
      OPEN(7,FILE='C:\DATA\BROILD3.PRN',STATUS='OLD')
      OPEN(8,FILE='C:\DATA\GQML.DAT',STATUS='UNKNOWN')
      OPEN(9,FILE='C:\DATA\HMAT.DAT',STATUS='UNKNOWN')
      OPEN(10,FILE='C:\DATA\PARM.DAT',STATUS='UNKNOWN')
      OPEN(11,FILE='C:\DATA\MLPARM.DAT',STATUS='UNKNOWN')
      OPEN(12,FILE='C:\DATA\RESID.PRN',STATUS='UNKNOWN')
      OPEN(13,FILE='C:\DATA\JUNK2.PRN',STATUS='UNKNOWN')
      OPEN(14,FILE='C:\DATA\WIBSX.PRN',STATUS='UNKNOWN')
      OPEN(15,FILE='C:\DATA\WTX.PRN',STATUS='UNKNOWN')
C
      CALL DFLT 
      NFILE=8
      NREAD=9
      NPUNCH=9
C
      READ(6,*,END=10)(X(I),I=1,NPX)
C
10    READ(7,*,END=11)(YEAR(I),CATF(I),CATS(I),HOGF(I),HOGS(I),CF(I),
     1  CS(I),SMF(I),SMS(I),WIBF(I),WIBS(I),WIBGD(I),SEAS1(I),
     2  SEAS2(I),SEAS3(I),SEAS4(I),SEAS5(I),HA(I),BP(I),TS(I),
     3  TRND(I),I=1,M)
C
11    WRITE(13,12)(YEAR(I),CATF(I),CATS(I),HOGF(I),HOGS(I),I=1,M)
12    FORMAT(F8.0,2X,4F12.4)
C
      DO I=1,175
        X1=DEXP(CS(I)/100.)
        X2=DEXP(SMS(I)/100.)
        FC(I)=DLOG((70.*X1/.56)+(30.*X2/20.))*100.
      END DO
C
      NOBS=0
C
C     SET PARAMETERS FOR OPT
C
      ACC=1.D-7
      ITERL=600
      IVER=3
      IST=2
      IFP=1
      FDFRAC=0.0001
      MX=1
      NP=NPX
      NQ=265000
      ITRFLG=1
C
C     FOR NEXT 5 CARDS SEE HANDBOOK
C
      CALL OPT(X,NP,FU0,DFP,ITERL,MX,IER,ACC,FUNC1,ALABEL)
      INFLG=1
      CALL PUNCH(X,NP,0)
      CALL OPTMOV(1,HES,NPSQ)
      K=0
      DO I=1,NP
        DO J=1,NP
          K=K+1
          HES2(I,J)=-HES(K)
        END DO
      END DO
      DO I=1,NP
        STD(I)=DSQRT(HES2(I,I))
        TRAT(I)=X(I)/STD(I)
      END DO
      WRITE(11,30)(ALABEL(I),X(I),STD(I),TRAT(I),I=1,NP)
30    FORMAT(' ',A8,2X,3F24.12)
      WRITE(11,35) IER
35    FORMAT('IER FROM OPT =',5X,I6)
      WRITE(11,40) FU0
40    FORMAT('THE FUNCTION VALUE IS=',2X,F24.12/)
      WRITE(10,45)(X(I),I=1,NP)
45    FORMAT(2X,F24.12)
      WRITE(12,50)(YEAR(I),EP1(I),EP2(I),EP3(I),EP4(I),EP5(I),EP6(I),
     1  EP7(I),EP8(I),I=11,175)
50    FORMAT(F6.1,2X,8F16.8)
C
      WRITE(14,60)(YEAR(I),WIBS(I),WIBSE(I),WT1(I),WT2(I),WT3(I),
     1  WT4(I),WIBSX(I),ALPHA(I),EP8(I),EEP8(I),UEP8(I),I=11,175)
60    FORMAT(F6.1,2X,11F16.4)
C
      WRITE(15,65)(WT1(I),WT2(I),WT3(I),WT4(I),I=11,175)
65    FORMAT(4F24.12)
C
      WRITE(14,70)
70    FORMAT(//'THE WEIGHTS IN THE Ex-Ante PRICE EXPECTATION'/)
      WRITE(14,80)(ALAB2(I),WGHT(I),I=1,NZ)
80    FORMAT(' ',A8,2X,F16.8)
C
      WRITE(14,90)
90    FORMAT(//'THE VARIANCE OF THE Ex-Ante PRICE EXPECTATION'/)
      WRITE(14,100)WIBSV
100   FORMAT(10X,'Variance WIBSX =',2X,F16.8)
C
      WRITE(14,110)
110   FORMAT(//'The Submatrix Sigma'/)
      WRITE(14,130)((ETA(I,J),J=1,NC),I=1,NC)
C
      WRITE(14,120)
120   FORMAT(//'The Submatrix S'/)
      WRITE(14,130)((SIGMA(I,J),J=1,NC),I=1,NC)
130   FORMAT(3(3F16.8/))
C
      WRITE(14,140)
140   FORMAT(//'The Variance-Covariance Matrix Omega'/)
      WRITE(14,150)((OMEGA(I,J),J=1,NX),I=1,NX)
150   FORMAT(8(8F16.8/))
C
C     NOW FIND THE DETERMINANT OF THE VAR-COV MATRIX OMEGA
C
      DO I=1,NX
        DO J=1,NX
          TEM1(I,J)=OMEGA(I,J)
        END DO
      END DO  
C
      CALL DETER(TEM1,NX,NXP,INDX,DETX)
C
      WRITE(14,151)
151   FORMAT(//'The det of the Variance-Covariance Matrix Omega'/)
      WRITE(14,152)DETX
152   FORMAT(F24.8/)
C
      CALL CHOLESK(OMEGA,CHOLO,NZ)
      WRITE(14,153)
153   FORMAT(//'The Cholesky Decompositon of Covariance Matrix Omega'/)
      WRITE(14,154)((CHOLO(I,J),J=1,NX),I=1,NX)
154   FORMAT(8(8F16.8/))
C
C     STORE THE WEIGHTS
C
      DO I=1,NZ
        STORWT(I)=WGHT(I)
      END DO
C
C     CREATE FIRST DIFFERENCED VARIABLES
C
      DO I=3,175
        DWIBS(I)=WIBS(I)-WIBS(I-1)
        DCATF(I)=CATF(I-1)-CATF(I-2)
        DHOGF(I)=HOGF(I-1)-HOGF(I-2)
        DCF(I)=CF(I-1)-CF(I-2)
        DCATS(I)=CATS(I)-CATS(I-1)
        DHOGS(I)=HOGS(I)-HOGS(I-1)
        DCS(I)=CS(I)-CS(I-1)
        WIBSE(I)=WIBSE(I)-WIBS(I-1)
      END DO
C
C     COMPUTE THE R-SQUARES AND OTHER MEASURES OF FIT FOR THE
C     ESTIMATED EQUATIONS
C
      CALL OLS(DWIBS,WIBSE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(DWIBS,WIBSE,M,PRMSE)
      CALL DURBIN(DWIBS,EP1,M,RA2,RHOA,DWA)
      CALL DWAT(EP1,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-10.375))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N11
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N12
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N13
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N14
	ENDIF
      END DO
C
      CALL BOX(EP1,RH11,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP1,RH12,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP1,RH13,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP1,RH14,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I)-2.)/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,160)
160   FORMAT(//'THE REGRESSION EQUATION FOR WHOLESALE BROILER PRICE'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,161)
161   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,162)
162   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(DCATF,CATFE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(DCATF,CATFE,M,PRMSE)
      CALL DURBIN(DCATF,EP2,M,RA2,RHOA,DWA)
      CALL DWAT(EP2,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-10.375))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N21
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N22
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N23
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N24
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EP2,RH21,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP2,RH22,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP2,RH23,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP2,RH24,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,170)
170   FORMAT('THE REGRESSION EQUATION FOR LIVE CATTLE f(t-1)'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,171)
171   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,172)
172   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(DHOGF,HOGFE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(DHOGF,HOGFE,M,PRMSE)
      CALL DURBIN(DHOGF,EP3,M,RA2,RHOA,DWA)
      CALL DWAT(EP3,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-10.375))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N21
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N22
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N23
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N24
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EP3,RH31,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP3,RH32,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP3,RH33,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP3,RH34,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,180)
180   FORMAT('THE REGRESSION EQUATION FOR LIVE HOGS f(t-1)'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,181)
181   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,182)
182   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(DCF,CFE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(DCF,CFE,M,PRMSE)
      CALL DURBIN(DCF,EP4,M,RA2,RHOA,DWA)
      CALL DWAT(EP4,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-17.6))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N41
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N42
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N43
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N44
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EP4,RH41,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP4,RH42,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP4,RH43,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP4,RH44,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,190)
190   FORMAT('THE REGRESSION EQUATION FOR CORN f(t-1)'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,191)
191   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,192)
192   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(DCATS,CATSE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(DCATS,CATSE,M,PRMSE)
      CALL DURBIN(DCATS,EP5,M,RA2,RHOA,DWA)
      CALL DWAT(EP5,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-10.375))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N51
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N52
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N53
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N54
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EP5,RH51,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP5,RH52,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP5,RH53,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP5,RH54,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,210)
210   FORMAT('THE REGRESSION EQUATION FOR CATTLE s(t)'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,211)
211   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,212)
212   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(DHOGS,HOGSE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(DHOGS,HOGSE,M,PRMSE)
      CALL DURBIN(DHOGS,EP6,M,RA2,RHOA,DWA)
      CALL DWAT(EP6,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-10.375))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N61
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N62
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N63
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N64
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EP6,RH61,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP6,RH62,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP6,RH63,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP6,RH64,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,220)
220   FORMAT('THE REGRESSION EQUATION FOR HOGS s(t)'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,221)
221   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,222)
222   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(DCS,CSE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(DCS,CSE,M,PRMSE)
      CALL DURBIN(DCS,EP7,M,RA2,RHOA,DWA)
      CALL DWAT(EP7,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-10.375))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N71
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N72
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N73
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N74
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EP7,RH71,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP7,RH72,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP7,RH73,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP7,RH74,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,230)
230   FORMAT('THE REGRESSION EQUATION FOR CORN s(t)'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,231)
231   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,232)
232   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(BP,BPE,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(BP,BPE,M,PRMSE)
      CALL DURBIN(BP,EP8,M,RA2,RHOA,DWA)
      CALL DWAT(EP8,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-10.375))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=N81
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=N82
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=N83
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=N84
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EP8,RH81,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EP8,RH82,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EP8,RH83,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EP8,RH84,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I)-8.)/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,250)
250   FORMAT('THE REGRESSION EQUATION FOR BROILER PRODUCTION'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,251)
251   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,252)
252   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      WRITE(14,400)
400   FORMAT('THE SAMPLE AUTOCORRELATIONS BY EQUAITON'/)
      WRITE(14,410)(I,RH14(I),RH24(I),RH34(I),RH44(I),RH54(I),RH64(I),
     1  RH74(I),RH84(I),I=1,24)
410   FORMAT(2X,I3,2X,8F12.3)
C
      QA1=X(86)
      QA2=X(87)
      QA3=X(88)
      QN=X(110)
C
C     RECOVER THE TRUE EXPECTATION ERROR a(t)
C
      DO I=11,175
        AT(I)=EP1(I)-ALPHA(I)
        ES1(I)=QA1*AT(I)
        ES2(I)=QA2*AT(I)
        ES3(I)=QA3*AT(I)
        EV1(I)=EP5(I)-ES1(I)
        EV2(I)=EP6(I)-ES2(I)
        EV3(I)=EP7(I)-ES3(I)
      END DO
C
C     COMPUTE R**2s OF SPOT FORECAST ERROR PROJECTIONS
C
      CALL OLS(EP5,ES1,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(EP5,ES1,M,PRMSE)
      CALL DURBIN(EP5,EV1,M,RA2,RHOA,DWA)
      CALL DWAT(EV1,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-1.))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=NS1
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=NS2
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=NS3
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=NS4
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EV1,RHS11,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EV1,RHS12,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EV1,RHS13,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EV1,RHS14,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,255)
255   FORMAT(//'THE PROJECTION FOR SPOT CATTLE PRICE'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,256)
256   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,257)
257   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(EP6,ES2,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(EP6,ES2,M,PRMSE)
      CALL DURBIN(EP6,EV2,M,RA2,RHOA,DWA)
      CALL DWAT(EV2,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-1.))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=NS1
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=NS2
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=NS3
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=NS4
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EV2,RHS21,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EV2,RHS22,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EV2,RHS23,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EV2,RHS24,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,260)
260   FORMAT('THE PROJECTION FOR SPOT HOG PRICE'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,261)
261   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,262)
262   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(EP7,ES3,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(EP7,ES3,M,PRMSE)
      CALL DURBIN(EP7,EV3,M,RA2,RHOA,DWA)
      CALL DWAT(EV3,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-1.))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=NS1
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=NS2
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=NS3
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=NS4
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EV3,RHS31,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EV3,RHS32,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EV3,RHS33,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EV3,RHS34,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,265)
265   FORMAT('THE PROJECTION FOR SPOT CORN PRICE'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,267)
267   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,268)
268   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      WRITE(14,450)
450   FORMAT('SAMPLE AUTOCORRELATIONS FOR SPOT PROJECTION EQNS'/)
      WRITE(14,460)(I,RHS14(I),RHS24(I),RHS34(I),I=1,24)
460   FORMAT(2X,I3,2X,3F12.3)
C
C     COMPUTE FUTURES ERROR PROJECTIONS
C
      DO I=11,175
        ES1(I)=QA1*ALPHA(I)
        ES2(I)=QA2*ALPHA(I)
        ES3(I)=QA3*ALPHA(I)
        EV1(I)=EP2(I)-ES1(I)
        EV2(I)=EP3(I)-ES2(I)
        EV3(I)=EP4(I)-ES3(I)
      END DO
C
C     COMPUTE R**2s OF FUTURES FORECAST ERROR PROJECTIONS
C
      CALL OLS(EP2,ES1,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(EP2,ES1,M,PRMSE)
      CALL DURBIN(EP2,EV1,M,RA2,RHOA,DWA)
      CALL DWAT(EV1,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-1.))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=NS1
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=NS2
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=NS3
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=NS4
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EV1,RHS11,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EV1,RHS12,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EV1,RHS13,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EV1,RHS14,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,275)
275   FORMAT(//'THE PROJECTION FOR FUTURES CATTLE PRICE'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,276)
276   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,277)
277   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(EP3,ES2,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(EP3,ES2,M,PRMSE)
      CALL DURBIN(EP3,EV2,M,RA2,RHOA,DWA)
      CALL DWAT(EV2,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-1.))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=NS1
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=NS2
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=NS3
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=NS4
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EV2,RHS21,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EV2,RHS22,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EV2,RHS23,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EV2,RHS24,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,280)
280   FORMAT('THE PROJECTION FOR FUTURES HOG PRICE'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,281)
281   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,282)
282   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      CALL OLS(EP4,ES3,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(EP4,ES3,M,PRMSE)
      CALL DURBIN(EP4,EV3,M,RA2,RHOA,DWA)
      CALL DWAT(EV3,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-1.))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=NS1
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=NS2
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=NS3
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=NS4
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EV3,RHS31,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EV3,RHS32,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EV3,RHS33,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EV3,RHS34,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,285)
285   FORMAT('THE PROJECTION FOR FUTURES CORN PRICE'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,287)
287   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,288)
288   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
C
      WRITE(14,470)
470   FORMAT('SAMPLE AUTOCORRELATIONS FOR FUTURES PROJECTION EQNS'/)
      WRITE(14,480)(I,RHS14(I),RHS24(I),RHS34(I),I=1,24)
480   FORMAT(2X,I3,2X,3F12.3)
C
C     COMPUTE PRODUCTION ERROR PROJECTIONS
C
      DO I=11,175
        ES1(I)=QN*AT(I)
        EV1(I)=EP8(I)-ES1(I)
      END DO
C
C     COMPUTE R**2s OF PRODUCTION FORECAST ERROR PROJECTIONS
C
      CALL OLS(EP8,ES1,M,AL,BETA,SEA,SEB,R2,DW,RHO)
      CALL RMSE(EP8,ES1,M,PRMSE)
      CALL DURBIN(EP8,EV1,M,RA2,RHOA,DWA)
      CALL DWAT(EV1,M,DWB)
C
C     COMPUTE THE ADJUSTED R**2.
C
      R2B=1.-((165.)/(165.-1.))*(1.-R2)
C
C     COMPUTE LJUNG-BOX AND BOX-PIERCE Q STATISTICS
C
      DO I=1,8
	IF (I .EQ. 1 .OR. I .EQ. 2) THEN
	  NDF(I)=NS1
	ENDIF
	IF (I .EQ. 3 .OR. I .EQ. 4) THEN
	  NDF(I)=NS2
	ENDIF
	IF (I .EQ. 5 .OR. I .EQ. 6) THEN
	  NDF(I)=NS3
	ENDIF
	IF (I .EQ. 7 .OR. I .EQ. 8) THEN
	  NDF(I)=NS4
	ENDIF
        TEST(I)=0.
      END DO
C
      CALL BOX(EV1,RHS11,M,NDF(1),TEST(1),TEST(2))
      CALL BOX(EV1,RHS12,M,NDF(3),TEST(3),TEST(4))
      CALL BOX(EV1,RHS13,M,NDF(5),TEST(5),TEST(6))
      CALL BOX(EV1,RHS14,M,NDF(7),TEST(7),TEST(8))
C
      DO I=1,8
	AA=(NDF(I))/2.
	CHI=TEST(I)/2.
	CALL GAMMP(AA,CHI,GVAL)
	PVAL(I)=1.-GVAL
      END DO
C
      WRITE(14,295)
295   FORMAT(//'THE PROJECTION FOR BROILER PRODUCTION'/)
      WRITE(14,300) AL,BETA,R2,R2B,DW,RHO,PRMSE
      WRITE(14,310) SEA,SEB,RA2,RHOA,DWB
C
      WRITE(14,296)
296   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(6) AND Q(12) STATISTICS'/)
      WRITE(14,320) TEST(1),TEST(2),TEST(3),TEST(4) 
      WRITE(14,330) PVAL(1),PVAL(2),PVAL(3),PVAL(4)
      WRITE(14,297)
297   FORMAT('  BOX-PIERCE AND LJUNG-BOX Q(18) AND Q(24) STATISTICS'/)
      WRITE(14,320) TEST(5),TEST(6),TEST(7),TEST(8)
      WRITE(14,340) PVAL(5),PVAL(6),PVAL(7),PVAL(8)
C
      WRITE(14,490)
490   FORMAT('SAMPLE AUTOCORRELATIONS FOR PRODUCTION PROJECTION EQN'/)
      WRITE(14,500)(I,RHS14(I),I=1,24)
500   FORMAT(2X,I3,2X,F12.3)
C
300   FORMAT(F12.5,2X,F12.5,2X,5F12.5)
310   FORMAT(F12.5,2X,F12.5,2X,3F12.5/)
320   FORMAT(5X,4F12.5)
330   FORMAT(5X,4F12.5/)
340   FORMAT(5X,4F12.5//)
C
C     COMPUTE THE GRADIENT OF THE WEIGHT PARAMETERS IN ALPHAt
C
C     FIRST, SET INITIALIZE THE GRADIENTS OF THE WEIGHTS TO ZERO
C
      DO I=1,NPX
        GRD1(I)=0.
        GRD2(I)=0.
        GRD3(I)=0.
        GRD4(I)=0.
        GRD5(I)=0.
        GRD6(I)=0.
        GRD7(I)=0.
        GRD8(I)=0.
      END DO
C
C       SAVE INITIAL PARAMTER VALUES
C
        XK(1)=X(72)
        XK(2)=X(84)
        XK(3)=X(85)
        XK(4)=X(86)
        XK(5)=X(87)
        XK(6)=X(88)
        XK(7)=X(89)
        XK(8)=X(90)
        XK(9)=X(91)
        XK(10)=X(92)
        XK(11)=X(93)
        XK(12)=X(94)
        XK(13)=X(95)
        XK(14)=X(96)
        XK(15)=X(97)
        XK(16)=X(98)
        XK(17)=X(99)
        XK(18)=X(100)
        XK(19)=X(101)
        XK(20)=X(102)
        XK(21)=X(103)
        XK(22)=X(104)
        XK(23)=X(105)
        XK(24)=X(106)
        XK(25)=X(107)
        XK(26)=X(108)
        XK(27)=X(109)
        XK(28)=X(110)
        XK(29)=X(111)
        XK(30)=X(112)
        XK(31)=X(113)
        XK(32)=X(114)
C
C     ITERATE THROUGH THE NUMBER OF PARAMETERS
C
      OUTSIDE: DO K=1,32
C
C     ITERATE ON THE 'NO CHANGE', 'UP' AND 'DOWN' DIFFERENCES
C
      INSIDE: DO JJ=1,3
C
        DO I=1,NZ
          WGHT(I)=0.
        END DO
C
        XX(1)=X(72)
        XX(2)=X(84)
        XX(3)=X(85)
        XX(4)=X(86)
        XX(5)=X(87)
        XX(6)=X(88)
        XX(7)=X(89)
        XX(8)=X(90)
        XX(9)=X(91)
        XX(10)=X(92)
        XX(11)=X(93)
        XX(12)=X(94)
        XX(13)=X(95)
        XX(14)=X(96)
        XX(15)=X(97)
        XX(16)=X(98)
        XX(17)=X(99)
        XX(18)=X(100)
        XX(19)=X(101)
        XX(20)=X(102)
        XX(21)=X(103)
        XX(22)=X(104)
        XX(23)=X(105)
        XX(24)=X(106)
        XX(25)=X(107)
        XX(26)=X(108)
        XX(27)=X(109)
        XX(28)=X(110)
        XX(29)=X(111)
        XX(30)=X(112)
        XX(31)=X(113)
        XX(32)=X(114)
C
        IF (JJ .EQ. 1) THEN
          DELT=0.0000
        ENDIF
        IF (JJ .EQ. 2) THEN
          DELT=DLT
        ENDIF
        IF (JJ .EQ. 3) THEN
          DELT=-DLT
        ENDIF
C
        XX(K)=(1.+DELT)*XX(K)          
C
        T6=XX(1)
        SIGA=XX(2)
        SIGB=XX(3)
        QA1=XX(4)
        QA2=XX(5)
        QA3=XX(6)
        CE1=XX(7)
        CE21=XX(8)
        CE2=XX(9)
        CE31=XX(10)
        CE32=XX(11)
        CE3=XX(12)
        CS1=XX(13)
        CS21=XX(14)
        CS2=XX(15)
        CS31=XX(16)
        CS32=XX(17)
        CS3=XX(18)
        P1=XX(19)
        P12=XX(20)
        P13=XX(21)
        P21=XX(22)
        P2=XX(23)
        P23=XX(24)
        P31=XX(25)
        P32=XX(26)
        P3=XX(27)
        QN=XX(28)
        CQS1=XX(29)
        CQS2=XX(30)
        CQS3=XX(31)
        SIGN=XX(32)
C
C     FILL IN THE CHOLESKY MATRICES THAT WILL BE USED TO ULTIMATELY
C     TO CONSTRUCT THE 'ETA' AND 'SIGMA' MATRICES
C
      DO I=1,NC
        DO J=1,NC
          CHOLE(I,J)=0.
          CHOLET(I,J)=0.
          ETA(I,J)=0.
          CHOLS(I,J)=0.
          CHOLST(I,J)=0.
          SIGMA(I,J)=0.
          BLK1(I,J)=0.
          BLK2(I,J)=0.
        END DO
      END DO
C
      CHOLE(1,1)=CE1
      CHOLS(1,1)=CS1
C
      CHOLE(2,1)=CE21
      CHOLE(2,2)=CE2
      CHOLS(2,1)=CS21
      CHOLS(2,2)=CS2
C
      CHOLE(3,1)=CE31
      CHOLE(3,2)=CE32
      CHOLE(3,3)=CE3
      CHOLS(3,1)=CS31
      CHOLS(3,2)=CS32
      CHOLS(3,3)=CS3
C
C     TAKE THE TRANSPOSE OF THE CHOLESKY MATRICES
C 
      CALL TRANS(CHOLE,CHOLET,NC)
      CALL TRANS(CHOLS,CHOLST,NC)
C
C     MULTIPLY 'CHOLx' WITH 'CHOLxT' TO OBTAIN THE 'ETA' AND 'SIGMA' 
C     SUBMATRICES OF THE MODEL'S VARIANCE-COVARIANCE MATRIX
C
      CALL MULT(CHOLE,CHOLET,ETA,NC)
      CALL MULT(CHOLS,CHOLST,SIGMA,NC)
C
C     CONSTRUCT nX1 VECTOR OF q's
C
      QA(1)=QA1
      QA(2)=QA2
      QA(3)=QA3
C
C     CONSTRUCT THE nxn MATRIX 'PSI' 
C
C     ROW 1.
C
      PSI(1,1)=P1
      PSI(1,2)=P12
      PSI(1,3)=P13
C
C     ROW 2.
C
      PSI(2,1)=P21
      PSI(2,2)=P2
      PSI(2,3)=P23
C
C     ROW 3.
C
      PSI(3,1)=P31
      PSI(3,2)=P32
      PSI(3,3)=P3
C
C     TAKE THE TRANSPOSE OF PSI
C
      CALL TRANS(PSI,PSIT,NC)
C
C     TAKE THE OUTTER PRODUCT OF THE q's AND PUT INTO MATRIX qq
C
      CALL OUTR(QA,QA,QQ,NC)
C
      DO I=1,NC
        DO J=1,NC
          BLK1(I,J)=(SIGB*SIGB)*QQ(I,J)+ETA(I,J)
          BLK2(I,J)=(SIGA*SIGA)*QQ(I,J)+SIGMA(I,J)
        END DO
      END DO
C
C     FILL IN THE ELEMENTS OF THE 8x8 VAR-COV MATRIX, OMEGA
C
C     ROW 1.
C
      OMEGA(1,1)=(SIGA*SIGA)+(SIGB*SIGB)
      OMEGA(1,2)=(SIGB*SIGB)*QA1
      OMEGA(1,3)=(SIGB*SIGB)*QA2
      OMEGA(1,4)=(SIGB*SIGB)*QA3
      OMEGA(1,5)=(SIGA*SIGA)*QA1
      OMEGA(1,6)=(SIGA*SIGA)*QA2
      OMEGA(1,7)=(SIGA*SIGA)*QA3
C
C     COLUMN 1.
C
      OMEGA(2,1)=(SIGB*SIGB)*QA1
      OMEGA(3,1)=(SIGB*SIGB)*QA2
      OMEGA(4,1)=(SIGB*SIGB)*QA3
      OMEGA(5,1)=(SIGA*SIGA)*QA1
      OMEGA(6,1)=(SIGA*SIGA)*QA2
      OMEGA(7,1)=(SIGA*SIGA)*QA3
C
      DO I=2,4
        DO J=2,4
          OMEGA(I,J)=BLK1(I-1,J-1)
        END DO
        DO J=5,7
          OMEGA(I,J)=PSI(I-1,J-4)
        END DO
      END DO
C
      DO I=5,7
        DO J=2,4
          OMEGA(I,J)=PSIT(I-4,J-1)
        END DO
        DO J=5,7
          OMEGA(I,J)=BLK2(I-4,J-4)
        END DO
      END DO
C
      OMEGA(1,8)=T6*(SIGB*SIGB)+QN*(SIGA*SIGA)
      OMEGA(2,8)=T6*(SIGB*SIGB)*QA1
      OMEGA(3,8)=T6*(SIGB*SIGB)*QA2
      OMEGA(4,8)=T6*(SIGB*SIGB)*QA3
      OMEGA(5,8)=(SIGA*SIGA)*QN*QA1+CQS1
      OMEGA(6,8)=(SIGA*SIGA)*QN*QA2+CQS2
      OMEGA(7,8)=(SIGA*SIGA)*QN*QA3+CQS3
      OMEGA(8,8)=(T6*T6)*(SIGB*SIGB)+(SIGN*SIGN)
C
      DO I=1,7
        OMEGA(8,I)=OMEGA(I,8)
      END DO
C
C     COMPUTE THE INVERSE OF THE VAR-COV MATRIX OMEGA
C
      DO I=1,NX
        DO J=1,NX
          STORE(I,J)=OMEGA(I,J)
        END DO
      END DO
C
C     CALL THE SUBROUTINE TO TAKE THE INVERSE OF STORE
C
      CALL INV(STORE,NX,NXP,INDX,OUT)
C
      DO I=1,NX
        DO J=1,NX
          OMEGAI(I,J)=OUT(I,J)
        END DO
      END DO
C
C     COMPUTE THE WEIGHTS TO USE IN THE OPTIMAL Ex Ante PRICE
C     EXPECTATIONS
C
      DO I=1,NZ
        DO J=1,NZ
          OMEGA2(I,J)=OMEGA(I,J)
        END DO
      END DO
C
      DO I=1,NZ
        DO J=1,NZ
          STOR2(I,J)=OMEGA2(I,J)
        END DO
      END DO
C
      CALL INV(STOR2,NZ,NZP,INDX2,OUT2) 
C
      DO I=1,NZ
        DO J=1,NZ
          OMEGA2I(I,J)=OUT2(I,J)
        END DO
      END DO
C
      DO I=1,NZ
        SIGZY(I)=0.
      END DO
C
      SIGZY(1)=SIGB*SIGB
      DO I=1,NC
        SIGZY(I+1)=(SIGB*SIGB)*QA(I)
      END DO      
      SIGZY(8)=T6*(SIGB*SIGB)
C
      CALL PSTMULT(SIGZY,OMEGA2I,WGHT,NZ)
C
      IF (JJ .EQ. 1) THEN
        W0T1=WGHT(1)
        W0T2=WGHT(2)
        W0T3=WGHT(3)
        W0T4=WGHT(4)
        W0T5=WGHT(5)
        W0T6=WGHT(6)
        W0T7=WGHT(7)
        W0T8=WGHT(8)
C
        DO I=1,NZ
          DIFWT(I)=WGHT(I)-STORWT(I)
        END DO
C
      ENDIF
C
      IF (JJ .EQ. 2) THEN
        WUT1=WGHT(1)
        WUT2=WGHT(2)
        WUT3=WGHT(3)
        WUT4=WGHT(4)
        WUT5=WGHT(5)
        WUT6=WGHT(6)
        WUT7=WGHT(7)
        WUT8=WGHT(8)
      ENDIF
C
      IF (JJ .EQ. 3) THEN
        WDT1=WGHT(1)
        WDT2=WGHT(2)
        WDT3=WGHT(3)
        WDT4=WGHT(4)
        WDT5=WGHT(5)
        WDT6=WGHT(6)
        WDT7=WGHT(7)
        WDT8=WGHT(8)
      ENDIF
C
      END DO INSIDE
C
      DRV1(K)=(((WUT1-W0T1)-(WDT1-W0T1))/(DLT*XK(K)))*0.5
      DRV2(K)=(((WUT2-W0T2)-(WDT2-W0T2))/(DLT*XK(K)))*0.5
      DRV3(K)=(((WUT3-W0T3)-(WDT3-W0T3))/(DLT*XK(K)))*0.5
      DRV4(K)=(((WUT4-W0T4)-(WDT4-W0T4))/(DLT*XK(K)))*0.5
      DRV5(K)=(((WUT5-W0T5)-(WDT5-W0T5))/(DLT*XK(K)))*0.5
      DRV6(K)=(((WUT6-W0T6)-(WDT6-W0T6))/(DLT*XK(K)))*0.5
      DRV7(K)=(((WUT7-W0T7)-(WDT7-W0T7))/(DLT*XK(K)))*0.5
      DRV8(K)=(((WUT8-W0T8)-(WDT8-W0T8))/(DLT*XK(K)))*0.5
C
      END DO OUTSIDE
C
      GRD1(72)=DRV1(1)
      GRD2(72)=DRV2(1)
      GRD3(72)=DRV3(1)
      GRD4(72)=DRV4(1)
      GRD5(72)=DRV5(1)
      GRD6(72)=DRV6(1)
      GRD7(72)=DRV7(1)
      GRD8(72)=DRV8(1)
C
      DO I=2,32
        GRD1(82+I)=DRV1(I)
        GRD2(82+I)=DRV2(I)
        GRD3(82+I)=DRV3(I)
        GRD4(82+I)=DRV4(I)
        GRD5(82+I)=DRV5(I)
        GRD6(82+I)=DRV6(I)
        GRD7(82+I)=DRV7(I)
        GRD8(82+I)=DRV8(I)
      END DO
C
      DO J=1,NPX
        GRDMT(1,J)=GRD1(J)
        GRDMT(2,J)=GRD2(J)
        GRDMT(3,J)=GRD3(J)
        GRDMT(4,J)=GRD4(J)
        GRDMT(5,J)=GRD5(J)
        GRDMT(6,J)=GRD6(J)
        GRDMT(7,J)=GRD7(J)
        GRDMT(8,J)=GRD8(J)
        GRDMT2(1,J)=GRD1(J)
        GRDMT2(2,J)=GRD5(J)
        GRDMT2(3,J)=GRD6(J)
        GRDMT2(4,J)=GRD7(J)
        GRDMT2(5,J)=GRD8(J)
        GRDMT3(1,J)=GRD2(J)
        GRDMT3(2,J)=GRD3(J)
        GRDMT3(3,J)=GRD4(J)
        GRDMT4(1,J)=GRD5(J)
        GRDMT4(2,J)=GRD6(J)
        GRDMT4(3,J)=GRD7(J)
        GRDMT5(1,J)=GRD1(J)
        GRDMT5(2,J)=GRD8(J)
      END DO
C
C     TAKE THE TRANSPOSE OF THE 8x114 MATRIX GRDMT
C
      CALL TRANS2(GRDMT,GRDM,NZ,NPX)
      CALL TRANS2(GRDMT2,GRDM2,4,NPX)
      CALL TRANS2(GRDMT3,GRDM3,3,NPX)
      CALL TRANS2(GRDMT4,GRDM4,3,NPX)
      CALL TRANS2(GRDMT5,GRDM5,2,NPX)
C
      CALL QUAD(GRD1,HES2,WGHTS(1),NPX)
      CALL QUAD(GRD2,HES2,WGHTS(2),NPX)
      CALL QUAD(GRD3,HES2,WGHTS(3),NPX)
      CALL QUAD(GRD4,HES2,WGHTS(4),NPX)
      CALL QUAD(GRD5,HES2,WGHTS(5),NPX)
      CALL QUAD(GRD6,HES2,WGHTS(6),NPX)
      CALL QUAD(GRD7,HES2,WGHTS(7),NPX)
      CALL QUAD(GRD8,HES2,WGHTS(8),NPX)
C
      DO I=1,NZ
        WGHT(I)=STORWT(I)
        WGHTS(I)=DSQRT(WGHTS(I))
        WGHTT(I)=WGHT(I)/WGHTS(I)
        WGHTW(I)=WGHT(I)*WGHT(I)/(WGHTS(I)*WGHTS(I))
        AA=1/2.
        CHIW(I)=WGHTW(I)/2.
        GVAL=0.
        CALL GAMMP(AA,CHIW(I),GVAL)
        PVW(I)=1.-GVAL
      END DO
C
      WGHT3(1)=WGHT(1)
      WGHT3(2)=WGHT(5)
      WGHT3(3)=WGHT(6)
      WGHT3(4)=WGHT(7)
      WGHT3(5)=WGHT(8)
      WGHT4(1)=WGHT(2)
      WGHT4(2)=WGHT(3)
      WGHT4(3)=WGHT(4)
      WGHT5(1)=WGHT(5)
      WGHT5(2)=WGHT(6)
      WGHT5(3)=WGHT(7)
      WGHT6(1)=WGHT(1)
      WGHT6(2)=WGHT(8)
C
      WRITE(14,700)
700   FORMAT(//'THE WEIGHTS IN THE Ex-Ante PRICE EXPECTATION'/)
      WRITE(14,710)(ALAB2(I),WGHT(I),WGHTS(I),WGHTT(I),DIFWT(I),
     1  WGHTW(I),PVW(I),I=1,NZ)
710   FORMAT(' ',A8,2X,5F16.8,2X,F6.4)
C
C     CONSTRUCT A WALD TEST OF THE RESTRICTIONS THAT ALL OF THE
C     WEIGHTS ARE JOINTLY EQUAL TO ZERO
C
      CALL MULT2(GRDMT,HES2,TEMP1,NZ,NPX,NPX)
      CALL MULT2(TEMP1,GRDM,TEMP2,NZ,NPX,NZ)
      CALL INV(TEMP2,NZ,NZP,INDX2,TEMP3) 
      CALL QUAD(WGHT,TEMP3,WALD1,NZ)
C
      AA=NZ/2.
      CHI=WALD1/2.
      CALL GAMMP(AA,CHI,GVAL)
      PV=1.-GVAL
C
      WRITE(14,720)
720   FORMAT(//'  WALD TEST OF THE JOINT RESTRICTIONS Wt=0 df=6'/)
      WRITE(14,730) WALD1,PV
730   FORMAT(5X,2F24.5)
C
C     CONSTRUCT A WALD TEST OF THE RESTRICTIONS THAT ALL OF THE
C     WEIGHTS ARE JOINTLY EQUAL TO ZERO EXCEPT FOR THOSE ON
C     FUTURES PRICES
C
      CALL MULT2(GRDMT2,HES2,TEMP4,4,NPX,NPX)
      CALL MULT2(TEMP4,GRDM2,TEMP5,4,NPX,4)
      CALL INV(TEMP5,4,4,INDX3,TEMP6) 
      CALL QUAD(WGHT3,TEMP6,WALD2,4)
C
      AA=4./2.
      CHI=WALD2/2.
      CALL GAMMP(AA,CHI,GVAL)
      PV=1.-GVAL
C
      WRITE(14,731)
731   FORMAT(//'  WALD TEST OF ALL BUT f(t-1) Wt=0 df=4'/)
      WRITE(14,732) WALD2,PV
732   FORMAT(5X,2F24.5)
C
C     CONSTRUCT A WALD TEST OF THE RESTRICTIONS THAT ALL OF THE
C     WEIGHTS ON FUTURES PRICES ARE JOINTLY EQUAL TO ZERO 
C
      CALL MULT2(GRDMT3,HES2,TEMP7,3,NPX,NPX)
      CALL MULT2(TEMP7,GRDM3,TEMP8,3,NPX,3)
      CALL INV(TEMP8,3,3,INDX4,TEMP9) 
      CALL QUAD(WGHT4,TEMP9,WALD3,3)
C
      AA=3./2.
      CHI=WALD3/2.
      CALL GAMMP(AA,CHI,GVAL)
      PV=1.-GVAL
C
      WRITE(14,733)
733   FORMAT(//'  WALD TEST OF f(t-1) Wt=0 df=3'/)
      WRITE(14,734) WALD3,PV
734   FORMAT(5X,2F24.5)
C
C     CONSTRUCT A WALD TEST OF THE RESTRICTIONS THAT ALL OF THE
C     WEIGHTS ON SPOT PRICES ARE JOINTLY EQUAL TO ZERO 
C
      CALL MULT2(GRDMT4,HES2,TEMP10,3,NPX,NPX)
      CALL MULT2(TEMP10,GRDM4,TEMP11,3,NPX,3)
      CALL INV(TEMP11,3,3,INDX4,TEMP12) 
      CALL QUAD(WGHT5,TEMP12,WALD4,3)
C
      AA=3./2.
      CHI=WALD4/2.
      CALL GAMMP(AA,CHI,GVAL)
      PV=1.-GVAL
C
      WRITE(14,735)
735   FORMAT(//'  WALD TEST OF s(t) Wt=0 df=3'/)
      WRITE(14,736) WALD4,PV
736   FORMAT(5X,2F24.5)
C
C     CONSTRUCT A WALD TEST OF THE RESTRICTIONS THAT BOTH OF THE
C     WEIGHTS ON PRICE EXPECTATION ERRORS AND PRODUCITON EQUATION
C     ERRORS ARE JOINTLY EQUAL TO ZERO 
C
      CALL MULT2(GRDMT5,HES2,TEMP13,2,NPX,NPX)
      CALL MULT2(TEMP13,GRDM5,TEMP14,2,NPX,2)
      CALL INV(TEMP14,2,2,INDX6,TEMP15) 
      CALL QUAD(WGHT6,TEMP15,WALD5,2)
C
      AA=2./2.
      CHI=WALD5/2.
      CALL GAMMP(AA,CHI,GVAL)
      PV=1.-GVAL
C
      WRITE(14,737)
737   FORMAT(//'  WALD TEST OF u(t) and v(t) Wt=0 df=2'/)
      WRITE(14,738) WALD5,PV
738   FORMAT(5X,2F24.5)
C
      CLOSE(6)
      CLOSE(7)
      CLOSE(8)
      CLOSE(9)
      CLOSE(10)
      CLOSE(11)
      CLOSE(13)
      CLOSE(14)
      CLOSE(15)
      STOP
      END
C
      SUBROUTINE FUNC1(X,NP,FU0,*)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(NX=8,NXP=8,NZ=8,NZP=8,NC=3,NCP=3,M=177,
     1  PI=3.1415926535898)
C
      DIMENSION X(NP),CHOLE(NC,NC),CHOLET(NC,NC),CHOLS(NC,NC),
     1  CHOLST(NC,NC),PSI(NC,NC),PSIT(NC,NC),BLK1(NC,NC),BLK2(NC,NC),
     2  QA(NC),QQ(NC,NC),STORE(NX,NX),INDX(NX),OUT(NX,NX),EPS(NX),
     3  TEM1(NX,NX),TEM2(NX,NX),OMEGA2(NZ,NZ),OMEGA2I(NZ,NZ),INDX2(NZ),
     4  SIGZY(NZ),STOR2(NZ,NZ),OUT2(NZ,NZ)
C
      DIMENSION YEAR(M),CATF(M),CATS(M),HOGF(M),HOGS(M),CF(M),CS(M),
     1  SMF(M),SMS(M),WIBF(M),WIBS(M),WIBGD(M),SEAS1(M),SEAS2(M),
     2  SEAS3(M),SEAS4(M),SEAS5(M),HA(M),BP(M),TS(M),TRND(M),FC(M),
     3  EP1(M),EP2(M),EP3(M),EP4(M),EP5(M),EP6(M),EP7(M),EP8(M),
     4  OMEGA(NX,NX),OMEGAI(NX,NX),WIBSE(M),CATFE(M),HOGFE(M),
     5  CFE(M),SMFE(M),CATSE(M),HOGSE(M),CSE(M),SMSE(M),BPE(M),
     6  WIBSX(M),WT1(M),WT2(M),WT3(M),WT4(M),WGHT(NZ),WT(NC),
     7  ETA(NC,NC),SIGMA(NC,NC),ALPHA(M),EEP8(M),UEP8(M)
C
C     NEXT CARD SHOWS HOW YOU COMMUNICATE DATA TO THE FUNCTION
C     SUBROUTINE
C
      COMMON/USER1/YEAR,CATF,CATS,HOGF,HOGS,CF,CS,SMF,SMS,WIBF,WIBS,
     1  WIBGD,SEAS1,SEAS2,SEAS3,SEAS4,SEAS5,HA,BP,TS,TRND,FC,EP1,EP2,
     2  EP3,EP4,EP5,EP6,EP7,EP8,OMEGA,OMEGAI,WIBSE,CATFE,HOGFE,CFE,
     3  SMFE,CATSE,HOGSE,CSE,SMSE,BPE,WIBSX,WT1,WT2,WT3,WT4,WGHT,WT,
     4  WIBSV,ETA,SIGMA,ALPHA,EEP8,UEP8
C
C     DECLARE THE PARAMETER NAMES
C
      A0=X(1)
      A1=X(2)
      A2=X(3)
      A3=X(4)
      A4=X(5)
      A5=X(6)
      A6=X(7)
      A7=X(8)
      A8=X(9)
      A9=X(10)
      A10=X(11)
      A11=X(12)
      A12=X(13)
      TAU1=X(14)
      B0=X(15)
      B1=X(16)
      B2=X(17)
      B3=X(18)
      B4=X(19)
      B5=X(20)
      B6=X(21)
      B7=X(22)
      B8=X(23)
      B9=X(24)
      B10=X(25)
      B11=X(26)
      B12=X(27)
      TAU2=X(28)
      C0=X(29)
      C1=X(30)
      C2=X(31)
      C3=X(32)
      C4=X(33)
      C5=X(34)
      C6=X(35)
      C7=X(36)
      C8=X(37)
      C9=X(38)
      C10=X(39)
      C11=X(40)
      C12=X(41)
      TAU3=X(42)
      D0=X(43)
      D1=X(44)
      D2=X(45)
      D3=X(46)
      D4=X(47)
      D5=X(48)
      D6=X(49)
      D7=X(50)
      D8=X(51)
      D9=X(52)
      D10=X(53)
      D11=X(54)
      D12=X(55)
      TAU4=X(56)
      F0=X(57)
      F1=X(58)
      TAU5=X(59)
      G0=X(60)
      G1=X(61)
      TAU6=X(62)
      H0=X(63)
      H1=X(64)
      TAU7=X(65)
      T0=X(66)*100.
      T1=X(67)
      T2=X(68)
      T3=X(69)
      T4=X(70)
      T5=X(71)
      T6=X(72)
      T7=X(73)
      T8=X(74)
      T9=X(75)
      T10=X(76)
      T11=X(77)
      T12=X(78)
      T13=X(79)
      T14=X(80)
      T15=X(81)
      T16=X(82)
      T17=X(83)
      SIGA=X(84)
      SIGB=X(85)
      QA1=X(86)
      QA2=X(87)
      QA3=X(88)
      CE1=X(89)
      CE21=X(90)
      CE2=X(91)
      CE31=X(92)
      CE32=X(93)
      CE3=X(94)
      CS1=X(95)
      CS21=X(96)
      CS2=X(97)
      CS31=X(98)
      CS32=X(99)
      CS3=X(100)
      P1=X(101)
      P12=X(102)
      P13=X(103)
      P21=X(104)
      P2=X(105)
      P23=X(106)
      P31=X(107)
      P32=X(108)
      P3=X(109)
      QN=X(110)
      CQS1=X(111)
      CQS2=X(112)
      CQS3=X(113)
      SIGN=X(114)
C
C     FILL IN THE CHOLESKY MATRICES THAT WILL BE USED TO ULTIMATELY
C     TO CONSTRUCT THE 'ETA' AND 'SIGMA' MATRICES
C
      DO I=1,NC
        DO J=1,NC
          CHOLE(I,J)=0.
          CHOLET(I,J)=0.
          ETA(I,J)=0.
          CHOLS(I,J)=0.
          CHOLST(I,J)=0.
          SIGMA(I,J)=0.
          BLK1(I,J)=0.
          BLK2(I,J)=0.
        END DO
      END DO
C
      CHOLE(1,1)=CE1
      CHOLS(1,1)=CS1
C
      CHOLE(2,1)=CE21
      CHOLE(2,2)=CE2
      CHOLS(2,1)=CS21
      CHOLS(2,2)=CS2
C
      CHOLE(3,1)=CE31
      CHOLE(3,2)=CE32
      CHOLE(3,3)=CE3
      CHOLS(3,1)=CS31
      CHOLS(3,2)=CS32
      CHOLS(3,3)=CS3
C
C     TAKE THE TRANSPOSE OF THE CHOLESKY MATRICES
C 
      CALL TRANS(CHOLE,CHOLET,NC)
      CALL TRANS(CHOLS,CHOLST,NC)
C
C     MULTIPLY 'CHOLx' WITH 'CHOLxT' TO OBTAIN THE 'ETA' AND 'SIGMA' 
C     SUBMATRICES OF THE MODEL'S VARIANCE-COVARIANCE MATRIX
C
      CALL MULT(CHOLE,CHOLET,ETA,NC)
      CALL MULT(CHOLS,CHOLST,SIGMA,NC)
C
C     CONSTRUCT nX1 VECTOR OF q's
C
      QA(1)=QA1
      QA(2)=QA2
      QA(3)=QA3
C
C     CONSTRUCT THE nxn MATRIX 'PSI' 
C
C     ROW 1.
C
      PSI(1,1)=P1
      PSI(1,2)=P12
      PSI(1,3)=P13
C
C     ROW 2.
C
      PSI(2,1)=P21
      PSI(2,2)=P2
      PSI(2,3)=P23
C
C     ROW 3.
C
      PSI(3,1)=P31
      PSI(3,2)=P32
      PSI(3,3)=P3
C
C     TAKE THE TRANSPOSE OF PSI
C
      CALL TRANS(PSI,PSIT,NC)
C
C     TAKE THE OUTTER PRODUCT OF THE q's AND PUT INTO MATRIX qq
C
      CALL OUTR(QA,QA,QQ,NC)
C
      DO I=1,NC
        DO J=1,NC
          BLK1(I,J)=(SIGB*SIGB)*QQ(I,J)+ETA(I,J)
          BLK2(I,J)=(SIGA*SIGA)*QQ(I,J)+SIGMA(I,J)
        END DO
      END DO
C
C     FILL IN THE ELEMENTS OF THE 8x8 VAR-COV MATRIX, OMEGA
C
C     ROW 1.
C
      OMEGA(1,1)=(SIGA*SIGA)+(SIGB*SIGB)
      OMEGA(1,2)=(SIGB*SIGB)*QA1
      OMEGA(1,3)=(SIGB*SIGB)*QA2
      OMEGA(1,4)=(SIGB*SIGB)*QA3
      OMEGA(1,5)=(SIGA*SIGA)*QA1
      OMEGA(1,6)=(SIGA*SIGA)*QA2
      OMEGA(1,7)=(SIGA*SIGA)*QA3
C
C     COLUMN 1.
C
      OMEGA(2,1)=(SIGB*SIGB)*QA1
      OMEGA(3,1)=(SIGB*SIGB)*QA2
      OMEGA(4,1)=(SIGB*SIGB)*QA3
      OMEGA(5,1)=(SIGA*SIGA)*QA1
      OMEGA(6,1)=(SIGA*SIGA)*QA2
      OMEGA(7,1)=(SIGA*SIGA)*QA3
C
      DO I=2,4
        DO J=2,4
          OMEGA(I,J)=BLK1(I-1,J-1)
        END DO
        DO J=5,7
          OMEGA(I,J)=PSI(I-1,J-4)
        END DO
      END DO
C
      DO I=5,7
        DO J=2,4
          OMEGA(I,J)=PSIT(I-4,J-1)
        END DO
        DO J=5,7
          OMEGA(I,J)=BLK2(I-4,J-4)
        END DO
      END DO
C
      OMEGA(1,8)=T6*(SIGB*SIGB)+QN*(SIGA*SIGA)
      OMEGA(2,8)=T6*(SIGB*SIGB)*QA1
      OMEGA(3,8)=T6*(SIGB*SIGB)*QA2
      OMEGA(4,8)=T6*(SIGB*SIGB)*QA3
      OMEGA(5,8)=(SIGA*SIGA)*QN*QA1+CQS1
      OMEGA(6,8)=(SIGA*SIGA)*QN*QA2+CQS2
      OMEGA(7,8)=(SIGA*SIGA)*QN*QA3+CQS3
      OMEGA(8,8)=(T6*T6)*(SIGB*SIGB)+(SIGN*SIGN)
C
      DO I=1,7
        OMEGA(8,I)=OMEGA(I,8)
      END DO
C
C     COMPUTE THE INVERSE OF THE VAR-COV MATRIX OMEGA
C
      DO I=1,NX
        DO J=1,NX
          STORE(I,J)=OMEGA(I,J)
        END DO
      END DO
C
C     CALL THE SUBROUTINE TO TAKE THE INVERSE OF STORE
C
      CALL INV(STORE,NX,NXP,INDX,OUT)
C
      DO I=1,NX
        DO J=1,NX
          OMEGAI(I,J)=OUT(I,J)
        END DO
      END DO
C
C     COMPUTE THE WEIGHTS TO USE IN THE OPTIMAL Ex Ante PRICE
C     EXPECTATIONS
C
      DO I=1,NZ
        DO J=1,NZ
          OMEGA2(I,J)=OMEGA(I,J)
        END DO
      END DO
C
      DO I=1,NZ
        DO J=1,NZ
          STOR2(I,J)=OMEGA2(I,J)
        END DO
      END DO
C
      CALL INV(STOR2,NZ,NZP,INDX2,OUT2) 
C
      DO I=1,NZ
        DO J=1,NZ
          OMEGA2I(I,J)=OUT2(I,J)
        END DO
      END DO
C
      DO I=1,NZ
        SIGZY(I)=0.
      END DO
C
      SIGZY(1)=SIGB*SIGB
      DO I=1,NC
        SIGZY(I+1)=(SIGB*SIGB)*QA(I)
      END DO      
      SIGZY(8)=T6*(SIGB*SIGB)
C
      CALL PSTMULT(SIGZY,OMEGA2I,WGHT,NZ)
C
C     COMPUTE THE VARIANCE OF THE PREDICTION
C
      CALL QUAD(SIGZY,OMEGAI,WIBSV,NZ)
C
C     SET SOME COUNTERS TO ZERO
C
      SUM1=0.
C
      ERROR: DO I=11,175
C
C     SPECIFY THE ERROR CORRECTION TERMs
C
C     1.  Error correction term for 12-city wholesale iced broilers
C
      ERC1=WIBS(I-1)-71.1069-0.404432*CATS(I-1)-0.218263*HOGS(I-1)
     1    -0.126369*CS(I-1)
C
C     2.  Error correction terms for futures prices
C
      ERC2=-(CATS(I-1)-CATF(I-2))
      ERC3=-(HOGS(I-1)-HOGF(I-2))
      ERC4=-(CS(I-1)-CF(I-2))
C
C     3.  Error correction terms for spot prices
C
      ERC5=CATS(I-1)-CATF(I-2)
      ERC6=HOGS(I-1)-HOGF(I-2)
      ERC7=CS(I-1)-CF(I-2)
C
C     SPECIFY THE STOCHASTIC EQUATIONS FOR THE MODEL
C
C     1.  Wholesale iced broiler spot price
C
      ER1=A0+A1*SEAS1(I)+A2*SEAS2(I)+A3*SEAS3(I)+A4*SEAS4(I)
     1  +A5*SEAS5(I)
      ER2=A6*(WIBS(I-1)-WIBS(I-2))+A7*(WIBS(I-2)-WIBS(I-3))
     1  +A8*(CATS(I-1)-CATS(I-2))
      ER3=A9*(HOGS(I-1)-HOGS(I-2))+A10*(CS(I-1)-CS(I-2))
      ER4=A11*(BP(I-1)-BP(I-2))+A12*(BP(I-2)-BP(I-3))
     1  +TAU1*ERC1
      WIBSE(I)=WIBS(I-1)+ER1+ER2+ER3+ER4
      EP1(I)=WIBS(I)-WIBSE(I)
C
C     2.  Futrues price f(t-1) for live cattle
C
      ER1=B0+B1*SEAS1(I)+B2*SEAS2(I)+B3*SEAS3(I)+B4*SEAS4(I)
     1  +B5*SEAS5(I)
      ER2=B6*(WIBS(I-1)-WIBS(I-2))+B7*(WIBS(I-2)-WIBS(I-3))
     1  +B8*(CATS(I-1)-CATS(I-2))
      ER3=B9*(HOGS(I-1)-HOGS(I-2))+B10*(CS(I-1)-CS(I-2))
      ER4=B11*(BP(I-1)-BP(I-2))+B12*(BP(I-2)-BP(I-3))
     1  +TAU2*ERC2
      CATFE(I)=ER1+ER2+ER3+ER4
      EP2(I)=CATF(I-1)-CATF(I-2)-CATFE(I)
C
C     3.  Futrues price f(t-1) for live hogs
C
      ER1=C0+C1*SEAS1(I)+C2*SEAS2(I)+C3*SEAS3(I)+C4*SEAS4(I)
     1  +C5*SEAS5(I)
      ER2=C6*(WIBS(I-1)-WIBS(I-2))+C7*(WIBS(I-2)-WIBS(I-3))
     1  +C8*(CATS(I-1)-CATS(I-2))
      ER3=C9*(HOGS(I-1)-HOGS(I-2))+C10*(CS(I-1)-CS(I-2))
      ER4=C11*(BP(I-1)-BP(I-2))+C12*(BP(I-2)-BP(I-3))
     1  +TAU3*ERC3
      HOGFE(I)=ER1+ER2+ER3+ER4
      EP3(I)=HOGF(I-1)-HOGF(I-2)-HOGFE(I)
C
C     4.  Futrues price f(t-1) for corn
C
      ER1=D0+D1*SEAS1(I)+D2*SEAS2(I)+D3*SEAS3(I)+D4*SEAS4(I)
     1  +D5*SEAS5(I)
      ER2=D6*(WIBS(I-1)-WIBS(I-2))+D7*(WIBS(I-2)-WIBS(I-3))
     1  +D8*(CATS(I-1)-CATS(I-2))
      ER3=D9*(HOGS(I-1)-HOGS(I-2))+D10*(CS(I-1)-CS(I-2))
      ER4=D11*(BP(I-1)-BP(I-2))+D12*(BP(I-2)-BP(I-3))
     1  +TAU4*ERC4
      CFE(I)=ER1+ER2+ER3+ER4
      EP4(I)=CF(I-1)-CF(I-2)-CFE(I)
C
C     5.  Spot price s(t) for live cattle
C
      ER1=F0+F1*(CATF(I-1)-CATF(I-2))+TAU5*ERC5
      CATSE(I)=ER1
      EP5(I)=CATS(I)-CATS(I-1)-CATSE(I)
C
C     6.  Spot price s(t) for live hogs
C
      ER1=G0+G1*(HOGF(I-1)-HOGF(I-2))+TAU6*ERC6
      HOGSE(I)=ER1
      EP6(I)=HOGS(I)-HOGS(I-1)-HOGSE(I)
C
C     7.  Spot price s(t) for corn
C
      ER1=H0+H1*(CF(I-1)-CF(I-2))+TAU7*ERC7
      CSE(I)=ER1
      EP7(I)=CS(I)-CS(I-1)-CSE(I)
C
C     8. Broiler production (supply) equation bp(t)
C
      ER1=T0+T1*SEAS1(I)+T2*SEAS2(I)+T3*SEAS3(I)+T4*SEAS4(I)
     1  +T5*SEAS5(I)
      ER2=T6*WIBSE(I)+T7*FC(I-1)+T8*FC(I-2)+T9*HA(I-1)+T10*BP(I-1)
     1  +T11*BP(I-2)
      ER3=+T12*BP(I-3)+T13*BP(I-4)+T14*BP(I-5)+T15*BP(I-6)+T16*BP(I-7)
     1  +T17*BP(I-8)
      BPE(I)=ER1+ER2+ER3
      EP8(I)=BP(I)-BPE(I)
C
C     COMPUTE THE OPTIMAL EXPECTED WHOLSALE ICED BROILER SPOT PRICE,
C     WIBSX
C
C     OPTIMAL PREDICTION OF PEOPLE'S Ex Ante PRICE EXPECTATIONS
C
      WTX1=WIBSE(I)+WGHT(1)*EP1(I)+WGHT(2)*EP2(I)+WGHT(3)*EP3(I)
     1  +WGHT(4)*EP4(I)
      WTX2=WGHT(5)*EP5(I)+WGHT(6)*EP6(I)+WGHT(7)*EP7(I)
     1  +WGHT(8)*EP8(I)
      WIBSX(I)=WTX1+WTX2
C
      ALPHA(I)=WIBSX(I)-WIBSE(I)
C
      EEP8(I)=T6*ALPHA(I)
      UEP8(I)=EP8(I)-EEP8(I)
C
      WT1(I)=WGHT(1)*EP1(I)
      WT2(I)=WGHT(2)*EP2(I)+WGHT(3)*EP3(I)+WGHT(4)*EP4(I)
      WT3(I)=WGHT(5)*EP5(I)+WGHT(6)*EP6(I)+WGHT(7)*EP7(I)
      WT4(I)=WGHT(8)*EP8(I)
C
C     STORE THE CURRENT VALUES OF THE RESIDUALS
C
      EPS(1)=EP1(I)
      EPS(2)=EP2(I)
      EPS(3)=EP3(I)
      EPS(4)=EP4(I)
      EPS(5)=EP5(I)
      EPS(6)=EP6(I)
      EPS(7)=EP7(I)
      EPS(8)=EP8(I)
C
C     EVALUATE THE QUADRATIC TERM (EPS)'*INV(OMEGA)*EPS
C
      DO K=1,NX
        DO L=1,NX
          TEM2(K,L)=OUT(K,L)
        END DO
      END DO
C
      CALL QUAD(EPS,TEM2,QAD,NX)
      SUM1=SUM1+QAD
C
      END DO ERROR
C
C     NOW FIND THE DETERMINANT OF THE VAR-COV MATRIX OMEGA
C
      DO I=1,NX
        DO J=1,NX
          TEM1(I,J)=OMEGA(I,J)
        END DO
      END DO  
C
      CALL DETER(TEM1,NX,NXP,INDX,DET)
      IF (DET .LE. 0.) RETURN 1
C
C     THE LIKELIHOOD FUNCTION
C    
      FU0=-(165.*8./2.)*DLOG(2.*PI)-(165./2.)*DLOG(DET)-0.5*SUM1
C
      RETURN
      END      
C
C  THE FOLLOWING SUBROUTINE PERFORMES L-U DECOMPOSITION
C
      SUBROUTINE LUDCMP(A,NN,NNP,INDX,D)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION A,VV,AAMAX,SUM,DUM,D
      PARAMETER (NMAX=100)
      DIMENSION A(NNP,NNP),INDX(NN),VV(NMAX)
      TINY=1.0D-20
      D=1.
      DO 12 I=1,NN
	AAMAX=0.
	DO 11 J=1,NN
	  IF (DABS(A(I,J)) .GT. AAMAX) AAMAX=DABS(A(I,J))
11      CONTINUE
	IF (AAMAX .EQ. 0.) PAUSE 'Singular matrix.'
	VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,NN
	IF (J .GT. 1) THEN
	  DO 14 I=1,J-1
	    SUM=A(I,J)
	    IF (I .GT. 1) THEN
	      DO 13 K=1,I-1
		SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
	      A(I,J)=SUM
	    ENDIF
14        CONTINUE
	ENDIF
	AAMAX=0.
	DO 16 I=J,NN
	  SUM=A(I,J)
	  IF (J .GT. 1) THEN
	    DO 15 K=1,J-1
	      SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
	    A(I,J)=SUM
	  ENDIF
	  DUM=VV(I)*DABS(SUM)
	  IF (DUM .GE. AAMAX) THEN
	    IMAX=I
	    AAMAX=DUM
	  ENDIF
16      CONTINUE
	IF (J .NE. IMAX) THEN
	  DO 17 K=1,NN
	    DUM=A(IMAX,K)
	    A(IMAX,K)=A(J,K)
	    A(J,K)=DUM
17        CONTINUE
	  D=-D
	  VV(IMAX)=VV(J)
	ENDIF
	INDX(J)=IMAX
	IF (J .NE. NN) THEN
	  IF (A(J,J) .EQ. 0.) A(J,J)=TINY
	  DUM=1./A(J,J)
	  DO 18 I=J+1,NN
	    A(I,J)=A(I,J)*DUM
18        CONTINUE
	ENDIF
19    CONTINUE
      IF (A(NN,NN) .EQ. 0.) A(NN,NN)=TINY
      RETURN
      END
C
C  THE NEXT SUBROUTINE DOES FORWARD AND BACKWARD 
C  SUBSTITUTION
C
      SUBROUTINE LUBKSB(A,NN,NNP,INDX,B)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION A,B,SUM
      DIMENSION A(NNP,NNP),INDX(NN),B(NN)
      II=0
      DO 12 I=1,NN
	LL=INDX(I)
	SUM=B(LL)
	B(LL)=B(I)
	IF (II .NE. 0) THEN
	  DO 11 J=II,I-1
	    SUM=SUM-A(I,J)*B(J)
11        CONTINUE
	ELSE IF (SUM .NE. 0.) THEN
	  II=I
	ENDIF
	B(I)=SUM
12    CONTINUE
      DO 14 I=NN,1,-1
	SUM=B(I)
	IF (I .LT. NN) THEN
	  DO 13 J=I+1,NN
	    SUM=SUM-A(I,J)*B(J)
13        CONTINUE
	ENDIF
	B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END
C
C  THE FOLLOWING SUBROUTINE CALCULATES THE
C  INVERSE OF AN NxN MATRIX.
C
      SUBROUTINE INV(A,NN,NNP,INDX,Y)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION A,Y
      DIMENSION A(NNP,NNP),Y(NNP,NNP),INDX(NNP) 
C
C  FIRST, SET UP THE IDENTITY MATRIX
C
      DO 140 I=1,NN
      DO 130 J=1,NN
      Y(I,J)=0.
130   CONTINUE
      Y(I,I)=1.
140   CONTINUE
C
C  DECOMPOSE THE MATRIX JUST ONCE
C
      CALL LUDCMP(A,NN,NNP,INDX,D)
C
C  FIND THE INVERSE BY COLUMNS
C
      DO 150 J=1,NN
      CALL LUBKSB(A,NN,NNP,INDX,Y(1,J))
150   CONTINUE
C
C  THE MATRIX Y WILL NOW CONTAIAN THE INVERSE
C  OF THE ORIGINAL MATRIX SIGMA, WHICH WILL
C  HAVE BEEN DESTROYED.
C
      RETURN
      END
C
C  THE FOLLOWING SUBROUTINE EVALUATES THE
C  DETERMINANT OF AN NxN MATRIX
C
      SUBROUTINE DETER(A,NN,NNP,INDX,D)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION A,D
      DIMENSION A(NNP,NNP),INDX(NNP)
      CALL LUDCMP(A,NN,NNP,INDX,D)
      DO 10 I=1,NN
      D=D*A(I,I)
10    CONTINUE
      RETURN
      END      
C
C     SUBROUTINE TO MULTIPLY TOGETHER TWO REAL NxN MATRICES.  
C     THE RESULTS ARE PUT IN MATRIX C
C
      SUBROUTINE MULT(A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N,N),B(N,N),C(N,N)
      DO 10 I=1,N
      DO 10 J=1,N
	C(I,J)=0.
10    CONTINUE      
      DO 20 I=1,N
      DO 20 J=1,N
      DO 20 K=1,N
	C(I,J)=C(I,J)+A(I,K)*B(K,J)
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO MULTIPLY TOGETHER A REAL NxM MATRIX AND A REAL
C     MxL MATRIX.  THE RESULTS ARE PUT IN NxL MATRIX C
C
      SUBROUTINE MULT2(A,B,C,N,M,L)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N,M),B(M,L),C(N,L)
      DO 10 I=1,N
      DO 10 J=1,L
	C(I,J)=0.
10    CONTINUE      
      DO 20 I=1,N
      DO 20 J=1,L
      DO 20 K=1,M
	C(I,J)=C(I,J)+A(I,K)*B(K,J)
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO PRE-MULTIPLY A (Nx1) REAL VECTOR B BY A (NxN) REAL
C     MATRIX A.  THE RESULTS ARE PUT IN (Nx1) MATRIX C
C
      SUBROUTINE PREMULT(A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N,N),B(N),C(N)
      DO 10 I=1,N
	C(I)=0.
10    CONTINUE      
      DO 20 I=1,N
      DO 20 J=1,N
	C(I)=C(I)+A(I,J)*B(J)
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO PRE-MULTIPLY A (NxN) REAL MATRIX B BY A (1xN) REAL
C     MATRIX A.  THE RESULTS ARE PUT IN (1xN) MATRIX C
C
      SUBROUTINE PSTMULT(A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N),B(N,N),C(N)
      DO 10 I=1,N
	C(I)=0.
10    CONTINUE      
      DO 20 I=1,N
      DO 20 J=1,N
	C(I)=C(I)+A(J)*B(J,I)
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO PREMULTIPLY AND POSTMULTIPLY THE NxN MATRIX B BY 
C     THE 1xN (resp. Nx1) VECTOR A
C
      SUBROUTINE QUAD(A,B,C,N)      
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NN=200)
      DIMENSION A(N),B(N,N),D(NN)
      SUM1: DO I=1,N
	D(I)=0.
	DO 10 J=1,N
	  D(I)=D(I)+A(J)*B(J,I)
10      CONTINUE
      END DO SUM1
      C=0.
      SUM2: DO I=1,N
	C=C+A(I)*D(I)
      END DO SUM2
      RETURN
      END
C
C     SUBROUTINE TO MULTIPLY TOGETHER A REAL 1xN VECTOR AND 
C     A REAL NxN MATRIX.  THE RESULTS ARE PUT IN 1xN VECTOR C
C
      SUBROUTINE VMULT(A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N),B(N,N),C(N)
      DO 10 I=1,N
	C(I)=0.
10    CONTINUE      
      DO 20 I=1,N
      DO 20 J=1,N
	C(I)=C(I)+A(J)*B(J,I)
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO MULTIPLY TOGETHER A REAL NxN MATRIX AND 
C     A REAL Nx1 VECTOR.  THE RESULTS ARE PUT IN Nx1 VECTOR C
C
      SUBROUTINE VMULT2(A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N),B(N,N),C(N)
      DO 10 I=1,N
	C(I)=0.
10    CONTINUE      
      DO 20 I=1,N
      DO 20 J=1,N
	C(I)=C(I)+B(I,J)*A(J)
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO OBTAIN THE DOT PRODUCT OF A REAL 1xN VECTOR AND 
C     A REAL Nx1 VECTOR.  THE RESULTS ARE PUT IN A SCALAR C
C
      SUBROUTINE VECR(A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N),B(N)
      C=0.
      DO 20 I=1,N
	C=C+A(I)*B(I)
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO OBTAIN THE OUTTER PRODUCT OF A REAL Nx1 VECTOR AND
C     A REAL 1xN VECTOR.  THE RESULTS ARE PUT IN A NxN MATRIX C
C
      SUBROUTINE OUTR(A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N),B(N),C(N,N)
      DO I=1,N
	DO J=1,N
	  C(I,J)=A(I)*B(J)
	END DO
      END DO
      RETURN
      END
C
C     SUBROUTINE TO TAKE THE TANSPOSE OF A NxN MATRIX
C     THE RESULTS ARE PUT IN MATRIX B
C
      SUBROUTINE TRANS(A,B,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N,N),B(N,N)
      DO 10 I=1,N
      DO 10 J=1,N
	B(I,J)=A(J,I)
10    CONTINUE      
      RETURN
      END
C
C     SUBROUTINE TO TAKE THE TANSPOSE OF AN MxK MATRIX
C     THE RESULTS ARE PUT IN A KxM MATRIX B
C
      SUBROUTINE TRANS2(A,B,M,K)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(M,K),B(K,M)
      DO 10 I=1,K
      DO 10 J=1,M
	B(I,J)=A(J,I)
10    CONTINUE      
      RETURN
      END

C
C     SUBROUTINE TO EVALUATE THE TRACE OF A NxN MATRIX
C     THE RESULTS ARE PUT IN SCALAR B
C
      SUBROUTINE TRACE(A,B,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N,N)
      B=0.0
      DO 10 I=1,N
	B=B+A(I,I)
10    CONTINUE      
      RETURN
      END
C
C     SUBROUTINE TO OBTAIN THE CHOLESKY DECOMPOSITION OF A REAL, 
C     SYMMETRIC, NxN MATRIX, A.  THE RESULTS PUT IN LOWER TRIANGULAR
C     MATRIX C
C
      SUBROUTINE CHOLESK(A,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N,N),C(N,N)
C
      DO I=1,N
	DO J=1,N
	  C(I,J)=0.0
	END DO
      END DO
C
      DO I=1,N
	DO J=1,N
	  IF (I .EQ. J) THEN
	    IF (I .EQ. 1) THEN
	      C(I,J)=DSQRT(A(I,J))
	    ENDIF
	    IF (I .GT. 1) THEN
	      SUM1=0.0
	      DO K=2,I
		K2=K-1
		SUM1=SUM1+C(I,K2)*C(I,K2)
	      END DO
	      C(I,J)=DSQRT(A(I,J)-SUM1)
	    ENDIF
	  ENDIF
	  IF (J .GT. I) THEN
	    IF (I .EQ. 1) THEN
	      C(I,J)=A(I,J)/C(I,I)
	      C(J,I)=C(I,J)
	    ENDIF
	    IF (I .GT. 1) THEN
	      SUM1=0.0
	      DO K=2,I
		K2=K-1
		SUM1=SUM1+C(I,K2)*C(J,K2)
	      END DO
	      C(I,J)=(A(I,J)-SUM1)/C(I,I)
	      C(J,I)=C(I,J)
	    ENDIF
	  ENDIF
	END DO
      END DO
C
C     CONVERT C BACK TO LOWER TRIANGULAR
C
      DO I=1,N
	DO J=I+1,N
	  C(I,J)=0.0
	END DO
      END DO
C
      RETURN
      END     
C
C     FUNCTION TO EVALUATE THE EXPONENTIAL OF A SQUARED VARIABLE
C
      DOUBLE PRECISION FUNCTION EX(Z)
      EX=DEXP(-.5*Z*Z)
      RETURN
      END
C
C     FUNCTION TO SQUARE A VARIABLE
C
      DOUBLE PRECISION FUNCTION SQ(Z)
      SQ=Z*Z
      RETURN
      END
C
C     FUNCTION TO CUBE A VARIABLE
C
      DOUBLE PRECISION FUNCTION CU(Z)
      CU=Z*Z*Z
      RETURN
      END
C
C     SUBROUTINE TO EVALUATE AVERAGES
C
      SUBROUTINE AVG(A,N,X)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N)
      SUM1=0.
      DO 10 I=11,175
	SUM1=SUM1+A(I)
10    CONTINUE
      X=SUM1/165.
      RETURN
      END
C
C     YET ANOTHER SUBROUTINE TO EVALUATE AVERAGES
C
      SUBROUTINE AVG2(A,N,X)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION SUM1,X,A
      DIMENSION A(N)
      SUM1=0.
      DO 10 I=8,N
	SUM1=SUM1+A(I)
10    CONTINUE
      X=SUM1/(177.-7.)
      RETURN
      END
C
C     YET EVEN ANOTHER SUBROUTINE TO EVALUATE AVERAGES
C
      SUBROUTINE AVG3(A,N,X,U)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N),U(N)
      SUM1=0.
      DO 10 I=1,N
	SUM1=SUM1+A(I)
10    CONTINUE
      X=SUM1/N
      DO 20 I=1,N
        U(I)=A(I)-X
20    CONTINUE
      RETURN
      END
C
C     SUBROUTINE TO EVALUATE MEAN AND STD. DEVIATIONS
C
      SUBROUTINE STAT(A,N,X,Y)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION T
      DIMENSION A(N)
      SUM1=0.
      SUM2=0.
      T=N-0.0
      DO 10 I=1,N
	SUM1=SUM1+A(I)
10    CONTINUE
      X=SUM1/T
      DO 20 I=1,N
	SUM2=SUM2+(A(I)-X)*(A(I)-X)
20    CONTINUE
      Y=DSQRT(SUM2/(T-1.))
      RETURN
      END
C
C     SUBROUTINE TO EVALUATE MEAN AND STD. DEVIATIONS AND MIN AND MAX
C
      SUBROUTINE STAT2(A,N,X,Y,RMIN,RMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION A,X,Y,RMIN,RMAX,SUM1,SUM2
      PARAMETER (T=697.0)
      DIMENSION A(N)
      SUM1=0.
      SUM2=0.
      DO 10 I=8,N
	SUM1=SUM1+A(I)
	IF (I .EQ. 8) THEN
	  RMIN=A(I)
	  RMAX=A(I)
	ENDIF
	IF (I .GE. 6) THEN
	  IF(A(I) .LE. RMIN) THEN
	    RMIN=A(I)
	  ENDIF
	  IF(A(I) .GE. RMAX) THEN
	    RMAX=A(I)
	  ENDIF
	ENDIF
10    CONTINUE
      X=SUM1/T
      DO 20 I=8,N
	SUM2=SUM2+(A(I)-X)*(A(I)-X)
20    CONTINUE
      Y=DSQRT(SUM2/(T-1.))
      RETURN
      END
C
C     SUBROUTINE TO RETURN LJUNG-BOX AND BOX-PIERCE Q-STATISTICS
C
      SUBROUTINE BOX(EP,RHO,N,M,Q1,Q2)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION SCR
      DIMENSION EP(N),RHO(M)
      NUP=175
      DENOM=0.
      DO I=11,NUP
	DENOM=DENOM+EP(I)*EP(I)
      END DO
      DO K=1,M
	RHO(K)=0.
	KBAR=K+11
	DO I=KBAR,NUP
	  RHO(K)=RHO(K)+EP(I)*EP(I-K)
	END DO
	RHO(K)=RHO(K)/DENOM
      END DO
      Q1=0.
      Q2=0.
      DO K=1,M
	Q1=Q1+RHO(K)*RHO(K)
	SCR=165.-K
	Q2=Q2+RHO(K)*RHO(K)/SCR
      END DO
      Q1=165.*Q1
      Q2=165.*(165.+2.)*Q2
      RETURN
      END
C
C     SUBROUTINE TO RETURN LJUNG-BOX AND BOX-PIERCE Q**2-STATISTICS
C
      SUBROUTINE BOXSQ(EP,RHO,N,M,Q1,Q2)
      IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION SCR
      DIMENSION EP(N),RHO(M)
      VAR=0.
      DO I=8,N
        VAR=VAR+EP(I)
      END DO
      VAR=VAR/697.
      DENOM=0.
      DO I=8,N
	DENOM=DENOM+(EP(I)-VAR)*(EP(I)-VAR)
      END DO
      DO K=1,M
	RHO(K)=0.
	KBAR=K+8
	DO I=KBAR,N
	  RHO(K)=RHO(K)+(EP(I)-VAR)*(EP(I-K)-VAR)
	END DO
	RHO(K)=RHO(K)/DENOM
      END DO
      Q1=0.
      Q2=0.
      DO K=1,M
	Q1=Q1+RHO(K)*RHO(K)
	SCR=697.-K
	Q2=Q2+RHO(K)*RHO(K)/SCR
      END DO
      Q1=697.*Q1
      Q2=697.*(697.+2.)*Q2
      RETURN
      END
C
C     SUBROUTINE TO RETURN TO COMPUTE HOSKING'S MULTIVARIATE 
C     PORTMANTEAU Q STATISTIC
C
      SUBROUTINE HOSK(E1,E2,E3,E4,N,M,PBARH1,PBARH2)
      PARAMETER (MM=4,NN=697,MMP=4)
      DOUBLE PRECISION T,TADJ
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION E1(N),E2(N),E3(N),E4(N),C0(MM,MM),CD(MM,MM),
     1  CDT(MM,MM),C(MM,MM),CI(MM,MM),STOR(MM,MM),STOR2(MM,MM),
     2  STOR3(MM,MM),EPS1(MM),EPS2(MM),INDX(MM)
C
      SUM1=0.0
      SUM2=0.0
      T=697.0
      PORT: DO K=0,M 
	TADJ=T-K
	DO 20 I=1,MM
	  DO 10 J=1,MM
	    C(I,J)=0.0
10        CONTINUE
20      CONTINUE
	JUD=K+1
	SWEEP: DO I=JUD,NN
	  EPS1(1)=E1(I+7)*1000.
	  EPS1(2)=E2(I+7)*1000.
	  EPS1(3)=E3(I+7)*1000.
          EPS1(4)=E4(I+7)*1000.
	  EPS2(1)=E1(I+7-K)*1000.
	  EPS2(2)=E2(I+7-K)*1000.
	  EPS2(3)=E3(I+7-K)*1000.
          EPS2(4)=E4(I+7-K)*1000.
	  CALL OUTR(EPS1,EPS2,STOR,MM)
	  DO 40 II=1,MM
	    DO 30 JJ=1,MM
	      C(II,JJ)=C(II,JJ)+STOR(II,JJ)/T
30          CONTINUE
40        CONTINUE
	END DO SWEEP
	IF (K .EQ. 0) THEN
	  DO I=1,MM
	    DO J=1,MM
	      C0(I,J)=C(I,J)
	    END DO
	  END DO 
	  CALL INV(C,MM,MMP,INDX,C0)
	ENDIF
	IF (K .GE. 1) THEN
	  DO 60 I=1,MM
	    DO 50 J=1,MM
	      CD(I,J)=C(I,J)
	      CI(I,J)=C0(I,J)
50          CONTINUE        
60        CONTINUE
	  CALL TRANS(CD,CDT,MM)
	  CALL MULT(CDT,CI,STOR,MM)
	  DO 80 I=1,MM
	    DO 70 J=1,MM
	      CI(I,J)=C0(I,J)
	      CD(I,J)=C(I,J)
70          CONTINUE        
80        CONTINUE
	  CALL MULT(CD,CI,STOR2,MM)
	  CALL MULT(STOR,STOR2,STOR3,MM)
	  CALL TRACE(STOR3,TRC,MM)
	  SUM1=SUM1+TRC/TADJ
	  SUM2=SUM2+TRC
	ENDIF
      END DO PORT
      PBARH1=T*SUM2
      PBARH2=(T*T)*SUM1
      RETURN
      END
C
C     SUBROUTINE TO EVALUATE PERCENT RMSE
C
      SUBROUTINE RMSE(Y,YHAT,N,PRMSE)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(N),YHAT(N)
      SUM1=0.
      DO I=8,N
	SUM1=SUM1+(Y(I)-YHAT(I))*(Y(I)-YHAT(I))
      END DO
      PRMSE=DSQRT(SUM1/(177.-7.))
      RETURN
      END
C
C     FUNCTION TO COMPUTE CHI-SQUARED VALUES
C
      SUBROUTINE GAMMP(A,X,OUT)
      IMPLICIT REAL*8 (A-H,O-Z)
      IF(X.LT.0..OR.A.LE.0.)PAUSE
      IF(X.LT.A+1.)THEN
	CALL GSER(GAMSER,A,X,GLN)
	GP=GAMSER
      ELSE
	CALL GCF(GAMMCF,A,X,GLN)
	GP=1.-GAMMCF
      ENDIF
      OUT=GP
      RETURN
      END
C
      SUBROUTINE GSER(GAMSER,A,X,GLN)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (ITMAX=100,EPS=3.E-7)
      CALL GAMMLN(A,GLN)
      IF(X.LE.0.)THEN
	IF(X.LT.0.)PAUSE
	GAMSER=0.
	RETURN
      ENDIF
      AP=A
      SUM=1./A
      DEL=SUM
      DO 11 N=1,ITMAX
	AP=AP+1.
	DEL=DEL*X/AP
	SUM=SUM+DEL
	IF(ABS(DEL).LT.ABS(SUM)*EPS)GO TO 1
11    CONTINUE
      PAUSE 'A too large, ITMAX too small'
1     GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
      RETURN
      END
C
      SUBROUTINE GCF(GAMMCF,A,X,GLN)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (ITMAX=100,EPS=3.E-7)
      CALL GAMMLN(A,GLN)
      GOLD=0.
      A0=1.
      A1=X
      B0=0.
      B1=1.
      FAC=1.
      DO 11 N=1,ITMAX
	AN=FLOAT(N)
	ANA=AN-A
	A0=(A1+A0*ANA)*FAC
	B0=(B1+B0*ANA)*FAC
	ANF=AN*FAC
	A1=X*A0+ANF*A1
	B1=X*B0+ANF*B1
	IF(A1.NE.0.)THEN
	  FAC=1./A1
	  G=B1*FAC
	  IF(ABS((G-GOLD)/G).LT.EPS)GO TO 1
	  GOLD=G
	ENDIF
11    CONTINUE
      PAUSE 'A too large, ITMAX too small'
1     GAMMCF=EXP(-X+A*ALOG(X)-GLN)*G
      RETURN
      END
C
      SUBROUTINE GAMMLN(CC,XX)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER
      DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0,
     *    -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
      DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/
      X=CC-ONE
      TMP=X+FPF
      TMP=(X+HALF)*LOG(TMP)-TMP
      SER=ONE
      DO 11 J=1,6
	X=X+ONE
	SER=SER+COF(J)/X
11    CONTINUE
      XX=TMP+LOG(STP*SER)
      RETURN
      END
C
C     SUBROUTINE TO COMPUTE DW TEST STATISTICS
C
      SUBROUTINE DWAT(E,N,DW)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NN=177,NN2=175)
      DIMENSION E(N)
      SUMTOP=0.
      DO I=12,NN2
        SUMTOP=SUMTOP+(E(I)-E(I-1))*(E(I)-E(I-1))
      END DO
      SUMBOT=0.
      DO I=11,NN2
        SUMBOT=SUMBOT+E(I)*E(I)
      END DO
      DW=SUMTOP/SUMBOT
      RETURN
      END
C
C     SUBROUTINE TO COMPUTE D-W AND 1'st-Order AR 
C
      SUBROUTINE DURBIN(YA,E,N,R2,RHO,DW)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NN=177,NN2=175)
      DIMENSION YA(N),E(N),EL(NN)
      CALL STAT(YA,N,YABAR,YASIG)
      BOT=0.
      TOP=0.
      DO 10 I=11,NN2
        TOP=TOP+E(I)*E(I)
        BOT=BOT+(YA(I)-YABAR)*(YA(I)-YABAR)
10    CONTINUE
      R2=1.-(TOP/BOT)
      DO 20 I=12,NN2
        EL(I)=E(I-1)
20    CONTINUE
      RXX=0.0
      RXY=0.0
      RSQ1=0.0
      DO 30 I=12,NN2
	RXY=RXY+EL(I)*E(I)
	RXX=RXX+EL(I)*EL(I)
        RSQ1=RSQ1+(E(I)-EL(I))*(E(I)-EL(I))
30    CONTINUE
      RSQ2=0.0
      DO 40 I=11,NN2
        RSQ2=RSQ2+E(I)*E(I)
40    CONTINUE
      RHO=RXY/RXX
      DW=RSQ1/RSQ2
      RETURN
      END
C
C     SUBROUTINE TO DO A SIMPLE REGRESSION OF ONE (N-1)-DIMENSION 
C     VARIABLE ON ANOTHER
C
      SUBROUTINE OLS(Y,X,N,ALPHA,BETA,SEA,SEB,R2,DW,RHO)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NN=177,NN2=175)
      DIMENSION Y(N),X(N),UHAT(NN),UHATL(NN)
      SUM1=0.0
      SUM2=0.0
      DO 10 I=11,NN2
	SUM1=SUM1+Y(I)
	SUM2=SUM2+X(I)
10    CONTINUE
      YBAR=SUM1/165.
      XBAR=SUM2/165.
      SXX=0.0
      SXY=0.0
      SYY=0.0
      DO 20 I=11,NN2
	SXX=SXX+(X(I)-XBAR)*(X(I)-XBAR)
	SXY=SXY+(X(I)-XBAR)*(Y(I)-YBAR)
	SYY=SYY+(Y(I)-YBAR)*(Y(I)-YBAR)
20    CONTINUE
      BETA=SXY/SXX
      ALPHA=YBAR-BETA*XBAR
      R2=BETA*SXY/SYY
      S2=(SYY*(1.-R2))/(165.-2.)
      SEA=SQRT(S2*((1./(165.))+((XBAR*XBAR)/SXX)))
      SEB=SQRT(S2/SXX)
      DO 30 I=11,NN2
	UHAT(I)=Y(I)-ALPHA-BETA*X(I)
	IF (I .GE. 12) THEN
	  UHATL(I)=UHAT(I-1)
	ENDIF
30    CONTINUE
      RXX=0.0
      RXY=0.0
      TOP=0.0
      DO 40 I=12,NN2
	RXY=RXY+UHATL(I)*UHAT(I)
	RXX=RXX+UHATL(I)*UHATL(I)
	TOP=TOP+(UHAT(I)-UHATL(I))*(UHAT(I)-UHATL(I))
40    CONTINUE
      BOT=0.0
      DO 50 I=11,NN2
	BOT=BOT+(UHAT(I)*UHAT(I))
50    CONTINUE
      RHO=RXY/RXX
      DW=TOP/BOT
      RETURN
      END
C
C     SUBROUTINE TO DO MULTIPLE REGRESSION
C
      SUBROUTINE OLSM(Y,X,M,K,BETA,SBETA,TRAT,COV,R2,R2B,RSS,SIG,U)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (L=1,KD=100,MD=300)
      DIMENSION Y(M),X(M,K),BETA(K),SBETA(K),TRAT(K),COV(K,K),U(M),
     1  XT(KD,MD),STOR1(KD,KD),STOR2(KD,L),STOR3(KD,L),STOR4(MD,L),
     2  YHAT(MD),IDX(KD)
      CALL TRANS2(X,XT,M,K)
      CALL MULT2(XT,X,STOR1,K,M,K)
      CALL INV(STOR1,K,K,IDX,COV)
      CALL MULT2(XT,Y,STOR2,K,M,L)
      CALL MULT2(COV,STOR2,STOR3,K,K,L)
      DO I=1,K
	BETA(I)=STOR3(I,1)
      END DO
      CALL MULT2(X,BETA,STOR4,M,K,L)
      SIG=0.0
      SUM1=0.0
      DO I=1,M
	U(I)=Y(I)-STOR4(I,1)
	YHAT(I)=STOR4(I,1)
	SIG=SIG+U(I)*U(I)
	SUM1=SUM1+U(I)
      END DO
      RSS=SIG
      SIG=SIG/(M-K)
      IF (DABS(SUM1) .GE. 1.0D-5) THEN
	WRITE(*,100) 
100     FORMAT(' NO WAY! OLS ERRORS DONT SUM TO ZERO!')
	PAUSE
      END IF
      DO I=1,K
	SBETA(I)=DSQRT(COV(I,I)*SIG)
	TRAT(I)=BETA(I)/SBETA(I)
      END DO
      CALL STAT(Y,M,YBAR,YSIG)
      SYY=0.0
      SXY=0.0
      SXX=0.0
      DO I=1,M
	SYY=SYY+(Y(I)-YBAR)*(Y(I)-YBAR)
	SXY=SXY+(YHAT(I)-YBAR)*(Y(I)-YBAR)
	SXX=SXX+(YHAT(I)-YBAR)*(YHAT(I)-YBAR)
      END DO
      R2=(SXY*SXY)/(SYY*SXX)
      R2B=1.-((M-1.)/(M-K))*(1-R2)
      RETURN
      END
