! Last change: DOS 3 Aug 2000 5:25 pm ! *** copyright 2000 *** ! *** filename mh2.f95 *** John F. Monahan ** ! ********************** program mh2 ! mh2 -- second trial of Metropolis-Hastings algorithm ! using log-series posterior ! and random walk, uniform trial distribution ! follow Geweke's approach for checking stationarity implicit none integer, parameter :: nn = 14 integer, parameter :: n = 2**nn integer, parameter :: na = 2**(nn-3) ! first 1/8 of series integer, parameter :: nb = 2**(nn-1) ! last 1/2 of series integer nd,i,j,kaccept complex, dimension(na) :: xa complex, dimension(nb) :: xb real pstarold,pstarnew,told,tnew,ltn,lp,alpha real fn,sxa,sxb,s0a,s0b,sxx,xi,t,ran ! 21 FORMAT(' M-H with series length',i12,' log(length)=',i4,' accepts',i12) 22 FORMAT(' First part of series, first',i10,' observations') 23 FORMAT(' Last part of series, last',i10,' observations') 24 FORMAT(' Mean and Variance estimates',2f12.6) 25 FORMAT(' bandwidth d=',i4,' spectral var',f14.8) 26 FORMAT(' Mean and spectral std err',2f12.6) 27 FORMAT(' Final t statistic',f12.6) ! OPEN( UNIT=6, FILE='mh2.out' ) ! initialize sxa = 0. sxb = 0. told = .5 pstarold = 0. ! wrong, but it shouldn't matter kaccept = 0 tnew = ran(5151917) ! main loop do i = 1,n tnew = told + (ran(i) - .5)/2.5 ! add unif( -.2, +.2 ) ! get pstar for tnew pstarnew = 0. if ( (tnew .GT. 0. ) .AND. (tnew .LT. 1.) ) then ltn = log(1.-tnew) lp = 16.*log(tnew) + ltn - 10.*log(-ltn) if ( lp > -60. ) pstarnew = exp(lp) end if ! get alpha if ( pstarnew .lt. pstarold ) then alpha = pstarnew / pstarold else alpha = 1. end if if ( ran(i) .lt. alpha ) then ! acceptance kaccept = kaccept + 1 told = tnew pstarold = pstarnew end if ! output xi = told ! first 1/8 of series in xa if ( i .le. na ) then xa(i) = xi sxa = sxa + xi end if ! last 1/2 of series in xb if ( i .gt. nb ) then xb(i-nb) = xi sxb = sxb + xi end if end do write(6,21) n,nn,kaccept ! calculations on first part write(6,22) na nd = int(sqrt(real(na))) fn = real(na) ! get mean, variance, autocorrelation sxa = sxa / fn sxx = (xa(1) - sxa)*(xa(1) - sxa) do i = 2,na xi = REAL(xa(i)) sxx = sxx + (xi - sxa)*(xi - sxa) end do sxx = sxx / (fn-1.) write(6,24) sxa,sxx ! get spectral density estimate at 0. call fft2n(xa,nn-3,1.) s0a = conjg(xa(2))*xa(2)/2. ! just to match SAS do j = 1,nd s0a = s0a + conjg(xa(j+1))*xa(j+1) end do s0a = (2./fn) * s0a/real(2*nd+1) ! (2/n) from periodogram write(6,25) nd,s0a s0a = sqrt(s0a/real(na)) ! get standard error write(6,26) sxa,s0a ! calculations on last half write(6,23) nb nd = int(sqrt(real(nb))) fn = real(nb) sxb = sxb / fn sxx = (xb(1) - sxb)*(xb(1) - sxb) do i = 2,nb xi = real(xb(i)) sxx = sxx + (xi - sxb)*(xi - sxb) end do sxx = sxx / (fn-1.) write(6,24) sxb,sxx ! get spectral density estimate at 0. call fft2n(xb,nn-1,1.) s0b = conjg(xb(2))*xb(2)/2. ! just to match SAS do j = 1,nd s0b = s0b + conjg(xb(j+1))*xb(j+1) end do s0b = (2./fn) * s0b/real(2*nd+1) ! (2/n) from periodogram write(6,25) nd,s0b s0b = sqrt(s0b/real(nb)) ! get standard error write(6,26) sxb,s0b ! results t = ( sxa - sxb ) / sqrt( s0a*s0a + s0b*s0b ) write(6,27) t end program mh2 REAL FUNCTION RAN(IDUM) ! UNIFORM PSEUDORANDOM NUMBER GENERATOR ! A LINEAR CONGRUENTIAL GENERATOR X(N+1) = MOD( A*X(N), 2**31 - 1 ) ! CODING FOLLOWS ALGORITHM 1 OF HORMANN AND DERFLINGER, WITH ! PERSONAL MODIFICATIONS TO AVOID OVERFLOWS ! ! W. HORMANN & G. DERFLINGER (1993) "A PORTABLE RANDOM NUMBER GENERATOR ! WELL SUITED FOR THE REJECTION METHOD," ACM TOMS VOL 19, PP.489-495. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM INTEGER, PARAMETER :: AHI = 12121 ! PART OF MULTIPLIER INTEGER, PARAMETER :: ALOW = 23166 ! PART OF MULTIPLIER INTEGER, PARAMETER :: ALOW2 = 46332 ! PART OF MULTIPLIER INTEGER, PARAMETER :: B15 = 32768 ! 2**15 INTEGER, PARAMETER :: B16 = 65536 ! 2**16 INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: XHI = 0 INTEGER, SAVE :: XLOW = 0 INTEGER MID1,MID2,MID,X,K ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( (XHI .EQ. 0) .AND. (XLOW .EQ. 0) ) THEN XHI = IDUM / 65536 XLOW = IDUM - XHI*65536 END IF ! ( FIRST CALL ) ! MULTIPLIER IS A = 397204094 = AHI*(2**15) + ALOW MID1 = AHI*XLOW MID2 = ALOW2*XHI ! TEST FOR OVERFLOW IF( MID1-1 .GT. P-MID2 ) THEN ! HERE MID IS > 2**31, SO WRITE AS MID - 2**31 MID = (MID1-1) - (P-MID2) K = 1 ELSE ! HERE MID IS < 2**31, SO WRITE AS POSITIVE MID = MID1 + MID2 K = 0 END IF ! ( MID1-1 .GT. P-MID2 ) ! NOW GET TO MAIN STEP ! SUBTRACT P = 2**31 - 1 IN ADVANCE MID2 = MID / B16 X = (ALOW*XLOW - P) + AHI*XHI + MID2 + K*B15 IF( X .LT. 0 ) X = X + P X = X + ( ( MID - MID2*B16 )*B15 - P ) IF( X .LT. 0 ) X = X + P ! WILL NEED TWO PARTS OF X FOR NEXT CALL XHI = X / B16 XLOW = X - XHI*B16 ! 2**16 2**15 RAN = ( REAL(XHI) + ( REAL(XLOW) / 65536. ) )/32768. RETURN END FUNCTION RAN SUBROUTINE FFT2N(A,N,SGN) ! ! SUBROUTINE FFT2N -- DISCRETE FOURIER TRANSFORM ! ! A COMPLEX VECTOR INPUT: VECTOR TO BE TRANSFORMED ! OUTPUT: TRANSFORMED VECTOR ! N INTEGER INPUT: LOG(BASE 2) OF LENGTH OF VECTOR A ! LENGTH OF A IS 2 ** N ! SGN REAL INPUT: INDICATES FORWARD OR INVERSE TRANSFO ! SGN = 1. THEN FORWARD TRANSFORM ! SGN = -1. THEN INVERSE TRANSFORM ! ! DISCRETE FOURIER TRANSFORM OF A VECTOR A(J) IS GIVEN BY ! NN-1 ! T( A ) = SUM A(J) EXP( - I SGN 2 PI J K / NN ) ! J K=0 ! ! WHERE I = SQRT( -1 ) ! NN = 2 ** N ! PI = 3.1415923565... ! ! NOTE THAT THE INDEXING OF A IS SHIFTED BY ONE -- ! A(0) STORED IN A(1), ... , A( NN-1 ) STORED IN A( NN ) ! ! J F MONAHAN (NOV,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695-8203 ! RECODED DECEMBER 1999, APRIL 2000 FOR FORTRAN 95 ! ! THIS IS A 'POWER OF 2' ALGORITHM BASED ON THE PRESENTATION IN ! ! A. V. AHO, J. E. HOPCROFT, AND J. D. ULLMAN (1974) THE DESIGN AND ! ANALYSIS OF COMPUTER ALGORITHMS, ADDISON-WESLEY, READING. ! ! ! BEGIN WITH DECLARATIONS IMPLICIT NONE INTEGER, INTENT(IN) :: N COMPLEX, DIMENSION(1:N), INTENT(IN OUT) :: A REAL, INTENT(IN) :: SGN REAL C,S COMPLEX, DIMENSION(25) :: WK COMPLEX WS,ALEFT,ARIGHT LOGICAL CARRY,NEW LOGICAL, DIMENSION(25) :: RBITS INTEGER I,J,L,LM,LL,M,MM,M2,NN ! ! COMPUTE WK(I) = EXP( - 2 PI IMAG SGN / 2**I ) ! ! SPECIAL VALUES FIRST WK(1) = CMPLX( -1., 0. ) WK(2) = CMPLX( 0., -SGN ) ! SKIP IF N IS ONLY 1 OR 2 IF( N .GT. 2 ) THEN ! AVOID SINES AND COSINES -- JUST USE SQUARE ROOTS C = 0. S = 1. DO I = 3,N C = SQRT( ( 1. + C ) / 2. ) S = S / ( 2.*C ) WK(I) = CMPLX( C, -SGN*S ) END DO ! LOOP ON I END IF ! ( N .GT. 2 ) ! ! THIS IS THE MAIN LOOP OF THE ALGORITHM ! DO MM = 1,N ! CONVERT TO REVERSE LOOP M = N - MM ! M STEPS DOWN FROM N-1 TO 0 ! ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I ! M2 M2 = 2**M ! LM CONTROLS THE NUMBER OF SUBLOOPS LM = 2**(MM-1) ! DO LL = 1,LM ! L IS THE MAJOR INDEX VARIABLE L = 2*M2*(LL-1) ! INITIALIZE TO GET W ** REV( L / M2 ) WS = CMPLX(1., 0.) CARRY = .TRUE. ! START AT 2 (SHIFT) SINCE L/M2 IS ALWAYS EVEN DO I = 2,N NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) WS = WS * WK(I) END DO ! LOOP ON I ! WS IS NOW W ** REV( L / M2 ) ! NOW DO THE REAL COMPUTATION DO J = 1,M2 ALEFT = A( L + J ) ARIGHT = WS*A( L + J + M2 ) A( L + J ) = ALEFT + ARIGHT A( L + J + M2 ) = ALEFT - ARIGHT END DO ! LOOP ON J ! DONE WITH LOOP ON L END DO ! LOOP ON LL ! DONE WITH MAIN LOOP ON M END DO ! LOOP ON MM ! ! SCRAMBLE THE NUMBERS -- SWITCH A(J) AND A( REV(J) ) ! BUT DON'T DO IT TWICE ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I NN = 2 ** N ! LOOP THROUGH ALL THE NUMBERS DO LL = 1,NN ! DO THE BIT REVERSAL LM = 0 M2 = NN CARRY = .TRUE. ! START AT 1 HERE -- USE ALL THE NUMBERS DO I = 1,N M2 = M2 / 2 NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) LM = LM + M2 END DO ! LOOP ON I ! DON'T SWITCH TWICE IF( LL-1 .LT. LM ) THEN ALEFT = A(LL) ARIGHT = A(LM+1) A(LL) = ARIGHT A(LM+1) = ALEFT END IF ! ( LL-1 .LT. LM ) END DO ! LOOP ON LL ! IF INVERSE TRANSFORM, REMEMBER TO DIVIDE BY NN IF( SGN .GT. 0. ) RETURN C = REAL(NN) DO I = 1,NN A(I) = A(I) / C END DO ! LOOP ON I RETURN END SUBROUTINE FFT2N ! *** end of filename mh2.f95 *********************