C     PROGRAM 6.1  ARMA
C
C  ...  This program analyses time series via ARMA modeling  ...
C
C     Inputs:
C        IDEV:  Input devie (=5: CONSOLE)
C        M:     AR order
C        L:     MA order
C        SIG2:  Innovation variance
C        A(I):  AR coefficients
C        B(I):  MA coefficients
C        K:     Maximum lag of autocovariance function (K<MJ+1)
C     Parameter:
C        MJ:    Adjustable dimension of A, B, etc..
C        NF:    Number of frequencies in evaluating spectrum
C     Programmed by Y.I and G.K.
C
      PARAMETER( NF=200,MJ=100,IDEV=1 )
      IMPLICIT REAL*8(A-H,O-Z )
      DIMENSION  A(MJ), B(MJ), PAR(MJ), ROOTA(MJ,2), ROOTB(MJ,2)
      DIMENSION  G(0:MJ), COV(0:MJ), SP(0:NF)
      DIMENSION  WRK1(0:MJ), WRK2(0:MJ), WRK3(MJ)
      DATA  PAR/MJ*0.0D0/
C
      WRITE( 6,* )  'K = ?'
      READ( 5,* )   K
      OPEN( IDEV,FILE='A:$LASERLIB$ARMA.DAT' )
      READ(IDEV,*)  M, L
      READ(IDEV,*)  SIG2
      READ(IDEV,*)  (A(I),I=1,M)
      READ(IDEV,*)  (B(I),I=1,L)
      CLOSE( IDEV )
C
      CALL  IMPULS( M,L,A,B,K,G )
      CALL  ARMCOV( M,L,A,B,SIG2,K,COV )
      CALL  PARCOR( A,M,PAR )
      CALL  ARCOEF( PAR,M,A )
      IF( L.GT.0 )  CALL  ARYULE( COV,1000,K,WRK1,WRK2,PAR,WRK3,MAR )
      CALL  ARMASP( A,M,B,L,SIG2,NF,SP )
      CALL  CHROOT( A,M,ROOTA,MJ )
      CALL  CHROOT( B,L,ROOTB,MJ )
      CALL  PRARMA( M,L,A,B,G,K,COV,K,PAR,SP,NF,ROOTA,ROOTB,MJ )
      CALL  PTARMA( G,K,COV,K,PAR,SP,NF,ROOTA,M,ROOTB,L,MJ )
      STOP
      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  PTARMA( G,K,COV,K2,PARCOR,SP,NF,ROOTA,M,
     *                    ROOTB,L,MJ )
C
C  ...  This subroutine draws impulse response function, PARCOR,
C       autocovariance function, PARCOR, power spectrum and
C       characteristics roots of the AR and MA operators  ...
C
C     Inputs:
C        G(I):    Impulse response function (I=0,K)
C        COV(I):  Autocovariance function (I=0,K2)
C        SP(I):   Power spectrum (I=0,NF)
C        PARCOR(I):  Partial autocorrelation function (I=1,K2)
C        ROOTA(I):   Characteristic roots of AR operator (I=1,M)
C        ROOTB(I):   Characteristic roots of MA operator (I=1,L)
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER  TITLE*72, VNAME*8
      DIMENSION  G(0:K), COV(0:K2), PARCOR(K)
      DIMENSION  ROOTA(MJ,2), ROOTB(MJ,2), SP(0:NF)
      DIMENSION  VNAME(5), YM(5)
      COMMON  /CMDATA/  TITLE
C
      WX = 6.0D0
      WY = 6.0D0
      TITLE    = '        '
      VNAME(1) = 'AR ORDER'
      VNAME(2) = 'MA ORDER'
      YM(1)    = M
      YM(2)    = L
C
      CALL  PLOTS
