C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************
C
C  THIS DRIVER PROGRAM SHOWS HOW TO USE THE NON-GAUSSIAN ARMA
C  COEFFICIENT ESTIMATION PROCESS
C
C***********************************************************************
C        NGMLE-----------IS THE NON-GAUSSIAN MLE SUBROUTINE CALL WHICH
C                        MUST BE COMPILED WITH THE FOLLOWING LIBRARIES
C        POPFNC.FOR------THIS PROVIDES THE MLE CALCULATIONS SUPPLIED TO
C                        THE SEARCH LIBRARY
C        SEARCH.FOR------USES A DETERMINISTIC GRID SEARCH FOR SMALL 
C                        PARAMETER SPACES OR STOCHASTIC TECHNIQUES FOR
C                        LARGE PARAMETER SPACES
C        COVAR.FOR-------CALCULATES THE COVARIANCE MATRIX OF THE 
C                        PARAMETER ESTIMATES
C        PUBLIC.FOR------PUBLIC/INTERNET AVAILABLE ROUTINES FOR RANDOM
C                        NUMBER GENERATION, MATRIX INVERSION,
C                        AND POLYNOMIAL ROOT FINDING ALGORITHMS
C        DNSITY.FOR------STATISTICAL DISTRIBUTION OF THE NOISE PROCESS
C***********************************************************************
        INTEGER MAXX,MAXETA,NPOP,MAXL
C  MAXX,MAXETA ARE THE DIMENSIONS OF IMPORTANT ARRAYS DESCRIBED LATER
C  MUST BE GREATER THAN OR EQUAL TO NX,NP
        PARAMETER(MAXX=2000,MAXETA=6)
C  NPOP,MAXL ARE CONSTANTS USED IN THE STOCHASTIC SEARCH ROUTINES
        PARAMETER(NPOP=100,MAXL=100)
C  DIMENSIONING/TYPING OF ALL VARIABLES USED IN THE NGMLE ROUTINE
        REAL X(MAXX),Z(MAXX),ETA(MAXETA),ETALIM(2,MAXETA),STDEV(MAXETA)
        REAL COVMAT(MAXETA,MAXETA),COVVEC(MAXETA*(MAXETA+1)/2)
        REAL RSPC(24*(MAXETA)+(MAXETA+1)*NPOP+
     1     2*MAXX+10*(MAXL+1)+10*MAXETA*MAXETA)
        COMPLEX CSPC(4*MAXETA)
        INTEGER ISPC(4*MAXETA+2*NPOP)
        INTEGER P,Q,SUPSAG,CONV2M,NX,I,NDIV,NREP,NZ,NNZC,J,K
        REAL C1,C2,XMLE

C  OPEN FILES TO RETRIEVE DATA AND STORE DATA
C  INPUT
C        X23.DAT---------THE INPUT SEQUENCE DATA
C  OUTPUT
C        MLE.X23C-------THE MLE COEFFICIENT DATA WITH ALL INITIAL
C                        CONDITIONS
C        Z.X23C---------MOVING AVERAGE DATA

        OPEN(UNIT=7,FILE='X23.DAT',STATUS='UNKNOWN')
        OPEN(UNIT=8,FILE='MLE.X23C',STATUS='UNKNOWN')
        OPEN(UNIT=9,FILE='Z.X23C',STATUS='UNKNOWN')

C  INITIALIZE THE DATA ASSOCIATED WITH THE MLE FUNCTION USING THE
C  FOLLOWING INFORMATION
C
C  INPUT TO THE MLE ROUTINE
C        NX------LENGTH OF DATA TO ANALYZE                     (INTEGER)
C        X-------THE TIME SERIES DATA SEQUENCE                    (REAL)
C        P-------FOR AN ARMA(P,Q) THE P ASSOCIATED WITH THIS   (INTEGER)
C        Q-------FOR AN ARMA(P,Q) THE Q ASSOCIATED WITH THIS   (INTEGER)
C        ETALIM--THE MIN AND MAX VALUES OF THE ETA'S              (REAL)
C                MATRIX WHERE FIRST ROW IS THE MINIMUM AND
C                SECOND ROW IS THE MAXIMUM. I.E.
C                ETALIM(1,J)---MINIMUM JTH COEFFICIENT VALUE
C                ETALIM(2,J)---MAXIMUM JTH COEFFICIENT VALUE
C                1 <= J <= (P+Q+1)
C                WHEN J=(P+Q+1) ==> STANDARD DEVIATION OF THE NOISE
C        C1------DISTRIBUTIONAL CONSTANT NUMBER 1                 (REAL)
C        C2------DISTRIBUTIONAL CONSTANT NUMBER 2                 (REAL)
C        NDIV----GRID DENSITY FOR SEARCH (NUMBER OF DIVISIONS) (INTEGER)
C        NREP----NUMBER OF GRID REDUCTIONS (NUMBER OF REPETITIONS) (INT)
C        SUPSAG--SUPPRESS THE SIMULATED ANNEALING AND GENETIC ALGORITHM
C                SUPSAG =1 FOR THE SUPPRESSION                 (INTEGER)
C        RSPC----REAL WORK SPACE VARIABLE WHOSE MIN DIMENSION IS  (REAL)
C                24*(NP)+(NP+1)*NPOP+2*NX+10*(MAXL+1)+10*NP*NP)
C                WHERE NP=P+Q+1, MAXL=100, NPOP=100
C        CSPC----COMPLEX WORK SPACE VARIABLE WHOSE MIN DIMENSION  (REAL)
C                IS 4*NP
C        ISPC----INTEGER WORK SPACE WHOSE MIN DIMENSION IS     (INTEGER)
C                4*(NP)+2*NPOP
C  OUTPUT OF THE MLE ROUTINE
C        Z-------MOVING AVERAGE SEQUENCE                          (REAL)
C        NZ------LENGTH OF Z SEQUENCE                          (INTEGER)
C        XMLE----THE MLE FUNCTION                                 (REAL)
C        ETA-----THE COEFFICIENTS OF THE MLE PROCESS          (REAL VEC)
C        COVVEC--THE COVARIANCE VECTOR ASSOCIATED WITH ETA    (REAL VEC)
C        STDEV---THE STANDARD DEVIATION OF THE COEFFICIENTS OF ETA(RVEC)
C        NNZC----THE NUMBER OF NON-ZERO COEFFICIENTS           (INTEGER)

