C	BMM.FOR ARTIFICAL DATA EXAMPLE
C	DOUBLE LOG FUNCTIONAL FORM ONLY
	IMPLICIT REAL*8(A-H,O-Z)
	DOUBLE PRECISION XD(50,4),XS(50,4),XB(50,5),P(50),Q(50)
	DOUBLE PRECISION ALABEL(20),X(20)
	DOUBLE PRECISION XX(2),FPD(2),UD(50),US(50),UQ(50),UP(50)
	COMMON/USER1/XD,XS,XB,P,Q
	COMMON/USER2/N,KD,KS,KB,K,XK1,XK2,XK3,FF
	COMMON/USER3/I,UD,US,UP,UQ
	COMMON/USER4/X,AA,BB,CC,DD,EE,FFF,PD,PS,TA
	COMMON/BSTACK/AINT(5000)
	COMMON/BSTAK/NQ,NTOP
	COMMON/BTRAT/ITRFLG
	COMMON/BINPUT/INFLG
	COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
	EXTERNAL BMM
	EXTERNAL DFP
	EXTERNAL DIST
	EXTERNAL DENS
	EXTERNAL BIVNOR
	EXTERNAL FQQQQ
	EXTERNAL FQQQQ1
	OPEN(8,FILE='c:\DATA\GQOPT\T1.DAT',STATUS='OLD')
	CALL DFLT
C	NUMBER OF VARIABLES
	KD=3
	KS=3
	KB=3
	K=KD+KS+KD
	NP=K+6
C	RESTRICTIONS FOR VARIANCES
	XK1=0.5
	XK2=1.0
	XK3=0.5
C	FUNCTIONAL FORM DOUBLE LOG EXAMPLE
	FF=4.  
C       GENERATE EXAMPLE DATA
	IPT=0
	JPT=0
	NPP=2
	MAX=1
        N=50
	NQ=5000
	ITRFLG=1
	ITERL=100
	ACC=1.D-8
	CALL LABEL(ALABEL,NPP)
	DO 10 I=1,N
C     	FOR DOUBLE-LOG XD(I,1)=XS(I,1)=1. MUST BE SET
	XD(I,1)=1.
	XS(I,1)=1.
	XB(I,1)=1.
	U2=RAND(102)
	XD(I,2)=90.*U2+20.*(1-U2)
	U3=RAND(11243)
	XD(I,3)=120.*U3+40.*(1-U3)
	U4=RAND(424)
	XS(I,2)=7.*U4+1.*(1-U4)
	U5=RAND(588)
	XS(I,3)=400.*U5+135.*(1.-U5)
	U6=RAND(424)
	XB(I,2)=10.*U6+20.*(1-U6)
	U7=RAND(88)
	XB(I,3)=10.*U7+20.*(1.-U7)
	UD(I)=DSQRT(.5)*GAUS(123)
	US(I)=DSQRT(1.)*GAUS(3243)
	UQ(I)=DSQRT(.5)*GAUS(1332)
	UP(I)=DSQRT(1.)*GAUS(14544)
10	CONTINUE
	DO 20 I=1,N
	XX(1)=10.0
	XX(2)=20.0
C	SOLVE FOR Q AND P
	CALL NEWRAP(XX,NPP,ACC,ITERL,ALABEL,FQQQQ,IER)
	Q(I)=XX(1)
	P(I)=XX(2)
20	CONTINUE
C
C	EXAMPLE DATA GENERATION COMPLETE   
C
	CALL DFLT
	NQ=5000
	MAX=1
	INFLG=0
	ITRFLG=1
C	SET STARTING VALUES
	X(1)=20.
	X(2)=-0.028
	X(3)=0.123
	X(4)=10.0
	X(5)=0.154
	X(6)=-0.057
	X(7)=0.0
	X(8)=0.2
	X(9)=-0.2
	X(10)=-0.077
	X(11)=0.286
	X(12)=.5
	X(13)=1.0
	X(14)=.5
	X(15)=1.
	CALL LABEL(ALABEL,NP)
	CALL OPT(X,NP,F,DFP,ITERL,MAX,IER,ACC,BMM,ALABEL)	
