C     MAIN PROGRAM ==============================================
      COMMON ML,NL
      DOUBLE PRECISION YO(200),XO(200,10),A(10),EPS,QD,CN
      DOUBLE PRECISION Y(200),X(200,10)
      INTEGER  IUX,IU,IV,IUU
      CHARACTER *40 OFIL,IFIL
C     ==============================================================
      WRITE(*,*)'THE GENERALIZED MAXIMUM ENTROPY LEUVEN ESTIMATION'
      WRITE(*,*)'COMPUTER PROGRAM AND METHOD BY PROF. SK MISHRA'
      WRITE(*,*)'DEPT. OF ECONOMICS, NEHU, SHILLONG'
C     PARAMETERS FOR DIMENSION: CHANGE THEM IF THE PROBLEM IS LARGER
      EPS=0.00001  ! CHANGE ACCURACY REQUIREMENT IF NEEDED
      NL=200
      ML=10
      NTRY=1
      WRITE(*,*)'------- HOW TO STORE DATA IN THE INPUT FILE --------'
      WRITE(*,*)' FOR DETAILS READ NOTES GIVEN BELOW THE PROGRAM'
      WRITE(*,*)'----------------------------------------------------'
      WRITE(*,*)'DATA ARE TO BE STORED SUCH THAT EACH VARIABLE HAS BEEN'
      WRITE(*,*)'PUT IN M COLUMNS (OF N ROWS) LIKE Y X1 X2 X3...XM'
      WRITE(*,*)'SO THAT IT IS A MATRIX IN N ROWS AND M+1 COLUMNS'
      WRITE(*,*)'THE LAST VARIABLE IS A UNIT COLUMN FOR INTERCEPT'
      WRITE(*,*)'----------------------------------------------------'
      WRITE(*,*)'------NAME THE FILE IN WHICH DATA Y AND X ARE STORED'
      WRITE(*,*)'-------------FILE NAME WITHIN SINGLE QUOTES'
      READ(*,*) IFIL
      WRITE(*,*)'OUTPUT FILE (WITHIN SINGLE QUOTES) ?'
      READ(*,*) OFIL
      NOFIL=7
C     CHANGE THE FORMAT IF REQUIRED
10    FORMAT(8F10.2)
C          -------------------------------------------------------------
      OPEN(UNIT=NOFIL,FILE=OFIL)
      WRITE(*,*)'FEED METRIC - SUM{ABS(U)**METRIC}'
      READ(*,*) NPOW
      WRITE(*,*)'FEED METRIC FOR GMEL: 1=ABS(MMEL), 2=EUCLID(MEL)'
      READ(*,*) NPR
      WRITE(*,*)' '
      WRITE(*,*)'MEL1 (0) OR MEL2 (1)'
      READ(*,*) MEL1
      WRITE(*,*)' '
      WRITE(*,*)'FEED NO. OF OBSERVATIONS & NO. OF INDEPENDET INCLUDING
     * THE LAST UNITARY VARIABLES'
      READ(*,*) N,M
      WRITE(*,*)' '
      WRITE(*,*)' --- SEED MUST LIE BETWEEN -32767 AND 32767'
      WRITE(*,*)'SEED FOR RANDOM WALK AND ITERATIONS'
      READ(*,*) IUX,NITR
      WRITE(*,*)' '
      WRITE(*,*)' '
      WRITE(*,*)' ---------------------- RESULTS FOLLOW ------------'
      NNH=1
      DO 100 ITR=1,NTRY
C      WRITE(*,*)'TRIAL NO. = ',ITR
      OPEN(11,FILE=IFIL)
      DO 1 I=1,N
      READ(11,*) Y(I),(X(I,J),J=1,M),QQQ
      YO(I)=Y(I)
      DO 2 J=1,M
      XO(I,J)=X(I,J)
2     CONTINUE
1     CONTINUE
      CLOSE(11)
11    FORMAT(6F8.2)
C     ---------------------------------------------------------
      CALL RWALK(A,M,N,YO,XO,IUX,EPS,NITR,QD,NPOW,NNH,NOFIL,NPR,CN,
     * ITR,MEL1)