C  NOW DEFINE THE INPUTS TO THE NGMLE ROUTINE FOR A QUICK LOOK
C  AT AN ARMA(2,3)

        NX=2000
        READ(7,*) (X(I),I=1,NX)

        P=2
        Q=3

       ETALIM(1,1)=0.5
        ETALIM(2,1)=1.5
        ETALIM(1,2)=0.25
        ETALIM(2,2)=0.65
        ETALIM(1,3)=0.25
        ETALIM(2,3)=0.75
        ETALIM(1,4)=-1.5
        ETALIM(2,4)=-1.00
        ETALIM(1,5)=0.25
        ETALIM(2,5)=0.50
        ETALIM(1,6)=0.5
        ETALIM(2,6)=2.0

        C1=2.548
        C2=9

        NDIV=100
        NREP=1
        SUPSAG=0

C***********************************************************************
C  AFTER DEFINING THE INPUTS TO THIS ROUTINE NOW CALL AND FIND THE
C  SIGNIFICANT PARAMETERS

       CALL NGMLE(NX,X,P,Q,ETALIM,C1,C2,NDIV,NREP,SUPSAG,RSPC,CSPC,ISPC,
     1            Z,NZ,XMLE,ETA,COVVEC,STDEV,NNZC)

C***********************************************************************
C  SAVE THE OUTPUT SEQUENCE IF DESIRED FOR FURTHER ANALYSIS
        DO 10 I=1,NZ
                WRITE(9,100) Z(I)
10      CONTINUE

C  PRINT THE INPUTS TO THE NON-GAUSSIAN ROUTINE
        WRITE(8,120)
        WRITE(8,130) NX,P,Q,C1,C2
        IF(SUPSAG.EQ.1) THEN
                WRITE(8,140)
        ELSE
                WRITE(8,150)
                WRITE(8,160)
        ENDIF
        WRITE(8,170) NDIV
        WRITE(8,180) NREP
        WRITE(8,100)
        WRITE(8,190)
        DO 20 I=1,P+Q+1
                WRITE(8,110) ETALIM(1,I),ETALIM(2,I)
20      CONTINUE

C  PRINT THE OUTPUTS FROM THIS ROUTINE (FIRST THE MLE)
        WRITE(8,100)
        WRITE(8,200)
        WRITE(8,110) XMLE

C  PRINT THE OUTPUTS FROM THIS ROUTINE (SECOND THE COVARIANCE MATRIX)
C  THE FORM OF WHICH IS IN LOWER TRIANGULAR SINCE IT IS A POSITIVE
C  DEFINITE SYMMETRIC MATRIX.  CONVERT THE VECTOR TO A MATRIX
        WRITE(8,100)
        WRITE(8,210)
        CONV2M=0
        DO 40 I=1,NNZC
                DO 30 J=I,NNZC
                        CONV2M=CONV2M+1
                        COVMAT(I,J)=COVVEC(CONV2M)
                        COVMAT(J,I)=COVMAT(I,J)
30              CONTINUE
                WRITE(8,100) (COVMAT(I,K),K=1,NNZC)
40      CONTINUE
        WRITE(8,100)

C  PRINT THE OUTPUTS FROM THIS ROUTINE (LAST THE ESTIMATION ACCURACY)
        WRITE(8,220)
        DO 50 I=1,NNZC
                WRITE(8,110) ETA(I),STDEV(I)
50      CONTINUE

100     FORMAT(6F13.6)
110     FORMAT(2F13.6)
120     FORMAT('INPUTS TO THE NON-GAUSSIAN MLE SUBROUTINE')
130     FORMAT('NX=',I6,'  P=',I3,'  Q=',I3,'  C1=',E13.6,'  C2=',E13.6)
140     FORMAT('SUPPRESS THE STOCHASTIC ROUTINES')
150     FORMAT('LET PROGRAM DETERMINE THE OPTIMAL')
160     FORMAT('METHOD TO SEARCH FOR THE SOLUTION')
170     FORMAT('NUMBER OF GRID DIVISIONS=',I4)
180     FORMAT('NUMBER OF REPETITIONS TO NARROW SEARCH=',I3)
190     FORMAT('COEFFICIENT/PARAMETER LIMITS')
200     FORMAT('THE MLE(LOG) FOR THE FOLLOWING OUTPUT')
210     FORMAT('THE COVARIANCE MATRIX')
220     FORMAT('THE PARAMETERS AND THEIR STANDARD DEVIATIONS')

        CLOSE(7)
        CLOSE(8)
        CLOSE(9)
        STOP
        END
C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