C	DETERMINE BARGAINING COEFFICIENT, SELLER PRO OF SURPLUS
C	DETERMINE PREDICTIONS FOR Q, PD, P, PS
C	RELEVANT FOR DOUBLE-LOG ONLY
	JPT=0
	IPT=0
	DO 50 I=1,N
	XX(1)=Q(I)
	XX(2)=P(I)
C	SOLVE FOR Q AND P
	CALL NEWRAP(XX,NPP,ACC,ITERL,ALABEL,FQQQQ1,IER)
C	PREDICTED BARGAINING LIMITS
	PD=AA*(XX(1)**BB)+EE
	PS=CC*(XX(1)**DD)-FFF
C	CONSUMER SURPLUS
	CSUR=(AA/(1+BB))*XX(1)**(BB+1)+EE*XX(1)-XX(2)*XX(1)
C	PRODCUER SURPLUS
	PSUR=(XX(2)*XX(1))-((CC/(DD+1))*XX(1)**(DD+1)-FFF*XX(1))
C	SELLER SURPLUS AS PROPORTION OF TOTAL
	SPS=PSUR*100/(CSUR+PSUR)
C	PRINT PREDICTED Q, TAU, SELLER PRO OF SURPLUS
C	PRINT PREDICTED PD,PREDICTED P, PREDICTED PS
50	WRITE(8,200)I,XX(1),TA,SPS,PD,XX(2),PS
200	FORMAT(I3,6(2X,F10.3))
	STOP
	END

	SUBROUTINE BMM(A,NP,F,*)
	IMPLICIT REAL*8(A-H,O-Z)
	DOUBLE PRECISION XD(50,4),XS(50,4),XB(50,5),P(50),Q(50)
	DOUBLE PRECISION A(20)
C
C     CALCULATES THE DENSITY FUNCTION FOR THE BILATERAL MONOPOLY MODEL 
C     OCZKOWSKI(1999, ECONOMIC MODELLING, 16, 53-69) 
C     MODEL DEFINED BY EQN (7)-(12) P 58.
C
C
C     NP   = TOTAL NUMBER OF PARAMETERS (NP=KD+KS+KB+6)                    
C     N    = NUMBER OF OBSERVATIONS                                     
C     KD   = NUMBER OF COEFFICIENTS IN DEMAND EQUATION (EXCLUDE QUNATITY COEFF)                  
C     KS   = NUMBER OF COEFFICIENTS IN SUPPLY EQUATION (EXCLUDE QUNATITY COEFF)                 
C     KB   = NUMBER OF COEFFICIENTS IN BARGAINING EQUATION
C     K    = KD+KS+KB
C     A( ) = COEFFICIENTS TO BE ESTIMATED; 
C	     A(1),...,A(KD) ARE DEMAND EQN COEFFICIENTS (EXCLUDE QUANTITY COEFF) 
C	     A(KD+1),...,A(KD+KS) ARE SUPPLY EQN COEFFICIENTS (EXCLUDE QUANTITY COEFF)    
C            A(KD+KS+1),...,A(KD+KS+KB) ARE BARGAINING EQN COEFFICIENTS             
C            A(K+1) IS THE DEMAND QUANTITY COEFFICENT
C	     A(K+2) IS THE SUPPLY QUNATITY COEFFICENT 
C            A(K+3) IS THE DEMAND ERROR VARIANCE
C            A(K+4) IS THE SUPPLY ERROR VARIANCE
C	     A(K+5) IS THE PRICE ERROR VARIANCE
C            A(K+6) IS THE QUANTITY ERROR VARIANCE
C     Q( ) = OBSERVATIONS ON TRANSACTED AMOUNT; PASSED THROUGH COMMON
C     XD( )= OBSERVATIONS ON RIGHT HAND VARIABLES IN DEMAND (EXCLUDING QUANTITY); 
C            MUST BE DIMENSIONED XD(N,KD) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C     XS( )= OBSERVATIONS ON RIGHT HAND VARIABLES IN SUPPLY (EXCLUDING QUANTITY); 
C            MUST BE DIMENSIONED XS(N,KS) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C     XB( )= OBSERVATIONS ON RIGHT HAND VARIABLES IN BARGAINING EQN; 
C            MUST BE DIMENSIONED XB(N,KB) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C     P( ) = OBSERVATIONS ON PRICE MUST BE    
C	     DIMENSIONED P(N) IN MAIN PROGRAM AND PASSED THROUGH COMMON
C
C	     XK1,XK2,XK3 ARE RESTRICTIONS FOR TYING ERROR VARIANCES (SET IN MAIN)
C	     FF CHOICE FOR FUNCTIONAL FORM OF DEMAND AND SUPPLY (SET IN MAIN)
C
C	     CALL EXTERNAL PROGRAMS: BIVNOR, DIST, DENS
C
C     F   = THE FUNCTION VALUE CALCULATED
C
	COMMON/USER1/XD,XS,XB,P,Q
	COMMON/USER2/N,KD,KS,KB,K,XK1,XK2,XK3,FF
