C....*...1.........2.........3.........4.........5.........6.........7.*
C     GMM     6/20/87, 3/18/02
C
C     PURPOSE
C     SUBROUTINE FOR GENERALIZED METHOD OF MOMENTS.  THE SYSTEM OF
C     EQUATIONS IS PRESUMED TO BE OF THE FORM E=Q(Y,X,THETA) WHERE E IS
C     M-VECTOR, Y A KY-VECTOR, X A KX-VECTOR, AND THETA IS AN LT-VECTOR
C     OF PARAMETERS TO BE ESTIMATED SUBJECT TO THETA=G(RHO) WHERE RHO
C     IS AN LR-VECTOR.
C
C     USAGE
C     CALL GMM(R,EPS,LIMIT)
C
C     USER SUPPLIED SUBROUTINES REQUIRED:
C     0.  PARMS
C         USAGE:
C         CALL PARMS(N,M,KY,KX,KZ,LT,LR,RHO,ITER1,ITER2,TOL)
C         N     - NUMBER OF OBSERVATIONS.  RETURNED BY PARMS.
C                 INTEGER*4
C         M     - NUMBER OF EQUATIONS IN THE SYSTEM.  RETURNED BY PARMS.
C                 INTEGER*4
C         KY    - NUMBER OF ENDOGENOUS VARIABLES IN THE SYSTEM.
C                 RETURNED BY PARMS.
C                 INTEGER*4
C         KX    - NUMBER OF EXOGENOUS VARIABLES IN THE SYSTEM.  RETURNED
C                 BY PARMS.
C                 INTEGER*4
C         KZ    - NUMBER OF INSTRUMENTAL VARIABLES.  RETURNED BY PARMS.
C                 INTEGER*4
C         LT    - LENGTH OF THETA.  RETURNED BY PARMS.
C                 INTEGER*4
C         LR    - LENGTH OF RHO.  RETURNED BY PARMS.
C                 INTEGER*4
C         RHO   - AN LR-VECTOR CONTAINING A START VALUE FOR RHO; THETA=
C                 G(RHO)  WHERE THETA IS THE PARAMETER VECTOR.  G IS A
C                 FUNCTION DEFINING THE CONSTRAINTS ON THE SYSTEM.
C                 RETURNED BY PARMS.
C                 REAL*8
C         ITER1 - A LIMIT ON THE NUMBER OF ITERATIONS FOR COMPUTING RHO
C                 ITER1=50 IS OFTEN REASONABLE.  RETURNED BY PARMS.
C                 INTEGER*4
C         ITER2 - THE NUMBER OF ITERATES ON THE ESTIMATE OF THE VARIANCE
C                 MATRIX OF THE MOMENT EQUATIONS.  ITER2=1 GIVES THE
C                 ONE-STEP ESTIMATE FROM THE DEFAULT START.  A LARGER
C                 VALUE GIVES A MULTI-STEP ESTIMATE.  THE DEFAULT START
C                 IS THE KRONECKER PRODUCT OF THE IDENTITY OF ORDER M
C                 WITH THE KZ BY KZ MATRIX OF SUM OF SQUARES AND CROSS
C                 PRODUCTS OF THE INSTRUMENTAL VARIABLES.  TO SUBSTITUTE
C                 A USER SUPPLIED START MATRIX, SET ITER2 TO A NEGATIVE
C                 VALUE, ADD THE FOLLOWING STATEMENT TO PARMS:
C                 COMMON /ZFIXS/ S
C                 AND SUPPLY THE MATRIX S, AN M*KZ BY M*KZ MATRIX
C                 CONTAINING THE INVERSE OF THE DESIRED VARIANCE MATRIX
C                 STORED COLUMNWISE (STORAGE MODE 0).  S IS REAL*8.
C                 IF ONE NEEDS AN INVERSION ROUTINE, THE USAGE:
C                 CALL DSWEEP(S,M*KZ,1.D-13,IER)
C                 WILL INVERT S IN PLACE.
C                 ITER2 IS RETURNED BY PARMS.
C                 INTEGER*4
C          TOL  - TOLERANCE TO DETERMINE IF CONVERGENCE HAS OBTAINED.
C                 A CONSERVATIVE VALUE IS 1.D-13. RETURNED BY PARMS.
C                 REAL*8
C     1.  MODEL
C         USAGE:
C         CALL MODEL(YT,XT,THETA,ET,DERQT,IN)
C         ARGUMENTS:
C         YT     - A KY-VECTOR OF ENDOGENOUS VARIABLES.  SUPPLIED BY THE
C                  CALLING PROGRAM.
C                  REAL*8
C         XT     - A KX-VECTOR OF EXOGENOUS VARIABLES.  SUPPLIED BY THE
C                  CALLING PROGRAM.
C                  REAL*8
C         THETA  - AN LT-VECTOR OF PARAMETER VALUES.  SUPPLIED BY THE
C                  CALLING PROGRAM.
C                  REAL*8
C         ET     - THE COMPUTED M-VECTOR OF RESIDUALS CORRESPONDING TO
C                  XT, YT, AND THETA.  RETURNED BY MODEL.
C                  REAL*8
C         DERET  - JACOBIAN OF ET WITH RESPECT TO THETA.  A MATRIX OF
C                  OF ORDER M BY LT STORED COLUMNWISE (STORAGE MODE 0).
C                  RETURNED BY MODEL.
C                  REAL*8
C         IN     - SWITCH.  IF IN=1 IT IS NOT NECESSARY TO COMPUTE
C                  DERQT.  SUPPLIED BY CALLING PROGRAM.
C                  INTEGER*4
C     2.  DATA
C         USAGE:
C         CALL DATA(IT,YT,XT,ZT)
C         ARGUMENTS:
C         IT - OBSERVATION NUMBER; IT=1,2,...,N.  SUPPLIED BY THE
C              CALLING PROGRAM.  DATA IS CALLED WITH IT OCCURRING IN
C              RANDOM, NOT SERIAL, ORDER.
C              INTEGER*4
C         YT - A KY-VECTOR OF ENDOGENOUS VARIABLES CORRESPONDING TO IT.
C              RETURNED BY DATA.
C              REAL*8
C         XT - A KX-VECTOR OF EXOGENOUS VARIABLES CORRESPONDING TO IT.
C              RETURNED BY DATA.
C              REAL*8
C         ZT - A KZ-VECTOR OF INSTRUMENTAL VARIABLES CORRESPONDING TO
C              IT.  RETURNED BY DATA.
C              REAL*8
C     3.  CONSTR
C         USAGE:
C         CALL CONSTR(RHO,THETA,DERG,IN)
C         ARGUMENTS:
C         RHO   - AN LR-VECTOR OF PARAMETERS.  SUPPLIED BY THE CALLING
C                 PROGRAM.
C                 REAL*8
C         THETA - A VECTOR OF LENGTH LT CONTAINING THE LEFT HAND SIDE
C                 OF THE CONSTRAINT EQUATIONS THETA=G(RHO).  RETURNED
C                 BY CONSTR.
C                 REAL*8
C         DERG -  A MATRIX OF ORDER LT BY LR CONTAINING THE JACOBIAN
C                 OF G(RHO).  RETURNED BY CONSTR.  STORED COLUMNWISE
C                 (STORAGE MODE 0).
C         IN    - SWITCH.  IF IN=1 IT IS NOT NECESSARY TO COMPUTE
C                 DERG.  SUPPLIED BY CALLING PROGRAM.
C                 INTEGER*4
C
C     LIBRARY SUBROUTINES CALLED
C     DSWEEP, DGMPRD, DGMAPB, DGMABC, DCOND2, SPACER, HEADER, AFC,
C     DGMCPY, DGMPNT, DGMABP
C
C     ARGUMENTS
C     R     - VECTOR OF WORKSPACE OF LENGTH LIMIT.
C             REAL*8
C             ON RETURN R CONTAINS:
C             -----------------------------------------------------
C             VARIABLE  LENGTH     COMMENT
C             -----------------------------------------------------
C             RHO1      LR         ESTIMATE OF RHO
C             S         M*KZ*M*KZ  INVERSE OF MOMENT EQS. VARIANCE
C             C         LR*LR      ESTIMATED VARIANCE OF RHO1
C             RHO2      LR         PARTIAL GAUSS-NEWTON STEP FROM RHO1
C             OBJ       1          VALUE OF THE GMM OBJECTIVE FUNCTION
C             THETA     LT         ESTIMATE OF THETA
C             GCG       LT*LT      ESTIMATED VARIANCE OF THETA
C             -----------------------------------------------------
C     EPS   - VALUE USED AS A RELATIVE TOLERANCE WHEN TESTING FOR
C             DEGENERATE RANK IN MATRIX INVERSIONS.  A REASONABLE VALUE
C             IS 1.D-13.  SUPPLIED BY THE CALLING PROGRAM.
C             REAL*8
C     LIMIT - LENGTH OF THE WORK VECTOR R.  BE SURE THAT R IS DIMENSION-
C             ED AT LEAST AS LARGE AS LIMIT IN THE CALLING PROGRAM.
C             SUPPLIED BY THE CALLING PROGRAM.  IF THE WORK SPACE IS
C             INADEQUATE, THE ROUTINE WILL TERMINATE WITH A MESSAGE
C             STATING THE REQUIRED LENGTH.
C             INTEGER*4
C
C     COMMENTS
C
C     THE USAGE:
C     COMMON /ZLNSIZ/ LNSIZE
C     LNSIZE=80
C     WILL SET THE LINESIZE TO 80.  LINESIZES BETWEEN 80 AND 133 ARE
C     PERMITTED.  THE DEFAULT LINESIZE IS 133, 132 PLUS CARRIAGE CONTROL
C     CHARACTER.
C
C     THE USAGE:
C     COMMON /ZPRINT/ ISHORT
C     ISHORT=IPNT
C     CONTROLS PRINTING.
C     IF IPNT=3 THERE IS NO OUTPUT.
C     IF IPNT=2 PRINTING OF RHO HAT AND ITS VARIANCE IS SUPPRESSED.
C     IF IPNT=1 PRINTING OF THETA HAT AND ITS VARIANCE IS SUPPRESSED.
C     IF IPNT=0 FULL OUTPUT, DEFAULT.
C
C     THE USAGE:
C     COMMON /ZLAG/ LAG,IWTS
C     LAG=NLAG
C     IWTS=1
C     ALLOWS THE USER TO SET THE NUMBER OF LAGS USED TO COMPUTE THE
C     ESTIMATE OF THE VARIANCE OF THE MOMENT EQUATIONS TO NLAGS AND
C     SUPPRESS THE USES OF PARZEN WEIGHTS.  SET IWTS=0 TO USE PARZEN
C     WEIGHTS WITH USER SPECIFIED LAG LENGTH.
C
C     WARNING
C     WITH THE MICROSOFT COMPILER, THESE STATEMENTS MUST BE PRESENT:
C     COMMON /ZPRINT/ ISHORT
C     COMMON /ZLNSIZ/ LNSIZE
C     COMMON /ZLAG/   LAG,IWTS
C     ISHORT=0
C     LNSIZE=133
C     LAG=-1
C     IWTS=1
C     THESE ARE DEFAULTS, SET AS ABOVE TO OBTAIN OPTIONS.
C
C     REFERENCE
C     A. RONALD GALLANT (1987).  NONLINEAR STATISTICAL MODELS.  NEW
C     YORK: WILEY.  CHAPTER 6, SECTION 3.
C
C     PROGRAMMER
C     A. RONALD GALLANT
C     DEPARTMENT OF STATISTICS
C     NORTH CAROLINA STATE UNIVERSITY
C     RALEIGH, NORTH CAROLINA 27695-8203
C     (919) 737-2531
C
C     EXAMPLE
C     REAL*8 R(4096)
C     CALL GMM(R,1.D-10,4096)
C     STOP
C     END
C     SUBROUTINE PARMS(N,M,KY,KX,KZ,LT,LR,RHO,ITER1,ITER2,TOL)
C     IMPLICIT REAL*8 (A-H,O-Z)
C     REAL*8 RHO(4)
C     COMMON /ZPRINT/ ISHORT
C     COMMON /ZLNSIZ/ LNSIZE
C     COMMON /ZLAG/   LAG,IWTS
C     ISHORT=0
C     LNSIZE=80
C     LAG=-1
C     IWTS=1
C     N=9
C     M=2
C     KY=2
C     KX=2
C     KZ=3
C     LT=4
C     LR=2
C     DO 10 I=1,LR
C10   RHO(I)=1.001D0
C     ITER1=50
C     ITER2=1
C     TOL=1.D-10
C     RETURN
C     END
C     SUBROUTINE DATA(IT,Y,X,Z)
C     IMPLICIT REAL*8 (A-H,O-Z)
C     REAL*8 Y(2),X(2),Z(3)
C     REAL*4 XX(2,9)
C    &/1.,1.,2.,4.,3.,9.,4.,16.,5.,25.,8.,3.,3.,.2,6.,7.,3.,9./
C     REAL*4 EE(2,9)
C    &/.10480,.22368,.24130,.42167,.37570,.65833,.06473,.28676,.06773,
C    & .77921,.99562,.96301,.89579,.85475,.07583,.28647,.76390,.38564/
C     X(1)=XX(1,IT)
C     X(2)=XX(2,IT)
C     Z(1)=X(1)
C     Z(2)=X(2)
C     Z(3)=DSQRT(X(1)**2+X(2)**2)
C     Y(1)=DLOG(X(1)+X(2))+EE(1,IT)
C     Y(2)=DEXP(X(1))+EE(2,IT)
C     RETURN
C     END
C     SUBROUTINE MODEL(Y,X,T,E,Q,IN)
C     IMPLICIT REAL*8 (A-H,O-Z)
C     REAL*8 Y(2),X(2),T(4),E(2),Q(2,4)
C     E(1)=Y(1)-DLOG(T(1)*X(1)+T(2)*X(2))
C     E(2)=Y(2)-T(3)*DEXP(T(4)*X(1))
C     IF(IN.EQ.1) RETURN
C     Q(1,1)=-X(1)/(T(1)*X(1)+T(2)*X(2))
C     Q(1,2)=-X(2)/(T(1)*X(1)+T(2)*X(2))
C     Q(1,3)=0.D0
C     Q(1,4)=0.D0
C     Q(2,1)=0.D0
C     Q(2,2)=0.D0
C     Q(2,3)=-DEXP(T(4)*X(1))
C     Q(2,4)=-T(3)*X(1)*DEXP(T(4)*X(1))
C     RETURN
C     END
C     SUBROUTINE CONSTR(RHO,G,GG,IN)
C     IMPLICIT REAL*8 (A-H,O-Z)
C     REAL*8 RHO(1),G(1),GG(1)
C     REAL*8 DG(4,2)/1.D0,0.D0,1.D0,0.D0,0.D0,1.D0,0.D0,1.D0/
C     LR=2
C     LT=4
C     CALL DGMPRD(DG,RHO,G,LT,LR,1)
C     IF(IN.EQ.1) RETURN
C     CALL DGMCPY(DG,GG,LT,LR)
C     RETURN
C     END
C     ANSWER
C     RHO=(1.00941,0.99888)
C     VAR(RHO)=| 1.0011D-05 -1.2632D-06|
C              |-1.2632D-06  1.5941D-07|
C
C
      SUBROUTINE GMM(R,EPS,LIMIT)
      IMPLICIT INTEGER*4 (A-Z)
      save
      REAL*8 R(1),EPS,TOL
      CALL PARMS(N,M,KY,KX,KZ,LT,LR,R,ITER1,ITER2,TOL)
      RHO1=1
      S   =RHO1+LR
      C   =S   +M*KZ*M*KZ
      RHO2=C   +LR*LR
      W   =RHO2+LR
      CALL Z1GMM(N,M,KY,KX,KZ,LT,R(S),LR,R(RHO1),R(RHO2),R(C),R(W),
     &ITER1,ITER2,LIMIT,TOL,EPS)
      RETURN
      END
      SUBROUTINE Z1GMM(N,M,KY,KX,KZ,LT,S,LR,RHO1,RHO2,C,W,
     &ITER1,ITER2,LIMIT,TOL,EPS)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      INTEGER*4 ALPHA,BETA
      REAL*8 FIXS(1)
      REAL*8 S(1),RHO1(LR),RHO2(LR),C(LR,LR),W(1)
      LOGICAL*4    DONE
      CHARACTER*1  LFC(8),COMMA,CPAREN,BLANK,DIGIT(10)
      CHARACTER*4  XFC(2),BLANKS
      CHARACTER*8  RFC,SAVEFC
      CHARACTER*8  FMT2(5),FMT3(5),FMT4(36)
      CHARACTER*8  FMT5(12),FMT6(12),FMT7(12),FMT8(12)
      CHARACTER*40 CFMT2,CFMT3
      CHARACTER*96 CFMT4A,CFMT4B,CFMT4C,CFMT5,CFMT6,CFMT7,CFMT8
      EQUIVALENCE (RFC,LFC(1),XFC(1))
      EQUIVALENCE (FMT2,CFMT2),(FMT3,CFMT3),(FMT5,CFMT5),(FMT6,CFMT6)
      EQUIVALENCE (FMT4(1),CFMT4A),(FMT4(13),CFMT4B),(FMT4(25),CFMT4C)
      EQUIVALENCE (FMT5,CFMT5),(FMT6,CFMT6),(FMT7,CFMT7),(FMT8,CFMT8)
      COMMON /ZFIXS/  FIXS
      COMMON /ZLNSIZ/ LNSIZE
      COMMON /ZPRINT/ ISHORT
      COMMON /ZLAG/   NLAGS,IWTS
      DATA COMMA/','/,CPAREN/')'/,BLANK/' '/,BLANKS/'    '/
      DATA DIGIT/'0','1','2','3','4','5','6','7','8','9'/
      DATA FMT2 /'(1H ,24X',',I4,4X, ','        ',' I13,8X,','        '/
      DATA FMT3 /'(1H ,24X',',4X,4X, ',' 25X,   ',' I13,8X,','        '/
      DATA FMT4 /
     &'(1H ,24X',',1X,4H  ','  ,8X,20','H       ','        ','     ,  ',
     &'8X,9HPAR','AMETER,8','X,21H   ','        ','        ','     )  ',
     &'(1H ,24X',',1X,4HST','EP,8X,20','H OBJECT','IVE FUNC','TION ,  ',
     &'8X,9H  N','UMBER ,8','X,21HINT','ERMEDIAT','E ESTIMA','TE   )  ',
     &'(1H ,24X',',1X,4H--','--,8X,20','H-------','--------','-----,  ',
     &'8X,9H---','------,8','X,21H---','--------','--------','--   )  '/
      DATA FMT5/
     &'(1H ,24X',',1X,51H-','---   IN','VERSION ','ERROR EN','COUNTERE',
     &'D, G-INV','ERSE USE','D  )    ','        ','        ','        '/
      DATA FMT6/
     &'(1H ,24X',',1X,51H-','---   UN','ABLE TO ','IMPROVE ','ON THIS ',
     &'ESTIMATE','        ','   )    ','        ','        ','        '/
      DATA FMT7/
     &'(1H ,24X',',1X,51H-','---   CH','ANGE WAS',' LESS TH','AN PRESE',
     &'T TOLERA','NCE     ','   )    ','        ','        ','        '/
      DATA FMT8/
     &'(1H ,24X',',1X,51H-','---   IT','ERATION ','LIMIT EX','CEEDED  ',
     &'        ','        ','   )    ','        ','        ','        '/
      LNSIZ=LNSIZE
      IF((LNSIZE.LT.72).OR.(LNSIZE.GT.133)) LNSIZ=133
      IOBJ  =1
      ITHETA=IOBJ  +1
      IGCG  =ITHETA+LT
      IG    =IGCG  +LT*LT
      IW    =IG    +LT*LR
      LW    =IW    +LT*LR
      MKZ=M*KZ
      IER3=0
      IOUT=3
      IF (NLAGS.LT.0) THEN
      ISWVAR=0
      ELSE
      LAG=NLAGS
      ISWVAR=3
      IF(IWTS.EQ.0) ISWVAR=2
      ENDIF
      SAVEFC=FMT2(1)
      IPAD=(LNSIZ-80)/2
      IF(IPAD.LT.0) IPAD=0
      IF(IPAD.LT.24) THEN
      IPAD10=IPAD/10
      IPAD1=IPAD-10*IPAD10
      RFC=SAVEFC
      LFC(6)=DIGIT(IPAD10+1)
      LFC(7)=DIGIT(IPAD1+1)
      IF(IPAD.EQ.0) XFC(2)=BLANKS
      SAVEFC  =RFC
      FMT2(1) =RFC
      FMT3(1) =RFC
      FMT4(1) =RFC
      FMT4(13)=RFC
      FMT4(25)=RFC
      FMT5(1) =RFC
      FMT6(1) =RFC
      FMT7(1) =RFC
      FMT8(1) =RFC
      ENDIF
      CALL Z2GMM(N,M,KY,KX,KZ,LT,LR,LIMIT,SAVEFC)
      IF(ITER2.GE.0) THEN
      ISTART=0
      CALL Z3GMM(N,M,KY,KX,KZ,LT,LR,S,RHO1,LAG,1,EPS,IER3,W(LW))
      ELSE
      ISTART=1
      CALL DGMCPY(FIXS,S,MKZ,MKZ)
      ENDIF
      DO 50 ICTR2=ISTART,IABS(ITER2)
      IF(ISHORT.LE.2) THEN
      CALL SPACER(0)
      CALL SPACER(4)
      CALL HEADER('// THE INVERSE OF THE VARIANCE-COVARIANCE MATRIX/
     &IS FIXED AS SHOWN HERE FOR THE COMPUTATIONS WHICH FOLLOW///_')
      CALL DGMPNT(S,MKZ,MKZ)
      CALL SPACER(0)
      CALL SPACER(2)
      CALL HEADER('//MODIFIED GAUSS-NEWTON ITERATIONS///_')
      CALL SPACER(2)
      CALL SPACER(2)
      WRITE(IOUT,CFMT4A)
      WRITE(IOUT,CFMT4B)
      WRITE(IOUT,CFMT4C)
      ENDIF
      CALL DGMCPY(RHO1,RHO2,LR,1)
      OBJLAG=0.D0
      DO 30 ICTR1=0,ITER1
      CALL DGMCPY(RHO2,RHO1,LR,1)
      OBJLAG=OBJ
      CALL Z4GMM(N,M,KY,KX,KZ,LT,LR,S,RHO1,C,OBJ,RHO2,EPS,IER,W(LW))
      IER=IER+1000*IER3
      W(IOBJ)=OBJ
      DONE=.FALSE.
      IF(ICTR1.NE.0) CALL Z5GMM(LR,RHO1,RHO2,OBJ,OBJLAG,TOL,DONE)
      IF(ISHORT.LE.2) THEN
      CALL AFC(OBJ,RFC,25,16)
      LFC(8)=COMMA
      FMT2(3)=RFC
      CALL AFC(RHO1(1),RFC,25,16)
      LFC(8)=CPAREN
      FMT2(5)=RFC
      CALL SPACER(1)
      I=1
      WRITE(IOUT,CFMT2) ICTR1,OBJ,I,RHO1(I)
      IF(LR.EQ.1) GO TO 16
      DO 15 I=2,LR
      CALL AFC(RHO1(I),RFC,25,16)
      LFC(8)=CPAREN
      FMT3(5)=RFC
15    WRITE(IOUT,CFMT3) I,RHO1(I)
16    CONTINUE
      IF(IER.GT.9) WRITE(IOUT,CFMT5)
      IF(MOD(IER,10).NE.0) WRITE(IOUT,CFMT6)
      IF(DONE) WRITE(IOUT,CFMT7)
      IF(ICTR1.EQ.ITER1) WRITE(IOUT,CFMT8)
      ENDIF
      IF(DONE) GO TO 31
      IF(MOD(IER,10).NE.0) GO TO 31
30    CONTINUE
31    CONTINUE
      IF(ISHORT.LE.2) THEN
      CALL SPACER(2)
      CALL HEADER('//CHECK THESE ITERATIONS TO SEE IF CONVERGENCE HAS/BE
     &EN ACHIEVED BEFORE USING THE FOLLOWING RESULTS///_')
      ENDIF
      if(ictr2.lt.iabs(iter2)) then
      CALL Z3GMM(N,M,KY,KX,KZ,LT,LR,S,RHO1,LAG,ISWVAR,EPS,IER3,W(LW))
      endif
50    CONTINUE
      IF(ISHORT.LE.2) THEN
      CALL SPACER(0)
      CALL SPACER(2)
      CALL HEADER('//ESTIMATED VALUE OF RHO/RHO-HAT///_')
      CALL DGMPNT(RHO1,LR,1)
      CALL SPACER(2)
      CALL HEADER('//ESTIMATED VARIANCE-COVARIANCE MATRIX OF RHO-HAT///_
     &')
      CALL DGMPNT(C,LR,LR)
      ENDIF
      IN=0
      CALL CONSTR(RHO1,W(ITHETA),W(IG),IN)
      CALL DGMABP(C,W(IG),W(IW),LR,LR,LT)
      CALL DGMPRD(W(IG),W(IW),W(IGCG),LT,LR,LT)
      IF(ISHORT.LE.1) THEN
      CALL SPACER(0)
      CALL SPACER(2)
      CALL HEADER('//ESTIMATED VALUE OF THETA/THETA-HAT///_')
      CALL DGMPNT(W(ITHETA),LT,1)
      CALL SPACER(2)
      CALL HEADER('//ESTIMATED VARIANCE-COVARIANCE MATRIX OF THETA-HAT//
     &/_')
      CALL DGMPNT(W(IGCG),LT,LT)
      ENDIF
      RETURN
      END
      SUBROUTINE Z2GMM(N,M,KY,KX,KZ,LT,LR,LIMIT,SAVEFC)
      IMPLICIT INTEGER*4 (A-Z)
      save
      CHARACTER*8 SAVEFC,FMT1(12)
      CHARACTER*96 CFMT1
      EQUIVALENCE (FMT1,CFMT1)
      COMMON /ZPRINT/ ISHORT
      DATA FMT1/
     &'(1H ,24X',',7X,32HD','IMENSION',' R AT LE','AST AS L','ARGE AS,',
     &'I6,29H A','ND SET L','IMIT TO ','CORRESPO','ND     )','        '/
      FMT1(1)=SAVEFC
      IOUT=3
C     WORKSPACE USED BY GMM
      RHO1  =1
      S     =RHO1  +LR
      C     =S     +M*KZ*M*KZ
      RHO2  =C     +LR*LR
      W     =RHO2  +LR
      IOBJ  =W
      ITHETA=IOBJ  +1
      IGCG  =ITHETA+LT
      IG    =IGCG  +LT*LT
      IW    =IG    +LT*LR
      LW1   =IW    +LT*LR
C     INCREMENTAL WORKSPACE USED BY Z3GMM
      ITHETA=LW1
      IG  =ITHETA+LT
      IYT =IG    +LT*LR
      IXT =IYT   +MAX0(KY,M)
      IZT =IXT   +KX
      IZTL=IZT   +KZ
      IE  =IZTL  +KZ
      IEL =IE    +M
      IQ  =IEL   +M
      LW2 =IQ    +M*LT
C     INCREMENTAL WORKSPACE USED BY Z4GMM
      ISTEP =LW1
      ITHETA=ISTEP +LR
      IQVQ  =ITHETA+LT
      IQVE  =IQVQ  +LT*LT
      IW    =IQVE  +LT
      IG    =IW    +LT*LT
      LW3   =IG    +LT*LR
C     INCREMENTAL WORKSPACE USED BY Z6GMM OF Z4GMM
C     NOTE THAT IT EXCEEDS THAT USED BY Z7GMM
      IYT   =LW1
      IXT   =IYT   +MAX0(KY,M)
      IZT   =IXT   +KX
      IE    =IZT   +KZ
      IQ    =IE    +M
      IZE   =IQ    +M*LT
      IZQ   =IZE   +KZ*M
      LW4   =IZQ   +KZ*M*LT
C     TOTAL WORKSPACE
      LW=MAX0(LW1,LW2,LW3,LW4)
      IF(LW.GT.LIMIT) GO TO 50
      IF(ISHORT.GE.3) RETURN
      CALL SPACER(0)
      CALL SPACER(4)
      CALL HEADER('//WORKSPACE REQUIREMENTS///_')
      CALL SPACER(4)
      WRITE(IOUT,CFMT1) LW
C     CALL SPACER(4)
C     CALL HEADER('//GMM VERSION 6-20-87.  FOR DOCUMENTATION PRINT//
C    &NB4139.GMM.FBA  ///_')
C     CALL HEADER('//PLEASE REPORT ANY PROBLEMS WITH THIS PROGRAM TO://
C    &A. RONALD GALLANT/DEPARTMENT OF STATISTICS/NORTH CAROLINA STATE UN
C     &IVERSITY/RALEIGH, NORTH CAROLINA 27695-8203/(919) 737-2531///_')
      RETURN
50    CALL SPACER(0)
      CALL SPACER(4)
      CALL HEADER('//INSUFFICIENT WORKSPACE///_')
      CALL SPACER(4)
      WRITE(IOUT,CFMT1) LW
      CALL SPACER(4)
      CALL HEADER('//COMPUTATIONS TERMINATED///_')
      STOP
      END
      SUBROUTINE Z5GMM(LR,RHO1,RHO2,OBJ,OBJLAG,TOL,DONE)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 RHO1(1),RHO2(1),NORM1,NORM12
      LOGICAL*4 DONE
      DONE=.TRUE.
      NORM1=0.D0
      NORM12=0.D0
      DO 10 I=1,LR
      NORM1=NORM1+RHO1(I)**2
10    NORM12=NORM12+(RHO1(I)-RHO2(I))**2
      NORM1=DSQRT(NORM1)
      NORM12=DSQRT(NORM12)
      IF(NORM12.GE.TOL*NORM1) DONE=.FALSE.
      IF(DABS(OBJ-OBJLAG).GE.TOL*OBJLAG) DONE=.FALSE.
      RETURN
      END
      SUBROUTINE Z3GMM(N,M,KY,KX,KZ,LT,LR,S,RHO,LAG,ISW,EPS,IER3,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 S(1),RHO(1),W(1),TMP(1)
      INTEGER*4 ALPHA,BETA,ROW,COL
C
C     SET UP POINTERS
C
      ITHETA=0
      IG    =ITHETA+LT
      IYT   =IG    +LT*LR
      IXT   =IYT   +MAX0(KY,M)
      IZT   =IXT   +KX
      IZTL  =IZT   +KZ
      IE    =IZTL  +KZ
      IEL   =IE    +M
      IQ    =IEL   +M
      LW    =IQ    +M*LT
C
      MKZ=M*KZ
      FN=N
      IN=1
      CALL CONSTR(RHO,W(ITHETA+1),W(IG+1),IN)
C
C     SET TO (I @ Z'Z) IF ISW.EQ.1
C
      IF (ISW.EQ.1) THEN
      DO 100 J=1,MKZ
      DO 100 I=1,MKZ
100   S(MKZ*(J-1)+I)=0.D0
      DO 140 IT=1,N
      CALL DATA(IT,W(IYT+1),W(IXT+1),W(IZT+1))
      DO 120 BETA=1,M
      DO 120 J=1,KZ
      COL=KZ*(BETA-1)+J
      ALPHA=BETA
      DO 120 I=1,KZ
      ROW=KZ*(ALPHA-1)+I
      LOC=MKZ*(COL-1)+ROW
120   S(LOC)=S(LOC)+W(IZT+I)*W(IZT+J)
140   CONTINUE
C
C     STANDARD CASE, ISW.NE.1
C
      ELSE
      DO 200 J=1,MKZ
      DO 200 I=1,MKZ
200   S(MKZ*(J-1)+I)=0.D0
      IF ((ISW.EQ.2).OR.(ISW.EQ.3)) GO TO 205
      LAG=DEXP(DLOG(FN)/5.D0)+0.5D0
205   CONTINUE
      FLAG=LAG
C     LAG ZERO CASE:
      DO 240 IT=1,N
      CALL DATA(IT,W(IYT+1),W(IXT+1),W(IZT+1))
      CALL MODEL(W(IYT+1),W(IXT+1),W(ITHETA+1),W(IE+1),W(IQ+1),IN)
      DO 220 BETA=1,M
      DO 220 J=1,KZ
      COL=KZ*(BETA-1)+J
      DO 220 ALPHA=1,M
      DO 220 I=1,KZ
      ROW=KZ*(ALPHA-1)+I
      LOC=MKZ*(COL-1)+ROW
220   S(LOC)=S(LOC)+W(IZT+I)*W(IE+ALPHA)*W(IZT+J)*W(IE+BETA)
240   CONTINUE
C     POSITIVE LAGS CASE:
      IF (LAG.LE.0) GO TO 360
      WEIGHT=1.D0
      DO 350 L=1,LAG
      IF (ISW.EQ.3) GO TO 300
      FL=L
      X=FL/FLAG
      IF (X.LE.0.5D0) WEIGHT=1.D0-6.D0*(X**2)+6.D0*(X**3)
      IF (X.GT.0.5D0) WEIGHT=2.D0*((1.D0-X)**3)
300   CONTINUE
      DO 340 IT=1+L,N
      CALL DATA(IT  ,W(IYT+1),W(IXT+1),W(IZT +1))
      CALL MODEL(W(IYT+1),W(IXT+1),W(ITHETA+1),W(IE +1),W(IQ+1),IN)
      CALL DATA(IT-L,W(IYT+1),W(IXT+1),W(IZTL+1))
      CALL MODEL(W(IYT+1),W(IXT+1),W(ITHETA+1),W(IEL+1),W(IQ+1),IN)
      DO 320 BETA=1,M
      DO 320 J=1,KZ
      COL=KZ*(BETA-1)+J
      DO 320 ALPHA=1,M
      DO 320 I=1,KZ
      ROW=KZ*(ALPHA-1)+I
      LOC=MKZ*(COL-1)+ROW
320   S(LOC)=S(LOC)
     &       +WEIGHT*(W(IZT +I)*W(IE +ALPHA)*W(IZTL+J)*W(IEL+BETA))
     &       +WEIGHT*(W(IZTL+I)*W(IEL+ALPHA)*W(IZT +J)*W(IE +BETA))
340   CONTINUE
350   CONTINUE
360   CONTINUE
      ENDIF
C
C     INVERT
C
      IF(LW.GT.MKZ) CALL DCOND2(S,MKZ,W(ITHETA+1),0)
      CALL DSWEEP(S,MKZ,EPS,IER3)
      IF(LW.GT.MKZ) CALL DCOND2(S,MKZ,W(ITHETA+1),1)
C
      RETURN
      END
      SUBROUTINE Z4GMM(N,M,KY,KX,KZ,LT,LR,S,RHO,C,OBJ,RHO2,EPS,IER,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 S(1),RHO(LR),RHO2(LR),C(LR,LR),W(1)
      INTEGER*4 ALPHA,BETA
C
C     SET UP WORKSPACE
C
      ISTEP =1
      ITHETA=ISTEP +LR
      IQVQ  =ITHETA+LT
      IQVE  =IQVQ  +LT*LT
      IW    =IQVE  +LT
      IG    =IW    +LT*LT
      LW    =IG    +LT*LR
C
C     COMPUTE A FULL GAUSS-NEWTON STEP
C
      IN=1
      CALL CONSTR(RHO,W(ITHETA),W(IG),IN)
      CALL Z6GMM(N,M,KY,KX,KZ,S,W(ITHETA),LT,W(IQVQ),W(IQVE),OBJ,W(IW))
      IN=0
      CALL CONSTR(RHO,W(ITHETA),W(IG),IN)
      CALL DGMABC(W(IG),W(IQVQ),W(IG),C,LT,LR,LT,LR)
      CALL DCOND2(C,LR,W(IW),0)
      CALL DSWEEP(C,LR,EPS,IER2)
      CALL DCOND2(C,LR,W(IW),1)
      CALL DGMAPB(W(IG),W(IQVE),RHO2,LT,LR,1)
      CALL DGMPRD(C,RHO2,W(ISTEP),LR,LR,1)
      DO 20 I=1,LR
20    W(ISTEP-1+I)=-W(ISTEP-1+I)
C
C     COMPUTE A MODIFIED GAUSS-NEWTON STEP
C
      OBJ1=OBJ
      IER1=0
      V=1.1D0
      DO 60 L=1,40
      IF(L.LE.6) V=V-.1D0
      IF(L.GT.6) V=V*.5D0
      DO 40 I=1,LR
40    RHO2(I)=RHO(I)+V*W(ISTEP-1+I)
      DO 45 I=1,LR
      IF(RHO(I).NE.RHO2(I)) GO TO 47
45    CONTINUE
      GO TO 65
47    CONTINUE
      IN=1
      CALL CONSTR(RHO2,W(ITHETA),W(IG),IN)
      CALL Z7GMM(N,M,KY,KX,KZ,S,W(ITHETA),LT,OBJ2,W(IW))
      IF(OBJ2.LT.OBJ1) GO TO 70
60    CONTINUE
65    CONTINUE
      IER1=1
70    CONTINUE
      IER=IER1+10*IER2
      RETURN
      END
      SUBROUTINE Z6GMM(N,M,KY,KX,KZ,S,THETA,LT,QVQ,QVE,OBJ,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 S(1),THETA(LT),QVQ(LT,LT),QVE(LT),W(1),TMP(1)
      INTEGER*4 ALPHA,BETA,ROW,COL
C
C     SET UP POINTERS
C
      IYT=0
      IXT=IYT+MAX0(KY,M)
      IZT=IXT+KX
      IE =IZT+KZ
      IQ =IE +M
      IEZ=IQ +M*LT
      IQZ=IEZ+KZ*M
      LW =IQZ+KZ*M*LT
C
      MKZ=M*KZ
      FN=N
C
C     COMPUTE SUM Q@Z & SUM E@Z FOR T = 1, ..., N
C
      DO 200 ALPHA=1,M
      DO 200 I=1,KZ
200   W(IEZ+KZ*(ALPHA-1)+I)=0.D0
      DO 210 COL=1,LT
      DO 210 ALPHA=1,M
      DO 210 I=1,KZ
210   W(IQZ+MKZ*(COL-1)+KZ*(ALPHA-1)+I)=0.D0
      IN=0
      DO 240 IT=1,N
      CALL DATA(IT,W(IYT+1),W(IXT+1),W(IZT+1))
      CALL MODEL(W(IYT+1),W(IXT+1),THETA,W(IE+1),W(IQ+1),IN)
      DO 220 ALPHA=1,M
      DO 220 I=1,KZ
      LOC=IEZ+KZ*(ALPHA-1)+I
220   W(LOC)=W(LOC)+W(IZT+I)*W(IE+ALPHA)
      DO 230 COL=1,LT
      DO 230 ALPHA=1,M
      DO 230 I=1,KZ
      LOC=IQZ+MKZ*(COL-1)+KZ*(ALPHA-1)+I
230   W(LOC)=W(LOC)+W(IZT+I)*W(IQ+M*(COL-1)+ALPHA)
240   CONTINUE
C
C     COMPUTE QVQ, QVE & OBJ
C
      CALL DGMABC(W(IQZ+1),S,W(IQZ+1),QVQ,MKZ,LT,MKZ,LT)
      CALL DGMABC(W(IQZ+1),S,W(IEZ+1),QVE,MKZ,LT,MKZ,1)
      CALL DGMABC(W(IEZ+1),S,W(IEZ+1),TMP,MKZ,1,MKZ,1)
      OBJ=TMP(1)
      RETURN
      END
      SUBROUTINE Z7GMM(N,M,KY,KX,KZ,S,THETA,LT,OBJ,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 S(1),THETA(LT),W(1),TMP(1)
      INTEGER*4 ALPHA,BETA
C
C     SET UP POINTERS
C
      IYT=0
      IXT=IYT+MAX0(KY,M)
      IZT=IXT+KX
      IE =IZT+KZ
      IQ =IE +M
      IEZ=IQ +M*LT
      IQZ=IEZ+KZ*M
      LW =IQZ+KZ*M*LT
C
      MKZ=M*KZ
      FN=N
C
C     COMPUTE SUM E@Z FOR T = 1, ..., N
C
      DO 200 ALPHA=1,M
      DO 200 I=1,KZ
200   W(IEZ+KZ*(ALPHA-1)+I)=0.D0
      IN=1
      DO 240 IT=1,N
      CALL DATA(IT,W(IYT+1),W(IXT+1),W(IZT+1))
      CALL MODEL(W(IYT+1),W(IXT+1),THETA,W(IE+1),W(IQ+1),IN)
      DO 220 ALPHA=1,M
      DO 220 I=1,KZ
      LOC=IEZ+KZ*(ALPHA-1)+I
220   W(LOC)=W(LOC)+W(IZT+I)*W(IE+ALPHA)
240   CONTINUE
C
C     COMPUTE OBJ
C
      CALL DGMABC(W(IEZ+1),S,W(IEZ+1),TMP,MKZ,1,MKZ,1)
      OBJ=TMP(1)
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*
C     THESE STATEMENTS HAVE NO EFFECT WITH THE MICROSOFT COMPILER/LINKER
C     AND DEFAULT OPTIONS MUST BE SET BY THE USER.  SEE WARNING ABOVE.
      BLOCK DATA
      COMMON /ZPRINT/ ISHORT
      COMMON /ZLAG/   NLAGS,IWTS
      DATA ISHORT/1/
      DATA NLAGS/-1/,IWTS/1/
      END