C     call  plots( 1,0,0,1,0 )
C     call  form( 1 )
C     call  factor( 10.0 )
      CALL  HEADER( 'PROGRAM 6.1:  ARMA MODEL',24,2,VNAME,YM )
      CALL  DRAWY ( 'IMPULSE RESPONSE',16,3.0D0,15.0-WY,G,K+1,
     *               WX,WY,0,1)
      CALL  DELX( 0.0,DBLE(K),DX )
      CALL  DRAWY ( 'AUTOCOVARIANCE',14,WX+2.0,0.0D0,COV,K2+1,
     *               WX,WY,0,1)
      CALL  DRAWY2( 'PARCOR',6,WX+2.0,0.0D0,PARCOR,K,WX,WY,1,1,
     *               0.0D0,DBLE(K),DX,-1.0D0,1.0D0,1.0D0 )
      CALL  PTCIRC( -2*(WX+2.0),-(WY+2.0),WX,1 )
      CALL  PLOTSC( ROOTA(1,1),ROOTA(1,2),M, -1.0D0,1.0D0,
     *              -1.0D0,1.0D0,WX,WX,0.3D0,'*' )
      CALL  PLOTSC( ROOTB(1,1),ROOTB(1,2),L,-1.0D0,1.0D0,
     *              -1.0D0,1.0D0,WX,WX,0.3D0,'+' )
      CALL  SYMBOL( 0.0,SNGL(WX)+0.2,0.25,'CHARACTERISTIC ROOTS',
     *              0.0,20 )
      CALL  MAXMIN( SP(0),NF+1,SMIN,SMAX,DY )
      CALL  DRAWY2( 'POWER SPECTRUM (IN LOG SCALE)',29,WX+2.0,0.0D0,
     *              SP,NF,WX,WY,0,0,0.0D0,0.5D0,0.1D0,SMIN,SMAX,DY )
      CALL  PLOTE
C     call  plot( 0.0,0.0,999 )
C
      RETURN
      E N D
      SUBROUTINE  PRARMA( M,L,A,B,G,K,COV,K2,PARCOR,SP,NF,ROOTA,
     *                    ROOTB,MJ )
C
C  ...  This subroutine prints out impulse response function,
C       autocovariance function, PARCOR, power spectrum and
C       characteristic roots of the AR and MA operators  ...
C
C     Inputs:
C        A(I):    AR coefficient (I=1,M)
C        B(I):    MA coefficient (I=1,L)
C        G(I):    Impulse response function (I=0,K)
C        COV(I):  Autocovariance function (I=0,K2)
C        SP(I):   Power spectrum (I=0,NF)
C        PARCOR(I):  Partial autocorrelation function (I=1,K2)
C        ROOTA(I):   Characteristic roots of AR operator (I=1,M)
C        ROOTB(I):   Characteristic roots of MA operator (I=1,L)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  G(0:K), COV(0:K2), PARCOR(K2), A(M), B(L)
      DIMENSION  ROOTA(MJ,2), ROOTB(MJ,2), SP(0:NF)
C
      WRITE(6,600)
      WRITE(6,610)  M, L
      IF(M.GT.0)  WRITE(6,620)  (A(I),I=1,M)
      IF(L.GT.0)  WRITE(6,630)  (B(I),I=1,L)
      WRITE(6,640)  (G(I),I=1,K)
      WRITE(6,650)  (COV(I),I=0,K2)
      WRITE(6,660)  (PARCOR(I),I=1,K2)
      WRITE(6,670)  (SP(I),I=0,NF)
      WRITE(6,700)
      WRITE(6,680)
      DO 10  I=1,M
      AMP = DSQRT ( ROOTA(I,1)**2 + ROOTA(I,2)**2 )
      TAN = DATAN2( ROOTA(I,2),ROOTA(I,1) )
      DEG = TAN * 0.5729577951D02
   10 WRITE(6,690)  ROOTA(I,1), ROOTA(I,2), AMP, TAN, DEG
      WRITE(6,710)
      WRITE(6,680)
      DO 20  I=1,L
      AMP = DSQRT ( ROOTB(I,1)**2 + ROOTB(I,2)**2 )
      TAN = DATAN2( ROOTB(I,2),ROOTB(I,1) )
      DEG = TAN * 0.5729577951D02
   20 WRITE(6,690)  ROOTB(I,1), ROOTB(I,2), AMP, TAN, DEG