C                                                                       
C	TEST FOR USER ERROR                                               
C                                                                       
      	IF(K+6.GE.NP.AND.KD.GT.0.AND.KS.GT.0.AND.KB.GT.0) GOTO 1
      	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 VARIANCES; IF YES, TAKE NONSTANDARD RETURN  
C                                                                      
1     	IF(A(K+3).LE.0.0) RETURN 1
      	IF(A(K+4).LE.0.0) RETURN 1
      	IF(A(K+5).LE.0.0) RETURN 1
      	IF(A(K+6).LE.0.0) RETURN 1
C
C     	CHECK FOR VARIANCE ASSUMPTION
      	VD=A(K+3)
      	VS=A(K+4)
      	VP=A(K+5)
      	VQ=A(K+6)
      	IF(XK1.GT.0.0)VS=(1./XK1)*VD
      	IF(XK2.GT.0.0)VP=(1./XK2)*VD
      	IF(XK3.GT.0.0)VQ=(1./XK3)*VD
      	SVD=DSQRT(VD)
      	SVS=DSQRT(VS)
     	SVP=DSQRT(VP)
      	SVQ=DSQRT(VQ)
c
c
C
C
      	SUM=0.0

C     	OBSERVATION LOOP BEGINS HERE
      	DO 200 I=1,N
      	S1=0.0
      	S2=0.0
C     	LOOP TO COMPUTE RIGHT HAND SIDE OF DEMAND EQUATION                
      	DO 6 J=1,KD                                                      
6     	S1=S1+A(J)*XD(I,J)
      	G=S1+Q(I)*A(K+1)
C     	LOOP TO COMPUTE RIGHT HAND SIDE OF SUPPLY EQUATION                
      	DO 7 J=1,KS                                                      
7     	S2=S2+A(KD+J)*XS(I,J)
      	H=S2+Q(I)*A(K+2)
C     	CHECK FUNCTIONAL FORM AND CHOOSE TRANSFORMATION
      	IF(FF-1.)10,10,20
C     	LINEAR
10    	DEL=1/(A(K+2)-A(K+1))
      	GOTO 100 
C
C
20    	IF(FF-2.)30,30,40
C     	LOG-LINEAR
30    	G=DEXP(G)
      	H=DEXP(H)	
      	DEL=1/(A(K+2)*H-A(K+1)*G)
      	GOTO 100 
C
C
40      IF(FF-3.)50,50,60
C     	LINEAR-LOG
50     	G=DLOG(G)
      	H=DLOG(H)	
      	DEL=1/((A(K+2)-A(K+1))/Q(I))
      	GOTO 100 
C
C
60      IF(FF-4.)70,70,100
C     	DOUBLE-LOG
C     	FOR DOUBLE-LOG XD(1,I)=XS(1,I)=1. MUST BE SET
70      S1=1.
      	S2=1.
C     	LOOP TO COMPUTE RIGHT HAND SIDE OF DEMAND EQUATION                
      	DO 71 J=2,KD                                                      
71     	S1=S1*(XD(I,J)**(A(J)))
      	G=S1*(Q(I)**(A(K+1)))*A(1)