100   CONTINUE
      WRITE(NOFIL,*)'THESE ARE MEL AND OLS COEFFICIENTS RESPECTIVELY'
      CLOSE(NOFIL)
      WRITE(*,*)'RESULTS STORED IN FILE ',OFIL
      END
C     ------------------------------------------------------------------
      SUBROUTINE GENU(IU,IV,U,M)
C     GENERATES A VECOR U(M) OF UNIFORMLY DISTRIBUTED RANDOM NUMBERS
C     BETWEEN -1 AND +1. THE EUCLIDEAN NORM OF U IS UNITY.
C     REF: RAO, SS (1978) OPTIMIZATION: THEORY AND APPLICATIONS.
C     WILEY EASTERN LIMITED, NEW DELHI. PP. 252-257.
C     COMMON /RNDM/IU,IV
      DOUBLE PRECISION U(M),R,RAND
    1 R=0.D0
      DO I=1,M
      CALL RANDOM(IU,IV,RAND)
      U(I)=DBLE(RAND-.5)*2
      R=R+U(I)**2
      ENDDO
      IF(R.LE.1.D0) THEN
      R=DSQRT(R)
      DO I=1,M
      U(I)=U(I)/R
      ENDDO
      ELSE
      GOTO 1  ! 'NORM NOT MET'
      ENDIF
      RETURN
      END
C     ----------------------------------------------------------------
      SUBROUTINE RANDOM(IU,IV,RAN)
      DOUBLE PRECISION  RAN
      INTEGER IU,IV
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAN=IV
       IU=IV
       RAN=RAN*0.4656613E-09
       RETURN
       END
C     =================================================================
C     --------------------------------------------------------------
      SUBROUTINE RWALK(A,M,NV,YO,XO,IU,EPS,NITR, VO,NPOW,NNH,NOFIL,NPR,
     * CN,ITR,MEL1)
C     RANDOM WALK METHOD TO FIND MIN OF A MULTIVARIABLE FUNCTION
      INTEGER IU,IV
      DOUBLE PRECISION A(10),R(10),AR(10),YO(200),XO(200,10),AA(10)
      DOUBLE PRECISION LAMBDA,EPS,F,VO,VR,RNORM,YH,VOO,HMEL,CN,RAN
      DOUBLE PRECISION CVX(10,10),CVY(10),EVEC(10,10),EVAL(10)
      DOUBLE PRECISION V(10,10),W(10,10),P(10),ER(200),RSQUARE
      DIMENSION MM(10)
      CHARACTER *11 OFIL
      ITMAX=100
      VOO=10.0**10
      DO 999 IIT=1,NITR
      LAMBDA=20.0
C     INITIALISATION OF DECISION VARIABLES
      DO 7 I=1,M
7     A(I)=1.0
C     --------------- FUNCTION EVALUATION ---------------------
      F=0.0
      DO 11 I=1,NV
      YH=0.0
      DO 12 J=1,M
12    YH=YH+XO(I,J)*A(J)
      F=F+DABS(YO(I)-YH)**NPOW
      ER(I)=YO(I)-YH
11    CONTINUE
      CALL MEL(ER,A,M,NV,NPR,HMEL,MEL1)
      F=HMEL+ F*ABS(MEL1-1)
      VO=F
C     ---------------------------------------------------------
      IT=0
200   IT=IT+1
      IF(IT.GT.1000) THEN
      WRITE(*,*)'NO CONVERGENCE IN 1000 ITERATIONS',LAMBDA,IMP
      GOTO 1000
      ENDIF
      LAMBDA=LAMBDA/2.0
      IMP=0
      DO 100 II=1,ITMAX
C     GENERATE M UNIFORMLY DISTRIBUTED RANDOM NUMBERS (-1,1)
150   CALL GENU(IU,IV,R,M)

2     CONTINUE
C     NORMALISE THE RANDOM NUMBERS
      RNORM=0.0
      DO 3 I=1,M
      RNORM=RNORM+R(I)**2
