C     PROGRAM 10.1  ARMAFT
C
C  ...  ARMA MODEL FITTING  ...
C
C     Inputs:
C        M:       AR Order (M <= 20)
C        L:       MA Order (M <= 20)
C        IPARAM:  =0    Use defalt initail values
C                 =1    Read intial values
C        AR(I):   AR coefficients (I=1,M3)
C     Parameters:
C        NMAX:    Adjustable dimension of Y
C        MJ:      Adjustable dimension of XF, VF, etc.
C     @TEST.FILTER2O    NOV.29,1990, SEP.02,1992
C
      PARAMETER( NMAX=1000,MJ=20 )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  AA(20), AR(20), PAR(20), CMA(20)
      COMMON  /C92825/  OUTMIN, OUTMAX
      COMMON  /C92826/  Y(NMAX)
      COMMON  /C92907/  ALIMIT
      COMMON  /C92908/  M, L, N
      COMMON  /C92909/  FLK, SIG2
      COMMON  / CCC /  ISW, IPR, ISMT, IDIF
      EXTERNAL  FFARMA
C
C  ...  Read Model Orders  ...
C
      READ( 5,* )  M, L, IPARAM
C
C  ...  Set Defalt Parameters  ...
C
      IPR  = 7
      ALIMIT = 0.95D0
      CALL  SPARA1( M,L,AR,CMA,OUTMIN,OUTMAX,IOPT )
      IF( IPARAM.EQ.1 )  THEN
         READ( 5,* )  (AR(I),I=1,M)
         READ( 5,* )  (CMA(I),I=1,L)
      END IF
C
C  ...  Read Time Series  ...
C
      CALL  READTS( 1,Y,N )
C
C  ...  Subtrac Mean Value  ...
C
      SUM = 0.0D0
      DO 10 I=1,N
   10 SUM = SUM + Y(I)
      YMEAN = SUM/N
      DO 20 I=1,N
   20 Y(I) = Y(I) - YMEAN
C
      WRITE(6,*) M, L
      WRITE(6,*) (AR(I),I=1,M)
C
      CALL  PARCOR( AR,M,PAR )
      DO 30 I=1,M
   30 AA(I) = DLOG( (ALIMIT+PAR(I))/(ALIMIT-PAR(I)) )
      CALL  PARCOR( CMA,L,PAR )
      DO 40 I=1,L
   40 AA(M+I) = DLOG( (ALIMIT+PAR(I))/(ALIMIT-PAR(I)) )
C
C  ...  Maximum Likelihood Method  ...
C
      IF( IOPT.EQ.1 )  CALL  DAVIDN( FFARMA,AA,M+L,2 )
C
      DO 50 I=1,M
   50 PAR(I)  = ALIMIT*(DEXP(AA(I))-1.0D0)/(DEXP(AA(I))+1.0D0)
      CALL  ARCOEF( PAR,M,AR )
      DO 60 I=1,L
   60 PAR(I)  = ALIMIT*(DEXP(AA(M+I))-1.0D0)/(DEXP(AA(M+I))+1.0D0)
      CALL  ARCOEF( PAR,L,CMA )
C
      AIC = -2*FLK + 2*(M+L+1)
C
C  ...  Print out the Maximum Likelihood Estimates  ...
C
      CALL  PRARMA( M,L,AR,CMA,SIG2,FLK,AIC )
C
      STOP
      E N D
      SUBROUTINE  SPARA1( M,L,AR,CMA,OUTMIN,OUTMAX,IOPT )
C
C  ...  Set Default Parameters  ...
C
C     Inputs:
C       M:       AR order
C       L:       MA order
C     Outputs:
C       TAU2:    System noise variance
C       AR:      AR coefficients
C       OUTMIN:  Lower bound for outliers
C       OUTMAX:  Upper bound for outliers
C       IOPT:    (=1  MLE by numerical optimization)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  AR(20), PAR(20), CMA(20)
C
      DO 10 I=1,M
   10 PAR(I) = -(-0.6D0)**I
      CALL  ARCOEF( PAR,M,AR )
      DO 20 I=1,L
   20 PAR(I) = -(-0.6D0)**I
      CALL  ARCOEF( PAR,L,CMA )
C
      OUTMIN = -1.0D30
      OUTMAX =  1.0D30
      IOPT = 1
C
      RETURN
      E N D
      SUBROUTINE  PRARMA( M,L,AR,CMA,SIG2,FF,AIC )
C
C  ...  Print out Estimates  ...
C
C     Inputs:
C        M:      AR Order
C        L:      MA Order
C        AR:     AR coefficients
C        CMA:    MA coefficients
C        SIG2:   Innovation variance
C        FF:     Log-Likelihood of the mode
C        AIC:    AIC of the model
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER  TITLE*72
      DIMENSION  AR(*), CMA(*)
      COMMON     /CMDATA/  TITLE
