C....*...1.........2.........3.........4.........5.........6.........7.*
C     NLSYS    2/25/85
C
C     PURPOSE
C     SUBROUTINE FOR UNIVARIATE NONLINEAR LEAST SQUARES, MULTIVARIATE
C     NONLINEAR LEAST SQUARES (SEEMINGLY UNRELATED NONLINEAR REGRES-
C     SIONS), TWO STAGE NONLINEAR LEAST SQUARES, AND THREE-STAGE NON-
C     LINEAR LEAST SQUARES.  THE MODEL IS PRESUMED TO BE OF THE FORM
C     E=Q(Y,X,THETA) WHERE E IS M-VECTOR, Y A KY-VECTOR, X A KX-VECTOR,
C     AND THETA IS AN LT-VECTOR OF PARAMETERS TO BE ESTIMATED SUBJECT
C     TO THETA=G(RHO) WHERE RHO IS AN LR-VECTOR.
C
C     USAGE
C     CALL NLSYS(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.  M MAY BE SET TO 1
C                 TO ESTIMATE A SINGLE EQUATION BY LEAST SQUARES OR TWO-
C                 STAGE LEAST-SQUARES.  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.  SET  KZ=0  FOR
C                 UNIVARIATE OR MULTIVARIATE LEAST SQUARES COMPUTATIONS.
C                 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 EASONABLE.  RETURNED BY PARMS.
C                 INTEGER*4
C         ITER2 - THE NUMBER OF ITERATES ON THE VARIANCE MATRIX OF THE
C                 M-VECTOR OF ERRORS (ACROSS EQUATIONS VARIANCE MATRIX)
C                 ITER2=1 GIVES THE ONE-STEP SEEMINGLY UNRELATED OR
C                 THREE-STAGE LEAST-SQUARES ESTIMATE.  A LARGER VALUE
C                 GIVES A MULTI-STEP ESTIMATE.  THE START MATRIX IS THE
C                 IDENTITY WITH THIS USAGE.  TO SUBSTITUTE A USER
C                 SUPPLIED START MATRIX, SET ITER2 TO A NEGATIVE VALUE,
C                 ADD THE FOLLOWING STATEMENT TO PARMS:
C                 COMMON /ZFIXS/ S
C                 AND SUPPLY THE MATRIX S, AN M BY M MATRIX CONTAINING
C                 THE INVERSE OF THE DESIRED VARIANCE MATRIX STORED
C                 COLUMNWISE (STORAGE MODE 0).  S IS REAL*8.
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         DERQT  - 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.
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.  IF THERE ARE NONE, KZ=0, INCLUDE THE ARGUMENT AND
C              THE STATEMENT  REAL*8 ZT(1).  IF THERE ARE INSTRUMENTS,
C              KZ>0, THEN THE N BY KZ MATRIZ Z WITH ROWS ZT IS ASSUMED
C              TO HAVE ORTHONORMAL COLUMNS.  LIBRARY ROUTINE DGRAM MAY
C              BE USED TO ORTHOGANALIZE Z.
C              USAGE:
C              CALL DGRAM(Z,N,KZ,1.D-13,IR,WORK)
C              WORK IS A REAL*8 N-VECTOR.
C              ZT IS 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
C                 LEFT HAND SIDE OF THE CONSTRAINT EQUATIONS
C                 THETA=G(RHO).  RETURNED 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     SYSMGN, DSWEEP, DGMPRD, DGMAPB, DCOND2, SPACER, HEADER, AFC,
C     DGMCPY, DGMPNT, DGMABP, DGRAM
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*M     INVERSE OF ERROR VARIANCE
C             C         LR*LR   ESTIMATED VARIANCE OF RHO1
C             RHO2      LR      PARTIAL GAUSS-NEWTON STEP FROM RHO1
C             S2        M*M     INVERSE OF RESIDUALS FROM RHO1
C             OBJ       1       VALUE OF 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 DIMINSION-
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 /ZSHORT/ 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
C     REFERENCES
C     A. RONALD GALLANT.  SEEMINGLY UNRELATED NONLINEAR REGRESSIONS.
C     JOURNAL OF ECONOMETRICS, 3. (1975) PP. 35-50.
C     A. RONALD GALLANT.  THREE-STAGE LEAST-SQUARES ESTIMATION FOR A
C     SYSTEM OF SIMULTANEOUS, NONLINEAR, IMPLICIT EQUATIONS.  JOURNAL OF
C     ECONOMETRICS, 5.  (1977)  PP. 71-88.
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     COMMON /ZLNSIZ/ LNSIZE
C     LNSIZE=80
C     CALL NLSYS(R,1.D-13,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     N=5
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,5)/1.,1.,2.,4.,3.,9.,4.,16.,5.,25./
C     REAL*4 ZZ(3,5)/1.,-2.,2.,1.,-1.,-1.,1.,0.,-2.,1.,1.,-1.,1.,2.,2./
C     REAL*4 EE(2,5)/.10480,.22368,.24130,.42167,.37570,
C    & .77921,.99562,.96301,.89579,.85475/
C     X(1)=XX(1,IT)
C     X(2)=XX(2,IT)
C     Z(1)=ZZ(1,IT)/DSQRT(5.D0)
C     Z(2)=ZZ(2,IT)/DSQRT(10.D0)
C     Z(3)=ZZ(3,IT)/DSQRT(14.D0)
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
C     ANSWER
C     RHO=(1.08407,0.98497)
C     VAR(RHO)=| 4.3726D-05 -8.3102D-06|
C              |-8.3102D-06  1.5902D-06|
C
      SUBROUTINE NLSYS(R,EPS,LIMIT)
      IMPLICIT INTEGER*4 (A-Z)
      save
      REAL*8 R(1),TOL,EPS
      CALL PARMS(N,M,KY,KX,KZ,LT,LR,R,ITER1,ITER2,TOL)
      RHO1=1
      S   =RHO1+LR
      C   =S   +M*M
      RHO2=C   +LR*LR
      S2  =RHO2+LR
      W   =S2  +M*M
      CALL Z1NLSYS(N,M,KY,KX,KZ,LT,R(S),LR,R(RHO1),R(RHO2),R(C),R(S2),
     &R(W),ITER1,ITER2,LIMIT,TOL,EPS)
      RETURN
      END
      SUBROUTINE Z1NLSYS(N,M,KY,KX,KZ,LT,S,LR,RHO1,RHO2,C,S2,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(M,M),RHO1(LR),RHO2(LR),C(LR,LR),S2(M,M),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 /ZSHORT/ ISHORT
      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
      IOUT=3
      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 Z2NLSYS(N,M,KY,KX,KZ,LT,LR,LIMIT,SAVEFC)
      IF(KZ.NE.0) CALL Z3NLSYS(N,M,KY,KX,KZ,W)
      IF(ITER2.GE.0) THEN
      ISTART=0
      DO 5 ALPHA=1,M
      DO 5 BETA=1,M
      S(ALPHA,BETA)=0.D0
5     IF(ALPHA.EQ.BETA) S(ALPHA,BETA)=1.D0
      ELSE
      ISTART=1
      CALL DGMCPY(FIXS,S,M,M)
      ENDIF
      CALL DGMCPY(S,S2,M,M)
      DO 50 ICTR2=ISTART,IABS(ITER2)
      CALL DGMCPY(S2,S,M,M)
      IF(ISHORT.LE.2) THEN
      CALL SPACER(0)
      CALL SPACER(4)
      CALL HEADER('// THE ACROSS EQUATIONS VARIANCE-COVARIANCE MATRIX/
     &IS FIXED AS SHOWN HERE FOR THE COMPUTAIONS WHICH FOLLOW///_')
      CALL DSWEEP(S2,M,EPS,IER)
      CALL DGMPNT(S2,M,M)
      CALL DGMCPY(S,S2,M,M)
      CALL HEADER('//INVERSE OF THE ACROSS EQUATIONS/VARIANCE-COVARIANCE
     & MATRIX///_')
      CALL DGMPNT(S,M,M)
      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 SYSMGN(N,M,KY,KX,KZ,LT,LR,S,RHO1,C,OBJ,RHO2,S2,
     &EPS,IER,W)
      W(IOBJ)=OBJ
      DONE=.FALSE.
      IF(ICTR1.NE.0) CALL Z4NLSYS(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
50    CONTINUE
      IF(ISHORT.LE.1) 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.0) 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 Z2NLSYS(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 /ZSHORT/ 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 NLSYS
      RHO1  =1
      S     =RHO1  +LR
      C     =S     +M*M
      RHO2  =C     +LR*LR
      S2    =RHO2  +LR
      W     =S2    +M*M
      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 ZCKZ
      IZZ   =W
      IYT   =IZZ   +KZ*KZ
      IXT   =IYT   +KY
      IZT   =IXT   +KX
      LW2   =IZT   +KZ
C     INCREMENTAL WORKSPACE USED BY SYSMGN
      ISTEP =W
      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 ZSTEP1 OF SYSMGN
C     NOTE THAT IT EXCEEDS THAT USED BY ZSTEP2
      IYT   =IW
      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('//NLSYS VERSION 2-25-85.  FOR DOCUMENTATION PRINT//
C    &NB4139.NLSYS.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 Z3NLSYS(N,M,KY,KX,KZ,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 W(1)
      IZZ=0
      IYT=IZZ+KZ*KZ
      IXT=IYT+KY
      IZT=IXT+KX
      LW =IZT+KZ
      DO 10 I=1,KZ
      DO 10 J=1,KZ
10    W(KZ*(J-1)+I)=0.D0
      DO 20 IT=1,N
      CALL DATA(IT,W(IYT+1),W(IXT+1),W(IZT+1))
      DO 20 I=1,KZ
      DO 20 J=1,KZ
20    W(KZ*(J-1)+I)=W(KZ*(J-1)+I)+W(IZT+I)*W(IZT+J)
      TOL=1.D-6
      DO 30 I=1,KZ
      IF(DABS(W(KZ*(I-1)+I)-1.D0).GT.TOL) GO TO 40
      DO 30 J=1,KZ
      IF(I.EQ.J) GO TO 30
      IF(DABS(W(KZ*(J-1)+I)).GT.TOL) GO TO 40
30    CONTINUE
      RETURN
40    CONTINUE
      CALL SPACER(0)
      CALL SPACER(4)
      CALL HEADER('//THE MATRIX OF INSTRUMENTAL VARIABLES Z/DOES NOT HAV
     &E ORTHONORMAL COLUMNS/SEE Z''Z BELOW///_')
      CALL DGMPNT(W,KZ,KZ)
      CALL SPACER(4)
      CALL HEADER('//COMPUTATIONS TERMINATED///_')
      STOP
      END
      BLOCK DATA
      COMMON /ZSHORT/ ISHORT
      DATA ISHORT/0/
      END
      SUBROUTINE Z4NLSYS(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
