C....*...1.........2.........3.........4.........5.........6.........7.*
C     SYSMGN   6/1/83
C
C     PURPOSE
C     COMPUTE A MODIFIED GAUSS-NEWTON ITERATIVE STEP FOR A CONSTRAINED
C     NONLINEAR SEEMINGLY UNRELATED REGRESSION OR A CONSTRAINED THREE-
C     STAGE LEAST SQUARES REGRESSION.  THE MODEL IS PRESUMED TO BE OF
C     THE FORM ET=Q(YT,XT,THETA) WHERE ET IS AN M-VECTOR, YT A KY-
C     VECTOR, XT A KX-VECTOR, AND THETA IS AN LT-VECTOR OF PARAMETERS
C     TO BE ESTIMATED SUBJECT TO THE CONSTRAINT THETA=G(RHO) WHERE RHO
C     IS AN LR-VECTOR.
C
C     USAGE
C     CALL SYSMGN(N,M,KY,KX,KZ,LT,LR,S,RHO,C,OBJ,RHO2,S2,EPS,IER,W)
C
C     USER SUPPLIED SUBROUTINES REQUIRED:
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 Q 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 WITH
C              USAGE  CALL DGRAM(Z,N,KZ,1.D-13,IR,WORK)  MAY BE USED TO
C              ORTHOGANALIZE Z.  WORK IS A REAL*8 N-VECTOR.  ZT IS
C              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                 REAL*8
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
C
C     ARGUMENTS
C     N    - NUMBER OF OBSERVATIONS.  SUPPLIED BY THE CALLING PROGRAM.
C            INTEGER*4
C     M    - NUMBER OF EQUATIONS IN THE SYSTEM.  M=1 FOR FITTING ONE 
C            EQUATION IS PERMISSIBLE.  SUPPLIED BY THE CALLING PROGRAM.
C            INTEGER*4
C     KY   - NUMBER OF ENDOGENOUS VARIABLES IN THE SYSTEM.  SUPPLIED BY
C            THE CALLING PROGRAM.
C            INTEGER*4
C     KX   - NUMBER OF EXOGENOUS VARIABLES IN THE SYSTEM.  SUPPLIED BY
C            THE CALLING PROGRAM.
C            INTEGER*4
C     KZ   - NUMBER OF INSTRUMENTAL VARIABLES.  SET  KZ=0  TO COMPUTE
C            A SEEMINGLY UNRELATED NONLINEAR REGRESSION.  SUPPLIED BY
C            THE CALLING PROGRAM.
C            INTEGER*4
C     LT   - LENGTH OF THETA.  SUPPLIED BY THE CALLING PROGRAM.
C            INTEGER*4
C     LR   - LENGTH OF RHO.  SUPPLIED BY THE CALLING PROGRAM.
C            INTEGER*4
C     S    - AN M BY M MATRIX CONTAINING THE INVERSE OF THE ESTIMATED
C            VARIANCE-COVARIANCE MATRIX OF THE ERRORS.  SET S TO THE
C            IDENTITY TO OBTAIN AN ORDINARY LEAST SQUARES STEP (KZ=0)
C            OR A TWO-STAGE LEAST SQUARES STEP (KZ>0).  C IS COMPUTED
C            ERRONEOUSLY WITH THIS USAGE.  SUPPLIED BY THE CALLING
C            PROGRAM.  STORED COLUMNWISE (STORAGE MODE 0)
C            REAL*8
C     RHO  - AN LR-VECTOR CONTAINING AN ESTIMATE OF RHO; THETA=G(RHO)
C            WHERE THETA IS THE PARAMETER VECTOR.  G IS A FUNCTION
C            DEFINING THE CONSTRAINTS ON THE SYSTEM.  SUPPLIED BY THE
C            CALLING PROGRAM.
C            REAL*8
C     C    - AN LR BY LR MATRIX CONTAINING THE ESTIMATED VARIANCE-
C            COVARIANCE MATRIX OF RHO PROVIDED THAT RHO IS THE SEEMINGLY
C            UNRELATED REGRESSION OR THE THREE-STAGE LEAST SQUARES
C            ESTIMATE AND THAT S IS A CONSISTENT ESTIMATE.  RETURNED BY
C            SYSMGN.  STORED COLUMNWISE (STORAGE MODE 0).
C            REAL*8
C     OBJ  - VALUE OF THE SEEMINGLY UNRELATED OBJECTIVE FUNCTION  OBJ=
C            E'(SXI)E  OR THE THREE-STAGE LEAST SQUARES OBJECTIVE FUNC-
C            TION  OBJ=E'(SXZZ')E  ACCORDING AS  KZ=0 OR KZ>0  WITH THE
C            STACKED RESIDUAL VECTOR E OF LENGTH N*M EVALUATED AT
C            THETA=G(RHO).  RETURNED BY SYSMGN.
C            REAL*8
C     RHO2 - AN LR-VECTOR CONTAINING THE NEW ESTIMATE OF THETA.  RETURN-
C            ED BY SYSMGN.
C            REAL*8
C     S2   - AN M BY M MATRIX CONTAINING THE INVERSE OF THE ESTIMATED
C            VARIANCE-COVARIANCE MATRIX OF THE ERRORS COMPUTED FROM
C            RESIDUALS WITH THETA=G(RHO).  RETURNED BY SYSMGN.  STORED
C            COLUMNWISE (STORAGE MODE 0).
C            REAL*8
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     IER  - ERROR PARAMETER CODED AS FOLLOWS
C            IER=0         NO ERROR
C            MOD(IER,10)=1 NO STEP LENGTH COULD BE FOUND WHICH WOULD
C                          YIELD A RHO2 TO IMPROVE OBJ.
C            IER>1         AN INVERSION ERROR OCCURRED.  EITHER  C OR S2
C                          IS RETURNED AS A G-INVERSE.
C
C     W    - WORK VECTOR DIMENSIONED AT LEAST AS LARGE AS LR+2*LT+
C            LT**2+MAX(LT+LT**2+LR*LT,2*M+KX+KZ+M*LT+KZ*M+KZ*M*LT)
C            PUT M=MAX0(M,KY) IN THIS COMPUTATION.
C            REAL*8
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
C
      SUBROUTINE SYSMGN(N,M,KY,KX,KZ,LT,LR,S,RHO,C,OBJ,RHO2,S2,
     &EPS,IER,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 S(M,M),RHO(LR),RHO2(LR),C(LR,LR),S2(M,M),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 Z1SYSMGN(N,M,KY,KX,KZ,S,W(ITHETA),LT,W(IQVQ),W(IQVE),OBJ,S2,
     &EPS,IER3,W(IW))
      IN=0
      CALL CONSTR(RHO,W(ITHETA),W(IG),IN)
      CALL DGMPRD(W(IQVQ),W(IG),W(IW),LT,LT,LR)
      CALL DGMAPB(W(IG),W(IW),C,LT,LR,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 Z2SYSMGN(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+1000*IER3
      RETURN
      END
      SUBROUTINE Z1SYSMGN(N,M,KY,KX,KZ,S,THETA,LT,QVQ,QVE,OBJ,S2,
     &EPS,IER,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 S(M,M),THETA(LT),QVQ(LT,LT),QVE(LT),S2(M,M),W(1)
      INTEGER*4 ALPHA,BETA,ROW,COL,APLUS1
C
C     SET UP WORKSPACE
C
      IYT=0
      IXT=IYT+MAX0(KY,M)
      IZT=IXT+KX
      IE =IZT+KZ
      IQ =IE +M
      IZE=IQ +M*LT
      IZQ=IZE+KZ*M
      LW =IZQ+KZ*M*LT
C
C     INITIALIZE ACCUMULATORS
C
      OBJ=0.D0
      DO 10 ROW=1,LT
      QVE(ROW)=0.D0
      DO 10 COL=1,LT
10    QVQ(ROW,COL)=0.D0
      DO 20 I=1,M
      DO 20 J=1,M
20    S2(I,J)=0.D0
      FN=N
      IF(KZ.GT.0) GO TO 190
C
C     COMPUTE QVQ, QVE, & OBJ FOR A SEEMINGLY UNRELATED REGRESSION
C
      DO 170 IT=1,N
      CALL DATA(IT,W(IYT+1),W(IXT+1),W(IZT+1))
      IN=0
      CALL MODEL(W(IYT+1),W(IXT+1),THETA,W(IE+1),W(IQ+1),IN)
      DO 110 ALPHA=1,M
      DO 110 BETA=1,M
110   S2(ALPHA,BETA)=S2(ALPHA,BETA)+W(IE+ALPHA)*W(IE+BETA)/FN
      DO 120 ALPHA=1,M
120   OBJ=OBJ+S(ALPHA,ALPHA)*(W(IE+ALPHA)**2)
      MLESS1=M-1
      DO 125 ALPHA=1,MLESS1
      APLUS1=ALPHA+1
      DO 125 BETA=APLUS1,M
125   OBJ=OBJ+2.D0*S(ALPHA,BETA)*W(IE+ALPHA)*W(IE+BETA)
      DO 160 ROW=1,LT
      DO 130 ALPHA=1,M
      DO 130 BETA=1,M
130   QVE(ROW)=QVE(ROW)+S(ALPHA,BETA)*W(IE+ALPHA)*W(IQ+M*(ROW-1)+BETA)
      DO 160 COL=ROW,LT
      DO 140 ALPHA=1,M
      DO 140 BETA=1,M
140   QVQ(ROW,COL)=QVQ(ROW,COL)
     &       +S(ALPHA,BETA)*W(IQ+M*(ROW-1)+ALPHA)*W(IQ+M*(COL-1)+BETA)
160   CONTINUE
170   CONTINUE
      CALL DCOND2(S2,M,W(IYT+1),0)
      CALL DSWEEP(S2,M,EPS,IER)
      CALL DCOND2(S2,M,W(IYT+1),1)
      DO 180 ROW=1,LT
      DO 180 COL=ROW,LT
180   QVQ(COL,ROW)=QVQ(ROW,COL)
      RETURN
C
C     COMPUTE QVQ, QVE & OBJ FOR A THREE-STAGE LEAST-SQUARES REGRESSION
C
190   CONTINUE
      DO 200 ALPHA=1,M
      DO 200 I=1,KZ
200   W(IZE+KZ*(ALPHA-1)+I)=0.D0
      DO 210 COL=1,LT
      DO 210 ALPHA=1,M
      DO 210 I=1,KZ
210   W(IZQ+M*KZ*(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 215 ALPHA=1,M
      DO 215 BETA=1,M
215   S2(ALPHA,BETA)=S2(ALPHA,BETA)+W(IE+ALPHA)*W(IE+BETA)/FN
      DO 220 ALPHA=1,M
      DO 220 I=1,KZ
      LOC=IZE+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=IZQ+M*KZ*(COL-1)+KZ*(ALPHA-1)+I
230   W(LOC)=W(LOC)+W(IZT+I)*W(IQ+M*(COL-1)+ALPHA)
240   CONTINUE
      DO 280 ALPHA=1,M
      DO 280 BETA=1,M
      ACC=0.D0
      DO 250 I=1,KZ
250   ACC=ACC+W(IZE+KZ*(ALPHA-1)+I)*W(IZE+KZ*(BETA-1)+I)
      OBJ=OBJ+S(ALPHA,BETA)*ACC
      DO 265 COL=1,LT
      ACC=0.D0
      DO 260 I=1,KZ
260   ACC=ACC+W(IZQ+M*KZ*(COL-1)+KZ*(ALPHA-1)+I)
     &       *W(IZE+KZ*(BETA-1)+I)
265   QVE(COL)=QVE(COL)+S(ALPHA,BETA)*ACC
      DO 275 ROW=1,LT
      DO 275 COL=ROW,LT
      ACC=0.D0
      DO 270 I=1,KZ
270   ACC=ACC+W(IZQ+M*KZ*(COL-1)+KZ*(ALPHA-1)+I)
     &       *W(IZQ+M*KZ*(ROW-1)+KZ*(BETA-1)+I)
275   QVQ(ROW,COL)=QVQ(ROW,COL)+S(ALPHA,BETA)*ACC
280   CONTINUE
      CALL DCOND2(S2,M,W(IYT+1),0)
      CALL DSWEEP(S2,M,EPS,IER)
      CALL DCOND2(S2,M,W(IYT+1),1)
      DO 290 ROW=1,LT
      DO 290 COL=ROW,LT
290   QVQ(COL,ROW)=QVQ(ROW,COL)
      RETURN
      END
      SUBROUTINE Z2SYSMGN(N,M,KY,KX,KZ,S,THETA,LT,OBJ,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 S(M,M),THETA(LT),W(1)
      INTEGER*4 ALPHA,BETA,APLUS1
C
C     SET UP WORKSPACE
C
      IYT=0
      IXT=IYT+MAX0(KY,M)
      IZT=IXT+KX
      IE =IZT+KZ
      IQ =IE +M
      IZE=IQ +M*LT
      LW =IZE+KZ*M
      IF (KZ.GT.0) GO TO 190
C
C     COMPUTE  OBJ FOR A SEEMINGLY UNRELATED REGRESSION
C
      OBJ=0.D0
      DO 170 IT=1,N
      CALL DATA(IT,W(IYT+1),W(IXT+1),W(IZT+1))
      IN=1
      CALL MODEL(W(IYT+1),W(IXT+1),THETA,W(IE+1),W(IQ+1),IN)
      DO 120 ALPHA=1,M
120   OBJ=OBJ+S(ALPHA,ALPHA)*(W(IE+ALPHA)**2)
      MLESS1=M-1
      DO 125 ALPHA=1,MLESS1
      APLUS1=ALPHA+1
      DO 125 BETA=APLUS1,M
125   OBJ=OBJ+2.D0*S(ALPHA,BETA)*W(IE+ALPHA)*W(IE+BETA)
170   CONTINUE
      RETURN
C
C     COMPUTE  OBJ FOR A THREE-STAGE LEAST-SQUARES REGRESSION
C
190   CONTINUE
      OBJ=0.D0
      DO 200 ALPHA=1,M
      DO 200 I=1,KZ
200   W(IZE+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=IZE+KZ*(ALPHA-1)+I
220   W(LOC)=W(LOC)+W(IZT+I)*W(IE+ALPHA)
240   CONTINUE
      DO 280 ALPHA=1,M
      ACC=0.D0
      DO 250 I=1,KZ
250   ACC=ACC+W(IZE+KZ*(ALPHA-1)+I)**2
      OBJ=OBJ+S(ALPHA,ALPHA)*ACC
280   CONTINUE
      MLESS1=M-1
      DO 285 ALPHA=1,MLESS1
      APLUS1=ALPHA+1
      DO 285 BETA=APLUS1,M
      ACC=0.D0
      DO 255 I=1,KZ
255   ACC=ACC+W(IZE+KZ*(ALPHA-1)+I)*W(IZE+KZ*(BETA-1)+I)
      OBJ=OBJ+2.D0*S(ALPHA,BETA)*ACC
285   CONTINUE
      RETURN
      END
