SUBROUTINE MNED(A,NP,F,*)
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      DIMENSION XD(N,K1),XS(N,K2),Q(N),QSE(N),QDE(N),A(NP)
C                                                                       
C     CALCULATES THE DENSITY FUNCTION FOR THE MN EFFECTIVE DEMAND 
C     DISEQUILIBRIUM MODEL, OCZKOWSKI(1997, AUSTRALIAN ECONOMIC PAPERS,
C     36, 283-307). MODEL DEFINED BY EQN (1)-(4) & (9) PP287-8.
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+4)                    
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(K1+K2+1) IS THE DEMAND ERROR TERM VARIANCE
C	     A(K1+K2+2) IS THE SUPPLY ERROR TERM VARINACE 
C            A(K1+K2+3) IS THE DEMAND MANIPULATION COEFFICIENT
C            A(K1+K2+4) 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, DIST, 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+4.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
      IF(A(K1+K2+2).LE.0.0) RETURN 1
C
C     DEFINE PARAMETERS
      SIGD=DSQRT(A(K1+K2+1))
      SIGS=DSQRT(A(K1+K2+2))
      A2=A(K1+K2+3)
      B2=A(K1+K2+4)
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
      H1D=(Q(I)-S1-A2*(S1-QSE(I)))/SIGD
      H1S=(Q(I)-S2-B2*(S2-QDE(I)))/SIGS
      H2D=H1D
      H2S=(Q(I)-S2)/SIGS
      H3D=(Q(I)-S1)/SIGD
      H3S=H1S
      H4D=H3D
      H4S=H2S 
C
C     DEFINE STEP FUNCTIONS 
      Y=S1-QSE(I)
      Z=QDE(I)-S2
C
C     CALCULATE APPROXIMATIONS
      YA=SPLINE(Y)
      ZA=SPLINE(Z)
C
C     CALCULATE SEPERATE MN DENSITIES
      H1=(DENS(H1D)/SIGD)*(1-DIST(H1S))+(DENS(H1S)/SIGS)*(1-DIST(H1D))
      H2=(DENS(H2D)/SIGD)*(1-DIST(H2S))+(DENS(H2S)/SIGS)*(1-DIST(H2D))
      H3=(DENS(H3D)/SIGD)*(1-DIST(H3S))+(DENS(H3S)/SIGS)*(1-DIST(H3D))
      H4=(DENS(H4D)/SIGD)*(1-DIST(H4S))+(DENS(H4S)/SIGS)*(1-DIST(H4D))
C
C
C
C     CALCULATE DENSITY
      FDENS=YA*(1-ZA)*H1+YA*ZA*H2+(1-YA)*(1-ZA)*H3+(1-YA)*ZA*H4
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 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