C
      WRITE(6,600)
      WRITE(6,610) TITLE
      WRITE(6,620) M, L
      IF( M.GT.0 )  THEN
        WRITE(6,690)
        WRITE(6,700) (AR(I),I=1,M)
      END IF
      IF( L.GT.0 )  THEN
        WRITE(6,690)
        WRITE(6,700) (CMA(I),I=1,L)
      END IF
      WRITE(6,680) SIG2, FF, AIC
      RETURN
  600 FORMAT( 1H1,'PROGRAM 13.1:  ARMA MODEL FITTING' )
  610 FORMAT( 1H ,A72 )
  620 FORMAT( 1H ,'M =',I2,3X,'L =',I2 )
  650 FORMAT( 1H ,5D13.6 )
  680 FORMAT( 1H ,'SIG2 =',D13.6,3X,'LOG-LIKELIHOOD =',F13.4,3X,
     *            'AIC =',F13.4 )
  690 FORMAT( 1H ,'AR COEFFICIENTS' )
  700 FORMAT( 1H ,5F10.6 )
      E N D
      SUBROUTINE  FILTR3( Y,XF,VF,A,B,C,M,MJ,NS,N,OUTMIN,OUTMAX,
     *                    FF,OVAR )
C
C  ...  Kalman filter  ...
C
C     Inputs:
C        Y:      time series
C        N:      data length
C        NS:     Start position of filtering
C        XF:     Initial state vector
C        VF:     Initial covariance matrix
C        A:      Parameters of the matrix F
C        B:      Parameters of the matrix G
C        C:      Parameters of the matrix H
C        K:      Dimension of the state vector
C        MJ:     Adjustable dimension of XF, VF
C        OUTMIN: Lower limit for detecting outliers
C        OUTMAX: Upper limit for detecting outliers
C     Outputs:
C        FF:     Log likelihood
C        SIG2:   Estimated variance
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  A(MJ), B(MJ), C(MJ), Y(N)
      DIMENSION  XF(MJ), VF(MJ,MJ), XP(40), VP(40,40)
      DIMENSION  WRK(40,40), VH(40), GAIN(40)
      DATA   PI  /3.1415926535D0/
C
      OVAR = 0.0D0
      SDET = 0.0D0
      NSUM = 0
C
      DO 300  II=NS,N
C
C  ...  ONE STEP AHEAD PREDICTION  ...
C
      XP(M) = A(M)*XF(1)
      DO 100  I=1,M-1
  100 XP(I) = A(I)*XF(1) + XF(I+1)
C
      DO 110  J=1,M
      WRK(M,J) = A(M)*VF(1,J)
      DO 110  I=1,M-1
  110 WRK(I,J) = A(I)*VF(1,J) + VF(I+1,J)
C
      DO 120  I=1,M
      VP(I,M) = WRK(I,1)*A(M)
      DO 120  J=1,M-1
  120 VP(I,J) = WRK(I,1)*A(J) + WRK(I,J+1)
C
      DO 140  I=1,M
      DO 140  J=1,M
  140 VP(I,J) = VP(I,J) + B(I)*B(J)
C
C  ...  FILTERING  ...
C
      IF( Y(II).GT.OUTMIN .AND. Y(II).LT.OUTMAX .AND. II.LE.N )  THEN
C
      DO 210  I=1,M
  210 VH(I) = VP(I,1)
C
      PVAR = VH(1)
      PERR = Y(II) - XP(1)
C
      DO 230  I=1,M
  230 GAIN(I) = VH(I)/PVAR
C
      DO 250  I=1,M
  250 XF(I) = XP(I) + GAIN(I)*PERR
C
      DO 260  J=1,M
      DO 260  I=1,M
  260 VF(I,J) = VP(I,J) - GAIN(I)*VH(J)
C
      OVAR = OVAR + PERR**2/PVAR
      SDET = SDET + DLOG(PVAR)
      NSUM = NSUM + 1
C
C  ...  MISSING OBSERVATION  ...
C
      ELSE
      DO 270  I=1,M
      XF(I) = XP(I)
      DO 270  J=1,M
  270 VF(I,J) = VP(I,J)
      END IF
C
  300 CONTINUE
      OVAR = OVAR/NSUM
      FF = -0.5D0*(NSUM*DLOG(PI*2*OVAR) + SDET + NSUM)
C
      RETURN
      E N D
      SUBROUTINE  ISTAT3( M,L,AR,CMA,MJ,XF,VF )
C
C  ...  Initial state ...
C
C     Inputs:
C        M:     AR order
C        L:     MA order
C        AR:    AR coefficient
C        CMA:   MA coefficient
C        MJ:    Adjustable dimension of F
C     Outputs:
C         XF:   State vector, X(0|0)
C         VF:   State covarance matrix, V(0|0)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  AR(*), CMA(*)
      DIMENSION  XF(MJ), VF(MJ,MJ), COV(0:20), G(0:20)