C
      RETURN
  600 FORMAT( 1H ,'PROGRAM 6.1:  ARMA MODEL' )
  610 FORMAT( 1H ,'M =',I3,5X,'L =',I3 )
  620 FORMAT( 1H ,'**  AR COEFFICIENTS  **',/(10F8.4) )
  630 FORMAT( 1H ,'**  MA COEFFICIENTS  **',/(10F8.4) )
  640 FORMAT( 1H ,'**  IMPULSE RESPONSE FUNCTION  **',/(5F15.8) )
  650 FORMAT( 1H ,'**  AUTOCOVARIANCE FUNCTION  **',/(5D15.7) )
  660 FORMAT( 1H ,'**  PARCOR  **',/(10F8.4) )
  670 FORMAT( 1H ,'**  POWER SPECTRUM  **',/(10F8.3) )
  680 FORMAT( 9X,'REAL',8X,'IMAGINARY',10X,'SQRT(R**2+I**2)  ARCTAN(I/R)
     1',5X,'DEGREE' )
  690 FORMAT( 1H ,2D15.5,7X,2D15.5,F12.3 )
  700 FORMAT( 1H ,'**  CHARACTERISTIC ROOTS OF AR OPERATOR  **' )
  710 FORMAT( 1H ,'**  CHARACTERISTIC ROOTS OF MA OPERATOR  **' )
      E N D
      SUBROUTINE  PTCIRC( X0,Y0,WX,ISW )
C
C  ...  This subroutine draws circle  ...
C
      REAL*8  X0, Y0, WX
      C = 3.141593/180
C
      R  = WX/2
      XX = X0 + R
      YY = Y0 + R
      CALL  PLOT( XX,YY,-3 )
C
      CALL  NEWPEN( 1 )
      CALL  PLOT( R,0.0,3 )
      DO 10 I=1,360
        THETA = C*I
        XX = COS( THETA )*R
        YY = SIN( THETA )*R
   10 CALL  PLOT( XX,YY,2 )
C
      IF( ISW.EQ.1 )  THEN
      CALL  PLOT(  R,0.0,3 )
      CALL  PLOT( -R,0.0,2 )
      CALL  PLOT( 0.0, R,3 )
      CALL  PLOT( 0.0,-R,2 )
      END IF
C
      CALL  PLOT( -R,-R,-3 )
      RETURN
      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 ARMASP( A,M,B,L,SIG2,NF,SP )
C
C  ...  Logarithm of the power spectrum of the 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        NF:    Number of frequencies
C     Output:
C        SP(I): Power spectrum (in log scale)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(M), B(L)
      DIMENSION SP(0:NF), H(0:500), FR(0:500), FI(0:500)
C
      H(0) = 1.0D0
      DO 10 I=1,M
   10 H(I) = -A(I)
C
      CALL  FOURIE( H,M+1,NF+1,FR,FI )
C
      DO 20 I=0,NF
   20 SP(I) = SIG2/( FR(I)**2 + FI(I)**2 )
C
      H(0) = 1.0D0
      DO 30 I=1,L
   30 H(I) = -B(I)
      CALL  FOURIE( H,L+1,NF+1,FR,FI )
      DO 40 I=0,NF
   40 SP(I) = SP(I)*( FR(I)**2 + FI(I)**2 )
C
      DO 50 I=0,NF
   50 SP(I) = DLOG10( SP(I) )
C
      RETURN
      E N D
      SUBROUTINE  CHROOT( A,M,ROOT,MJ )
C
C  ...  Characteristic roots of the AR or MA operator  ...
C
C     Inputs:
C        A:     Vector of AR or MA coefficients
C        M:     Order of the model
C        MJ:    Adjustable dimension of ROOT
C     Output:
C        ROOT:  Characteristic roots (real part,imaginary part)
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  A(M), ROOT(MJ,2)
      DIMENSION  C(50), CW(50)
C
      IF( M.EQ.0 )  RETURN
      DO 10  I=1,M
   10 C(I) = -A(M-I+1)
      C(M+1) = 1.0D0
      MMAX = M
C
C  ... Characteristic roots of operator 1-A(1)*S- ... -A(M)*S**M=0  ...
C
      CALL  POLYRT( C,CW,MMAX,ROOT(1,1),ROOT(1,2),IER )
      IF( IER.NE.0 )   WRITE(6,600)
C
      RETURN
  600 FORMAT( 1H0,'*****  NON-CONVERGENCE AT POLYRT  *****' )
      E N D
      SUBROUTINE  POLYRT( A,B,M,ROOTR,ROOTI,IER )
C
C  ...  This subroutine finds the roots of the equation
C            A(1) + A(2)*S + ... + A(M)*(S**M) = 0
C       by Newton-Raphson method
C
C     Inputs:
C        A:     Coefficients of the equation
C        B:     Working area
C        M:     Degree of the polynomial
C     Outputs:
C        ROOTR:   Real parts of the roots
C        ROOTI:   Imaginary parts of the roots
C        IER:     Error code to indicate non-convergence
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION A(1), B(1), ROOTR(1), ROOTI(1)
C
      IFIT = 0
      ISW = 0
      JSW = 0
      K = M
      IER = 0
      KX = K
      KXX = K+1
      KJI = KXX
      K2 = 1
      DO 10  I=1,KXX
      A(I) = A(I)/A(K+1)
      J = KXX - I+1
   10 B(J) = A(I)
   20 XO = 0.5D-02
      YO = 0.1D-01
      IN = 0
   30 X = XO
C
      XO = -10.D0*YO
      YO = -10.D0*X
      X = XO
      Y = YO
      ICT = 0
      IN = IN+1
      GO TO 50
C
   40 IFIT = 1
      XPR = X
      YPR = Y
C
   50 UX = 0.D0
      UY = 0.D0
      V  = 0.D0
      YT = 0.D0
      XT = 1.D0
      U = B(K+1)
      IF( U .EQ. 0.D0 )  GO TO 140
      DO 60  I=1,K
      L = K-I+1
      XT2 = X*XT - Y*YT
      YT2 = X*YT+Y*XT
      U = U + B(L)*XT2
      V = V + B(L)*YT2
      FI = I
      UX = UX + FI*XT*B(L)
      UY = UY - FI*YT*B(L)
      XT = XT2
      YT = YT2
   60 CONTINUE
      SUM = UX**2 + UY**2
      IF( SUM .EQ. 0.D0 )  GO TO 100
      DX = (V*UY - U*UX)/SUM
      X = X + DX
      DY = -(U*UY + V*UX)/SUM
      Y = Y + DY
      IF( DABS(DY)+DABS(DX) .LT. 1.0D-10 )  GO TO 80
C
      ICT = ICT+1
      IF( ICT .LT. 500 )  GO TO 50
      ISW = 1
      IF( IN .GE. 5 )  GO TO 70
      IF( IFIT .NE. 0 )  GO TO 80
   65 ISW = 0
      IFIT = 0
      GO TO 30
C
   70 IF( IFIT .EQ. 0 )  GO TO 300
      JSW = 1
   80 DO 90  L=1,KXX
      J = KJI-L+1
      TEM = A(J)
      A(J) = B(L)
   90 B(L) = TEM
      ITEM = K
      K = KX
      KX = ITEM
  100 IF( IFIT .EQ. 0 )  GO TO 40
      IF( JSW .EQ. 1 )   GO TO 110
      IF( ISW-1 )  120,65,120
  110 X = XPR
      Y = YPR
      ISW = 0
      JSW = 0
  120 IFIT = 0
      IF( X .EQ. 0.D0 )  GO TO 130
      IF( DABS( Y/X ) .LT. 1.0D-08 )  GO TO 150
  130 ALPH = 2*X
      SUM = X*X + Y*Y
      K = K-2
      GO TO 160
C
  140 X = 0.D0
      KX = KX-1
      KXX = KXX-1
  150 Y = 0.D0
      SUM = 0.D0
      ALPH = X
      K = K-1
  160 B(2) = B(2) + ALPH*B(1)
  170 DO 180  L=2,K
  180 B(L+1) = B(L+1) + ALPH*B(L) - SUM*B(L-1)
  190 ROOTI(K2) = Y
      ROOTR(K2) = X
      K2 = K2+1
      IF( SUM .EQ. 0.D0 )  GO TO 200
      Y = -Y
      SUM = 0.D0
      GO TO 190
  200 IF (K .LE. 0 )  RETURN
      GO TO 20
C
  300 IER = 1
      RETURN
      E N D
      SUBROUTINE  AXISXY(X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC,IBOX,IX,IY)
C
C  ...  This subroutine draws X and Y axes.
C
C     Inputs:
C        X,Y:     location of the left bottom of the figure
C        WX,WY:   width and height of the figure
C        X0,X1:   lower and upper bounds of the X axis
C        Y0,Y1:   lower and upper bounds of the Y axis
C        DX,DY:   step width in X and Y axes
C        DWC:     size of characters
C        IBOX:    = 0   draw X and Y axes only
C                 = 1   draw window doundary
C        IX,IY:   number of additional click in each step
C     Modified:  8/31/90, 11/21/90
C
      REAL*8  X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC
C
      FX  = X
      FY  = Y
      FWX = WX
      FWY = WY
      FWC = DWC
      CALL  NEWPEN( 2 )
      CALL  PLOT( FX,FY,-3 )
      CALL  PLOT( 0.0,FWY,3 )
      CALL  PLOT( 0.0,0.0,2 )
      CALL  PLOT( FWX,0.0,2 )
      IF ( IBOX.EQ.1 )  THEN
        CALL  PLOT( FWX,FWY,2 )
        CALL  PLOT( 0.0,FWY,2 )
      END IF
C
C  ...  draw X axis ...
C
      WC = DWC
      NX = (X1-X0)/DX + 0.001
      NY = (Y1-Y0)/DY + 0.001
      MX = -DLOG10( DX ) + 0.001
      IF(DX.LT.1.0D0)  MX = -DLOG10( DX ) + 0.999
      IF( MX.LE.0 )  MX = -1
      D = WX*DX/(X1-X0)
      DO 10 I=1,NX+1
      XN = X0 + DX*(I-1)
      XX = D*(I-1)
      CALL  PLOT( XX,-FWC/2,3 )
      CALL  PLOT( XX,0.0,2 )
      IF( IBOX.EQ.1 )  THEN
        CALL  PLOT( XX,FWY-FWC/2,3 )
        CALL  PLOT( XX,FWY,2 )
      END IF
      LX = 0
      IF(XN.GT.0.0D0)  LX = ALOG10( XN ) + 0.001
      IF( MX.GT.0 )   XXX = WC*(LX+MX)*0.5 + WC*4/7.0
      IF( MX.LT.0 )   XXX = WC*LX*0.5      + WC*4/7.0
      CALL  NUMBER( XX-XXX,-(1.5*FWC+0.1),FWC,XN,0.0,MX )
   10 CONTINUE
      IF( IX.GT.1 )  THEN
      CALL  NEWPEN(1)
      DO 20 I=1,NX+1
      DO 20 J=1,IX
      XX = D*(I-1) + (D*J)/IX
      IF( XX.LE.WX )  THEN
        CALL  PLOT( XX,FWC/2,3 )
        CALL  PLOT( XX,0.0,2 )
        IF( IBOX.EQ.1 )  THEN
          CALL  PLOT( XX,FWY-FWC/2,3 )
          CALL  PLOT( XX,FWY,2 )
        END IF
      END IF
   20 CONTINUE
      CALL  NEWPEN( 2 )
      END IF
C
C  ...  draw Y axis  ...
C
      D = WY*DY/(Y1-Y0)
      MY = -DLOG10( DY ) + 0.99999
      IF( MY.LE.0 )  MY = -1
      DO 30 I=1,NY+1
      YN = Y0 + DY*(I-1)
      YY = D*(I-1)
      CALL  PLOT( -FWC/2,YY,3 )
      CALL  PLOT( 0.0,YY,2 )
      IF( IBOX.EQ.1 )  THEN
        CALL  PLOT( FWX-FWC/2,YY,3 )
        CALL  PLOT( FWX,YY,2 )
      END IF
      LY = 1
      IF(YN.LT. 0.0)  LY = 2
      IF(YN.GE. 1.0)  LY = ALOG10( YN )+1.00001
      IF(YN.LE.-1.0)  LY = ALOG10(-YN )+2.00001
      IF( MY.GT.0 )  YYY = WC*(LY+MY+1) + 0.20
      IF( MY.LT.0 )  YYY = WC*LY + 0.20
      CALL  NUMBER( -YYY,YY-FWC/2,FWC,YN,0.0,MY )
   30 CONTINUE
      IF( IX.GT.1 )  THEN
      CALL  NEWPEN(1)
      DO 40 I=1,NY+1
      DO 40 J=1,IY
      YY = D*(I-1) + (D*J)/IY
      IF( YY.LE.WY )  THEN
        CALL  PLOT( FWC/2,YY,3 )
        CALL  PLOT( 0.0,YY,2 )
        IF( IBOX.EQ.1 )  THEN
          CALL  PLOT( FWX-FWC/2,YY,3 )
          CALL  PLOT( FWX,YY,2 )
        END IF
      END IF
   40 CONTINUE
      CALL  NEWPEN( 2 )
      END IF
C
      RETURN
      E N D
      SUBROUTINE  HEADER( PNAME,NP,M,VNAME,VALUE )
C
C  ...  This subroutine draw title of data, program name, names
C       and values of parameters  ...
C
C     Inputs:
C        PNAME:    program name
C        NP:       number of characters of program name
C        M:        number of parameters
C        VNAME(I): name  of the I-th parameter
C        VALUE(I): value of the I-th parameter
C     Input via common:
C        TITLE:    title of the data set
C
      CHARACTER     PNAME*80, TITLE*72
      CHARACTER*8   VNAME(M)
      REAL*8        VALUE
      DIMENSION     VALUE(M), DAY(2), TIME(3)
      COMMON  /CMDATA/  TITLE
C
      CALL  NEWPEN( 2 )
      CALL  DATE( DAY )
      CALL  CLOCK( TIME,1 )
 
C
      CALL  SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 )
      CALL  SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP )
      CALL  NEWPEN( 1 )
      CALL  SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 )
      CALL  SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 )
      YY = 17.9
      DO 10 I=1,M
      IF( MOD(I,4).EQ.1 )  THEN
        YY = YY - 0.40
        X1 = 2.0
        X2 = 3.5
      ELSE
        X1 = X1 + 5.0
        X2 = X2 + 5.0
      END IF
      CALL  SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 )
   10 CALL  NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 )
