!     Last change:  J     8 Oct 2000   11:35 am
!                                       *** copyright 2000 ***
! *** filename chex45.f95               *** John F. Monahan **
!                                       **********************
      program chex45
!   checking Example 4.5
!  test of Durbin-Levinson algorithm for solving Toeplitz
!   systems of equations and the Yule-Walker Equations
!
      implicit none
      INTEGER, PARAMETER :: n = 4
      REAL, DIMENSION(n) :: r,b,x,one,xcopy,rshft
      REAL, DIMENSION(n,n) :: a
      REAL ph,th,det,s11,soz,szz,muhat,qform
      integer i,j,idet
!
!                   INTERFACE BLOCK
      INTERFACE
      SUBROUTINE CHLSKY(A,N,DET,IDET)
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: A
      REAL, INTENT(OUT) :: DET
      INTEGER, INTENT(IN) :: N
      INTEGER, INTENT(OUT) :: IDET
      END SUBROUTINE CHLSKY
      SUBROUTINE CHLSHI(A,N,Y)
      REAL, DIMENSION(:,:), INTENT(IN) :: A
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION(:), INTENT(IN OUT) :: Y
      END SUBROUTINE CHLSHI
      SUBROUTINE CHLSIH(A,N,Y)
      REAL, DIMENSION(:,:), INTENT(IN) :: A
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION(:), INTENT(IN OUT) :: Y
      END SUBROUTINE CHLSIH
      END INTERFACE
!
  21  format(' n =',i4,'  phi, theta are ',2f5.2)
  22  format(10f8.4)
  23  format(' det ',f12.6,'   idet',i6)
  24  format(' muhat ',f12.6,'  qform ',f12.6)
  25  format(i4,2f5.2)
  26  format(/' two vectors:  b, inv(A)y, then 1''inv(A) 1')
  27  format(/' correlation matrix  -- A')
  28  format(/'  Cholesky factor of A')
!
      DATA ph/ .9/, th/ .5/
      data x/ -2., 1., 0., 1. /
      data one/ 4*1. /
      data xcopy/ -2., 1., 0., 1. /
!
      OPEN( UNIT=6, FILE='chex45.out' )
!
!                   sample size and parameters
      WRITE(6,21) n,ph,th
!                   first compute correlations
      r(1) = (1.-ph*th)*(ph-th)/(1.+th*th-2.*ph*th)
!                   this is arma(1,1)
      do i = 2,n
      r(i) = ph*r(i-1)
      end do   ! loop on i
      write(6,22) (r(i),i=1,n)
!                   compute using levdrb - Levinson-Durbin algorithm
      call levdrb(n,r,b,x,s11,det,idet)
      write(6,23) det,idet
      write(6,26) 
      write(6,22) (b(i),i=1,n)
      write(6,22) (x(i),i=1,n)
      write(6,22) s11
!                   compute mle mean and quadratic form
      soz = 0.
      szz = 0.
      do i = 1,n
      soz = soz + x(i)
      szz = szz + xcopy(i)*x(i)
      END DO  ! loop on i
!                   muhat is ML mean estimate
      muhat = soz/s11
!                   qform is quadratic form for likelihood -- since
!                    we've used correlations, rescale by gam(0)
      qform = szz - muhat*soz
      write(6,24) muhat,qform
!
!                   compute the straightforward way
!
!                form the correlation matrix
      write(6,27)
      do i = 1,n
      do j = 1,n
      a(i,j) = 1.
      if( i .ne. j ) a(i,j) = r( abs(i-j) )
      end do  ! loop on j
      rshft(i) = - r(i)
!                   write out matrix
      write(6,22) (a(i,j),j=1,n)
      END DO  ! loop on i
!                   Cholesky factorization
      call chlsky(a,n,det,idet)
      write(6,23) det,idet
!                   write out Cholesky factor
      write(6,28)
      do i = 1,n
      write(6,22) (a(i,j),j=1,i)
      END DO  ! loop on i
!                   solve three systems in rshft, y and one
      write(6,26)
      call chlshi(a,n,rshft)
      call chlsih(a,n,rshft)
!                   write solution to Yule-Walker eqns
      write(6,22) (rshft(i),i=1,n)
      call chlshi(a,n,xcopy)
      call chlsih(a,n,xcopy)
!                   solution with x as rhs
      write(6,22) (xcopy(i),i=1,n)
!                   quadratic form in one
      call chlshi(a,n,one)
      s11 = 0.
      do i = 1,n
      s11 = s11 + one(i)*one(i)
      end do
      write(6,22) s11
      stop
      end program chex45
      subroutine levdrb(n,r,b,y,a11,det,idet)