C
      MM = MAX0( M,L+1 )
      DO 10  I=1,MJ
      XF(I) = 0.0D0
      DO 10  J=1,MJ
   10 VF(I,J) = 0.0D0
C
      CALL  ARMCOV( M,L,AR,CMA,1.0D0,MM,COV )
      CALL  IMPULS( M,L,AR,CMA,MM,G )
      VF(1,1) = COV(0)
      DO 50 I=2,MM
      SUM = 0.0D0
      DO 30 J=I,M
   30 SUM = SUM + AR(J)*COV(J-I+1)
      DO 40 J=I-1,L
   40 SUM = SUM - CMA(J)*G(J-I+1)
      VF(1,I) = SUM
   50 VF(I,1) = SUM
      DO 100 I=2,MM
      DO 100 J=I,MM
      SUM = 0.0D0
      DO 60 I1=I,M
      DO 60 J1=J,M
   60 SUM = SUM + AR(I1)*AR(J1)*COV(IABS(J1-J-I1+I))
      DO 70 I1=I,M
      DO 70 J1=J-I+I1,L
   70 SUM = SUM - AR(I1)*CMA(J1)*G(IABS(J1-J-I1+I))
      DO 80 I1=J,M
      DO 80 J1=I-J+I1,L
   80 SUM = SUM - AR(I1)*CMA(J1)*G(IABS(J1-I-I1+J))
      DO 90 I1=I-1,L
   90 SUM = SUM + CMA(I1)*CMA(J-I+I1)
      VF(I,J) = SUM
  100 VF(J,I) = SUM
C
      RETURN
      E N D
      SUBROUTINE  SETABC( M,L,AR,CMA,A,B,C )
C
C  ...  State space model for Seasonal Adjustment  ...
C
C     Input:
C       M:     AR Order
C       L:     MA Order
C       AR:    AR coefficient
C       CMA:   MA coefficient
C     Outputs:
C       A,B,C: Parameters of F, G, H
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  A(*), B(*), C(*)
      DIMENSION  AR(*), CMA(*)
C
      MM = MAX0( M,L+1 )
      DO 10  I=1,MM
      C(I) = 0.0D0
      A(I) = 0.0D0
   10 B(I) = 0.0D0
C
      DO 20  I=1,M
   20 A(I) = AR(I)
      B(1) = 1.0D0
      DO 30 I=1,L
   30 B(I+1) =-CMA(I)
      C(1) = 1.0D0
C
      RETURN
      E N D
      SUBROUTINE  FFARMA( K,AA,FF,GDUMMY,IFG )
C
C  ... Function for Maximum likelihood estimation (seasonal model) ...
C
C     Inputs:
C        K:      Number of parameters
C        AA(I):  Parameter vector
C     Outputs:
C        FF:     -(Log likelihood)
C        IFG:    =1 if some conditions are violated
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  PAR(20), GDUMMY(20)
      DIMENSION  AA(20), AR(20), CMA(20)
      DIMENSION  A(20), B(20), C(20)
      DIMENSION  XF(40), VF(40,40)
      COMMON  /C92826/   Y(1000)
      COMMON  /C92825/   OUTMIN, OUTMAX
      COMMON  /C92907/  ALIMIT
      COMMON  /C92908/   M, L, N
      COMMON  /C92909/   FLK, SIG2
      MJ  = 20
C
      MM = MAX0( M,L+1 )
      DO 10 I=1,M
   10 PAR(I)  = ALIMIT*(DEXP(AA(I))-1.0D0)/(DEXP(AA(I))+1.0D0)
      CALL  ARCOEF( PAR,M,AR )
      DO 20 I=1,L
   20 PAR(I)  = ALIMIT*(DEXP(AA(M+I))-1.0D0)/(DEXP(AA(M+I))+1.0D0)
      CALL  ARCOEF( PAR,L,CMA )
      IFG = 0
C
      CALL  SETABC( M,L,AR,CMA,A,B,C )
      CALL  ISTAT3( M,L,AR,CMA,MJ,XF,VF )
      CALL  FILTR3( Y,XF,VF,A,B,C,MM,MJ,1,N,OUTMIN,OUTMAX,FLK,SIG2 )
      FF = -FLK
      RETURN
C
  100 IFG = 1
      FF = 1.0D20
      RETURN
      E N D
      SUBROUTINE  READTS( IDEV,Y,N )
      REAL*8        Y(N)
      CHARACTER*72  TITLE
      COMMON  /CMDATA/  TITLE
C
C     OPEN( IDEV,FILE='temp.dat' )
      READ( IDEV,1 )  TITLE
      READ( IDEV,* )  N
      READ( IDEV,* )  (Y(I),I=1,N)
C     CLOSE( IDEV )
      RETURN
    1 FORMAT( A72 )
      E N D
      SUBROUTINE IMPULS( M,L,A,B,K,G )
