C     PROGRAM 13.2  TVAR
C
C  ...  Time Varying AR model  ...
C
C     Inputs:
C        M:       AR Order (M =< 10)
C        K:       Trend Order
C        NOBS:    Local Stationary Span (N/NOBS =< NDIM)
C        IOPT:    Search Method
C        NOUT:    Number of Outliers
C          LOUT(I):  Position of I-th Outlier
C     Parameters:
C        NMAX:    Adjustable dimension of Y
C        MJ:      Adjustable dimension of XF, VF, etc.
C        NDIM:    Adjustable dimension OF VFS, VSS, etc.
C     @TEST.FILTER2:  SEP.08,1990
C
      PARAMETER( NMAX=3000,MJ=20,NDIM=200 )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Y(NMAX)
      DIMENSION  A(2,MJ), MM(MJ), LOUT(10), LLOUT(NDIM)
      DIMENSION  XPS(MJ,NDIM), XFS(MJ,NDIM), XSS(MJ,NDIM)
      DIMENSION  VPS(MJ,MJ,NDIM), VFS(MJ,MJ,NDIM), VSS(MJ,MJ,NDIM)
      DIMENSION  XF(MJ), VF(MJ,MJ)
      DATA  OUTMIN, OUTMAX /-1.0D30, 1.0D30/
      DATA  LLOUT /NDIM*0/
C
      READ( 5,* )  M, K, NOBS, IOPT, NOUT
      IF( IOPT.EQ.1 )  READ(5,*)  TAU20, DELTA
      IF( NOUT.GT.0 )  THEN
        READ(5,*)  (LOUT(I),I=1,NOUT)
        DO 10 I=1,NOUT
        J = LOUT(I)/NOBS
        IF( J*NOBS-LOUT(I).GT.NOBS/2 )  J = J+1
   10   LLOUT(J) = 1
      END IF
      INUM = 19
      IF( IOPT.EQ.0 )  INUM = 9
C
C  ...  Read Time Series  ...
C
      CALL  READTS( 1,Y,N )
C
      FMAX  = -1.0D30
C
      CALL  SETCAR( M,K,A,MM )
      DO 100  II=1,INUM
      IF( IOPT.NE.0 )  TAU2 = TAU20 + DELTA*(II-9)
      IF( IOPT.EQ.0 .AND. K.EQ.1 )  TAU2 =10.0D0**(-II)
      IF( IOPT.EQ.0 .AND. K.GT.1 )  TAU2 =10.0D0**(-II-1)
      CALL  ISTCAR( M,K,MJ,XF,VF )
C
C  ...  Log-Likelihood Computation  ...
C
      CALL  FILTR2( Y,XF,VF,TAU2,M,K,N,NOBS,MJ,1,LLOUT,
     *              OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FLK,SIG2 )
      IF( FLK.GT.FMAX )  THEN
         FMAX  = FLK
         TAUMAX = TAU2
         SIG2M  = SIG2
      END IF
  100 WRITE(6,600)  TAU2, SIG2, FLK
      AIC = -2*FMAX + 2*(M+2)
      WRITE(6,*)  M, K
      WRITE(6,610)  TAUMAX, SIG2M, FMAX, AIC
      WRITE(6,*)  FF,SIG2
C     STOP
C
C  ... Fixed Interval Smoother  ...
C
      CALL  ISTCAR( M,K,MJ,XF,VF )
      CALL  FILTR2( Y,XF,VF,TAUMAX,M,K,N,NOBS,MJ,NDIM,LLOUT,
     *              OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 )
      NN = N/NOBS
      CALL  SMOTH1( A,MM,K,M,1,NN,NN,MJ,VFS,VPS,VSS,XFS,XPS,XSS )
C
C  ...  Plot PARCOR and print AR Coefficients  ...
C
C      CALL  PTCAR( XSS,XPS,N,NOBS,MJ,M,K,TAUMAX,SIG2,FF,AIC )
      CALL  PRCAR( XSS,N,NOBS,MJ,M,K,TAUMAX,SIG2,FF,AIC )