C
      RETURN
      E N D
      SUBROUTINE  DELX( XMIN,XMAX,DX )
C
C  ...  This subroutine determines step width in drawing axis ...
C
C     Inputs:
C        XMIN:   the minimum value in the axis
C        XMAX:   the maximum value in the axis
C     Output:
C        DX:     step size
C
      IMPLICIT REAL*8(A-H,O-Z)
      DDX = XMAX - XMIN
      DX  = INT( DLOG10( DDX ) )
      IF( DDX.LT.1.0D0 )  DX = DX - 1
      DX = 10.0D0**DX
      IF( DDX/DX.LT.3.0D0 )  DX = DX/2.0
      IF( DDX/DX.LT.3.0D0 )  DX = DX/2.0
      RETURN
      E N D
      SUBROUTINE  DRAWY( FNAME,NC,XO,YO,Y,N,WX,WY,IPOS,ISW )
C
C  ...  This subroutine draws time series data  ...
C
C     Inputs:
C        FNAME:  title of the data set
C        NC:     number of charcters of the name
C        XO,YO:  origin of the figure
C        Y(I):   time series
C        N:      data length
C        WX,WY:  width and height of the figure
C        IPOS:   starting position of time series
C        ISW:    = 1  connect data by line
C                = 2  connect data by dashed line
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*80
      DIMENSION  Y(N)
