SUBROUTINE BMM(A,NP,F,*) IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION XD(N,KD),XS(N,KS),XB(N,KB),P(N),Q(N) DOUBLE PRECISION A(NP) 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