3     CONTINUE
      RNORM=DSQRT(RNORM)
      IF(RNORM.GT.1.0) GOTO 150
      DO 4 I=1,M
      R(I)=R(I)/RNORM
4     CONTINUE
C     ADD RANDOM NUMBERS TIMES LAMBDA TO A VECTOR
      DO 5 I=1,M
      AR(I)=A(I)+LAMBDA*R(I)
5     CONTINUE
C     --------------- FUNCTION EVALUATION ---------------------
      F=0.0
      DO 13 I=1,NV
      YH=0.0
      DO 14 J=1,M
14    YH=YH+XO(I,J)*AR(J)
      F=F+DABS(YO(I)-YH)**NPOW
      ER(I)=YO(I)-YH
13    CONTINUE
      CALL MEL(ER,AR,M,NV,NPR,HMEL,MEL1)
      F=HMEL+ F*ABS(MEL1-1)
      VR=F
C     ---------------------------------------------------------
      IF(VR.LT.VO) THEN
      VO=VR
      DO 6 I=1,M
      A(I)=AR(I)
6     CONTINUE
      IMP=1
      ENDIF
100   CONTINUE
      IF((IMP.EQ.0).OR.(LAMBDA.GT.EPS)) GOTO 200
1000  CONTINUE
      IF(VOO.GT.VO) THEN
      VOO=VO
      DO 998 I=1,M
      AA(I)=A(I)
998   CONTINUE
      ENDIF
999   CONTINUE
      DO 997 I=1,M
      A(I)=AA(I)
997   CONTINUE
      VO=VOO
      WRITE(NOFIL,*)'MEL ESTIMATES OF REGRESSION COEFFICIENTS'
      CALL RSQR(NV,M,XO,YO,A,RSQUARE)
      WRITE(NOFIL,*)(A(I),I=1,M),'  RSQUARE=',RSQUARE
      DO 31 J=1,M
      CVY(J)=0.0
      DO 31 JJ=1,M
      CVX(J,JJ)=0.0
      DO 32 I=1,NV
32    CVX(J,JJ)=CVX(J,JJ)+XO(I,J)*XO(I,JJ)
      DO 33 I=1,NV
33    CVY(J)=CVY(J)+YO(I)*XO(I,J)
31    CONTINUE
      IF(ITR.EQ.1) THEN
      OPEN(9,FILE='XTX')
      WRITE(9,*) M
      DO 38 J=1,M
      WRITE(9,*)(CVX(J,JJ),JJ=1,M)
38    CONTINUE
      CLOSE(9)
      ENDIF
      NN=1
      CALL EIGEN(CVX,M,NN,V,W,P,MM)
C     =====================================
      DO 50 I=1,M
      F=0.0
      DO 51 J=1,M
51    F=F+V(J,I)*V(J,I)
      F=DSQRT(F)
      DO 52 J=1,M
52    V(J,I)=V(J,I)/F
50    CONTINUE
      DO 34 J=1,M
      DO 341 JJ=1,M
341   W(J,JJ)=0.0
      IF(DABS(CVX(J,J)).GT.1.0D-99) W(J,J)=1.0/CVX(J,J)
34    CONTINUE
      DO 35 J=1,M
      DO 35 JJ=1,M
      CVX(J,JJ)=0.0
      DO 35 I=1,M
      CVX(J,JJ)=CVX(J,JJ)+V(J,I)*W(I,JJ)
35    CONTINUE
      DO 36 J=1,M
      DO 36 JJ=1,M
      W(J,JJ)=0.0
      DO 36 I=1,M
      W(J,JJ)=W(J,JJ)+CVX(J,I)*V(JJ,I)
36    CONTINUE
      DO 37 J=1,M
      P(J)=0.0
      DO 37 JJ=1,M
      P(J)=P(J)+W(J,JJ)*CVY(JJ)
37    CONTINUE
      DO 53 J=1,M
      P(J)=P(J)/FLOAT(M)