C
      CALL  DELX( 0.0D0,DFLOAT(N),DX )
      CALL  MAXMIN( Y,N,YMIN,YMAX,DY )
      CALL  AXISXY( XO,YO,WX,WY,0.0D0,DBLE(N-1+IPOS),YMIN,YMAX,
     *              DX,DY,0.25D0,1,10,1 )
      IF( ISW.LE.1 ) CALL PLOTY (Y,N,YMIN,YMAX,WX,WY,IPOS,1 )
      IF( ISW.EQ.1 ) CALL PLOTY2(Y,N,YMIN,YMAX,WX,WY,IPOS,ISW,0.D0)
      IF( ISW.EQ.2 ) CALL PLOTY2(Y,N,YMIN,YMAX,WX,WY,IPOS,ISW,YMIN)
      CALL  NEWPEN( 1 )
      CALL  SYMBOL( 0.0,SNGL(WY)+0.2,0.3,FNAME,0.0,NC )
      RETURN
      END
      SUBROUTINE  DRAWY2( FNAME,NC,XO,YO,Y,N,WX,WY,IPOS,ISW,XMIN,
     *                    XMAX,DX,YMIN,YMAX,DY )
C
C  ...  This subroutine draws time series data  ...
C
C     Inputs:
C        FNAME:  title of the data set
C        NC:     number of charcters of the name
C        XO,YO:  origin of the figure
C        Y(I):   time series
C        N:      data length
C        WX,WY:  width and height of the figure
C        IPOS:   starting position of time series
C        ISW:    = 1  connect data by line
C                = 2  connect data by dashed line
C        XMIN,XMAX:  minimum and maximum of the X axis
C        DX:         step width of the X axis
C        YMIN,YMAX:  minimum and maximum of the Y axis
C        DY:         step width of the Y axis
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*80
      DIMENSION  Y(N)
