PROGRAM MULTICOL C DIMENSIONS ARE SPECIFIED DATA ARE READ AND PARAMETERS ARE SET C IN THE SUBROUTINE GMELE - THE LAST SUBROUTINE OF THIS PROGRAM C ------------------------------------------------------------------ C PROGRAM MULTICOL IS FOR RESOLVING THE PROBLEM OF MULTICOLLINRARITY C THIS PROGRAM ESTIMATES THE PARAMETERS OF A LINEAR (MULTIPLE) C REGRESSION MODEL (Y - XA + U) WHEN DATA (X) ARE MULTICOLLINEAR C PERFACT MULTICOLLINEARITY PROBLEMS ALSO ARE HANDLED. OPTIMIZATION C IS DONE BY THE DIFFERENTIAL EVOLUTION METHOD OF GLOBAL OPTIMATION. C ----------------------------------------------------------------- C THE "DIFFERENTIAL EVOLUTION ALGORITHM" OF GLOBAL OPTIMIZATION C WAS PROPOSED BY R. STORN AND KENNETH. V. PRICE IN 1995. REFERENCE: C "DIFFERENTIAL EVOLUTION - A SIMPLE AND EFFICIENT ADAPTIVE SCHEME C FOR GLOBAL OPTIMIZATION OVER CONTINUOUS SPACES" : TECHNICAL REPORT C INTERNATIONAL COMPUTER SCIENCE INSTITUTE, BERKLEY, 1995. C ----------------------------------------------------------------- C PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA) C ----------------------------------------------------------------- C DIFFERENTIAL EVOLUTION METAHEURISTIC FOR GLOBAL OPTIMIZATION IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! TYPE DECLARATION PARAMETER(NMAX=500,MMAX=50) ! MAXIMUM DIMENSION PARAMETERS PARAMETER (RX1=0.0, RX2=0.5) ! TO BE ADJUSTED SUITABLY, IF NEEDED C RX1 AND RX2 CONTROL THE SCHEME OF CROSSOVER. (0 <= RX1 <= RX2) <=1 C RX1 DETERMINES THE UPPER LIMIT OF SCHEME 1 (AND LOWER LIMIT OF C SCHEME 2; RX2 IS THE UPPER LIMIT OF SCHEME 2 AND LOWER LIMIT OF C SCHEME 3. THUS RX1 = .2 AND RX2 = .8 MEANS 0-20% SCHEME1, 20 TO 80 C PERCENT SCHEME 2 AND THE REST (80 TO 100 %) SCHEME 3. C ------------------------------ NOTE ----------------------------- C [RX1=0,RX2=1] (PURE EXPONENTIAL CROSSOVER) IS BEST IN MOST CASES C ------------------------------------------------------------------ C NOTE:(NCROSS=2) ! CROSS-OVER SCHEME (NCROSS <=0 OR =1 OR =>2) PARAMETER(IPRINT=500,EPS=1.D-08)!FOR WATCHING INTERMEDIATE RESULTS C IT PRINTS THE INTERMEDIATE RESULTS AFTER EACH IPRINT ITERATION AND C EPS DETERMINES ACCURACY FOR TERMINATION. IF EPS= 0, ALL ITERATIONS C WOULD BE UNDERGONE EVEN IF NO IMPROVEMENT IN RESULTS IS THERE. C ----------------------------------------------------------------- COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED) INTEGER IU,IV ! FOR RANDOM NUMBER GENERATION COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS * TITLE CHARACTER *70 FTIT ! TITLE OF THE FUNCTION COMMON /RIDGER/YY,Z,DSUM,DDSUM,QSUM,NVISIT DIMENSION YY(1000),Z(1000,20) C ----------------------------------------------------------------- C THE PROGRAM REQUIRES INPUTS FROM THE USER ON THE FOLLOWING ------ C (1) FUNCTION CODE (KF), (2) NO. OF VARIABLES IN THE FUNCTION (M); C (3) N=POPULATION SIZE (SUGGESTED 10 TIMES OF NO. OF VARIABLES, M, C FOR SMALLER PROBLEMS N=100 WORKS VERY WELL); C (4) PCROS = PROB. OF CROSS-OVER (SUGGESTED : ABOUT 0.85 TO .99); C (5) FACT = SCALE (SUGGESTED 0.5 TO .95 OR 1, ETC); C (6) ITER = MAXIMUM NUMBER OF ITERATIONS PERMITTED (5000 OR MORE) C (7) RANDOM NUMBER SEED (4 DIGITS INTEGER) C ---------------------------------------------------------------- DIMENSION X(NMAX,MMAX),Y(NMAX,MMAX),A(MMAX),FV(NMAX),IR(3) NVISIT=0 !SET TO 0 MEANS THAT DATA Y, X ARE NOT YED READ FROM FILE C ---------------------------------------------------------------- C ------- SELECT THE FUNCTION TO MINIMIZE AND ITS DIMENSION ------- CALL FSELECT(KF,M,FTIT) C SPECIFY OTHER PARAMETERS --------------------------------------- C WRITE(*,*)'POPULATION SIZE AND NO. OF ITERATIONS [ITER] ?' C WRITE(*,*)'SUGGESTED: POPULATION SIZE AT LEAST 10 TIMES M (NO. OF C & VARIABLES)' C WRITE(*,*) 'AND ITER 10000 OR LARGER' C READ(*,*) NP,ITER C ---------------------------------------------------------------- NP=M*10 ! YOU MAY RESET IT IF REQUIRED ITER=100000 ! YOU MAY RESET IT IF REQUIRED C ---------------------------------------------------------------- C WRITE(*,*)'CROSSOVER PROBABILITY [PCROS] AND SCALE [FACT] ?' C WRITE(*,*)'SUGGESTED : PCROS ABOUT 0.9; FACT=.5 OR LARGER BUT <=1' C READ(*,*) PCROS,FACT C ------------------------------------------------------------------ PCROS=0.9D0 ! YOU MAY RESET IT IF REQUIRED FACT=0.5D0 ! YOU MAY RESET IT IF REQUIRED C ------------------------------------------------------------------ WRITE(*,*)'RANDOM NUMBER SEED ?' WRITE(*,*)'A FOUR-DIGIT POSITIVE ODD INTEGER, SAY, 1171' READ(*,*) IU !SEED OF RANDOM NUMBER (4-DIGIT ODD NATURAL NUMBER) NFCALL=0 ! INITIALIZE COUNTER FOR FUNCTION CALLS GBEST=1.D30 ! TO BE USED FOR TERMINATION CRITERION C INITIALIZATION : GENERATE X(N,M) RANDOMLY DO I=1,NP DO J=1,M CALL RANDOM(RAND) ! GENERATES INITION X WITHIN X(I,J)=(RAND-.5D00)*2000 ! GENERATES INITION X WITHIN C RANDOM NUMBERS BETWEEN -1000 AND 1000 (BOTH EXCLUSIVE) ENDDO ENDDO IPCOUNT=0 DO 100 ITR=1,ITER ! ITERATION BEGINS C ----------------------------------------------------------------- C EVALUATE ALL X FOR THE GIVEN FUNCTION DO I=1,NP DO J=1,M A(J)=X(I,J) ENDDO CALL FUNC(A,M,F) C STORE FUNCTION VALUES IN FV VECTOR FV(I)=F ENDDO C ---------------------------------------------------------------- C FIND THE FITTEST (BEST) INDIVIDUAL AT THIS ITERATION FBEST=FV(1) KB=1 DO IB=2,NP IF(FV(IB).LT.FBEST) THEN FBEST=FV(IB) KB=IB ENDIF ENDDO C BEST FITNESS VALUE = FBEST : INDIVIDUAL X(KB) C ----------------------------------------------------------------- C GENERATE OFFSPRINGS DO I=1,NP ! I LOOP BEGINS C INITIALIZE CHILDREN IDENTICAL TO PARENTS; THEY WILL CHANGE LATER DO J=1,M Y(I,J)=X(I,J) ENDDO C SELECT RANDOMLY THREE OTHER INDIVIDUALS 20 DO IRI=1,3 ! IRI LOOP BEGINS IR(IRI)=0 CALL RANDOM(RAND) IRJ=INT(RAND*NP)+1 C CHECK THAT THESE THREE INDIVIDUALS ARE DISTICT AND OTHER THAN I IF(IRI.EQ.1.AND.IRJ.NE.I) THEN IR(IRI)=IRJ ENDIF IF(IRI.EQ.2.AND.IRJ.NE.I.AND.IRJ.NE.IR(1)) THEN IR(IRI)=IRJ ENDIF IF(IRI.EQ.3.AND.IRJ.NE.I.AND.IRJ.NE.IR(1).AND.IRJ.NE.IR(2)) THEN IR(IRI)=IRJ ENDIF ENDDO ! IRI LOOP ENDS C CHECK IF ALL THE THREE IR ARE POSITIVE (INTEGERS) DO IX=1,3 IF(IR(IX).LE.0) THEN GOTO 20 ! IF NOT THEN REGENERATE ENDIF ENDDO C THREE RANDOMLY CHOSEN INDIVIDUALS DIFFERENT FROM I AND DIFFERENT C FROM EACH OTHER ARE IR(1),IR(2) AND IR(3) C ===================== RANDOMIZATION OF NCROSS =================== C RANDOMIZES NCROSS NCROSS=0 CALL RANDOM(RAND) IF(RAND.GT.RX1) NCROSS=1 ! IF RX1=>1, SCHEME 2 NEVER IMPLEMENTED IF(RAND.GT.RX2) NCROSS=2 ! IF RX2=>1, SCHEME 3 NEVER IMPLEMENTED C ---------------------- SCHEME 1 ---------------------------------- C NO CROSS OVER, ONLY REPLACEMENT THAT IS PROBABILISTIC IF(NCROSS.LE.0) THEN DO J=1,M ! J LOOP BEGINS CALL RANDOM(RAND) IF(RAND.LE.PCROS) THEN ! REPLACE IF RAND < PCROS A(J)=X(IR(1),J)+(X(IR(2),J)-X(IR(3),J))*FACT ! CANDIDATE CHILD ENDIF ENDDO ! J LOOP ENDS ENDIF C ----------------------- SCHEME 2 --------------------------------- C THE STANDARD CROSSOVER SCHEME C CROSSOVER SCHEME (EXPONENTIAL) SUGGESTED BY KENNETH PRICE IN HIS C PERSONAL LETTER TO THE AUTHOR (DATED SEPTEMBER 29, 2006) IF(NCROSS.EQ.1) THEN CALL RANDOM(RAND) 1 JR=INT(RAND*M)+1 J=JR 2 A(J)=X(IR(1),J)+FACT*(X(IR(2),J)-X(IR(3),J)) 3 J=J+1 IF(J.GT.M) J=1 4 IF(J.EQ.JR) GOTO 10 5 CALL RANDOM(RAND) IF(PCROS.LE.RAND) GOTO 2 6 A(J)=X(I,J) 7 J=J+1 IF(J.GT.M) J=1 8 IF (J.EQ.JR) GOTO 10 9 GOTO 6 10 CONTINUE ENDIF C ------------------------ SCHEME 3 -------------------------------- C ESPECIALLY SUITABLE TO NON-DECOMPOSABLE (NON-SEPERABLE) FUNCTIONS C CROSSOVER SCHEME (NEW) SUGGESTED BY KENNETH PRICE IN HIS C PERSONAL LETTER TO THE AUTHOR (DATED OCTOBER 18, 2006) IF(NCROSS.GE.2) THEN CALL RANDOM(RAND) IF(RAND.LE.PCROS) THEN CALL NORMAL(RN) DO J=1,M A(J)=X(I,J)+(X(IR(1),J)+ X(IR(2),J)-2*X(I,J))*RN ENDDO ELSE DO J=1,M A(J)=X(I,J)+(X(IR(1),J)- X(IR(2),J))! FACT ASSUMED TO BE 1 ENDDO ENDIF ENDIF C ----------------------------------------------------------------- CALL FUNC(A,M,F) ! EVALUATE THE OFFSPRING IF(F.LT.FV(I)) THEN ! IF BETTER, REPLACE PARENTS BY THE CHILD FV(I)=F DO J=1,M Y(I,J)=A(J) ENDDO ENDIF ENDDO ! I LOOP ENDS DO I=1,NP DO J=1,M X(I,J)=Y(I,J) ! NEW GENERATION IS A MIX OF BETTER PARENTS AND C BETTER CHILDREN ENDDO ENDDO IPCOUNT=IPCOUNT+1 IF(IPCOUNT.EQ.IPRINT) THEN DO J=1,M A(J)=X(KB,J) ENDDO WRITE(*,*)' ' WRITE(*,*)(X(KB,J),J=1,M),' FBEST UPTO NOW = ',FBEST C WRITE(*,*)'TOTAL NUMBER OF FUNCTION CALLS =',NFCALL IF(DABS(FBEST-GBEST).LT.EPS) THEN WRITE(*,*) FTIT WRITE(*,*)'NO. OF VARIABLES =', M WRITE(*,*)' ' WRITE(*,*)'ESTIMATED PARAMETERS ARE AS FOLLOWS -----------------' WRITE(*,*)(X(KB,J),J=1,M) WRITE(*,*)'-----------------------------------------------------' C WRITE(*,*)'FINAL (TOTAL) PENALTY = ',FBEST IF(KF.EQ.1) THEN DO J=1,M A(J)=X(KB,J) ENDDO CALL FUNC(A,M,F) WRITE(*,*)'ESTIMATED PARAMETERS ARE AS FOLLOWS -----------------' WRITE(*,*)(A(J),J=1,M) WRITE(*,*)'FINAL (TOTAL) PENALTY = ',F WRITE(*,*)'PENALTY AS SUM ABS(Y-YHAT)**P = ',DSUM WRITE(*,*)'PENALTY AS PARMATER ENTROPY = ',DDSUM WRITE(*,*)'PENALTY AS RESIDUAL ENTROPY = ',QSUM WRITE(*,*)'COMPUTATION OVER. THANK YOU' ENDIF STOP ELSE GBEST=FBEST ENDIF IPCOUNT=0 ENDIF C ---------------------------------------------------------------- 100 ENDDO ! ITERATION ENDS : GO FOR NEXT ITERATION, IF APPLICABLE C ---------------------------------------------------------------- WRITE(*,*)'DID NOT CONVERGE. REDUCE EPS OR RAISE ITER OR DO BOTH' WRITE(*,*)'INCREASE N, PCROS, OR SCALE FACTOR (FACT)' END C ----------------------------------------------------------------- SUBROUTINE NORMAL(R) C PROGRAM TO GENERATE N(0,1) FROM RECTANGULAR RANDOM NUMBERS C IT USES BOX-MULLER VARIATE TRANSFORMATION FOR THIS PURPOSE. C ----------------------------------------------------------------- C ----- BOX-MULLER METHOD BY GEP BOX AND ME MULLER (1958) --------- C BOX, G. E. P. AND MULLER, M. E. "A NOTE ON THE GENERATION OF C RANDOM NORMAL DEVIATES." ANN. MATH. STAT. 29, 610-611, 1958. C IF U1 AND U2 ARE UNIFORMLY DISTRIBUTED RANDOM NUMBERS (0,1), C THEN X=[(-2*LN(U1))**.5]*(COS(2*PI*U2) IS N(0,1) C ALSO, X=[(-2*LN(U1))**.5]*(SIN(2*PI*U2) IS N(0,1) C PI = 4*ARCTAN(1.0)= 3.1415926535897932384626433832795 C 2*PI = 6.283185307179586476925286766559 C ----------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /RNDM/IU,IV INTEGER IU,IV C ----------------------------------------------------------------- CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1] U1=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1] CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1] U2=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1] R=DSQRT(-2.D0*DLOG(U1)) R=R*DCOS(U2*6.283185307179586476925286766559D00) RETURN END C ----------------------------------------------------------------- C RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE) SUBROUTINE RANDOM(RAND) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /RNDM/IU,IV INTEGER IU,IV IV=IU*65539 IF(IV.LT.0) THEN IV=IV+2147483647+1 ENDIF RAND=DBLE(IV+.0D0) IU=IV RAND=RAND*0.4656613E-09 RETURN END C ----------------------------------------------------------------- SUBROUTINE FSELECT(KF,M,FTIT) ! FUNCTION SELECTION C THE PROGRAM REQUIRES INPUTS FROM THE USER ON THE FOLLOWING ------ C (1) FUNCTION CODE (KF), (2) NO. OF VARIABLES IN THE FUNCTION (M); CHARACTER *70 TIT(10),FTIT WRITE(*,*)'----------------------------------------------------' DATA TIT(1)/'KF=1 GMELE REGRESSION IN N OBSERVATIONS M-VARIABLES'/ C ----------------------------------------------------------------- I=1 KF=1 !FUNCTION CODE WRITE(*,*)TIT(I) WRITE(*,*)'----------------------------------------------------' WRITE(*,*)'NO.OF VARIABLES/PARAMETERS IN REGRESSION MODEL [M] ?' READ(*,*) M FTIT=TIT(KF) ! STORE THE NAME OF THE CHOSEN FUNCTION IN FTIT RETURN END C ----------------------------------------------------------------- SUBROUTINE FUNC(X,M,F) C TEST FUNCTIONS FOR GLOBAL OPTIMIZATION PROGRAM IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /RNDM/IU,IV COMMON /KFF/KF,NFCALL,FTIT INTEGER IU,IV DIMENSION X(*) CHARACTER *70 FTIT NFCALL=NFCALL+1 ! INCREMENT TO NUMBER OF FUNCTION CALLS C KF IS THE CODE OF THE TEST FUNCTION C ------------------------------------------------------------------ IF(KF.EQ.1) THEN C GMELE METHOD TO RERSOLVE MULTICOLLINEARITY PROBLEM IN REGRESSION CALL GMELE(M,X,F) ! CALLS THE REGRESSION FUNCTION RETURN ENDIF C ----------------------------------------------------------------- C IF(KF.EQ.99) THEN C ***** FUNCTION C CALL SUBROUTINE(M,X,F) C RETURN C ENDIF C ================================================================= WRITE(*,*)'FUNCTION NOT DEFINED. PROGRAM ABORTED' STOP END C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE GMELE(M,X,F) ! GMELE REGRESSION (REF. PARIS, MISHRA) C ----------------------------------------------------------------- PARAMETER(N=13,MM=5,NP=1,NPP=2,NQ=1,NQQ=1)!SET THESE FOR YOURSELF PARAMETER (NGMELS=1) ! 0 IS OLS AND NONZERO IS GMELE ESTIMATOR C ------------------ SETTING GUIDELINES ------------------------- C N=NO.OF OBSERVATIONS(CASES); MM=M=NO.OF VARIABLES (= PARAMETERS) C NP=NORM FOR OBTAINING PROBABILITY 1= MMELE(MISHRA), 2= MELE(PARIS) C NPP=NORM FOR LEAST SQUARES (NPP=2); FOR LEAST ABSOLUTES (NPP=1) C NQQ=0 RESIDUAL PROBABILITY ARE NOT USED (MELE1 OF PARIS/MISHRA) C NQQ <> 0 RESIDUAL PROBABILITY ARE USED (MELE2 OF PARIS/MISHRA) C NQ IS NORM USED FOR PROBABILTY OF RESIDUALS (2 SQUARE; 1 ABSOLUTE) C (NOTE : NQ APPLIES ONLY IF NQQ IS NOT ZERO). C ----------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /RIDGER/YY,Z,DSUM,DDSUM,QSUM,NVISIT DIMENSION X(*),YY(1000),YHAT(N),Z(1000,20),P(MM) COMMON /RNDM/IU,IV CHARACTER *50 FIL C EXAMPLE OF THE FUNCTION TO MINIMIZE C ---------------- READ DATA Y X FROM FILE ------------------------ IF(NVISIT.EQ.0) THEN WRITE(*,*)'FILE NAME THAT STORES YOUR [Y, X] DATA ?' READ(*,*) FIL OPEN(8,FILE=FIL) DO I=1,N ! READS N OBSERVATIONS (CASES) READ(8,*) YY(I),(Z(I,J),J=1,MM) WRITE(*,1) YY(I),(Z(I,J),J=1,MM) ENDDO 1 FORMAT(8F10.2) CLOSE(8) NVISIT=1 ! DATA WAS READ; NVISIT = 1 (DATA ALREADY RED) WRITE(*,*)'COMPUTING --- PLEASE WAIT ' ENDIF C ---------------------------------------------------------------- DSUM=0.D0 ! SUM OF ABS(Y-YHAT)**NORM DDSUM=0.D0 ! ENTROPY PENALTY FOR ESTIMATED PARAMETERS QSUM=0.D0 ! ENTROPY PENALTY FOR ESTIMATED RESIDUALS C ---------------------------------------------------------------- QX=0.D0 SL=0.D0 DO I=1,M SL=SL+DABS(X(I))**NP ENDDO DO I=1,M P(I)=DABS(X(I))**NP/SL ENDDO C SUM OF ABS(Y-F(Z))NPP (2 FOR LEAST SQUARES, 1 FOR LEAST ABSOLUTE) S=0.D0 DO I=1,N SS=0.D0 DO J=1,MM SS=SS+X(J)*Z(I,J) ENDDO YHAT(I)=DABS(YY(I)-SS) S=S+YHAT(I)**NPP ENDDO C --------------------- APPLIES ONLY IF NQQ NOT ZERO --------------- IF(NQQ.NE.0) THEN SYHAT=0.D0 DO I=1,N SYHAT=SYHAT+YHAT(I)**NQ ENDDO SQ=0.D0 DO I=1,N QL=YHAT(I)**NQ IF(QL.GT.0.D0) SQ=SQ+QL*DLOG(QL) ENDDO QX=SQ+SYHAT*DLOG(SYHAT) ENDIF C ----------------------------------------------------------------- SS=0.D0 DO J=1,MM IF(P(I).GT.0.D0) SS=SS+P(I)*DLOG(P(I)) ENDDO SS=SS+SL*DLOG(SL) DSUM=S IF(NGMELS.EQ.0) THEN F=S ! FOR OLS - NO PENALTY OTHEN THAN SUM(ABS(Y-YHAT))**P ELSE DDSUM=SS QSUM=QX F=S+SS+QX ! GMELE ENDIF RETURN END