C This code is written in Vax Fortran. The notation generally follows c the notation in my paper ``Maximum Likelihood Estimation of Stationary c Univariate Fractionally Integrate Time Series Models.'' Revised Dec. 1990 c c Fallaw Sowell c c C YOU WILL NEED THE IMSL ROUTINE "DZPLRC" AND THE FUNCTION "DGAMMA" C c c This program is written to take advantage of the numerical optimization c library called GQOPT. c c ``GQOPT is a general purpose numerical optimization package designed c to run on IBM and VAX equipment, although some versions have been c adapted for Burroughs and Control Data computers.'' Information c concerning GQOPT is available from: Richard E. Quandt c Department of Economics c Princeton University c Princeton, NJ 08544 c c c c c Copyright 1993 Fallaw Sowell c Permission to copy and distribute granted provided no fee is charged c and this copyright notice is attached. c IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER P,Q COMMON/BTRAT/ITRFLG COMMON/BMYNOBS/NOBS COMMON/BSTACK/AINT(500) COMMON/BSTAK/NQ,NTOP COMMON/BSTOP/NVAR1,ISTOP(3) COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE COMMON/VARDAT/VARDAT COMMON/PI/DPI COMMON/DATA/X(500) COMMON/PARMS/MYNOBS,P,Q COMMON/COVAR/SIG(500),V2(500),CHL(124750) DIMENSION A(10),B(10),THETA(20),aaa(10),bbb(10),tratio(20) DIMENSION ALABEL(20) EXTERNAL FRAC,DFP OPEN(UNIT=10,FILE='MAX.DAT',STATUS='OLD') OPEN(UNIT=7 ,FILE='MAX.TXT',STATUS='NEW') OPEN(UNIT=12,FILE='TAB.DAT',STATUS='OLD',ACCESS='APPEND') OPEN(UNIT=13,FILE='SELECT.DAT',STATUS='OLD',ACCESS='APPEND') DPI = 4.0D0*ATAN(1.0D0) C C READ IN THE DATA FROM THE FILE "MAX.DAT" C THE MEAN AND STANDARD DEVIATION ARE CALCULATED AND THE DATA C ARE CENTERED AND STANDARDIZED. THE FINAL VALUE OF THE C LIKELIHOOD REPORTED IS CORRECTED FOR THE STANDARDIZATION. C READ(10,*) MYNOBS DO 10 I = 1,MYNOBS READ(10,*) X(I) SUM = SUM + X(I) 10 SQ = SQ + X(I)*X(I) DMEAN = SUM/MYNOBS SCMOM = SQ/MYNOBS VARDAT = SCMOM - DMEAN*DMEAN SDDAT = SQRT(VARDAT) WRITE(7,*) 'MEAN OF THE DATA = ',DMEAN WRITE(7,*) 'VARIANCE OF THE DATA = ',VARDAT DO 821 I = 1,MYNOBS 821 X(I) = (X(I) - DMEAN)/SDDAT WRITE(7,*) 'MEAN REMOVED FROM DATA' do 1001 ij = 1,4 p = ij - 1 IF((P.EQ.0)) THEN VAR = 1.0 GOTO 89 ENDIF RELERR = 1.0d-10 MAXIT = 2**28 CALL DNSPE(myNOBS,X,1,0,0.0,P,0,RELERR,MAXIT,MU,AaA,BbB,aVAR) var = avar IF(P.EQ.0) GOTO 89 DO 94 I = 1,P 94 A(I+1) = -AaA(I) 89 THETA(1) = .01 IF (P.EQ.0) GOTO 411 DO 111 I = 1,P 111 THETA(1+I) = A(I+1) 411 CONTINUE theta(2+p) = VAR WRITE(7,*) ' START VALUES ARE ' WRITE(7,*) 'D =',D WRITE(7,*) 'VAR = ',VAR WRITE(7,*) 'P = ',P WRITE(7,*) 'Q = ',Q WRITE(7,*) 'AUTO REGRESSIVE COEFFICIENTS' DO 50 I = 1,P+1 50 WRITE(7,5) I,A(I) WRITE(7,*) 'MOVING AVERAGE COEFFICIENTS' DO 60 I = 1,Q+1 60 WRITE(7,6) I,B(I) 5 FORMAT(' A(',I2,') = ',F7.4) 6 FORMAT(' B(',I2,') = ',F7.4) do 1000 ji = 1,4 q = ji - 1 if(q.eq.0) goto 213 theta(p+q+2) = theta(p+q+1) 212 theta(p+q+1) = -.01 213 continue C C START THE MAXIMIZATION ROUTINE. C NP = 2+P+Q IPT=2 JPT=2 MFILE=7 ACC=1.0D-13 ITERL=50 MX=1 NQ=1000 ISTOP(1) = 1 ISTOP(2) = 1 ITRFLG = 1 NOBS = 0 CALL LABEL(ALABEL,NP) CALL OPT(THETA,NP,F,DFP,ITERL,MX,IER,ACC,FRAC,ALABEL) CALL OPTOUT(0) WRITE(7,*) ' ' WRITE(7,*) ' LOG LIKELIHOOD AT OPTIMUM = ',F WRITE(7,*) ' ' WRITE(7,*) ' MODEL SELECTION CRITERIA' D2F = 2*F AIC = 2*F - 2.0D0*(NP-1) SIC = 2*F - (NP-1)*LOG( REAL(MYNOBS) ) WRITE(7,*) ' 2*LOGLIKE ',' AIC ',' SIC ' WRITE(7,65) D2F,AIC,SIC WRITE(13,65) D2F,AIC,SIC 65 FORMAT(3(1X,F11.5)) WRITE(7,*) ' ' do 300 i = 1,np-1 300 tratio(I) = theta(i)/ sqrt( -aint(np*(i-1) + i) ) write(12,801) ier, ( theta(i),i=1,(np-1) ) write(12,802) ( tratio(i),i=1,(np-1) ) 801 format(1x,i3,1x,7(f6.2,' ')) 802 format(1x,7(' ',f6.2)) 1000 continue 1001 continue STOP END C C THIS IS A SUBROUTINE TO EVALUATE A HYPERGEOMETRIC FUNCTION OF C THE FORM F(A,1;C,Z). THE APPROACH IS TO SIMPLY SUM THE SERIES C UNTIL THE TERMS ARE BELOW 1.0D-14 OR UNTIL THE FIRST 10,000 TERMS C HAVE BE SUMMED. C C NOTE: THE HYPERGEOMETRIC FUNCTION CONVERGES VERY SLOWLY WHEN Z IS C NEAR THE UNIT CIRCLE. C SUBROUTINE DHYPER(A,C,Z,HYP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMPLEX*16 Z,HYP,INC ACC = 1.0D-14 HYP = (1.0D0,0.0D0) INC = (1.0D0,0.0D0) DO 10 M = 1,10000 INC = INC*Z*(A+M-1.0D0)/(C+M-1.0D0) HYP = HYP + INC W = ABS(INC) 10 IF (W.LT.ACC) GOTO 20 WRITE(7,*) 'ALL THROUGH THE HYPER LOOP !!! 10,000 ITERATIONS' WRITE(7,*) 'A,C,Z,HYP ARE :',A,C,Z,HYP WRITE(*,*) 'ALL THROUGH THE HYPER LOOP !!! 10,000 ITERATIONS' WRITE(*,*) 'A,C,Z,HYP ARE :',A,C,Z,HYP 20 RETURN END C C THIS SUBROUTINE CALCULATES THE RATIO OF TWO GAMMA FUNCTIONS OF THE FORM C GAMMA(D+K)/GAMMA(1-D+K). THE PROBLEM IS THAT GAMMA(X) IS MACHINE C INFINITE ON SYSTEMS TYPICALLY AT X=35 OR X=56 OR X=157. C NOTE GAMMA(D+K)/GAMMA(1-D+K) = GAMMA(D-K)/GAMMA(1-D-K). C SUBROUTINE RATIO(D,K,R) IMPLICIT DOUBLE PRECISION (A-H,O-Z) R = DGAMMA(D)/DGAMMA(1-D) J = ABS(K) DO 10 I = 1,J 10 R = R*(D+I-1.)/(1.-D+I-1.) RETURN END C C THIS SUBROUTINE CALCULATES "THE LOG OF THE DETERMINANT" AND THE C CHOLESKY DECOMPOSITION OF "THE INVERSE" OF A TOEPLITZ MATRIX. C c the log determinant is of the toeplitz matrix not the log c determinant of the inverse of the toeplitz matrix c SUBROUTINE LEVIN(DLGDET,ERR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/COVAR/R(500),V2(500),A(124750) COMMON/PARMS/MYNOBS,P,Q ERR = 0 V2(1) = R(1) IF (V2(1).LE.(0.0)) THEN WRITE(*,*) 'ERROR IN LEVIN, R(1) IS NOT POSITIVE' WRITE(7,*) 'ERROR IN LEVIN, R(1) IS NOT POSITIVE' ERR = 1 GOTO 50 ENDIF A(1) = R(2)/R(1) V2(2) = V2(1)*(1. - A(1)**2) IF (V2(2).LE.0.0) THEN WRITE(*,*) 'ERROR IN LEVIN, V2(2) IS NOT POSITIVE' WRITE(7,*) 'ERROR IN LEVIN, V2(2) IS NOT POSITIVE' ERR = 2 GOTO 50 ENDIF DO 30 I = 2,MYNOBS-1 TEMP = 0.0D0 DO 10 K = 1,I-1 10 TEMP = TEMP + A(((I-1)*(I-2)/2)+K)*R(I-K+1) A((I*(I+1)/2)) = (R(I+1) - TEMP)/V2(I) DO 20 J = 1,I-1 20 A((I*(I-1)/2)+J) = A(((I-1)*(I-2)/2)+J) - & A((I*(I+1)/2))*A((I*(I-1)/2)+1-J) V2(I+1) = V2(I)*(1. - A((I*(I+1)/2))**2) IF (V2(I+1).LE.0.0) THEN WRITE(*,*) 'ERROR IN LEVIN, COVARIANCE NOT POSITIVE DEFINITE' WRITE(7,*) 'ERROR IN LEVIN, COVARIANCE NOT POSITIVE DEFINITE' ERR = 3 GOTO 50 ENDIF 30 CONTINUE DLGDET = 0.0D0 DO 40 I = 1,MYNOBS 40 DLGDET = DLGDET + LOG(V2(I)) 50 RETURN END C C SUBROUTINE TO EVALUATE THE LIKELIHOOD FUNCTION FOR A FRACTIONALLY C INTEGRATED TIME SERIES MODEL. ARIMA(P,D,Q) D<.499 C SUBROUTINE FRAC(THETA,NP,F,*) C C INPUTS TO THE ROUTINE ARE C D = THE DIFFERENCING PARAMETER ( D < .5, BUT NOT 0) C VAR = THE VARIANCE OF THE ERROR PROCESS. (VAR>0) C P = THE ORDER OF THE AUTOREGRESSIVE (AR) POLYNOMIAL (P<10) C Q = THE ORDER OF THE MOVING AVERAGE (MA) POLYNOMIAL (P<10) C A = ARRAY OF DIMENSION 10 WHERE THE FIRST P+1 ELEMENTS ARE C THE COEFFICIENTS OF THE AR POLYNOMIAL, (A(1) = 1) C A(1) + A(2)*Z + A(3)*Z**2 + A(4)*Z**3 + ... + A(P+1)*Z**P C THE EVALUATION ASSUMES THE PROCESS IS STATIONARY HENCE C THE ROOTS OF THE AR POLYNOMIAL MUST BE OUTSIDE THE UNIT C CIRCLE. C B = ARRAY OF DIMENSION 10 WHERE THE FIRST Q+1 ELEMENTS ARE C THE COEFFICIENTS OF THE MA POLYNOMIAL, (B(1) = 1) C B(1) + B(2)*Z + B(3)*Z**2 + B(4)*Z**3 + ... + B(Q+1)*Z**Q C C OUTPUTS FROM THE ROUTINE ARE C IERR = ERROR FLAG IF IERR = 0 AFTER THE ROUTINE THEN THE C EVALUATION WAS SUCCESSFUL OTHERWISE A C PROBLEM WAS ENCOUNTERED. THE VALUE OF C IER WILL IDENTIFY WHERE THE PROBLEM WAS. C 1 = D > .499 C 2 = VAR < .0001 C 3 = EITHER P>9 OR Q>9 C 4 = AR POLYNOMIAL HAS A ROOT INSIDE THE UNIT CIRCLE C ( ABS(RHO(I))>.99 ) C 5 = INVERSION OF THE COVARIANCE MATRIX FAILED. C F = VALUE OF THE LIKELIHOOD FUNCTION C C THE DATA IS PASSED TO THE SUBROUTINE THROUGH A COMMON STATEMENT C THE NUMBER OF OBSERVATIONS IS ALSO PASSED THROUGH COMMON STATEMENT C THIS PROGRAM IS WRITTEN TO NOT ALLOW ANY MORE THAN 500 OBSERVATIONS C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMPLEX*16 ZED(9),RHO(9),HYP(1060,9),RHO2P(9),SIG1 INTEGER S,P,Q COMMON/VARDAT/VARDAT COMMON/PI/DPI COMMON/DATA/X(500) COMMON/PARMS/MYNOBS,P,Q COMMON/COVAR/SIG(500),V2(500),CHL(124750) DIMENSION A(10),B(10),THETA(20),RAT(1040),C(530,9) C C MATCH THE THETA PARAMETERS TO THE PARAMETERS IN THE PAPER C D = THETA(1) A(1) = 1.0D0 IF (P.EQ.0) GOTO 2 DO 1 I = 2,P+1 1 A(I) = THETA(I) 2 B(1) = 1.0D0 IF (Q.EQ.0) GOTO 4 DO 3 I = 2,Q+1 3 B(I) = THETA(P+I) 4 VAR = THETA(1+P+Q+1) IERR = 0 C C CHECK TO SEE IF -5.0 < D < .499, VAR > .0001 AND P,Q<10. C IF ( (D.GE.(.499)) .OR. (D.LE.(-5.0)) ) THEN IERR = 1 GOTO 999 ENDIF IF (VAR.LE.(.0001)) THEN IERR = 2 GOTO 999 ENDIF IF ( (P.GE.10) .OR. (Q.GE.10)) THEN IERR = 3 GOTO 999 ENDIF IF (P.NE.0) GOTO 6 C C CODE TO CALCULATE THE COVARIANCE WHEN THERE ARE NO AR PARAMETERS C CONS = VAR*DGAMMA(1-2.0D0*D)/((DGAMMA(1.-D))*(DGAMMA(D))) ICHGVR = Q+MYNOBS+1 RAT(Q+ICHGVR) = DGAMMA(D+Q)/DGAMMA(1.-D+Q) DO 300 I = Q+ICHGVR-1,1,-1 300 RAT(I) = RAT(I+1)*(1.-D+I-ICHGVR)/(D+I-ICHGVR) DO 370 S = 1,MYNOBS SIG(S) = 0.0D0 DO 366 N = 1,Q+1 DO 365 M = 1,Q+1 365 SIG(S) = SIG(S) + B(N)*B(M)*RAT(N-M-(S-1)+ICHGVR) 366 CONTINUE 370 SIG(S) = CONS*SIG(S) GOTO 900 C C CODE TO CALCULATE THE COVARIANCE WHEN THERE ARE AR C C C CALCULATE THE ROOTS OF THE AR POLYNOMIAL USE IMSL ROUTINE "DZPLRC". C NOTE THAT THE SUBROUTINE CALCULATES THE ROOTS OF THE POLYNOMIAL, C BUT THE PROGRAM USES THE INVERSE OF THE ROOTS. C 6 CALL DZPLRC(P,A,RHO) C C CHECK TO SEE THAT NONE OF THE ROOTS OF THE AR POLYNOMIAL ARE TOO CLOSE C TO THE UNIT CIRCLE. C DO 10 I = 1,P RHO(I) = 1/RHO(I) W = ABS(RHO(I)) IF (W.GE.(.99)) THEN IERR = 4 GOTO 999 ENDIF 10 CONTINUE C C CALCULATE THE ZEDS C (CALCULATE 1/ZED THEN INVERT TO GET ZED) C DO 45 I = 1,P ZED(I) = RHO(I) DO 20 J = 1,P 20 ZED(I) = ZED(I)*( 1.0D0 - RHO(J)*RHO(I) ) IF (I.EQ.1) GOTO 35 DO 30 K = 1,I-1 30 ZED(I) = ZED(I)*(RHO(I) - RHO(K)) 35 IF (I.EQ.P) GOTO 43 DO 40 K = I+1,P 40 ZED(I) = ZED(I)*(RHO(I) - RHO(K)) 43 ZED(I) = 1.0D0/ZED(I) 45 CONTINUE C TO CALCULATE THE COVARIANCE TERMS WE NEED TO EVALUATE THE C FUNCTIONS ZED(I)*C(D,D,P+N-M-S,RHO(I)). IN THIS PROGRAM THESE VALUES C WILL BE DENOTED BY THE DOUBLE ARRAY C(*,*). TO EVALUATE C(*,*) WE C NEED TO CALCULATE HYPERGEOMETRIC FUNCTIONS F(D+H,1,1-D+H,RHO(I)). C THESE VALUES WILL BE DENOTED BY HYP(H-ICHGVR,I). ICHGVR IS A C CHANGE OF VARIABLES CONSTANT SO THAT THE INDICES ON THE ARRAYS ARE C ALWAYS POSITIVE. ICHGVR = P-Q-MYNOBS. C C FIRST WE NEED TO EVALUATE ALL THE HYPERGEOMETRIC C FUNCTIONS NEEDED. THIS WILL BE DONE RECURSIVELY C HENCE ONLY ONE HYPERGEOMETRIC FUNCTION NEEDS TO C BE EVALUATED FOR EACH ROOT. C ICHGVR = P-Q-MYNOBS DO 47 I = 1,P CALL & DHYPER(D-ICHGVR-1,1-D-ICHGVR-1,RHO(I),HYP((-1-2*ICHGVR),I)) J = (-1-2*ICHGVR) DO 46 J = -2-2*ICHGVR,1,-1 46 HYP(J,I) = DCMPLX(1.0D0) + & RHO(I)*HYP(J+1,I)*DCMPLX((D+J+ICHGVR)/(1.-D+J+ICHGVR)) 47 CONTINUE C C NOW CALCULATE THE C(*,*) FUNCTIONS C DO 49 I = 1,P 49 RHO2P(I) = RHO(I)**(2*P) COEF = DGAMMA(1. - 2*D)/(DGAMMA(1.-D) * DGAMMA(D)) CALL RATIO(D,ICHGVR,R) R = COEF*R DO 50 N = 1,2*Q+MYNOBS R = R*(D+ICHGVR+N-1.)/(1-D+ICHGVR+N-1.) DO 52 I = 1,P C(N,I) = ZED(I)*R*( RHO2P(I)*HYP(N,I) + & HYP(-N-2*ICHGVR,I) - 1.0D0 ) 52 CONTINUE 50 CONTINUE C C NOW CALCULATE THE AUTOCOVARIANCES C DO 70 S = 1,MYNOBS SIG(S) = 0.0D0 DO 68 I = 1,P SIG1 = 0.0D0 DO 66 N = 1,Q+1 DO 65 M = 1,Q+1 65 SIG1 = SIG1 + B(N)*B(M)*C(P+N-M-(S-1)-ICHGVR,I) 66 CONTINUE SIG(S) = SIG(S) + SIG1 68 CONTINUE 70 SIG(S) = VAR*SIG(S) C C THE INVERSE AND DETERMINANT OF THE COVARIANCE MATRIX IS CALCULATED C BY FIRST CALCULATING THE INVERSE AND DETERMINANT OF THE CORRELATION C MATRIX. C 900 SIGTMP = SIG(1) DO 80 S = 1,MYNOBS 80 SIG(S) = SIG(S)/SIGTMP C C LEVIN IS A SUBROUTINE THAT CALCULATES THE LOG OF THE DETERMINANT AND C THE CHOLESKY DECOMPOSITION OF THE INVERSE OF A TOEPLITZ MATRIX. C CALL LEVIN(DLGDET,IER) C C CHECK TO SEE IF THE LEVIN SUBROUTINE WAS SUCCESSFUL C IF (IER.NE.0) THEN IERR = 5 GOTO 999 ENDIF C C LETS EVALUATE THE LIKELIHOOD FUNCTION C DM = X(1)*X(1)/V2(1) + ( (X(2) - CHL(1)*X(1))**2)/V2(2) DO 90 I = 3,MYNOBS Z = X(I) DO 85 L = 1,I-1 85 Z = Z - CHL( (I*(I-1)/2)+1-L)*X(L) 90 DM = DM + Z*Z/V2(I) F=-(MYNOBS*LOG(2.0D0*DPI)+MYNOBS*LOG(SIGTMP)+DLGDET+DM/SIGTMP) & /(2.0D0) C CORRECTIONS FOR STANDARDIZING THE DATA. NOTE: DIVIDING BY THE SD C OF THE DATA ONLY CHANGES THE LIKELIHOOD FUNCTION BY A CONSTANT. C F = F - MYNOBS*LOG(VARDAT)/(2.0D0) GOTO 1000 999 CONTINUE WRITE(*,*) 'IERR =',IERR RETURN 1 1000 RETURN END