C
      CALL AXISXY(XO,YO,WX,WY,XMIN,XMAX,YMIN,YMAX,DX,DY,0.25D0,1,10,1)
      IF( ISW.LE.1 )  CALL  PLOTY ( Y,N,YMIN,YMAX,WX,WY,IPOS,1 )
      IF( ISW.EQ.1 )  CALL  PLOTY2( Y,N,YMIN,YMAX,WX,WY,IPOS,1,0.0D0 )
      IF( ISW.EQ.2 )  CALL  PLOTY2( Y,N,YMIN,YMAX,WX,WY,IPOS,1,YMIN )
      CALL  NEWPEN( 1 )
      CALL  SYMBOL( 0.0,SNGL(WY)+0.2,0.3,FNAME,0.0,NC )
      RETURN
      E N D
      SUBROUTINE  MAXMIN( X,N,XMIN0,XMAX0,DXL )
C
C  ...  This subroutine determines the minimum, the maximum and
C       the step width  ...
C
C     Inputs:
C        X(I):    data
C        N:       data length
C     Outputs:
C        XMIN0:   the minimum bound for the figure
C        XMAX0:   the maximum bound for the figure
C        DXL:     step width
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(N)
C
      XMIN = 1.0D30
      XMAX =-1.0D30
C
      DO 10 I=1,N
      IF( X(I) .LT. XMIN )   XMIN = X(I)
   10 IF( X(I) .GT. XMAX )  XMAX = X(I)
      DX = XMAX-XMIN
      IF( DLOG10(DX) .GE. 0.0D0 )  DXL = INT( DLOG10(DX) )
      IF( DLOG10(DX) .LT. 0.0D0 )  DXL = INT( DLOG10(DX) )-1.0
      DXL = 10.0**DXL
      IF( DX/DXL.GT.6.0D0 )  DXL = DXL*2.0
      DIF = INT( DX/DXL )
      XMIN0 = INT( XMIN/DXL )*DXL
      XMAX0 = XMIN0 + DIF*DXL
      IF( XMIN0 .GT. XMIN )  XMIN0 = XMIN0 - DXL
   30 IF( XMAX0 .GE. XMAX )  GO TO 40
        XMAX0 = XMAX0 + DXL
        GO TO 30
   40 CONTINUE