53    CONTINUE
      WRITE(*,*)' '
      WRITE(NOFIL,*)'OLS ESTIMATES OF REGRESSION COEFFICIENTS'
      CALL RSQR(NV,M,XO,YO,P,RSQUARE)
      WRITE(NOFIL,*)(P(J),J=1,M),'  RSQUARE=',RSQUARE
      RETURN
      END
      SUBROUTINE MEL(ER,A,M,N,NPR,HMEL,MEL1)
      DOUBLE PRECISION HMEL,S,A(10),P(10),ER(200),T
      S=0.0
      DO 1 I=1,M
      P(I)=DABS(A(I))**NPR
      S=S+P(I)
1     CONTINUE
      DO 2 I=1,M
      P(I)=P(I)/S
2     CONTINUE
      HMEL=S*DLOG(S)
      DO 3 I=1,M
      T=0.0
      IF(P(I).NE.0.0) T=P(I)*DLOG(P(I))
      HMEL=HMEL+T
3     CONTINUE
      IF(MEL1.EQ.0) RETURN
      S=0
      DO 4 I=1,N
      ER(I)=DABS(ER(I))**NPR
      S=S+ER(I)
4     CONTINUE
      DO 5 I=1,N
      ER(I)=ER(I)/S
5     CONTINUE
      S=S*DLOG(S)
      DO 6 I=1,N
      S=S+ER(I)*DLOG(ER(I))
6     CONTINUE
      HMEL=HMEL+S
      RETURN
      END
      SUBROUTINE EIGEN(A,N,NN,V,W,P,MM)
C     THIS SUBROUTINE IS ADAPTED FROM KRISHNAMURTHY AND SEN (1976)
      DOUBLE PRECISION A(10,10),V(10,10),W(10,10),P(10)
      DOUBLE PRECISION PMAX,EPLN,TAN,SIN,COS,AI,TT,TA,TB
      DIMENSION MM(10)
C     ------------ INITIALISATION -------------------------
      DO 50 I=1,N
      DO 51 J=1,N
      V(I,J)=0.0
51    W(I,J)=0.0
      P(I)=0.0
50    CONTINUE
      PMAX=0
      EPLN=0
      TAN=0
      SIN=0
      COS=0
      AI=0
      TT=0
      EPLN=1.0D-310
C     ------------------------------------------------------
      IF(NN.NE.0) THEN
      DO 3 I=1,N
      DO 3 J=1,N
      V(I,J)=0.0
      IF(I.EQ.J) V(I,J)=1.0
3     CONTINUE
      ENDIF
2     NR=0
5     MI=N-1
      DO 6 I=1,MI
      P(I)=0.0
      MJ=I+1
      DO 6 J=MJ,N
      IF(P(I).GT.DABS(A(I,J))) GO TO 6
      P(I)=DABS(A(I,J))
      MM(I)=J
6     CONTINUE
7     DO 8 I=1,MI
      IF(I.LE.1) GOTO 10
      IF(PMAX.GT.P(I)) GOTO 8
10    PMAX=P(I)
      IP=I
      JP=MM(I)
8     CONTINUE
      IF (PMAX.LE.EPLN) THEN
      GO TO 12
      ENDIF
      NR=NR+1
13    TA=2.0*A(IP,JP)
      TB=(DABS(A(IP,IP)-A(JP,JP))+
     * DSQRT((A(IP,IP)-A(JP,JP))**2+4.0*A(IP,JP)**2))
      TAN=TA/TB
      IF(A(IP,IP).LT.A(JP,JP)) TAN=-TAN
14    COS=1.0/DSQRT(1.0+TAN**2)
      SIN=TAN*COS
      AI=A(IP,IP)
      A(IP,IP)=(COS**2)*(AI+TAN*(2.0*A(IP,JP)+TAN*A(JP,JP)))
      A(JP,JP)=(COS**2)*(A(JP,JP)-TAN*(2.0*A(IP,JP)-TAN*AI))
      A(IP,JP)=0.0
      IF(A(IP,IP).GE.A(JP,JP)) GO TO 15
      TT=A(IP,IP)
      A(IP,IP)=A(JP,JP)
      A(JP,JP)=TT
      IF(SIN.GE.0) GO TO 16
      TT=COS
      GO TO 17