C     	LOOP TO COMPUTE RIGHT HAND SIDE OF SUPPLY EQUATION                
      	DO 72 J=2,KS                                                      
72     	S2=S2*(XS(I,J)**(A(J+KD)))
      	H=S2*(Q(I)**(A(K+2)))*A(KD+1)
      	DEL=1/((A(K+2)*H-A(K+1)*G)/Q(I))
C
C  	END DATA TRANSFORMATION
C
C     	DEFINE TAU AS LOGISTIC FN
100   	S3=0.
      	DO 110 J=1,KB
110   	S3=S3+XB(I,J)*A(KD+KS+J)         
      	TA=1/(1+DEXP(-S3))
C
C     DEFINE VARIOUS TERMS
      	D2=DEL*DEL
      	XA=1.-TA
	TA2=TA*TA
	XA2=(1.-TA)**2.
	TRI=VQ*(VP+VD*TA2+VS*XA2)
	TRI=TRI+D2*(VD*VS+VD*VP+VS*VP)
C	COMPUTE RHO
	R1=-SVS*SVD*(VQ*TA*XA-VP*D2)
	R2=DSQRT(VP*VQ+VD*VQ*TA2+VD*VP*D2)
	R3=DSQRT(VP*VQ+VS*VQ+XA2+VS*VP*D2)
	RHO=R1/(R2*R3)
C	COMPUTE STANDARD ERRORS AND MEANS
	SD2=(VD*(VP*VQ+VS*VQ*XA2+VS*VP*D2))/TRI
	SS2=(VS*(VP*VQ+VD*VQ*TA2+VD*VP*D2))/TRI
	XM1=G*(VP*VQ+VS*VQ*XA2+VS*VP*D2)
	XM2=H*VD*(VP*D2-VQ*TA*XA)
	XM3=P(I)*VD*(VQ*TA+VS*D2)
	XM4=Q(I)*VD*DEL*(VP+VS*XA)
	XMD=(XM1+XM2+XM3+XM4)/TRI
	XM1=H*(VP*VQ+VD*VQ*TA2+D2*VD*VP)
	XM2=G*VS*(VP*D2-VQ*TA*XA)
	XM3=P(I)*VS*(VD*D2+VQ*XA)
	XM4=Q(I)*VS*DEL*(VP+VD*TA)
	XMS=(XM1+XM2+XM3-XM4)/TRI
C	COMPUTE AA TERM
	A1=G*G*VP*VQ*VS*(VP*VQ+VQ*VS*XA2+D2*VP*VS)
	A2=H*H*VP*VQ*VD*(VP*VQ+VQ*VD*TA2+D2*VP*VD)
	A3=2.*G*H*VD*VS*VP*VQ*(TA*XA*VQ-D2*VP)
	A4=2.*G*VD*VS*VP*VQ*(P(I)*TA*VQ+P(I)*D2*VS+
     1 	DEL*Q(I)*XA*VS+DEL*Q(I)*VP)
	A5=2.*H*VD*VS*VP*VQ*(P(I)*XA*VQ+P(I)*D2*VD-
     1	DEL*Q(I)*TA*VD-DEL*Q(I)*VP)
	A6=P(I)*VD*VS*VQ*(P(I)*D2*VD*VS+P(I)*TA2*VD*VQ+
     1	P(I)*XA2*VS*VQ+2.*DEL*Q(I)*TA*VD*VP-
     1	2.*DEL*Q(I)*XA*VS*VP)
	A7=Q(I)*Q(I)*D2*VD*VS*VP*(VD*VS+VD*VP+VS*VP)
	A8=VD*VS*VP*VQ*TRI
	A9=(G*G/VD)+(H*H/VS)+(P(I)*P(I)/VP)+(Q(I)*Q(I)/VQ)
	AA=((A1+A2-A3+A4+A5+A6+A7)/A8)-A9
C
C	COMPUTE NORMALISED VARIBALES
	XLD=(P(I)-XMD)/DSQRT(SD2)
	XLS=(P(I)-XMS)/DSQRT(SS2)