C
      STOP
  600 FORMAT( 1H ,5X,F15.10,F12.6,F13.5 )
  610 FORMAT( 1H ,'TAUMAX =',F15.10,3X,'SIG2 =',F15.10,3X,
     *            'FF =',F13.5,3X,'AIC =',F13.4 )
      E N D
      SUBROUTINE  PRCAR( XSS,N,NOBS,MJ,M,K,TAU2,SIG2,FF,AIC )
C
C  ...  Print out Estimates  ...
C
C     Inputs:
C        XSS:    Smoothed state
C        N:      Data Length
C        NOBS:   Local Stationary Span
C        MJ:     Adjustable dimension of XSS
C        M,K:    AR Order and Trend Order
C        TAU2:   Variance of the system noise
C        SIG2:   Variance of the observational noise
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  XSS(MJ,*)
      COMMON  /CMDATA/  TITLE
C
      NN = N/NOBS
      WRITE(6,600)
      WRITE(6,610)  TITLE
      WRITE(6,620)  M, K, TAU2, SIG2, FF, AIC
      WRITE(6,630)  NN, N, NOBS
      WRITE(6,640)
      DO 10 I=1,NN
   10 WRITE(6,650)  I, (XSS(K*(J-1)+1,I),J=1,M)
      RETURN
C
  600 FORMAT( 1H1,'PROGRAM 9.1:  PREDICTION AND INTERPOLATION' )
  610 FORMAT( 1H ,'DATA: ',A72 )
  620 FORMAT( 1H ,'M =',I3,3X,'K =',I2,3X,'TAU2 =',D12.5,3X,'SIG2 =',
     *            D12.5,3X,'LOG-L =',F12.4,3X,'AIC =',F12.4 )
  630 FORMAT( 1H ,'NN =',I3,3X,'N =',I4,3X,'NOBS =',I3 )
  640 FORMAT( 1H ,'***  TIME VARYING AR COEFFOCOENTS ***' )
  650 FORMAT( 1H ,I5,(6F12.6) )
      E N D
      SUBROUTINE  SETCAR( M,K,A,MM )
C
C  ...  State space model for trend estimation  ...
C
C     Input:
C       M:   AR Order
C       K:   Trend Order
C     Outputs:
C       A:   Parameter of the Transition Matrix F
C       MM:  Dimension of the Component Model
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  A(2,*), MM(*)
C
      IF( K.EQ.1 )  THEN
        DO 10 I=1,M
   10   A(1,I) = 1.0D0
      END IF
      IF( K.EQ.2 )  THEN
        DO 20 I=1,M
        A(1,I) = 2.0D0
   20   A(2,I) =-1.0D0
      END IF
      DO 30 I=1,M
   30 MM(I) = K
C
      RETURN
      E N D
      SUBROUTINE  FILTR2( Y,XF,VF,Q,M,K,N,NOBS,MJ,NDIM,LLOUT,
     *                    OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 )
C
C  ...  Kalman Filter (for Time Varying AR model)  ...
C
C     Inputs:
C        Y:      time series
C        XF:     Initial state vector
C        VF:     Initial covariance matrix
C        Q:      K*K matrix, system noise covariance
C        M:      Dimension of the state vector
C        K:      Dimension of the system noise
C        N:      Data length
C        NOBS:   Local Stationary Span ( N/NOBS =< NDIM )
C        MJ:     Adjustable dimension of XF, VF
C        NDIM:   Adjustable dimension of VPS, VFS etc.
C        OUTMIN: Lower limit for detecting outliers
C        OUTMAX: Upper limit for detecting outliers
C     Outputs:
C        VFS:    Covariance matrices of the filter
C        VPS:    Covariance matrices of the predictor
C        XFS:    Mean vectors of the filter
C        XPS:    Mean vectors of the predictor
C        FF:     Log likelihood
C        SIG2:   Estimated observational noise variance
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Y(N), LLOUT(NDIM)
      DIMENSION  XF(MJ), VF(MJ,MJ), XP(40), VP(40,40)
      DIMENSION  XFS(MJ,NDIM),    XPS(MJ,NDIM)
      DIMENSION  VFS(MJ,MJ,NDIM), VPS(MJ,MJ,NDIM)
      DIMENSION  VH(40), GAIN(40)
      DATA   PI  /3.1415926535D0/
C
      MK = M*K
      SIG2 = 0.0D0
      SDET = 0.0D0
      NSUM = 0
      NS = M/NOBS+1
      NE = N/NOBS