C
      RETURN
      E N D
      SUBROUTINE  PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW )
C
C  ...  This subroutine draws time vs data plot of the data  ...
C
C     Inputs:
C        Y(I):   data
C        N:      data length
C        YMIN:   minimum bound for the Y axis
C        YMAX:   maximum bound for the Y axis
C        WX,WY:  width and height of the figure
C        IPOS:   starting position of the data
C        ISW:    = 1  connect data by straight line
C                = 2  connect data by dashed line
C                = 3  express data by character (character code=ISW)
C     Modified:  8/31/90
C
      REAL*8  Y(N), YMIN, YMAX, WX, WY
C
      DX = WX/(N-1+IPOS)
      DY = WY/(YMAX - YMIN)
C
      CALL  NEWPEN( 1 )
      DO 100  I=1,N
      XX = DX*(I-1+IPOS)
      YY = (Y(I) - YMIN)*DY
      IF( I.EQ.1 .AND. ISW.LT.3 ) THEN
         CALL PLOT( XX,YY,3 )
      ELSE
         IF(ISW.EQ.1)  CALL  PLOT( XX,YY,2 )
         IF(ISW.EQ.2)  CALL  DASHP( XX,YY,0.2 )
         IF(ISW.GE.3)  CALL  SYMBOL( XX,YY,0.2,ISW,0.0,-1 )
      END IF
  100 CONTINUE
C
      RETURN
      E N D
      SUBROUTINE  PLOTY2( Y,N,YMIN,YMAX,WX,WY,IOFF,ISW,YLEVEL )
C
C  ...  This subroutine draws box car graph of the data  ...
C
C     Inputs:
C        Y(I):    data
C        N:       data length
C        YMIN:    minimum bound for the Y axis
C        YMAX:    maximum bound for the Y axis
C        WX,WY:   width and height of the figure
C        IPOS:    starting position of the data
C        ISW:     = 1   draw box with thick line
C                 = 2   draw box with thin line
C        YLEVEL:  ground level for the box car
C
      REAL*8      Y(N), YMIN, YMAX, WX, WY, YLEVEL
C
      DX = WX/(N-1+IOFF)
      DY = WY/(YMAX - YMIN)