C	COMPUTE COMPONENTS OF DENSITY
	F1=DSQRT(SD2)*DSQRT(SS2)*DSQRT(1.-RHO*RHO)
	F1=F1/(3.141592654*VD*VS*VQ*VP)
	F2=DEXP(AA/2.)
	CALL BIVNOR(XLD,-XLS,-RHO,F3)
	FDENS=F1*F2*F3
C
C     	CHECK IF DENSITY POSITIVE; OTHERWISE TAKE NONSTANDARD RETURN
      	IF(FDENS.LE.0.0) RETURN 1
200   	SUM=SUM+DLOG(FDENS)
      	F=SUM
C	RESETING VARIANCES FOR RESTRICTED CASE
	A(K+4)=SVS*SVS
	A(K+5)=SVP*SVP
        A(K+6)=SVQ*SVQ
	RETURN 
	END


	SUBROUTINE BIVNOR(AH,AK,R,B)
	IMPLICIT REAL*8(A-H,O-Z)
C
C     	CALCULATES THE BIVARIATE NORMAL DISTRIBUTION FUNCTION
C
	TWOPI=6.283185307179587D0
	B=0.0D0
	IDIG=15
	GH=DIST(-AH)
	GH=GH/2.0D0
	GK=DIST(-AK)
	GK=GK/2.0
	IF(R)10,30,10
10	RR=1.0-R*R
	IF(RR)20,40,100
20	STOP
30	B=4.0*GH*GK
	GOTO 350
40	IF(R)50,70,70
50	IF(AH+AK)60,350,350
60	B=2.0*(GH+GK)-1.0
	GOTO 350
70	IF(AH-AK)80,90,90
80	B=2.0*GK
	GOTO 350
90	B=2.0*GH
	GOTO 350
100	SQR=DSQRT(RR)
	IF(IDIG-15)120,110,120
110	CON=TWOPI*1.D-15/2.0
	GOTO 140
120	CON=TWOPI/2.0
	DO 130 I=1,IDIG
	CON=CON/10.0
130	CONTINUE
140	IF(AH)170,150,170
150	IF(AK)190,160,190
160	B=DATAN(R/SQR)/TWOPI+0.25
	GOTO 350
170	B=GH
	IF(AH*AK)180,200,190
180	B=B-0.5
190	B=B+GK
	IF(AH)200,340,200
200	WH=-AH
	WK=(AK/AH-R)/SQR
	GW=2.0*GH
	IS=-1
210	SGN=-1.0
	T=0.0
	IF(WK)220,320,220
220	IF(DABS(WK)-1.0)270,230,240
230	T=WK*GW*(1.0-GW)/2.0
	GOTO 310
240	SGN=-SGN
	WH=WH*WK
	G2=DIST(WH)
	WK=1.0/WK
	IF(WK)250,260,260
250	B=B+0.5
260	B=B-(GW+G2)/2.0+GW*G2
270	H2=WH*WH
	A2=WK*WK
	H4=H2/2.0
	EX=DEXP(-H4)
	W2=H4*EX
	AP=1.0
	S2=AP-EX
	SP=AP
	S1=0.0
	SN=S1
	CONEX=DABS(CON/WK)
	GOTO 290
280	SN=SP
	SP=SP+1.0
	S2=S2-W2
	W2=W2*H4/SP
	AP=-AP*A2
290	CN=AP*S2/(SN+SP)
	S1=S1+CN
	IF(DABS(CN)-CONEX)300,300,280
300	T=(DATAN(WK)-WK*S1)/TWOPI
310	B=B+SGN*T
320	IF(IS)330,350,350
330	IF(AK)340,350,340
340	WH=-AK
	WK=(AH/AK-R)/SQR
	GW=2.0*GK
	IS=1
	GOTO 210
350	IF(B)360,370,370
360	B=0.0
370	IF(B-1.0)390,390,380
380	B=1.0
390	B=B
	RETURN 
	END
	

	FUNCTION DIST(X)
	IMPLICIT REAL*8(A-H,O-Z)