C
      DO 500  II=NS,NE
C
C  ...  ONE STEP AHEAD PREDICTION  ...
C
      IF( K.EQ.1 )  THEN
         DO 10 J=1,M
         XP(J) = XF(J)
         DO 10 I=1,M
   10    VP(I,J) = VF(I,J)
      END IF
      IF( K.EQ.2 )  THEN
         DO 20 J=1,M
         J2 = 2*J
         XP(J2-1) = 2*XF(J2-1) + XF(J2)
         XP(J2)   =  -XF(J2-1)
         DO 20 I=1,M
         I2 = 2*I
         VP(I2-1,J2-1) = 4*VF(I2-1,J2-1) + 2*VF(I2,J2-1)
     *                 + 2*VF(I2-1,J2)   + VF(I2,J2)
         VP(I2-1,J2)   =-2*VF(I2-1,J2-1) - VF(I2,J2-1)
         VP(I2,J2-1)   =-2*VF(I2-1,J2-1) - VF(I2-1,J2)
   20    VP(I2,J2)     =   VF(I2-1,J2-1)
      END IF
C
      DO 30 I=1,M
   30 VP(K*(I-1)+1,K*(I-1)+1) = VP(K*(I-1)+1,K*(I-1)+1) + Q
      IF( LLOUT(II).EQ.1 )  THEN
        DO 35 I=1,MK
   35   VP(I,I) = 1.0D3
      END IF
C
C  ...  SAVE MEAN AND COVARIANCE  ...
C
      IF( NDIM.GT.1 )  THEN
      DO 40  I=1,MK
      XPS(I,II) = XP(I)
      DO 40  J=1,MK
   40 VPS(I,J,II) = VP(I,J)
      END IF
C
C
C  ...  FILTERING  ...
C
      DO 400  JJ=1,NOBS
      I1 = NOBS*(II-1) + JJ
      IF( I1.LE.M )  GO TO 400
      IF( Y(I1).GT.OUTMIN.AND.Y(I1).LT.OUTMAX.AND. I1.LE.N ) THEN
C
      DO 210  I=1,MK
      SUM = 0.0D0
      DO 200  J=1,M
  200 SUM = SUM + VP(I,K*(J-1)+1)*Y(I1-J)
  210 VH(I) = SUM
C
      PERR = Y(I1)
      PVAR = 1.0D0
      DO 220  I=1,M
      PERR = PERR - Y(I1-I)*XP(K*(I-1)+1)
  220 PVAR = PVAR + Y(I1-I)*VH(K*(I-1)+1)
C
      DO 250  I=1,MK
  250 GAIN(I) = VH(I)/PVAR
C
      DO 290  I=1,MK
  290 XF(I) = XP(I) + GAIN(I)*PERR
C
      DO 310  I=1,MK
      DO 310  J=1,MK
  310 VF(I,J) = VP(I,J) - GAIN(I)*VH(J)
C
      IF( I1.GT.MJ )  THEN
        SIG2 = SIG2 + PERR**2/PVAR
        SDET = SDET + DLOG(PVAR)
        NSUM = NSUM + 1
      END IF
C
C  ...  MISSING OBSERVATION  ...
C
      ELSE
      DO 350  I=1,MK
      XF(I) = XP(I)
      DO 350  J=1,MK
  350 VF(I,J) = VP(I,J)
      END IF
      IF( JJ.NE.NOBS )  THEN
        DO 360 J=1,MK
        XP(J) = XF(J)
        DO 360 I=1,MK
  360   VP(I,J) = VF(I,J)
      END IF
  400 CONTINUE
C
C  ...  SAVE MEAN AND COVARIANCE  ...
C
      IF( NDIM.GT.1 )  THEN
      DO 370  I=1,MK
      XFS(I,II) = XF(I)
      DO 370  J=1,MK
  370 VFS(I,J,II) = VF(I,J)
      END IF
C
  500 CONTINUE
      SIG2 = SIG2/NSUM
      FF = -0.5D0*(NSUM*(DLOG(PI*2*SIG2) + 1) + SDET)
C
      DO 510 II=1,NS-1
      DO 510 J=1,MK
      XPS(J,II) = XPS(J,NS)
      XFS(J,II) = XFS(J,NS)
      DO 510 I=1,MK
      VPS(I,J,II) = VPS(I,J,NS)
  510 VFS(I,J,II) = VFS(I,J,NS)