C
C ...  Impulse Response Function  ...
C
C     Inputs:
C        M:     AR order
C        L:     MA order
C        A(I):  AR coefficient
C        B(I):  MA coefficient
C        K:     Required maximum lag of impulse respose
C     Output:
C        G(I):  Impulse response function
C     Y.I.
      IMPLICIT REAL*8( A-H,O-Z )
      DIMENSION A(*), B(*), G(0:K)
C
      G(0) = 1.0
      DO  20 I=1,K
      SUM = 0.0D0
      IF(I.LE.L) SUM = -B(I)
      DO  10 J=1,I
   10 IF(J.LE.M) SUM = SUM + A(J)*G(I-J)
   20 G(I) = SUM
C
      RETURN
      E N D
      SUBROUTINE ARMCOV( M,L,A,B,SIG2,K,COV )
C
C ...  Autocovariance Function of ARMA model  ...
C
C     Inputs:
C        M:     AR order
C        L:     MA order
C        A(I):  AR coefficient
C        B(I):  MA coefficient
C        SIG2:  innovation variance
C        K:     Required maximum lag of autocovariance
C     Output:
C        COV(I):  Autocovariance function
C     Y.I.
      IMPLICIT REAL*8( A-H,O-Z )
      DIMENSION  A(*), B(*), COV(0:K), G(0:100), X(30,30)
      DIMENSION  Z(100), UL(30,30), IPS(100)
C
      KMAX = MAX(M,L,K)
      CALL  IMPULS( M,L,A,B,KMAX,G)
C
      DO 10 I=1,M+1
      DO 10 J=1,M+1
   10 X(I,J) = 0.0D0
      DO 20 I=1,M+1
   20 X(I,I) = 1.0D0
      DO 30 I=1,M
      DO 30 J=2,M-I+2
   30 X(I,J) = X(I,J) - A(I+J-2)
      DO 40 I=2,M+1
      DO 40 J=1,I-1
   40 X(I,J) = X(I,J) - A(I-J)
C
      CALL  DECOM( M+1,X,30,UL,IPS )
C
      SUM = 1.0D0
      DO 50 J=1,L
   50 SUM = SUM - B(J)*G(J)
      Z(1)= SIG2*SUM
      DO 70 I=2,M+1
      SUM = 0.0D0
      DO 60 J=I-1,L
   60 SUM = SUM - B(J)*G(J-I+1)
   70 Z(I) = SIG2*SUM
C
      CALL  SOLVE( M+1,UL,30,Z,COV,IPS)
C
      DO 100 J=M+1,K
      SUM = 0.0D0
      DO 80 I=1,M
   80 SUM = SUM + A(I)*COV(J-I)
      DO 90 I=J,L
   90 SUM = SUM - B(I)*G(I-J)*SIG2
  100 COV(J) = SUM
C
      RETURN
      E N D
      SUBROUTINE DECOM( N,A,MJ,UL,IPS )
C
C  ...  UL decomposition:  A = L*U  ...
C
C     Inputs:
C        N:      Dimension of the matrix A
C        A(I,J): N*N positive definite matrix
C        MJ:     Adjustable dimension of A and UL
C     Outputs:
C        UL(I,J):  L-I and U
C        IPS:    Index vector
C     Y.I.
      IMPLICIT REAL*8(A-H,O-Z )
      DIMENSION A(MJ,*),UL(MJ,*),SCALES(100),IPS(100)
C
      DO 20 I=1,N
      IPS(I) = I
      RNORM = 0.0D0
      DO 10 J=1,N
      UL(I,J) = A(I,J)
      IF(RNORM.LT.ABS(UL(I,J)) ) RNORM = ABS( UL(I,J) )
   10 CONTINUE
      IF( RNORM .NE. 0.0D0 )  THEN
          SCALES(I) = 1/RNORM
      ELSE
          SCALES(I) = 0.0D0
          CALL  SING(0)
      END IF
   20 CONTINUE
C
      DO 60 K=1,N-1
      BIG = 0.0D0
      DO 30 I=K,N
      SIZE = ABS( UL(IPS(I),K) )*SCALES( IPS(I) )
      IF( BIG.LT.SIZE ) THEN
          BIG = SIZE
          INDEX = I
      END IF
   30 CONTINUE
      IF( BIG.EQ. 0.0D0 )  THEN
          CALL  SING(1)
          GO TO 60
      END IF
      IF( INDEX.NE.K ) THEN
      J = IPS(K)
      IPS(K) = IPS(INDEX)
      IPS(INDEX) = J
      END IF
C
      PIVOT = UL(IPS(K),K)
      DO 50 I=K+1,N
      TM = UL( IPS(I),K)/PIVOT
      UL( IPS(I),K) = TM
      IF( TM.NE. 0.0D0 )  THEN
      DO 40 J = K+1,N
   40 UL( IPS(I),J ) = UL( IPS(I),J)-TM*UL( IPS(K),J)
C     WRITE(6,*) (UL(IPS(I),J),J=1,N)
      END IF
   50 CONTINUE
   60 CONTINUE