C
C     	CALCULATES THE STANDARD NORMAL DISTRIBUTION FUNCTION
C
	Y=X
	IF(Y.LT. 0.0) Y = -Y
	P=0.2316419D0
	B1=0.319381530D0
	B2=-0.356563782D0
	B3=1.781477937D0
	B4=-1.821255978D0
	B5=1.330274429D0
	T=1./(1.+P*Y)
	DIST=1.-DENS(Y)*(B1*T+B2*T**2+B3*T**3+B4*T**4+B5*T**5)
	IF(X.GT.0.0) GOTO 10
	DIST = 1. - DIST
10	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


	SUBROUTINE FQQQQ(XX,NPP,FPD,*)
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION Q(50),XD(50,4),XS(50,4),P(50),XB(50,5)
	DOUBLE PRECISION XX(2),FPD(2),UD(50),US(50),UQ(50),UP(50)
	COMMON/USER1/XD,XS,XB,P,Q
	COMMON/USER2/N,KD,KS,KB,K,XK1,XK2,XK3,FF
	COMMON/USER3/I,UD,US,UP,UQ	
C
C	EMPLOYED ONLY TO GENERATE ARTIFICAL DATA (DOUBLE LOG ONLY)
C
	IF(XX(1).LE.0.0)RETURN 1
	IF(XX(2).LE.0.0)RETURN 1
	PD=20.0*(XD(I,2)**(-0.028))*(XD(I,3)**0.123)*(XX(1)**(-0.077))
	PD=PD+UD(I)
	PS=10.0*(XS(I,2)**0.154)*(XS(I,3)**(-0.057))*(XX(1)**0.286)
	PS=PS+US(I)
	PB=0.0+0.2*XB(I,2)-0.2*XB(I,3)
	TA=1./(1+DEXP(-PB))
	DEL=XX(1)/(0.286*PS+0.077*PD)
	FPD(1)=DEL*(PD-PS)+UQ(I)-XX(1)
	FPD(2)=TA*PD+(1-TA)*PS+UP(I)-XX(2)
	RETURN
	END 


	SUBROUTINE FQQQQ1(XX,NPP,FPD,*)
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION Q(50),XD(50,4),XS(50,4),P(50),XB(50,5),X(20)
	DOUBLE PRECISION XX(2),FPD(2),UD(50),US(50),UQ(50),UP(50)
	COMMON/USER1/XD,XS,XB,P,Q
	COMMON/USER2/N,KD,KS,KB,K,XK1,XK2,XK3,FF
	COMMON/USER3/I,UD,US,UP,UQ
	COMMON/USER4/X,AA,BB,CC,DD,EE,FFF,PD,PS,TA
C
C	EMPLOYED TO PREDICT Q AND P (DOUBLE LOG ONLY)
C
	IF(XX(1).LE.0.0)RETURN 1
	IF(XX(2).LE.0.0)RETURN 1
	S1=1.
	S2=1.
       	DO 10 J=2,KD
10     	S1=S1*(XD(I,J)**(X(J)))
       	AA=S1*X(1)
       	DO 11 J=2,KS
11     	S2=S2*(XS(I,J)**(X(J+KD)))
       	CC=S2*X(KD+1)
       	BB=X(K+1)
       	DD=X(K+2)
	ZD=(XX(2)-AA*(XX(1)**BB))/DSQRT(X(K+3))
	ZS=(XX(2)-CC*(XX(1)**DD))/DSQRT(X(k+4))
	EE=DSQRT(X(K+3))*(DENS(ZD)/DIST(-ZD))
	FFF=DSQRT(X(K+4))*(DENS(ZS)/DIST(ZS))
	S3=0.
	DO 12 J=1,KB
12	S3=S3+XB(I,J)*X(KD+KS+J)
	TA=1/(1+DEXP(-S3))
	A1=TA*AA
	A2=(1-TA)*CC
	A3=TA*EE-(1-TA)*FFF
	FPD(1)=XX(2)-A1*(XX(1)**BB)-A2*(XX(1)**DD)-A3	
	FPD(2)=CC*(XX(1)**DD)*(1+DD)-AA*(XX(1)**BB)*(1+BB)-EE+FFF
	RETURN
	END