! Last change: J 8 Oct 2000 11:25 am ! *** copyright 2000 *** ! *** filename chex44.f95 *** John F. Monahan ** ! ********************** program chex44 ! Check Examples 4.4 ! Compute likelihood information for arma(1,1) time series model ! using banded factorization of MA cov matrix ! implicit none INTEGER, PARAMETER :: n = 4 REAL, DIMENSION(n) :: z,bz,bo REAL, DIMENSION(n,n) :: bab REAL det,soo,soz,szz,muhat,qform REAL, PARAMETER :: phi = .9 REAL, PARAMETER :: th = .5 integer i,j,idet ! data is as short as can be data z/ -2., 1., 0., 1./ ! 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(/' three vectors, z, Bz, B1 ') 27 format(/' covariance matrix -- B A BT ') 28 format(/' Cholesky factor of B A BT ') 29 format(' first inv(L) B z, then inv(L) B one ') ! 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 END INTERFACE ! OPEN( UNIT=6, FILE='chex44.out') ! write(6,27) ! compute bz=B*z and bo=B*1 bz(1) = z(1) bo(1) = 1. do i = 2,n bz(i) = z(i) - phi*z(i-1) bo(i) = 1. - phi end do ! loop on i ! fill covariance matrix do i = 1,n do j = 1,n bab(i,j) = 0. ! fill with zeros END do ! loop on j ! except bab(i,i) = 1. + th*th ! main diagonal IF( i .gt. 1 ) bab(i,i - 1) = - th ! super and IF( i .gt. 1 ) bab(i - 1,i) = - th ! subdiagonals end do ! loop on i ! fix 1,1 element bab(1,1) = (1. + th*th - 2.*phi*th)/(1.-phi*phi) ! write out matrix DO i = 1,n write(6,22) (bab(i,j),j=1,i) END do ! loop on i ! ! write z, Bz, B1 write(6,26) write(6,22) (z(i),i=1,n) write(6,22) (bz(i),i=1,n) write(6,22) (bo(i),i=1,n) ! cholesky factorization call chlsky(bab,n,det,idet) write(6,23) det,idet ! write out cholesky factor write(6,28) do i = 1,n write(6,22) (bab(i,j),j=1,i) end do ! loop on i ! solve two systems in bz and bo write(6,29) call chlshi(bab,n,bz) write(6,22) (bz(i),i=1,n) call chlshi(bab,n,bo) write(6,22) (bo(i),i=1,n) ! get muhat and quadratic form soo = 0. soz = 0. szz = 0. do i = 1,n soo = soo + bo(i)*bo(i) soz = soz + bo(i)*bz(i) szz = szz + bz(i)*bz(i) end do ! loop on i ! form two quantities of interest muhat = soz / soo qform = szz - soz*muhat write(6,24) muhat,qform ! stop END program chex44 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 ! *** end of filename chex44.f95 *********************