!  solve Toeplitz systems A*b = -r and A*x = y using Levinson-Durbin
!  algorithm -- first one is Yule-Walker equations
!  A is correlation matrix -- A(i,j) = r( |i-j| ) -- Toeplitz
!
!  n    in   integer       size of problem/matrix A
!  r    in   real vector   correlations r(1),...,r(n)    r(0)=1
!  b    out  real vector   solution to Yule-Walker equations A*b = -r
!  y  in/out real vector   right hand side (in) / solution (out) A*x=y
!  a11  out  real          quadratic form in inverse(A), vector of 1's
!  det  out  real          determinant in det*(2**idet) form
!  idet out  integer
!
! J F Monahan (June,1992) Dept of Statistics, NCSU, Raleigh, 27695-8203
!             (Sept 1999, April 2000) rewritten for Fortran 95
!
      implicit none
      INTEGER, INTENT(IN) :: n
      REAL, DIMENSION( n ), INTENT(IN) :: r
      REAL, DIMENSION( n ), INTENT(IN OUT) :: b,y
      REAL, INTENT(OUT) :: a11,det
      INTEGER, INTENT(OUT) :: idet
      REAL d,e,f,bi,bic
      integer k,kh,km1,i,ic
!
      det = 1.
      idet = 0
!                   begin with k=1
      b(1) = - r(1)
      a11 = 1.
      if( n .eq. 1 ) return
      d = 1. - b(1)*b(1)
!                   loop on k
      do k = 2,n
      km1 = k - 1
      e = r(k)
      f = y(k)
      bi = 1.
      do i = 1,km1
      ic = k-i
      e = e + r(i)*b(ic)
      f = f - y(i)*r(ic)
      bi = bi + b(i)
      end do  ! loop on i
!                   compute the new elements at the bottom
      det = d * det
      call adjust(det,idet)
      b(k)= - e/d
      y(k) = f/d
      a11 = a11 + bi*bi/d
      d = d * (1. - b(k)*b(k))
!                   update remaining elements - but treat b special
      do i = 1,km1
      ic = k-i
      y(i) = y(i) + y(k)*b(ic)
      end do  ! loop on i
!                   careful with b so do not overwrite incorrectly
      kh = k / 2
      do i = 1,kh
      ic = k - i
      bi = b(i) + b(k)*b(ic)
      bic = b(ic) + b(k)*b(i)
      b(i) = bi
      b(ic) = bic
      END do  ! loop on i
!             end of loop on k
      END DO  ! loop on k
      return
      end subroutine levdrb
      SUBROUTINE ADJUST(D,I)
!  NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I)
      INTEGER, PARAMETER :: IBIG = -2147483644
      REAL, INTENT (IN OUT) :: D
      INTEGER, INTENT(IN OUT) :: I
      IF ( D .GT. 0.0 ) THEN
        DO WHILE ( D .LT. 1.0 )
        D = D*16.
        I = I-4
        ENDDO
        DO WHILE ( D .GT. 16.0 )
        D = D/16.
        I = I+4
        ENDDO
                        ELSE
      I = IBIG      ! IF D < 0 THEN I = - BIG
                        ENDIF
      RETURN
      END SUBROUTINE ADJUST
      SUBROUTINE CHLSKY(A,N,DET,IDET)
!  CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT)  FACTORIZATION
!   A = L * ( L - TRANSPOSE )   WHERE L IS LOWER TRIANGULAR
!     A IS OVERWRITTEN BY L IN ITS LOWER TRIANGULAR PART
!  SUBROUTINE ADJUST KEEPS DETERMINANT FROM EXPLODING USING
!  1 LE DET LE 16   AND   DETERMINANT OF A IS DET*2**IDET
!
!  J F MONAHAN (DEC,1983) DEPT OF STAT, N C S U, RALEIGH, N C 27650 USA
!      REVISED JANUARY, 2000 FOR FORTRAN 95
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: A
      REAL, INTENT(OUT) :: DET
      INTEGER, INTENT(IN) :: N
      INTEGER, INTENT(OUT) :: IDET
      REAL T,S
      INTEGER I,J,K,IM1,JM1
!                                    INITIALIZE FOR DETERMINANT
      DET = 1.                                                        
      IDET = 0                                                          
!                                      DO I-TH ROW
      DO I = 1,N
      T = A(I,I)                                                        