C
      IF( UL(IPS(N),N) .EQ. 0.0D0 )   CALL  SING(2)
      RETURN
      E N D
      SUBROUTINE  SOLVE( N,UL,MJ,B,X,IPS )
C
C  ...  Solve Ax=b using UL obtained by DECOM  ...
C
C     Inputs:
C        N:     Dimension of UL and B
C        UL:    LU decomposition of A
C        MJ:    Adjustable dimension of A
C        B:
C        IPS:   index vector
C     Output:
C        X:     Solution
C     Y.I.
      IMPLICIT REAL*8( A-H,O-Z )
      DIMENSION UL(MJ,*),B(*),X(*),IPS(100)
C
      DO 20 I=1,N
      SUM = 0.0D0
      DO 10 J=1,I-1
   10 SUM = SUM + UL(IPS(I),J)*X(J)
   20 X(I) = B(IPS(I)) - SUM
C
      DO 40 I=N,1,-1
      SUM = 0.0D0
      DO 30 J=I+1,N
   30 SUM = SUM + UL(IPS(I),J)*X(J)
   40 X(I) = ( X(I)-SUM )/UL(IPS(I),I)
      RETURN
  600 FORMAT(1H ,'N=',I10,/,(5X,'IPS=',I10 ) )
      E N D
      SUBROUTINE  SING( IFG )
C
C  ...  ERROR MESSAGE  ...
C
      IF( IFG.EQ.0) THEN
          WRITE(6,100)
      ELSE
          IF( IFG.EQ.1 ) THEN
              WRITE(6,200)
          ELSE
              WRITE(6,300)
          END IF
      END IF
      RETURN
  100 FORMAT(' MATRIX WITH ZERO ROW IN DECOMPOSE.')
  200 FORMAT(' SINGULAR MATRIX IN DECOMPOSE.ZERO DIDIVIDE IN SOLVE.')
  300 FORMAT(' CONVERGENCE IN IMPRUV.MATRIX IS NEARLY SINGULAR.')
      E N D
      SUBROUTINE  ARCOEF( PAR,K,A )
C
C  ...  This subroutine computes AR coefficients from PARCOR  ...
C
C     Inputs:
C        PAR(I):   PARCOR
C        K:        Order of AR model
C     Output:
C        A(I):     AR coefficient
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  PAR(K), A(K), AA(50)
C
      DO 30  II=1,K
      A(II)  = PAR(II)
      AA(II) = PAR(II)
      IF( II-1.LE.0 )  GO TO 30
      DO 10  J=1,II-1
   10 A(J) = AA(J) - PAR(II)*AA(II-J)
      IF( II.LT.K )  THEN
        DO 20  J=1,II-1
   20   AA(J) = A(J)
      END IF
   30 CONTINUE
      RETURN
      E N D
      SUBROUTINE  PARCOR( A,K,PAR )
C
C  ...  This subriutine computes PARCOR for AR coefficients  ...
C
C       Inputs:
C          A:   Vector of AR coefficients
C          K:   Order of the model
C
C       Output:
C          PAR: Vector of partial autocorrelations (PARCOR)
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  A(K), PAR(K), G(50)
C
      DO 10  I=1,K
   10 PAR(I) = A(I)
C
      DO 40  II=K-1,1,-1
      S = 1.D0 - PAR(II+1)**2
      DO 20  I=1,II
   20 G(I) = (PAR(I) + PAR(II+1)*PAR(II-I+1))/S
      I2 = (II+1)/2
      IF( MOD( II,2 ).EQ.1 )  G(I2) = PAR(I2)/(1.D0 - PAR(II+1))
      DO 30  I=1,II
   30 PAR(I) = G(I)
   40 CONTINUE
C
      RETURN
      E N D
      SUBROUTINE  DAVIDN( FUNCT,X,N,NDIF )
C
C  ...  6/20/83, 12/19/92
C
      IMPLICIT  REAL*8( A-H,O-Z )
      DIMENSION  X(40), DX(40), G(40), G0(40), Y(40)
      DIMENSION  H(40,40), WRK(40), S(40)
      COMMON     / CCC /  ISW, IPR, ISMT, IDIF
      COMMON     / DDD /  XM , AIC , SD
      EXTERNAL  FUNCT
      DATA        TAU2  /          1.0D-6  /
      DATA  EPS1, EPS2  / 1.0D-6 , 1.0D-6  /
      RAMDA  = 0.5D0
      CONST1 = 1.0D-30
      IPR = 0
C
C          INITIAL ESTIMATE OF INVERSE OF HESSIAN
C
      ICOUNT = 0
 1000 DO 20  I=1,N
      DO 10  J=1,N
   10 H(I,J) = 0.0D0
      S(I)   = 0.0D0
      DX(I)  = 0.0D0
   20 H(I,I) = 1.0D0
      ISW = 0