C
      CALL  NEWPEN( 3 )
      IF( N.GT.75 )  CALL  NEWPEN( 2 )
      IF( N.GT.150)  CALL  NEWPEN( 1 )
      IF( ISW.EQ.2 )  CALL  NEWPEN( 1 )
      DO 100  I=1,N
      XX = DX*(I-1+IOFF)
      YY = (Y(I) - YMIN)*DY
      Y0 = (YLEVEL-YMIN)*DY
      CALL  PLOT( XX,Y0,3 )
      CALL  PLOT( XX,YY,2 )
  100 CONTINUE
      CALL  NEWPEN( 2 )
C
      RETURN
      E N D
      SUBROUTINE  PLOTSC( X,Y,N,XMIN,XMAX,YMIN,YMAX,WX,WY,SIZE,CH )
C
C ... This subroutine draws scatter plot ...
C
C     Inputs:
C        X(I):      data for x-axis
C        Y(I):      data for y-axis
C        N:         data length
C        XMIN,XMAX: edges of x-axis
C        YMIN,YMAX: edges of y-axis
C        WX,WY:     width andhight of the figure
C        SIZE:      size of the symbol
C        CH:        character code
C
      REAL*8     X, Y, XMIN, XMAX, YMIN, YMAX, WX, WY, SIZE
      CHARACTER  CH*1
      DIMENSION  X(N), Y(N)
C
      DX = WX/(XMAX-XMIN)
      DY = WY/(YMAX-YMIN)
C
      DO 10 I=1,N
      XX = (X(I)-XMIN)*DX-SIZE/2
      YY = (Y(I)-YMIN)*DY-SIZE/2
      IF( XX.GE.0.0.AND.XX.LE.WX )  THEN
        IF( YY.GE.0.0.AND.YY.LE.WY )  THEN
          CALL SYMBOL( XX,YY,SIZE,CH,0.0,1 )
        END IF
      END IF
   10 CONTINUE
C
      RETURN
      E N D
      SUBROUTINE FOURIE( X,N,M,FC,FS )
C
C  ...  Discrete Fourier transformation by Goertzel method  ...
C
C     Inputs:
C        X(I):   data (I=1,N)
C        N:      data length
C        M:      number of Fourier components
C        FC(J):  Fourier cosine transform (J=1,M)
C        FS(J):  Fourier sine transform   (J=1,M)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION  X(N), FC(M), FS(M)
      DATA  PI/3.14159265358979D0/
C
      W = PI/(M-1)
      DO 20 I=1,M
      CI = DCOS(W*(I-1))
      SI = DSIN(W*(I-1))
      T1 = 0.0
      T2 = 0.0
      DO 10 J=N,2,-1
        T0 = 2*CI*T1 - T2 + X(J)
        T2 = T1
   10   T1 = T0
      FC(I) = CI*T1 - T2 + X(1)
   20 FS(I) = SI*T1
C
      RETURN
      E N D
      SUBROUTINE ARYULE( C,N,MAXM,SIG2,AIC,PARCOR,A,MAR )
C
C  ...  Yule-Walker method  ...
C
C     Inputs:
C        C(I):    Autocovariance function
C        N:       Data length
C        MAXM:    Highest AR order
C     Outputs:
C        SIG2(I): Innovation variance
C        AIC(I):  AIC
C        PARCOR(I):  PARCOR
C        AMIN:     AR coefficients of the best model
C        MAR:      Selected order of the model
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION  C(0:MAXM), SIG2(0:MAXM), AIC(0:MAXM)
      DIMENSION  PARCOR(MAXM), A(MAXM,MAXM)
      CONST = N*(DLOG(2*3.1415926535D0) + 1)
C
      SIG2(0) = C(0)
      AIC(0) = CONST + N*DLOG(SIG2(0)) + 2
      AICMIN = AIC(0)
      MAR = 0
C
      DO 50 M=1,MAXM
      SUM  = C(M)
      DO 10 I=1,M-1
   10 SUM  = SUM - A(I,M-1)*C(M-I)
      A(M,M) = SUM /SIG2(M-1)
      DO 20 J=1,M-1
   20 A(J,M) = A(J,M-1)-A(M,M)*A(M-J,M-1)
      SIG2(M) = SIG2(M-1)*(1.0D0-A(M,M)**2)
      AIC(M) = CONST + N*DLOG(SIG2(M)) + 2*(M+1)
      PARCOR(M) = A(M,M)
      IF( AIC(M).LT.AICMIN )  THEN
         AICMIN = AIC(M)
         MAR = M
      END IF
   50 CONTINUE
      RETURN
      E N D