!                                      FIRST ROW IS TRIVIAL
      IF(I .GT. 1) THEN
        IM1 = I-1
        DO J=1,IM1
        A(J,I) = 0.
!                                      WORK ON (I,J)-TH ELEMENT
        S = A(I,J)
        IF(J .GT. 1) THEN
          JM1 = J-1
          DO K = 1,JM1
          S = S - A(J,K) * A(I,K)
          ENDDO     ! LOOP ON K
                     ENDIF
        A(I,J)= S / A(J,J)
!                                 WORK ON DIAGONAL ELEMENT
        T = T - A(I,J) * A(I,J)
        ENDDO       ! LOOP ON J
!                                 FINISHED WITH LOOP ON J
!                                 UPDATE DETERMINANT WITH DIAGONAL
                   ENDIF
      DET = DET * T
!                                 KEEP DET FROM EXPLODING
      CALL ADJUST(DET,IDET)                                             
!                                 ABORT IF DET IS NEGATIVE OR ZERO
      IF(IDET.LT.-2000000000) RETURN                                    
!                                 FINISH A POSITIVE DIAGONAL
      A(I,I) = SQRT(T)                                                  
      ENDDO         ! LOOP ON I
      RETURN                                                            
      END SUBROUTINE CHLSKY
      SUBROUTINE CHLSHI(A,N,Y)                                          
!
!   OVERWRITES   Y   WITH      INVERSE( L ) * Y
!    WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A
!
!   TO COMPUTE QUADRATIC FORM IN THE INVERSE OF A POSITIVE DEFINITE
!   MATRIX A, THAT IS   Y' INVERSE( A ) Y THEN
!     FIRST:  FACTOR A BY CHLSKY(A,N,DET,IDET)
!             THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L'
!     SECOND:  COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y)
!     THIRD:  COMPUTE THE SUM OF SQUARES OF THE ELEMENTS OF Y
!
!  J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA
!      REVISED JUNE, 1999 FOR FORTRAN 95
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN) :: A
      REAL, DIMENSION(:), INTENT(IN OUT) :: Y
      INTEGER, INTENT(IN) :: N
      INTEGER I,J,IM1
!
!  SINCE THE MATRIX IS LOWER TRIANGULAR, SOLVE FROM THE TOP DOWN
!
      DO I = 1,N
!                            BRANCH AROUND ON THE FIRST ONE
      IF( I .GT. 1 ) THEN
        IM1 = I - 1
        DO J = 1,IM1
        Y(I) = Y(I) - A(I,J)*Y(J)
        ENDDO
                     ENDIF
      Y(I) = Y(I) / A(I,I)
      ENDDO         ! LOOP ON I
      RETURN                                                            
      END SUBROUTINE CHLSHI
      SUBROUTINE CHLSIH(A,N,Y)                                          
!                              ***  NOTICE TRANSPOSE  ***
!   OVERWRITES   Y   WITH      INVERSE( L' ) * Y
!    WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A
!
!   TO COMPUTE THE PRODUCT OF THE INVERSE OF A POSITIVE DEFINITE MATRIX
!   AND A VECTOR Y, THAT IS  INVERSE( A ) * Y  THEN
!     FIRST:  FACTOR A BY CHLSKY(A,N,DET,IDET)
!             THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L'
!     SECOND:  COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y)
!     THIRD:  COMPUTE INVERSE(L') * ( INVERSE(L)*Y ) BY CHLSIH(A,N,Y)
!
!  J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA
!      REVISED JUNE, 1999 FOR FORTRAN 95
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN) :: A
      REAL, DIMENSION(:), INTENT(IN OUT) :: Y
      INTEGER, INTENT(IN) :: N
      INTEGER I,II,IP1,J
!
!  SINCE THE MATRIX IS UPPER TRIANGULAR, SOLVE FROM THE BOTTOM UP
!
!                      DO THE LAST ONE SEPARATELY
      Y(N) = Y(N) / A(N,N)                                              
!                      IF N = 1 THEN QUIT ELSE START LOOPING
      IF( N .EQ. 1 ) RETURN                                             
      DO II = 2,N
!                      I BELOW NOW GOES FROM N-1 DOWN TO 1
        I = N + 1 - II
        IP1 = I + 1
        DO J = IP1,N
        Y(I) = Y(I) - A(J,I)*Y(J)
        ENDDO       ! LOOP ON J
        Y(I) = Y(I) / A(I,I)
      ENDDO         ! LOOP ON II,I
      RETURN                                                            
      END SUBROUTINE CHLSIH
! *** end of filename chex45.f95        *********************