C
      IF( NDIF.EQ.0 )  CALL  FUNCT( N,X,XM,G,IG )
      IF( NDIF.GE.1 )  CALL  FUNCND( FUNCT,N,X,XM,G,IG )
C
      IF( IPR .GE. 2 )   WRITE( 6,640 )   XM, SD, AIC
      WRITE(6,650) (X(I),I=1,N), XM, (G(I),I=1,N)
C
      ICC = 0
C      ITERATION
 2000 CONTINUE
      ICC = ICC + 1
      IF( ICC .EQ. 1 )  GO TO 120
C
      DO 40  I=1,N
   40 Y(I) = G(I) - G0(I)
      DO 60  I=1,N
      SUM = 0.0D0
      DO 50  J=1,N
   50 SUM = SUM + Y(J)*H(I,J)
   60 WRK(I) = SUM
      S1 = 0.0D0
      S2 = 0.0D0
      DO 70  I=1,N
      S1 = S1 + WRK(I)*Y(I)
   70 S2 = S2 + DX(I) *Y(I)
      IF( S1.LE.CONST1 .OR. S2.LE.CONST1 )  GO TO 900
C     IF( S1 .LE. S2 )   GO TO 100
C
C          UPDATE THE INVERSE OF HESSIAN MATRIX
C
C               ---  BROYDEN-FLETCHER-GOLDFARB-SHANNO TYPE CORRECTION  -
C
  100 CONTINUE
      STEM = S1 / S2 + 1.0D0
      DO 110  I=1,N
      DO 110  J=I,N
      H(I,J) = H(I,J)-(DX(I)*WRK(J)+WRK(I)*DX(J)-DX(I)*DX(J)*STEM)/S2
  110 H(J,I) = H(I,J)
C
  120 CONTINUE
      SS = 0.0D0
      DO 140  I=1,N
      SUM = 0.0D0
      DO 130  J=1,N
  130 SUM = SUM + H(I,J)*G(J)
      SS  = SS + SUM**2
  140 S(I) = -SUM
C
      S1 = 0.0D0
      S2 = 0.0D0
      DO 150  I=1,N
      S1 = S1 + S(I)*G(I)
  150 S2 = S2 + G(I)*G(I)
C     DS2 = DSQRT(S2)
C     GTEM = DABS(S1)/DS2
C     IF( GTEM.LE.TAU1 .AND. DS2.LE.TAU2 )  GO TO 900
      IF( S1.LT.0.0D0 )  GO TO 200
      DO 170  I=1,N
      DO 160  J=1,N
  160 H(I,J) = 0.0D0
      H(I,I) = 1.0D0
  170 S(I) = -S(I)
  200 CONTINUE
C
      ED = XM
C
C          LINEAR  SEARCH
C
      CALL  LINEAR( FUNCT,X,S,RAMDA,ED,N,IG )
C
      IF( IPR .GE. 2 )  WRITE( 6,630 )  RAMDA, ED, SD, AIC
C
      S1 = 0.0D0
      DO 210  I=1,N
      DX(I) = S(I)*RAMDA
      S1 = S1 + DX(I)*DX(I)
      G0(I) = G(I)
  210 X(I) = X(I) + DX(I)
      XMB = XM
      ISW = 0
C
      IF( NDIF.EQ.0 )  CALL  FUNCT( N,X,XM,G,IG )
      IF( NDIF.GE.1 )  CALL  FUNCND( FUNCT,N,X,XM,G,IG )
      WRITE(6,650) (X(I),I=1,N), XM, (G(I),I=1,N)
C
      S2 = 0.D0
      DO 220  I=1,N
  220 S2 = S2 + G(I)*G(I)
      IF( DSQRT(S2) .LT. TAU2 )  GO TO 900
      IF( XMB/XM-1.D0.LT.EPS1 .AND. DSQRT(S1).LT.EPS2 )  GO TO 900
C     IF( ICC .GE. 5 )  GO TO 900
      GO TO 2000
  900 CONTINUE
      IF( IPR .LE. 0 )  RETURN
      WRITE( 6,600 )
      WRITE( 6,610 )  (X(I),I=1,N)
      WRITE( 6,620 )
      WRITE( 6,610 )  (G(I),I=1,N)
C
      ICOUNT  = ICOUNT + 1
      S2 = 0.D0
      DO 910  I=1,N
  910 S2 = S2 + G(I)*G(I)
      IF( S2.GT.1.0D0.AND.ICOUNT.LE.1 )   GO TO 1000
      RETURN
  600 FORMAT( 1H0,'-----  X  -----' )
  610 FORMAT( 1H ,10D13.5 )
  620 FORMAT( 1H0,'***  GRADIENT  ***' )
  630 FORMAT( 1H ,'LAMBDA =',D15.7,3X,'(-1)LOG LIKELIHOOD =',D23.15,
     *        3X,'SD =',D22.15,5X,'AIC =',D23.15 )
  640 FORMAT( 1H ,26X,'(-1)LOG-LIKELIHOOD =',D23.15,3X,'SD =',D22.15,
     *        5X,'AIC =',D23.15 )
  650 FORMAT( 10X,5F12.7 )
      E N D
      SUBROUTINE  FUNCND( FUNCT,M,A,F,G,IFG )