16    TT=-COS
17    COS=DABS(SIN)
      SIN=TT
15    DO 18 I=1,MI
      IF(I-IP) 19, 18, 20
20    IF(I.EQ.JP)GO TO 18
19    IF(MM(I).EQ.IP) GO TO 21
      IF(MM(I).NE.JP) GO TO 18
21    K=MM(I)
      TT=A(I,K)
      A(I,K)=0.0
      MJ=I+1
      P(I)=0.0
      DO 22 J=MJ,N
      IF(P(I).GT.DABS(A(I,J))) GO TO 22
      P(I)=DABS(A(I,J))
      MM(I)=J
22    CONTINUE
      A(I,K)=TT
18    CONTINUE
      P(IP)=0.0
      P(JP)=0.0
      DO 23 I=1,N
      IF(I-IP) 24, 23, 25
24    TT=A(I,IP)
      A(I,IP)=COS*TT+SIN*A(I,JP)
      IF(P(I).GE.DABS(A(I,IP))) GO TO 26
      P(I)=DABS(A(I,IP))
      MM(I)=IP
26    A(I,JP)=-SIN*TT+COS*A(I,JP)
      IF(P(I).GE.DABS(A(I,JP))) GO TO 23
30    P(I)=DABS(A(I,JP))
      MM(I)=JP
      GO TO 23
25    IF(I.LT.JP) GO TO 27
      IF(I.GT.JP) GO TO 28
      IF(I.EQ.JP) GO TO 23
27    TT=A(IP,I)
      A(IP,I)=COS*TT+SIN*A(I,JP)
      IF(P(IP).GE.DABS(A(IP,I))) GO TO 29
      P(IP)=DABS(A(IP,I))
      MM(IP)=I
29    A(I,JP)=-TT*SIN+COS*A(I,JP)
      IF(P(I).GE.DABS(A(I,JP))) GO TO 23
      GO TO 30
28    TT=A(IP,I)
      A(IP,I)=TT*COS+SIN*A(JP,I)
      IF(P(IP).GE.DABS(A(IP,I))) GO TO 31
      P(IP)=DABS(A(IP,I))
      MM(IP)=I
31    A(JP,I)=-TT*SIN+COS*A(JP,I)
      IF(P(JP).GE.DABS(A(JP,I))) GO TO 23
      P(JP)=DABS(A(JP,I))
      MM(JP)=I
23    CONTINUE
      IF(NN.EQ.0) GOTO 7
      DO 32 I=1,N
      TT=V(I,IP)
      V(I,IP)=TT*COS+SIN*V(I,JP)
      V(I,JP)=-TT*SIN+COS*V(I,JP)
32    CONTINUE
      GO TO 7
12    RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE RSQR(NV,M,XO,YO,A,R2)
C     COMPUTING RSQUARE (BETWEEN OBSERVED & EXPECTED Y) GOODNESS OF FIT
      DOUBLE PRECISION YO(200),XO(200,10),A(10),R2,V1,V2,V3,C,S1,S2
      DOUBLE PRECISION YH,YSH,YHS,YHSS
C     -----------------------------------------------------------------
      N=NV
C     -----------------------------------------------------------------
      V1=0.0
      V2=0.0
      C=0.0
      S1=0.0
      S2=0.0
      DO I=1,N
      YH=0.0
      DO J=1,M
      YH=YH+XO(I,J)*A(J)
      ENDDO
      S1=S1+YO(I)
      V1=V1+YO(I)**2
      S2=S2+YH
      V2=V2+YH**2
      C=C+YO(I)*YH
      ENDDO
      V1=V1/N-(S1/N)**2    ! VARIANCE OF Y
      V2=V2/N-(S2/N)**2    ! VARIANCE OF YHAT
      V3=C/N-(S1/N)*(S2/N) ! COVARIANCE BETWEEN Y AND YHAT
      R2=V3**2/(V1*V2)     ! RSQUARE
      RETURN
      END