! Last change: DOS 3 Aug 2000 5:11 pm ! *** copyright 2000 *** ! *** filename chex131.f95 *** John F. Monahan ** ! ********************** program chex131 ! Example 13.1 -- Gibbs sampling from posterior distribution of ! a variance components problem -- see chex122a, chex122b, mixed2 ! ! exchange variance components problem ! data are 100*(rate-1.6) for USD/GBP rate in May 1997 ! ! f error/within variance component ! g trader/tmt variance component ! beta mean parameter ! implicit none integer, parameter :: nobs = 2500 ! number of observations integer, parameter :: p = 5 ! number of groups integer, parameter :: pp3 = p+3 ! total number of parameters integer, dimension(p) :: ni integer, parameter :: nbig = 15 ! total sample size real, dimension(p) :: yb,theta real, dimension(pp3) :: th,thmean real, dimension(pp3,pp3) :: thcov real, parameter :: w = 13.1 ! within ss real b,f,g real a1,b1,a2,b2,beta0,phi0,ft,sumth real a1new,b1new,a2new,b2new,beta0new,phi0new,gaminew,thetanw real ran,gnroul,gchirv ! integer i,j,t ! data for problem ! p number of groups ! ni number of obs in group i ! nbig sum of ni ! w within/error sum of squares ! a1,b1 prior for error component is inversegamma( a1, b1 ) ! a2,b2 prior for treatment component is inversegamma( a2, b2 ) ! b0,phi0 prior for mean in Normal( b0, phi0 ) ! yb mean for group i ! 21 format(i6,8f9.4) 28 format(/'Gibbs sampling -- Example 13.1 with',i8,' observations',& & e12.6/10x,'Posterior mean',10x,'Covariance Matrix') 29 format(' theta(1)',f12.6,6x,f12.6/ & & ' theta(2)',f12.6,6x,2f12.6/ & & ' theta(3)',f12.6,6x,3f12.6) ! data ni/ 2, 1, 7, 3, 2 / data yb/ 6.8, 10.3, 5.9, 6.6, 8.5/ ! a1 = 1. b1 = 1. a2 = 4. b2 = 4. beta0 = 8. phi0 = 16. ! ! output unit for most results open( unit=6, file='chex131.out' ) ! output unit for diagnostics open( unit=8, file='chex131.dgn' ) ! ! ! initialize sumth = ran(5151917) thmean = 0. thcov = 0. ! starting values b = 7. ! mean parameter f = 1. ! error variance component g = 1. ! trader variance component theta = yb ! trader effects ! ! main iteration loop do t = 1,nobs ft = real(t) ! ! first the distribution of phi given others a1new = a1 + nbig/2. b1new = 0. do i = 1,p b1new = b1new + ni(i)*(yb(i)-theta(i))*(yb(i)-theta(i)) end do ! loop on i b1new = b1 + w/2. + b1new/2. ! generate new phi f = gchirv(a1new*2.) ! is chi(a1new*2) f = b1new / (f*f/2.) ! f*f/2 is gamma(a1new) th(1) = f ! copy into param vector ! ! next the distribution of gamma given others a2new = a2 + p/2. b2new = 0. do i = 1,p b2new = b2new + (theta(i)-b)*(theta(i)-b) end do ! loop on i b2new = b2 + b2new/2. ! generate new gamma g = gchirv(a2new*2.) ! is chi(a2new*2) g = b2new / (g*g/2.) ! g*g/2. is gamma(a2new) th(2) = g ! copy into param vector ! ! now distribution of beta given others phi0new = 1./( 1./phi0 + float(p)/g ) sumth = 0. do i = 1,p sumth = sumth + theta(i) end do ! loop on i beta0new = (beta0/phi0 + sumth/g)* phi0new ! generate new beta b = beta0new + SQRT(phi0new)*gnroul(t) th(3) = b ! copy into param vector ! ! finally generate new theta's do i = 1,p gaminew = 1./( real(ni(i))/f + 1./g ) thetanw = ( real(ni(i))*yb(i)/f + b/g )*gaminew ! generate new theta theta(i) = thetanw + SQRT(gaminew)*gnroul(i) th(i+3) = theta(i) ! copy into param vector end do ! loop on i ! ! write out results to file for later analysis write(8,21) t,th ! ! compute statistics do i = 1,pp3 thmean(i) = thmean(i) + th(i) if( t .gt. 1 ) then do j = 1,i thcov(i,j) = thcov(i,j) + & & (ft*th(i)-thmean(i))*(ft*th(j)-thmean(j))/(ft*(ft-1.)) end do ! loop on j end if ! ( t .gt. 1 ) end do ! loop on i end do ! loop on t -- iterations ft = real(nobs) thmean = thmean / ft thcov = thcov / (ft-1.) write(6,21) nobs,thmean do i = 1,pp3 write(6,21) i,(thcov(i,j),j=1,i) end do ! loop on i write(6,28) nobs write(6,29) (thmean(i),(thcov(i,j),j=1,i),i=1,3) stop end program chex131 REAL FUNCTION GCHIRV(A) ! ! ALGORITHM TO GENERATE RANDOM VARIABLES WITH THE CHI DISTRIBUTION ! WITH A DEGREES OF FREEDOM, FOR A GREATER THAN OR EQUAL TO ONE ! ! J F MONAHAN (MAY, 1986) DEPT OF STATISTICS, NCSU, RALEIGH, NC USA ! RECODED (FEB, 2000) FOR FORTRAN 95 ! IMPLICIT NONE REAL, INTENT(IN) :: A REAL U,V,Z,ZZ,RNUM,W,S,VMAX REAL, SAVE :: ALPHM1,BETA,VMIN,VDIF ! REAL, PARAMETER :: EMHLF = .6065307 ! EXP(-1/2) REAL, PARAMETER :: RSQRT2 = .7071068 ! 1/SQRT(2) REAL, PARAMETER :: EMHLF4 = .1516327 ! EXP(-1/2)/4 REAL, PARAMETER :: EQTRT2 = 2.568051 ! 2*EXP(1/4) REAL, PARAMETER :: C = 1.036961 ! 4*EXP(-1.35) ! REAL, SAVE :: ALPHA = 0. ! GIVE INITIAL VALUE & SAVE REAL RAN ! IS THIS ALPHA THE SAME AS THE LAST ONE? IF( A .NE. ALPHA ) THEN ! DO A LITTLE SETUP ALPHA = A ALPHM1 = ALPHA - 1. BETA = SQRT( ALPHM1 ) ! GET DIMENSIONS OF BOX VMAX = EMHLF * ( RSQRT2 + BETA )/( .5 + BETA ) VMIN = -BETA IF( BETA .GT. 0.483643 ) VMIN = EMHLF4/ALPHA - EMHLF VDIF = VMAX - VMIN END IF ! ( A .NE. ALPHA) ! START ( AND RESTART ) ALGORITHM HERE DO ! NOTE UNRESTRICTED DO U = RAN(1) V = VDIF*RAN(2) + VMIN Z = V / U GCHIRV = Z + BETA ! RETURN ON SUCCESS ! DO SOME QUICK REJECT CHECKS FIRST IF( Z .LE. - BETA ) CYCLE ! REJECT ZZ = Z * Z RNUM = 2.5 - ZZ IF( V .LT. 0. ) RNUM = RNUM + ZZ * Z / ( 3. * (Z + BETA ) ) ! DO QUICK INNER CHECK IF( U .LT. RNUM/EQTRT2 ) RETURN ! ACCEPT IF( ZZ .GT. C / U + 1.4 ) CYCLE ! REJECT ! ABOVE WAS KNUTH'S NORMAL OUTER CHECK W = 2. * LOG( U ) ! NOW THE REAL CHECK S = - ( ZZ / 2. + Z * BETA ) IF( BETA .GT. 0. ) S = ALPHM1*LOG(1.+Z/BETA) + S IF( W .LE. S ) RETURN ! ACCEPT END DO ! UNRESTRICTED DO RETURN END FUNCTION GCHIRV REAL FUNCTION GNROUL(IX) ! RATIO OF UNIFORMS ALGORITHM FOR GENERATING NORMAL(0,1) RV'S ! USING LEVA'S QUADRATIC INNER AND OUTER BOUNDS ! ! A J KINDERMAN AND J F MONAHAN (1977) "COMPUTER GENERATION OF RANDOM ! VARIABLES USING THE RATIO OF UNIFORM DEVIATES," ACM TRANSACTIONS ON ! MATHEMATICAL SOFTWARE, VOLUME 3, PP.257-260 ! ! J L LEVA (1992) "A FAST NORMAL RANDOM NUMBER GENERATOR," ACM ! TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 18, PP. 449-453. ! IMPLICIT NONE INTEGER, INTENT(IN) :: IX ! DUMMY ARGUMENT REAL, PARAMETER :: R = 1.7155277 ! FIRST CONSTANT R = SQRT(2/E) ! CONSTANTS S,T CENTER OF ELLIPSES REAL, PARAMETER :: S = .449871 REAL, PARAMETER :: T = -.386595 ! SHAPE PARAMETERS OF ELLIPSES REAL, PARAMETER :: A = .19600 REAL, PARAMETER :: B = .25472 ! RADII OF ELLIPSES REAL, PARAMETER :: R1 = .27597 REAL, PARAMETER :: R2 = .27846 REAL RAN REAL U,V,X,Y,Q ! START / RESTART HERE DO ! NOTE UNRESTRICTED DO ! GENERATE (U,V) UNIFORMLY IN BOX U = RAN(1) V = R*(RAN(2) - 0.5) X = U - S Y = ABS(V) - T ! COMPUTE QUADRATIC PIECE Q = X*X + Y*(A*Y - B*X) ! COMPUTE DEVIATE BEFORE TESTS GNROUL = V / U ! INNER BOUND -- QUICK ACCEPT CHECK IF( Q .LE. R1 ) RETURN ! OUTER BOUND -- QUICK REJECT CHECK IF( A .GT. R2 ) CYCLE ! FINAL RATIO OF UNIFORMS CHECK IF( GNROUL*GNROUL .LE. -4.*LOG(U) ) RETURN ! REJECT -- START OVER END DO ! END OF UNRESTRICTED DO END FUNCTION GNROUL REAL FUNCTION RAN(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RAN = REAL(IX)*RP RETURN END FUNCTION RAN ! *** end of filename chex131.f95 *********************