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