C	EXAMPLE OF GTZED.FOR USING ARTIFICAL DATA DESCRIBED IN
C	OCZKOWSKI (1990, JQE, P196).

	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION XD(50,4),XS(50,4),Q(50),QSE(50),QDE(50)
	DIMENSION ALABEL(11),AINT(5000),QD(50),QS(50),X(11)
	COMMON/USER1/Q,XD,XS,QSE,QDE
	COMMON/USER2/N,K1,K2
	COMMON/USER3/APP
	COMMON/BSTACK/AINT
	COMMON/BSTAK/NQ,NTOP
	COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
	COMMON/BTRAT/ITRFLG
	COMMON/BINPUT/INFLG
        EXTERNAL GTZED
	EXTERNAL DFP
	EXTERNAL DENS
	EXTERNAL SPLINE
	OPEN(8,FILE='c:\DATA\GQOPT\T1.DAT',STATUS='OLD')
	OPEN(9,FILE='c:\DATA\GQOPT\T2.DAT',STATUS='OLD')
	CALL DFLT
        K1=4
        K2=4
	NP=K1+K2+3
	MAX=1
	N=50
	NQ=5000
	ITRFLG=1
	ITERL=100

C	GENERATE EXAMPLE DATA
	DO 10 I=1,N
	XD(I,1)=1.
	XS(I,1)=1.
	U1=RAND(110)
	XD(I,2)=10.*U1+2.*(1-U1)
	XS(I,2)=XD(I,2)
	U2=RAND(102)
	XD(I,3)=90.*U2+20.*(1-U2)
	U3=RAND(11243)
	XD(I,4)=120.*U3+40.*(1-U3)
	U4=RAND(424)
	XS(I,3)=7.*U4+1.*(1-U4)
	U5=RAND(55388)
10	XS(I,4)=400.*U5+135.*(1.-U5)
	DO 15 I=1,N
	QD(I)=18.0*XD(I,1)-13.0*XD(I,2)-0.37*XD(I,3)+1.6*XD(I,4)
15	QS(I)=0.0*XS(I,1)+3.5*XS(I,2)-0.54*XS(I,3)+0.2*XS(I,4)
	DO 20 I=2,N
        QSE(I)=QS(I-1)
20      QDE(I)=QD(I-1)
        QSE(1)=QS(1)
        QDE(1)=QD(1)
	DO 30 I=1,N
        QD(I)=QD(I)+0.2*(DMAX1(0,QD(I)-QSE(I)))
	QS(I)=QS(I)-0.35*(DMIN1(0,QDE(I)-QS(I)))
	UU=8.*GAUS(123)
30	Q(I)=DMIN1(QD(I),QS(I))+UU

C	STARTING VALUES FOR PARAMETERS
	DO 40 J=1,NP
40	X(J)=0.0
	X(K1+K2+1)=10.
	CALL LABEL(ALABEL,NP)

C	SET ACCUARCY AND APPROXIMATION FOR 1ST RUN
	ACC=1.D-4
	APP=70.0
	INFLG=0
        CALL OPT(X,NP,F,DFP,ITERL,MAX,IER,ACC,GTZED,ALABEL)

C	SET ACCURACY AND APPROXIMATION FOR 2ND RUN
	ACC=1.D-12
	APP=1.D-4
        CALL PUNCH(X,NP)
        REWIND 9
	INFLG=1
        CALL OPT(X,NP,F,DFP,ITERL,MAX,IER,ACC,GTZED,ALABEL)
        STOP
	END

      SUBROUTINE GTZED(A,NP,F,*)
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      DIMENSION XD(50,4),XS(50,4),Q(50),QSE(50),QDE(50),A(11)
C                                                                       
C     CALCULATES THE DENSITY FUNCTION FOR THE GTZ EFFECTIVE DEMAND 
C     DISEQUILIBRIUM MODEL, OCZKOWSKI(1990, JOURNAL OF QUNATITATIVE
C     ECONOMICS,9, 185-201). MODEL DEFINED BY EQN (14)-(19) PP 102-3.
C
C     SUBROUTINE SHOULD BE CALLED AT LEAST TWICE BY CHANGING THE SPLINE
C     APPROXIMATION COEFFICIENT (APP) TO BE SET IN MAIN PROGRAM AND PASSED
C     THROUGH AS COMMON/USER3/APP
C
C     NP   = TOTAL NUMBER OF PARAMETERS (NP=K1+K2+3)                    
C     N    = NUMBER OF OBSERVATIONS                                     
C     K1   = NUMBER OF COEFFICIENTS IN DEMAND EQUATION                  
C     K2   = NUMBER OF COEFFICIENTS IN SUPPLY EQUATION                  
C     A( ) = COEFFICIENTS TO BE ESTIMATED; A(1),...,A(K1) ARE DEMAND    
C            EQUATION COEFFICIENTS, A(K1+1),...,A(K1+K2) ARE SUPPLY     
C            EQUATION COEFFICIENTS, 
C            A(NP-2) IS THE QUANTITY ERROR TERM VARIANCE 
C            A(NP-1) IS THE DEMAND MANIPULATION COEFFICIENT
C            A(NP) IS THE SUPPLY MANIPULATION COEFFICIENT
C     Q( ) = OBSERVATIONS ON TRANSACTED AMOUNT; PASSED THROUGH COMMON
C     XD( )= OBSERVATIONS ON RIGHT HAND VARIABLES IN DEMAND; MUST BE
C	     DIMENSIONED XD(N,K1) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C     XS( )= OBSERVATIONS ON RIGHT HAND VARIABLES IN SUPPLY; MUST BE
C	     DIMENSIONED XS(N,K2) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C     QSE( )= SUPPLY EXPECTATIONS BY DEMANDERS MUST BE    
C	     DIMENSIONED QSE(N) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C     QDE( )= DEMAND EXPECTATIONS BY SUPPILERS MUST BE    
C	     DIMENSIONED QDE(N) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C
C	     CALL EXTERNAL PROGRAMS: DENS, SPLINE
C
C     F   = THE FUNCTION VALUE CALCULATED
C
      COMMON/USER1/Q,XD,XS,QSE,QDE
      COMMON/USER2/N,K1,K2
      COMMON/USER3/APP
      COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE                     