C
C  ...  FUNCTION EVALUATION AND NUMERICAL DIFFERENCING  ...
C
      IMPLICIT   REAL*8( A-H,O-Z )
      DIMENSION  A(M) , G(M) , B(20),GD(5)
      COMMON  / CCC /  ISW , IPR, ISMT, IDIF
      COMMON  /CMFUNC/ DJACOB,FC,SIG2,AIC,FI,SIG2I,AICI,GI(20),GC(20)
C     DATA       ICNT /0/
      CONST = 0.00001D0
C
      CALL  FUNCT( M,A,F,GD,IFG )
      FB = F
      IF( ISW .GE. 1 )   RETURN
C
C     WRITE( 6,600 )   (A(I),I=1,M)
      DO 10  I=1,M
   10 B(I) = A(I)
C
      DO 30  II=1,M
      B(II) = A(II) + CONST
      CALL  FUNCT( M,B,FF,GD,IFG )
      IF( IDIF .EQ. 1 )  GO TO 20
      B(II) = A(II) - CONST
      CALL  FUNCT( M,B,FB,GD,IFG )
   20 G(II) = (FF-FB)/(CONST*IDIF)
      IF( G(II) .GT. 1.0D20 )  G(II) = (F-FB)/CONST
      IF( G(II) .LT.-1.0D20 )  G(II) = (FF-F)/CONST
      IF( FB.GT.F .AND. FF.GT.F )  G(II) = 0.0D0
   30 B(II) = A(II)
C
      RETURN
C
  600 FORMAT( 3X,'---  PARAMETER  ---',(/,3X,5D13.5) )
  610 FORMAT( 3X,'---  GRADIENT  ---',(/,3X,5D13.5) )
      E N D
      SUBROUTINE  LINEAR( FUNCT,X,H,RAM,EE,K,IG )
C
C  ...  LINEAR SEARCH  ...
C
      IMPLICIT  REAL*8( A-H,O-Z )
      INTEGER  RETURN,SUB
      DIMENSION  X(1), H(1), X1(40)
      DIMENSION  G(40)
      COMMON     / CCC /  ISW , IPR, ISMT, IDIF
C     IPR = 7
C     C1 = 1.0D-7
C     C2 = 1.0D-30
C
      ISW = 1
      IG  = 0
      RAM = 0.1D0
      CONST2 = 1.0D-60
   14 DO 15 I=1,K
   15 IF( DABS(H(I)) .GT. 1.0D10 )  GO TO 16
      GO TO 18
   16 DO 17 I=1,K
   17 H(I) = H(I)*1.0D-10
      GO TO 14
   18 CONTINUE
      HNORM = 0.D0
      DO 10  I=1,K
   10 HNORM = HNORM + H(I)**2
      HNORM = DSQRT( HNORM )
C
      RAM2 = RAM
      E1 =EE
      RAM1 = 0.D0
C
      DO 20  I=1,K
   20 X1(I) = X(I) + RAM2*H(I)
      CALL  FUNCT( K,X1,E2,G,IG )
      IF(IPR.GE.7)  WRITE(6,2)  RAM2,E2
C
      IF( IG .EQ. 1  )  GO TO 50
      IF( E2 .GT. E1 )  GO TO 50
   30 RAM3 = RAM2*4.D0
      DO 40  I=1,K
   40 X1(I) = X(I) + RAM3*H(I)
      CALL  FUNCT( K,X1,E3,G,IG )
      IF( IG.EQ.1 )  GO TO  500
      IF( IPR.GE.7 )  WRITE(6,3)  RAM3,E3
      IF( E3 .GT. E2 )  GO TO 70
      IF(RAM3.GT.1.0D10 .AND. E3.LT.E1)  GO TO 45
      IF(RAM3.GT.1.0D10 .AND. E3.GE.E1)  GO TO 46
      RAM1 = RAM2
      RAM2 = RAM3
      E1 = E2
      E2 = E3
      GO TO 30
C
   45 RAM = RAM3
      EE = E3
      RETURN
C
   46 RAM = 0.0D0
      RETURN
C
   50 RAM3 = RAM2
      E3 = E2
      RAM2 = RAM3*0.1D0
      IF( RAM2*HNORM .LT. CONST2 )  GO TO  400
      DO 60  I=1,K
   60 X1(I) = X(I) + RAM2*H(I)
      CALL  FUNCT( K,X1,E2,G,IG )
      IF(IPR.GE.7)  WRITE(6,4)  RAM2,E2
      IF( E2.GT.E1 )  GO TO 50
C
   70 ASSIGN 80 TO RETURN
      GO TO 200