C
      RETURN
      E N D
      SUBROUTINE  ISTCAR( M,K,MJ,XF,VF )
C
C  ...  Initial State (for Time-Varying AR Model)  ...
C
C     Inputs:
C        M:     AR order
C        K:     Trend order
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  XF(MJ), VF(MJ,MJ)
C
      MK = M*K
      DO 10  J=1,MK
      XF(J) = 0.0D0
      DO 10  I=1,MK
   10 VF(I,J) = 0.0D0
C
      DO 20  I=1,MK
   20 VF(I,I) = 1.0D2
C
      RETURN
      E N D
      SUBROUTINE  SMOTH1( A,M,MMAX,NC,NS,N,NE,MJ,VFS,VPS,VSS,
     *                    XFS,XPS,XSS )
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  A(MMAX,NC), M(NC)
      DIMENSION  XS(40), VS(40,40), VP(40,40)
      DIMENSION  XFS(MJ,N), XPS(MJ,N), XSS(MJ,N)
      DIMENSION  VFS(MJ,MJ,N), VPS(MJ,MJ,N), VSS(MJ,MJ,N)
      DIMENSION  WRK(40,40), SGAIN(40,40)
      DIMENSION  I0(10)
C
      I0(1) = 0
      DO 10 I=2,NC
   10 I0(I) = I0(I-1) + M(I-1)
      MM = I0(NC) + M(NC)
C
C  ...  SMOOTHING  ...
C
      NSS = MIN0( N,NE )
      DO 20  I=1,MM
      XS(I)      = XFS(I,NSS)
      XSS(I,NSS) = XFS(I,NSS)
      DO 20  J=1,MM
      VS(I,J)      = VFS(I,J,NSS)
   20 VSS(I,J,NSS) = VFS(I,J,NSS)
      DO 30 II=NSS+1,NE
      DO 30 I=1,MM
      XSS(I,II) = XFS(I,II)
      DO 30 J=1,MM
   30 VSS(I,J,II) = VFS(I,J,II)
C
      DO 500  II=NSS-1,NS,-1
C
      NZERO = 0
      DO 100 I=1,MM
  100 IF( VFS(I,I,II).GT.1.0D-12 )  NZERO = NZERO + 1
C
      IF( NZERO.EQ.0 )  THEN
         DO 110  I=1,MM
         XS(I)     = XFS(I,II)
         XSS(I,II) = XFS(I,II)
         DO 110  J=1,MM
         VS(I,J)     = VFS(I,J,II)
  110    VSS(I,J,II) = VFS(I,J,II)
C
      ELSE
      DO 410  I=1,MM
      DO 410  J=1,MM
  410 VP(I,J) = VPS(I,J,II+1)
C
      CALL  GINVRS( VP,VDET,MM,40 )
C
      DO 420  I=1,MM
      DO 420  L=1,NC
      WRK(I,I0(L)+M(L)) = VFS(I,I0(L)+1,II)*A(M(L),L)
      DO 420  J=1,M(L)-1
  420 WRK(I,I0(L)+J) = VFS(I,I0(L)+1,II)*A(J,L) + VFS(I,I0(L)+J+1,II)
C
      DO 440  I=1,MM
      DO 440  J=1,MM
      SUM = 0.0D0
      DO 430 IJ=1,MM
  430 SUM = SUM + WRK(I,IJ)*VP(IJ,J)
  440 SGAIN(I,J) = SUM
C
      DO 450  I=1,MM
      XS(I) = XFS(I,II)
      DO 450  J=1,MM
      WRK(I,J) = 0.0D0
  450 VS (I,J) = VFS(I,J,II)
C
      DO 460  J=1,MM
      DO 460  I=1,MM
  460 XS(I) = XS(I) + SGAIN(I,J)*(XSS(J,II+1) - XPS(J,II+1))
C
      DO 470  J=1,MM
      DO 470 IJ=1,MM
      DO 470  I=1,MM
  470 WRK(I,J) = WRK(I,J) + SGAIN(I,IJ)*(VSS(IJ,J,II+1)-VPS(IJ,J,II+1))