C                                                                       
C     TEST FOR USER ERROR                                               
C                                                                       
      IF(K1+K2+3.GE.NP.AND.K1.GT.0.AND.K2.GT.0 ) GOTO 5
      IF(IPT.GT.0) WRITE (NFILE,1000)                                   
      IF(JPT.GT.0) WRITE (MFILE,1000)                                   
1000  FORMAT(' ERROR IN K1 OR K2 OR NP...EXECUTION TERMINATED')         
      STOP                                                              
C                                                                       
C     CHECK FOR NONPOSITIVE VARIANCE; IF YES, TAKE NONSTANDARD RETURN  
C                                                                      
5     IF(A(K1+K2+1).LE.0.0) RETURN 1              
C
C     DEFINE PARAMETERS
      SIG=DSQRT(A(K1+K2+1))
      A2=A(K1+K2+2)
      B2=A(K1+K2+3)
C
      SUM=0.0

C     OBSERVATION LOOP BEGINS HERE
      DO 100 I=1,N
      S1=0.0
      S2=0.0
C     LOOP TO COMPUTE RIGHT HAND SIDE OF DEMAND EQUATION                
      DO 10 J=1,K1                                                      
10    S1=S1+A(J)*XD(I,J)
C     LOOP TO COMPUTE RIGHT HAND SIDE OF SUPPLY EQUATION                
      DO 20 J=1,K2                                                      
20    S2=S2+A(K1+J)*XS(I,J)
C
C     DEFINE NORMALISED RESIDUALS
      F1=(Q(I)-S1-A2*(S1-QSE(I)))/SIG
      F2=(Q(I)-S1)/SIG
      G1=(Q(I)-S2-B2*(QDE(I)-S2))/SIG
      G2=(Q(I)-S2)/SIG
C
C     DEFINE STEP FUNCTIONS 
      X1=S1-S2*(1-B2)-B2*QDE(I)
      X2=S1*(1+A2)-B2*QDE(I)-(1-B2)*S2-A2*QSE(I)
      X3=S1-S2
      X4=S1*(1+A2)-A2*QSE(I)-S2
      Y=S1-QSE(I)
      Z=QDE(I)-S2
C
C     CALCULATE APPROXIMATIONS
      X1A=SPLINE(X1)
      X2A=SPLINE(X2)
      X3A=SPLINE(X3)
      X4A=SPLINE(X4)
      YA=SPLINE(Y)
      ZA=SPLINE(Z)
C
C     CALCULATE PRODUCT OF DENISTIES AND SWITCHES
      D1=(DENS(F1)/SIG)*YA*((1-ZA)*(1-X2A)+ZA*(1-X4A))
      D2=(DENS(F2)/SIG)*(1-YA)*((1-ZA)*(1-X1A)+ZA*(1-X3A))
      D3=(DENS(G1)/SIG)*(1-ZA)*((1-YA)*X1A+YA*X2A)
      D4=(DENS(G2)/SIG)*ZA*((1-YA)*X3A+YA*X4A)
C
C     CALCULATE DENSITY
      FDENS=D1+D2+D3+D4
C     CHECK IF DENSITY POSITIVE; OTHERWISE TAKE NONSTANDARD RETURN
      IF(FDENS.LE.0.0) RETURN 1
100   SUM=SUM+DLOG(FDENS)
      F=SUM
      RETURN
      END                                                               


      FUNCTION SPLINE(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/USER3/APP
C
C     CALCULATES THE TISHLER & ZANG 5TH ORDER SPLINE APPROXIMATION
C
	APPM=-APP
	IF(X.LE.APPM) SPLINE=0.
	IF(APPM.LT.X.AND.X.LT.APP)SPLINE=(3/(16*APP**5))*(X**5)
     1	-(5/(8*APP**3))*(X**3)+(15/(16*APP))*X+0.5
	IF(APP.LE.X)SPLINE=1.
	RETURN
	END




	FUNCTION DENS(X)
	IMPLICIT REAL*8(A-H,O-Z)
C
C      CALCULATES THE STANDARD NORMAL DENISTY FUNCTION
C
	PI=3.141592653589D0*2.D0
	Z1=1./(PI**0.5)
	Z2=DEXP(-(X*X)/2.)
	DENS=Z1*Z2
	RETURN
	END