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