C
      DO 480  J=1,MM
      DO 480 IJ=1,MM
      DO 480  I=1,MM
  480 VS(I,J) = VS(I,J) + WRK(I,IJ)*SGAIN(J,IJ)
      DO 485 I=1,MM
  485 IF( VS(I,I).LT.0.0D0 )  VS(I,I) = 0.0D0
C
      DO 490  I=1,MM
      XSS(I,II) = XS(I)
      DO 490  J=1,MM
  490 VSS(I,J,II) = VS(I,J)
C     WRITE(6,*) (XS(I),I=1,MM)
C     DO 88 I=1,MM
C  88 WRITE(6,89) (VS(I,J),J=1,MM)
C  89 FORMAT( 1H ,11D12.4 )
      END IF
C
  500 CONTINUE
C
      RETURN
      E N D
      SUBROUTINE  READTS( IDEV,Y,N )
      REAL*8        Y(N)
      CHARACTER*72  TITLE
      COMMON  /CMDATA/  TITLE
C
      OPEN( IDEV,FILE='temp.dat' )
      READ( IDEV,1 )  TITLE
      READ( IDEV,* )  N
      READ( IDEV,* )  (Y(I),I=1,N)
      CLOSE( IDEV )
      RETURN
    1 FORMAT( A72 )
      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  GINVRS( A,DET,M,MJ )
C
C  ...  Generalized inverse of a square matrix A  ...
C
C     Inputs:
C        A:     M*M matrix
C        M:     Dimension of A
C        MJ:    Adjustable dimension of A
C     Outputs:
C        A:     Generalize inverse of A
C        DET:   Determinant of A
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  A(MJ,MJ), IND(50)
      EPS = 1.0D-10
C
      DO 10  I=1,M
   10 IND(I) = I
      DO 60  L=1,M
      AMAX = 0.0D0
      DO 20  I=L,M
      IF( A(IND(I),IND(I)).GT.AMAX )  THEN
        AMAX = A(IND(I),IND(I))
        I0 = I
      END IF
   20 CONTINUE
      IF( AMAX.GT.EPS*A(IND(1),IND(1)) )  THEN
         IMAX = IND(I0)
         DO 30  I=I0,L+1,-1
   30    IND(I) = IND(I-1)
         IND(L) = IMAX
         LMAX   = L
         DO 40  I=L+1,M
         A(IND(I),IMAX) = -A(IND(I),IMAX)/A(IMAX,IMAX)
         DO 40  J=L+1,M
   40    A(IND(I),IND(J)) = A(IND(I),IND(J))
     *                    + A(IND(I),IMAX)*A(IMAX,IND(J))
      ELSE
         DO 50  I=L,M
         DO 50  J=L,M
   50    A(IND(I),IND(J)) = 0.0D0
         GO TO 70
      END IF
   60 CONTINUE
C
   70 DET = 1.0D0
      DO 80  I=1,M
C     DET = DET*A(IND(I),IND(I))
C     IF( A(IND(I),IND(I)).GT.EPS*AMAX )  THEN
      IF( A(IND(I),IND(I)).GT.0.0D0 )  THEN
        A(IND(I),IND(I)) = 1.0D0/A(IND(I),IND(I))
      ELSE
        A(IND(I),IND(I)) = 0.0D0
      END IF
   80 CONTINUE
C
      MS = MIN0( M-1,LMAX )
      DO 200  L=MS,1,-1
      DO 100 J=L+1,M
      SUM = 0.0D0
      DO 90  I=L+1,M
   90 SUM = SUM + A(IND(I),IND(L))*A(IND(I),IND(J))
  100 A(IND(L),IND(J)) = SUM
      SUM = A(IND(L),IND(L))
      DO 110  I=L+1,M
  110 SUM = SUM + A(IND(L),IND(I))*A(IND(I),IND(L))
      A(IND(L),IND(L)) = SUM
      DO 120  I=L+1,M
  120 A(IND(I),IND(L)) = A(IND(L),IND(I))
      IMAX = IND(L)
      DO 130 I=L+1,M
      IF( IMAX.GT.IND(I) )  THEN
         IND(I-1) = IND(I)
         IND(I)   = IMAX
      END IF
  130 CONTINUE
  200 CONTINUE
C
      RETURN
      E N D