C
   80 DO 90  I=1,K
   90 X1(I) = X(I) + RAM*H(I)
      CALL  FUNCT( K,X1,EE,G,IG )
      IF(IPR.GE.7)  WRITE(6,5)  RAM,EE
C
      IFG = 0
      ASSIGN  300 TO  SUB
      ASSIGN 200 TO SUB
   95 ASSIGN 130 TO RETURN
      IF( RAM .GT. RAM2 )  GO TO 110
      IF( EE .GE. E2 )  GO TO 100
      RAM3 = RAM2
      RAM2 = RAM
      E3 =E2
      E2 =EE
      GO TO  SUB,( 200,300 )
C
  100 RAM1 = RAM
      E1 = EE
      GO TO  SUB,( 200,300 )
C
  110 IF( EE .LE. E2 )  GO TO 120
      RAM3 = RAM
      E3 = EE
      GO TO  SUB,( 200,300 )
C
  120 RAM1 = RAM2
      RAM2 = RAM
      E1 = E2
      E2 = EE
      GO TO  SUB,( 200,300 )
C
  130 DO 140  I=1,K
  140 X1(I) = X(I) + RAM*H(I)
      CALL  FUNCT( K,X1,EE,G,IG )
      IF( IPR.GE.7 )  WRITE(6,6)  RAM,EE
      ASSIGN 200 TO SUB
      IFG = IFG+1
  600 FORMAT( 1H ,6D20.13 )
C     IF( HNORM*(RAM3-RAM1) .GT. 1.0D-12 )  GO TO 95
      IF( RAM3-RAM1 .GT. RAM2*1.0D-1 )  GO TO 95
C     IF( E1-EE .GT. 1.0D-15 )  GO TO 95
C     IF( E3-EE .GT. 1.0D-15 )  GO TO 95
C     IF(DX.LE.C1*XN+C2 .AND. DF.LE.C1*DABS(EE)+C2)  RETURN
C     IF( IFG .LE. 2 )  GO TO 95
C
      IF( E2 .LT. EE )  RAM = RAM2
      RETURN
C
C      -------  INTERNAL SUBROUTINE SUB1  -------
  200 IF( RAM3-RAM2 .GT. 5.0D0*(RAM2-RAM1) )  GO TO 202
      IF( RAM2-RAM1 .GT. 5.0D0*(RAM3-RAM2) )  GO TO 204
      A1 = (RAM3-RAM2)*E1
      A2 = (RAM1-RAM3)*E2
      A3 = (RAM2-RAM1)*E3
      B2 = (A1+A2+A3)*2.D0
      B1 = A1*(RAM3+RAM2) + A2*(RAM1+RAM3) + A3*(RAM2+RAM1)
      IF( B2 .EQ. 0.D0 )  GO TO 210
      RAM = B1 /B2
C
      IF( RAM .LE. RAM1 )  RAM = (RAM1 + RAM2)/2.0D0
      IF( RAM .GE. RAM3 )  RAM = (RAM2 + RAM3)/2.0D0
      IF( DABS(RAM-RAM2) .LE. 1.0D-15 )  RAM = (RAM2*4.D0 + RAM3)/5.0D0
      GO TO RETURN ,( 80,130 )
  202 RAM = (4.0D0*RAM2 + RAM3)/5.0D0
      GO TO RETURN, (80,130)
  204 RAM = (RAM1 + 4.0D0*RAM2)/5.0D0
      GO TO RETURN, (80,130)
C
  210 IG = 1
      RAM = RAM2
      RETURN
C
C      -------  INTERNAL SUBROUTINE SUB2  -------
C
  300 IF( RAM3-RAM2 .GT. RAM2-RAM1 )  GO TO 310
      RAM = (RAM1+RAM2)*0.5D0
      GO TO RETURN ,( 80,130 )
C
  310 RAM = (RAM2+RAM3)*0.5D0
      GO TO RETURN ,( 80,130 )
C
  400 RAM = 0.D0
      RETURN
C
  500 RAM = (RAM2+RAM3)*0.5D0
  510 DO 520  I=1,K
  520 X1(I) = X(I) + RAM*H(I)
      CALL  FUNCT( K,X1,E3,G,IG )
      IF( IPR.GE.7 )  WRITE(6,7)  RAM,E3
      IF( IG.EQ.1 )  GO TO 540
      IF( E3.GT.E2 )  GO TO 530
      RAM1 = RAM2
      RAM2 = RAM
      E1 = E2
      E2 = E3
      GO TO 500
C
  530 RAM3 = RAM
      GO TO 70
C
  540 RAM = (RAM2+RAM)*0.5D0
      GO TO 510
C
C ------------------------------------------------------------
    1 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E1 =',D25.17 )
    2 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E2 =',D25.17 )
    3 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E3 =',D25.17 )
    4 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E4 =',D25.17 )
    5 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E5 =',D25.17 )
    6 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E6 =',D25.17 )
    7 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E7 =',D25.17 )
      E N D