SUBROUTINE ACFAR(NP,RX,BETA,SIGMA) C This routine finds the AR(p) representation corresponding to a positive C definite autocovariance function for lags 0 through p, by solving the C Yule-Walker equations. It is the inverse of ARACF. C Positive definiteness of the autocovariance function is not checked! C Inputs: C NP Parameter p C RX Autocovariance function C Outputs: C BETA(j) Coefficient on lag j in autoregression C SIGMA Standard deviation of innovation IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION RX(0:NP),BETA(NP) CALL DCOPY(NP,RX(0),1,V1,1) CALL DCOPY(NP-1,RX(1),1,V1(NP+1),1) CALL DLSLTO(NP,V1,RX(1),1,BETA) SIGMA=DSQRT(RX(0)-DDOT(NP,RX(1),1,BETA,1)) RETURN END SUBROUTINE ARACF(NP,BETA,SIGMA,RX) C This routine finds the autocovariance function for lags 0 through p, C corresponding to an AR(p) representation, by solving the system of p+1 C linear equations. It is the inverse of ACFAR. C Invertibility of the AR(p) representation is not checked! C Inputs: C NP Parameter p C BETA(j) Coefficient on lag j in autoregression C Outputs: C RX(j) Autocovariance at lag j C SIGMA Standard deviation of innovation IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION RX(0:NP),BETA(NP) NP1=NP+1 CALL UMISET(NP1,A1,LD) DO 10 I=1,NP1 DO 10 K=1,NP J=1+IABS(K+1-I) 10 A1(I,J)=A1(I,J)-BETA(K) CALL DSET(NP,0.0D0,V1(2),1) V1(1)=SIGMA**2 CALL DLSLRG(NP1,A1,LD,V1,1,RX) RETURN END SUBROUTINE ARPACF(NP,BETA,R) C This routine computes the p partial autocorrelation function coefficients, C given an autoregression of order p, y(t)=Beta(1)y(t-1)+...+Beta(p)y(t-p) C + Epsilon(t). C Inputs: C NP Order of autoregression p C BETA AR coefficients, BETA(j)=Beta(j) C Output: C R Vector of partial autocorrelation coefficients, in order. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION BETA(NP),R(NP) COMMON/SCRA/V1(10100),V2(10100) CALL DCOPY(NP,BETA,1,V2,1) IF(NP.EQ.1)GO TO 30 DO 20 K=NP,2,-1 RR=V2(K) AA=1.0D0/(1.0D0-RR**2) K1=K-1 DO 10 I=1,K1 10 V1(I)=(V2(I)+RR*V2(K-I))*AA 20 CALL DCOPY(K1,V1,1,V2,1) 30 CALL DCOPY(NP,V2,1,R,1) RETURN END C SUBROUTINE PACFAR(NP,R,BETA) C This routine computes the autoregression of order p corresponding to the C first p autocorrrelation function coefficients. It is the inverse of C ARPACF. C Inputs: C NP Order of autoregression p C R Vector of partial autocorrelation coefficients, in order. C Output: C BETA AR coefficients, in order. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(NP),BETA(NP) COMMON/SCRA/V1(10100),V2(10100) V1(1)=R(1) IF(NP.EQ.1)GO TO 30 DO 20 K=2,NP RR=R(K) DO 10 I=1,K-1 10 V2(I)=V1(I)-RR*V1(K-I) V2(K)=RR 20 CALL DCOPY(K,V2,1,V1,1) 30 CALL DCOPY(NP,V1,1,BETA,1) RETURN END C REAL*8 FUNCTION ARPALJ(NP,R) C This function finds the log determinant Jacobian of transformation from p C partial autocorrelation function coefficients to the p coefficients in an C AR(p) -- i.e., log{det[partial(beta)/partial(r)]}. C Inputs: C NP Order of autoregression p C R Vector of partial autocorrelation coefficients, in order. C Output: C ARPALJ Log determinant Jacobian of transformation IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(NP) AA=0.0D0 IF(NP.EQ.1)GO TO 30 AA=DLOG(1.0D0-R(2)) IF(NP.EQ.2)GO TO 30 DO 10 J=3,NP,2 10 AA=AA+((J-1)/2)*DLOG(1.0D0-R(J)**2) C 10 AA=AA+DLOG(1.0D0-R(J)**(J-1)) IF(NP.EQ.3)GO TO 30 DO 20 J=4,NP,2 20 AA=AA+((J-2)/2)*DLOG(1.0D0-R(J)**2)+DLOG(1.0D0-R(J)) C 20 AA=AA+DLOG((1.0D0-R(J))*(1.0D0-R(J)**(J-2))) 30 ARPALJ=AA RETURN END REAL*8 FUNCTION ARSDI1(ALAM1,ALAM2) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This function integrates the product of the spectral density function C of a univariate autogregressive process, and the function C cos(lambda*tau), between the frequencies lambda1 and lambda2, and C -lambda2 and -lambda1. C If lambda1=0 and lambda2=pi, this provides the autocovariance at lag C tau. C If tau=0, this provides the power of the spectrum over the indicated C frequencies. C Inputs: C ALAM1 Lambda1 C ALAM2 Lambda2 C Inputs (/AUAR/): C B Coefficients of autoregression C H Precision of shock to autoregression C TAU tau C NP Number of coefficients in autoregression C Output: C ARSDI1 Integrated value of product described above IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL FARSD1 PARAMETER(LP=100) COMMON/AUAR/B(0:LP),H,TAU,NP COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII CALL DQDAG(FARSD1,ALAM1,ALAM2,1.0D-4,1.0D-4,2,AA,ERR) ARSDI1=PII*AA RETURN END REAL*8 FUNCTION BERN(P) C This function generates a Bernoulli random variable. C Input: C P Probability parameter C Output: C BERN = 1 with probability p, 0 with probability 1-p IMPLICIT REAL*8 (A-H,O-Z) AA=0.0D0 IF(DRNUNF().LT.P)AA=1.0D0 BERN=AA RETURN END C This is a type logical version of the above, returning .TRUE. with C probability P. LOGICAL FUNCTION LBERN(P) IMPLICIT REAL*8 (A-H,O-Z) LBERN=DRNUNF().LT.P RETURN END SUBROUTINE BET(N,R,P) C This routine generates one realization from a mutlivariate beta, or C Dirichlet, distribution. The kernel density of this distribution is C [p1^(r1-1)]*[p2^(r2-1)]*...*[pn^(rn-1)]. C Inputs: C N Order of distributions C R Entry i contains ri. C Outputs: C P Entries (1,...,n) contain a realization of (p1,...,pn) from the C distribution. C Reference: C Johnson and Kotz (1972), Section 40.5 IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(N),P(N) AA=0.0D0 DO 10 I=1,N CALL DRNCHI(1,2.0D0*R(I),P(I)) 10 AA=AA+P(I) AA=1.0D0/AA DO 20 I=1,N 20 P(I)=AA*P(I) RETURN END REAL*8 FUNCTION BETK(N,R,P) C This routine evaluates the log density kernel of a vector p with a C multivariate Beta(n,r) distribution. The density kernel is C [p1^(r1-1)]*[p2^(r2-1)]*...*[pn^(rn-1)]. C Inputs: C N Order of distributions C R Entry i contains ri. C P Entry i contains pi. C Outputs: C BETK Value of log density kernel IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(N),P(N) AA=0.0D0 DO 10 I=1,N 10 AA=AA+(R(I)-1.0D0)*DLOG(P(I)) BETK=AA RETURN END REAL*8 FUNCTION BETN(N,R) C This routine evaluates the log normalizing constant for the log C density kernel of a random variable computed in BETN. This C normalizing constant, when multiplied by the kernel in BETK, C yields the p.d.f. of a random vector p ~ Beta(n,r). C Inputs: C N Order of distributions C R Entry i contains ri. C Outputs: C BETK Value of log density kernel IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(N) AA=0.0D0 BB=0.0D0 DO 10 I=1,N AA=AA+R(I) 10 BB=BB+DLNGAM(R(I)) BETN=DLNGAM(AA)-BB RETURN END SUBROUTINE BET0(NFILE,N,TITLE,M,R,LDR,ANORM,P,LDP) C This routine reads in the parameters for multiple multivariate beta C distributions. It prints out informataion about the distributions, C computes and returns normalizing constants for the densities, C and one draw from each distribution. C Inputs: C NFILE File number for reading input information C N Number of beta distributions C TITLE Information printed at start C M Number of categories for each multivariate beta distribution C LDR Leading dimension of R matrix C LDP Leading dimension of P matrix C Outputs: C R Row i contains the M parameters of distribution i. C ANORM Log normalizing constant for joint multivariate beta distribution C P Row i has one drawing from distribution i. IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*)TITLE PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION R(LDR,M),P(LDR,M) WRITE(6,10)TITLE 10 FORMAT(/,' Multivariate beta prior distribution for ',A) AA=0.0D0 DO 50 I=1,N READ(NFILE,*,END=110,ERR=120)(V1(J),J=1,M) DO 20 J=1,M 20 IF(V1(J).LE.0.0D0)GO TO 130 AA=AA+BETN(M,V1) CALL BET(M,V1,V2) WRITE(6,30)(V1(J),J=1,M) 30 FORMAT(' Parameters:',4(1PE16.7),/,(12X,4(1PE16.7))) WRITE(6,40)(V2(J),J=1,M) 40 FORMAT(' Drawing:',4(1PE16.7),/,(12X,4(1PE16.7))) DO 50 J=1,M R(I,J)=V1(J) 50 P(I,J)=V2(J) ANORM=AA RETURN 110 CALL PEND('UNANTICIPATED END-OF-FILE READING BETA PRIOR') 120 CALL PEND('ERROR READING BETA PRIOR') 130 CALL PINTR(' BETA PRIOR PARAMETER',V1(J),0.0D0,1.0D1000) END SUBROUTINE CENMOM(RAW,CEN,SKEW,EKURT) C This routine computes the first four central moments, skewness and C excess kurtosis of a distribution, given the first four raw moments. C Inputs: C RAW(4) First four raw moments C Outputs: C CEN(4) First four central moments C SKEW Skewness coefficient C EKURT Excess kurtosis IMPLICIT REAL*8 (A-H,O-Z) DIMENSION RAW(4),CEN(4) AMU=RAW(1) CEN(1)=0.0D0 CEN(2)=RAW(2)-RAW(1)**2 CEN(3)=RAW(3)-3.0D0*AMU*RAW(2)+2.0D0*AMU**3 CEN(4)=RAW(4)-4.0D0*AMU*RAW(3)+6.0D0*(AMU**2)*RAW(2)-3.0D0*AMU**4 SKEW=CEN(3)/CEN(2)**1.5D0 EKURT=CEN(4)/CEN(2)**2-3.0D0 RETURN END REAL*8 FUNCTION CHI(S2,ANU) C This routine generates a random variable z, where C S2*z~chi-square(nu). C Inputs: C S2 Numerator parameter s2 C ANU Degrees of freedom parameter nu C Output: C CHI Random variable z IMPLICIT REAL*8 (A-H,O-Z) CALL DRNCHI(1,ANU,CHISQ) CHI=CHISQ/S2 RETURN END REAL*8 FUNCTION CHIK(S2,ANU,X) C This routine evaluates the log density kernel of a random variable C z for which S2*z~chi-square(ANU), at the point x. C The log density kernel is C ((nu-2)/2)*log(x) - s2*x/2 C Inputs: C S2 Numerator parameter s2 C ANU Degrees of freedom parameter nu C X Point of evaluation x C Output C CHIK Value of log density kernel IMPLICIT REAL*8 (A-H,O-Z) CHIK=0.5D0*((ANU-2.0D0)*DLOG(X)-S2*X) RETURN END REAL*8 FUNCTION CHIN(S2,ANU) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This routine evaluates the log normalizing constant for the log C density kernel of a random variable computed in CHIK. This C normalizing constant, when multiplied by the kernel in CHIK, C yields the p.d.f. of a random variable z for which S2*z~chi- C square(ANU) C Inputs: C S2 Numerator parameter s2 C ANU Degrees of freedom parameter nu C Output C CHIN Value of log normalizing constant IMPLICIT REAL*8 (A-H,O-Z) COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII AA=0.5D0*ANU BB=DLNGAM(AA) CHIN=-AA*DLOG2-DLNGAM(AA)+AA*DLOG(S2) RETURN END SUBROUTINE CHI0(NFILE,N,TITLE,LORD,ANU0,S20,ANORM,ACHI) C This routine reads in the constants and degrees-of-freedom parameters C for multiple chi-square distributions. It prints out informataion C about the distributions, computes and returns normalizing constants C for the densities, and one draw from each distribution. C Inputs: C NFILE File number for reading input information C N Number of chi-square distributions C TITLE Information printed at start C LORD If .TRUE., drawings must be in ascending order C Outputs: C ANU0 Degrees of freedom parameters C S20 Constant parametrs C ANORM Log normalizing constant for joint chi-square distribution C ACHI Drawings from distribution IMPLICIT REAL*8 (A-H,O-Z) REAL*4 CPSEC DATA TIME,SD,JUMP/60.0D0,1.0D-3,2500/ CHARACTER*(*)TITLE PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION ANU0(N),S20(N),ACHI(N) IF(LORD)GO TO 100 WRITE(6,10)TITLE 10 FORMAT(/,' Chi-square prior for ',A,/, 1 ' Degrees of freedom Scale parameter Initial draw') AA=0.0D0 DO 20 I=1,N READ(NFILE,*,END=40,ERR=50)ANU0(I),S20(I) CALL PINTR('Degrees of freedom',ANU0(I),0.0D0,1.D1000) CALL PINTR('Constant parameter',S20(I),0.0D0,1.D1000) AA=AA+CHIN(S20(I),ANU0(I)) ACHI(I)=CHI(S20(I),ANU0(I)) 20 WRITE(6,30)ANU0(I),S20(I),ACHI(I) 30 FORMAT(1PE19.7,1PE20.7,1PE17.7) ANORM=AA RETURN 40 CALL PEND('UNANTICIPATED END-OF-FILE READING CHISQUARE PRIOR') 50 CALL PEND('ERROR READING CHISQUARE PRIOR') C Ordered distribution 100 WRITE(6,110)TITLE 110 FORMAT(/,' Ordered chi-square prior for ',A,/, 1 ' Degrees of freedom Scale parameter Initial draw') AA=0.0D0 DO 120 I=1,N READ(NFILE,*,END=40,ERR=50)ANU0(I),S20(I) CALL PINTR('Degrees of freedom',ANU0(I),0.0D0,1.D1000) CALL PINTR('Constant parameter',S20(I),0.0D0,1.D1000) ACHI(I)=CHI(S20(I),ANU0(I)) 120 AA=AA+CHIN(S20(I),ANU0(I)) CALL CHIO(N,S20,ANU0,20,ACHI) DO 130 I=1,N 130 WRITE(6,30)ANU0(I),S20(I),ACHI(I) C Adjust normalization constant for ordering NITER=0 NTOT=0 CLOCK=CPSEC() 140 DO 160 IREP=1,JUMP DO 150 I=1,N 150 V1(I)=CHI(S20(I),ANU0(I)) 160 IF(LORDER(N,V1))NTOT=NTOT+1 NITER=NITER+JUMP P=DREAL(NTOT)/NITER SE=DSQRT(P*(1.0D0-P)/NITER) CLOCK1=CPSEC()-CLOCK IF((SE.GT.SD).AND.(CLOCK1.LT.TIME))GO TO 140 IF(NTOT.EQ.0)CALL PEND 1 ('Ordering of precisions inconsistent with prior parameters') ANORM=AA-DLOG(P) WRITE(6,170)SE,CLOCK1 170 FORMAT( 1 ' ...Standard deviation of approximation of log prior kernel:' 2 ,1PE15.5,/, 3 ' Computation time:' 4 ,0PF8.2) RETURN END SUBROUTINE CHI0P(NFILE,N,TITLE,LORD,JFIX,HFIX,ANU0,S20,ANORM,ACHI) C This routine reads in the constants and degrees-of-freedom parameters C for multiple chi-square distributions. It prints out informataion C about the distributions, computes and returns normalizing constants C for the densities, and one draw from each distribution. One of the C "chi-squares" is fixed at a pre-set value. C Inputs: C NFILE File number for reading input information C N Number of chi-square distributions C TITLE Information printed at start C LORD If .TRUE., drawings must be in ascending order C JFIX Precision to be fixed C HFIX Value of fixed precision C Outputs: C ANU0 Degrees of freedom parameters C S20 Constant parametrs C ANORM Log normalizing constant for joint chi-square distribution C ACHI Drawings from distribution IMPLICIT REAL*8 (A-H,O-Z) REAL*4 CPSEC DATA TIME,SD,JUMP/60.0D0,1.0D-3,2500/ CHARACTER*(*)TITLE PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION ANU0(N),S20(N),ACHI(N) IF(LORD)GO TO 100 WRITE(6,10)TITLE 10 FORMAT(/,' Chi-square prior for ',A,/, 1 ' Degrees of freedom Scale parameter Initial draw') AA=0.0D0 DO 30 I=1,N IF(I.NE.JFIX)GO TO 20 ACHI(I)=HFIX WRITE(6,15)ACHI(I) 15 FORMAT(15X,'----',16X,'----',1PE17.7) GO TO 30 20 READ(NFILE,*,END=40,ERR=50)ANU0(I),S20(I) CALL PINTR('Degrees of freedom',ANU0(I),0.0D0,1.D1000) CALL PINTR('Constant parameter',S20(I),0.0D0,1.D1000) AA=AA+CHIN(S20(I),ANU0(I)) ACHI(I)=CHI(S20(I),ANU0(I)) WRITE(6,25)ANU0(I),S20(I),ACHI(I) 25 FORMAT(1PE19.7,1PE20.7,1PE17.7) 30 CONTINUE ANORM=AA RETURN 40 CALL PEND('UNANTICIPATED END-OF-FILE READING CHISQUARE PRIOR') 50 CALL PEND('ERROR READING CHISQUARE PRIOR') C Ordered distribution 100 WRITE(6,110)TITLE 110 FORMAT(/,' Ordered chi-square prior for ',A,/, 1 ' Degrees of freedom Scale parameter Initial draw') AA=0.0D0 DO 120 I=1,N ACHI(I)=HFIX IF(I.EQ.JFIX)GO TO 120 READ(NFILE,*,END=40,ERR=50)ANU0(I),S20(I) CALL PINTR('Degrees of freedom',ANU0(I),0.0D0,1.D1000) CALL PINTR('Constant parameter',S20(I),0.0D0,1.D1000) ACHI(I)=CHI(S20(I),ANU0(I)) AA=AA+CHIN(S20(I),ANU0(I)) 120 CONTINUE CALL CHIOP(N,S20,ANU0,20,JFIX,ACHI) DO 130 I=1,N IF(I.EQ.JFIX)WRITE(6,15)ACHI(I) 130 IF(I.NE.JFIX)WRITE(6,25)ANU0(I),S20(I),ACHI(I) C Adjust normalization constant for ordering NITER=0 NTOT=0 CLOCK=CPSEC() V1(JFIX)=HFIX 140 DO 160 IREP=1,JUMP DO 150 I=1,N 150 IF(I.NE.JFIX)V1(I)=CHI(S20(I),ANU0(I)) 160 IF(LORDER(N,V1))NTOT=NTOT+1 NITER=NITER+JUMP P=DREAL(NTOT)/NITER SE=DSQRT(P*(1.0D0-P)/NITER) CLOCK1=CPSEC()-CLOCK IF((SE.GT.SD).AND.(CLOCK1.LT.TIME))GO TO 140 IF(NTOT.EQ.0)CALL PEND 1 ('Ordering of precisions inconsistent with prior parameters') ANORM=AA-DLOG(P) WRITE(6,170)SE,CLOCK1 170 FORMAT( 1 ' ...Standard deviation of approximation of log prior kernel:' 2 ,1PE15.5,/, 3 ' Computation time:' 4 ,0PF8.2) RETURN END SUBROUTINE CHIO(N,S2,ANU,NITER,H) C This routine draws from N independent chi square distributions, subject C to the restriction that the realizations are monotonically nondescending. C Inputs: C N Number of distributions C S2(Nx1) Constant parameters for distributions C ANU(Nx1) Degrees of freedom parameters for distributions C NITER Number of Gibbs iterations to perform C Input/outputs: C H Random vector IMPLICIT REAL*8 (A-H,O-Z) DIMENSION S2(N),ANU(N),H(N) DIMENSION V1(10) V1(1)=CHI(S2(1),ANU(1)) DO 10 I=2,N V1(I)=CHI(S2(I),ANU(I)) 10 IF(V1(I).LE.V1(I-1))RETURN CALL DCOPY(N,V1,1,H,1) RETURN END SUBROUTINE CHIO1(N,S2,ANU,NITER,CHI) C This routine draws from N independent chi square distributions, subject C to the restriction that the realizations are monotonically nondescending. C Inputs: C N Number of distributions C S2(Nx1) Constant parameters for distributions C ANU(Nx1) Degrees of freedom parameters for distributions C NITER Number of Gibbs iterations to perform C Input/outputs: C CHI Random vector IMPLICIT REAL*8 (A-H,O-Z) DIMENSION S2(N),ANU(N),CHI(N) N1=N-1 DO 20 ITER=1,NITER CHI(1)=XCHIT(S2(1),ANU(1),0.0D0,CHI(2),.FALSE.) IF(N.EQ.2)GO TO 20 DO 10 I=2,N1 10 CHI(I)=XCHIT(S2(I),ANU(I),CHI(I-1),CHI(I+1),.FALSE.) 20 CHI(N)=XCHIT(S2(N),ANU(N),CHI(N1),-1.0D0,.TRUE.) RETURN END SUBROUTINE CHIOP(N,S2,ANU,NITER,JFIX,CHI) C This routine draws from N independent chi square distributions, subject C to the restriction that the realizations are monotonically nondescending. C One of the "chi-squares" is fixed at a preset value. C Inputs: C N Number of distributions C S2(Nx1) Constant parameters for distributions C ANU(Nx1) Degrees of freedom parameters for distributions C NITER Number of Gibbs iterations to perform C JFIX Value of (pre-)fixed precision C Input/outputs: C CHI Random vector IMPLICIT REAL*8 (A-H,O-Z) DIMENSION S2(N),ANU(N),CHI(N) N1=N-1 DO 20 ITER=1,NITER IF(JFIX.NE.1)CHI(1)=XCHIT(S2(1),ANU(1),0.0D0,CHI(2),.FALSE.) IF(N.EQ.2)GO TO 20 DO 10 I=2,N1 10 IF(I.NE.JFIX)CHI(I)=XCHIT(S2(I),ANU(I),CHI(I-1),CHI(I+1),.FALSE.) 20 IF(N.NE.JFIX)CHI(N)=XCHIT(S2(N),ANU(N),CHI(N1),-1.0D0,.TRUE.) RETURN END SUBROUTINE CPROD(N1,N2,LIST1,LIST2,ZZ,LDZ,XY,LDXY) C This routine extracts the sums of cross products of a subset of C the variables in a sums-of-squares-and-cross-products matrix. C It is useful for reordering variables. C Inuts: C N1 Number of rows in output sums-of-cross-products matrix C N2 Number of columns in output sums-of-cross-products matrix C LIST1 N1x1 vector of ZZ entries corresponding to rows of XY C LIST1 N2x1 vector of ZZ entries corresponding to columns of XY C ZZ Sums-of-squares-and-cross-products matrix C LDZ Leading dimension of Z C LDXY Leading dimension of XY C Output: C XY N1xN2 sums-of-cross-products matrix IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION LIST1(N1),LIST2(N2),ZZ(LDZ,LDZ),XY(LDXY,LDXY) DO 10 I1=1,N1 J1=LIST1(I1) DO 10 I2=1,N2 J2=LIST2(I2) 10 A1(I1,I2)=ZZ(J1,J2) DO 20 I1=1,N1 DO 20 I2=1,N2 20 XY(I1,I2)=A1(I1,I2) RETURN END SUBROUTINE DAT0(NFILE,NT1,NT2,NVARS,ZZ,LDZ) C This routine reads a data file, writes some summary statistics to C output, and returns the sums-of-squares-and-cross-products matrix. C Inputs: C NFILE File number for read C NT1 First record (observation) to include C NT1 Last record (observation) to include C LDZ Leading dimension of ZZ C Output: C NVARS Number of variables available on data file C ZZ Sums-of-squares-and-cross-products matrix IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION ZZ(LDZ,LDZ) CALL FOPEN('Data',NFILE,.FALSE.) IT=0 READ(NFILE,*,END=110,ERR=120)NOBS,NVARS NK=NVARS WRITE(6,10)NK,NOBS,NT1,NT2 10 FORMAT(/,' Number of variables read from file:',I6,/, 2 ' Number of observations in file:',I6,/, 3 ' Observations used for sample:',I6,' -',I6) CALL PINTI('Columns to be read',NK,1,MIN0(LD,LDZ)) CALL PINTI('t1',NT1,1,NT2) CALL PINTI('t2',NT2,NT1,NOBS) NT=NT2-NT1+1 CALL UM0SET(NK,NK,ZZ,LDZ) DO 50 IT=1,NT2 READ(NFILE,*,END=110,ERR=120)(V1(IVARS),IVARS=1,NVARS) IF(IT.LT.NT1)GO TO 50 IF(IT.GT.NT1)GO TO 30 DO 20 IK=1,NK A1(1,IK)=V1(IK) A1(2,IK)=V1(IK) A1(3,IK)=0.0D0 20 A1(4,IK)=V1(IK) 30 DO 40 IK=1,NK IF(V1(IK).LT.A1(1,IK))A1(1,IK)=V1(IK) IF(V1(IK).GT.A1(2,IK))A1(2,IK)=V1(IK) A1(3,IK)=A1(3,IK)+V1(IK) DO 40 JK=IK,NK 40 ZZ(IK,JK)=ZZ(IK,JK)+V1(IK)*V1(JK) 50 CONTINUE CALL DFULL(NK,ZZ,LDZ) WRITE(6,60) 60 FORMAT(/,' Col',5X,'Minimum',5X,'Maximum',8X,'Mean',2X, 1 'Stan. dev.',7X,'First',8X,'Last') DO 70 IK=1,NK AMEAN=A1(3,IK)/NT AVAR=ZZ(IK,IK)/NT-AMEAN**2 SDEV=0.0D0 IF(AVAR.GT.0.0D0)SDEV=DSQRT(AVAR) 70 WRITE(6,80)IK,A1(1,IK),A1(2,IK),AMEAN,SDEV,A1(4,IK),V1(IK) 80 FORMAT(I4,6(1PE12.4)) RETURN 110 CALL PEND('UNANTICIPATED END-OF-FILE READING DATA') 120 CALL PEND('ERROR READING DATA') END SUBROUTINE DAT1(NFILE,NT1,NT2,NVARS,D,LDD,NCD) C This routine reads a data file, writes some summary statistics to C output, and returns records t1 through t2 in the corresponding rows of D. C Inputs: C NFILE File number for read C NT1 First record (observation) to include, t1 C NT2 Last record (observation) to include, t2 C LDD Leading dimension of D C NCD Number of columns allocated to D in calling routine C Output: C NVARS Number of variables available on data file C D Records t1 through t2 of data file IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION D(LDD,NCD) CALL FOPEN('Data',NFILE,.FALSE.) IT=0 READ(NFILE,*,END=130,ERR=140)NOBS,NVARS NK=NVARS WRITE(6,10)NK,NOBS,NT1,NT2 10 FORMAT(' Number of variables read from file:',I6,/, 2 ' Number of observations in file:',I6,/, 3 ' Observations used for sample:',I6,' -',I6) NKMAX=MIN0(LD,NCD) IF(NK.GT.NKMAX)GO TO 110 IF((NT1.GT.NT2).OR.(NT1.LE.0).OR.(NT2.GT.NOBS))GO TO 120 IF(NT2.GT.LDD)GO TO 150 NT=NT2-NT1+1 DO 50 IT=1,NT2 READ(NFILE,*,END=130,ERR=140)(V1(IVARS),IVARS=1,NVARS) IF(IT.LT.NT1)GO TO 50 IF(IT.GT.NT1)GO TO 30 DO 20 IK=1,NK A1(1,IK)=V1(IK) A1(2,IK)=V1(IK) A1(3,IK)=0.0D0 A1(4,IK)=V1(IK) 20 A1(5,IK)=0.0D0 30 DO 40 IK=1,NK IF(V1(IK).LT.A1(1,IK))A1(1,IK)=V1(IK) IF(V1(IK).GT.A1(2,IK))A1(2,IK)=V1(IK) A1(3,IK)=A1(3,IK)+V1(IK) A1(5,IK)=A1(5,IK)+V1(IK)**2 40 D(IT,IK)=V1(IK) 50 CONTINUE WRITE(6,60) 60 FORMAT(/,' Col',5X,'Minimum',5X,'Maximum',8X,'Mean',2X, 1 'Stan. dev.',7X,'First',8X,'Last') DO 70 IK=1,NK AMEAN=A1(3,IK)/NT AVAR=A1(5,IK)/NT-AMEAN**2 SDEV=0.0D0 IF(AVAR.GT.0.0D0)SDEV=DSQRT(AVAR) 70 WRITE(6,80)IK,A1(1,IK),A1(2,IK),AMEAN,SDEV,A1(4,IK),V1(IK) 80 FORMAT(I4,6(1PE12.4)) RETURN 110 WRITE(6,115)NVARS,NKMAX 115 FORMAT(/,' *** FILE CONTAINS',I4,' VARIABLES, EXCEEDS CAPACITY OF' 1 ,I4) CALL TERM 120 WRITE(6,125)NT1,NT2,NOBS 125 FORMAT(/,' *** RECORDS',I6,' -',I6' REQUESTED, ONLY 1-',I6, 1 ' AVAILABLE') CALL TERM 130 WRITE(6,135)IT 135 FORMAT(/,' *** UNANTICIPATED END-OF-FILE, OBSERVATION',I8) CALL TERM 140 WRITE(6,125)IT 145 FORMAT(/,' **** READ ERROR, OBSERVATION',I8) CALL TERM 150 WRITE(6,155)NT2,LDD 155 FORMAT(/,' *** FILE CONTAINS',I6,' OBSERVATIONS, EXCEEDS', 1 ' CAPACITY OF',I6) CALL TERM END SUBROUTINE DAT2(NFILE,NT1,NT2,NLAG,NVARS,ZZ,LDZ) C This routine reads a data file, writes some summary statistics to C output, and returns the sums-of-squares-and-cross-products matrix C for variables and their lags. C Inputs: C NFILE File number for read C NT1 First record (observation) to include C NT2 Last record (observation) to include C NLAG Number of lags of each variable C LDZ Leading dimension of ZZ C Output: C NVARS Number of variables available on data file C ZZ Sums-of-squares-and-cross-products matrix IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION ZZ(LDZ,LDZ) CALL FOPEN('Data',NFILE,.FALSE.) IT=0 READ(NFILE,*,END=110,ERR=120)NOBS,NVARS NK=NVARS NZ=(NLAG+1)*NK WRITE(6,10)NK,NOBS,NT1,NT2 10 FORMAT(/,' Number of variables read from file:',I6,/, 2 ' Number of observations in file:',I6,/, 3 ' Observations used for sample:',I6,' -',I6) CALL PINTI('Columns to be read',NK,1,MIN0(LD,LDZ)) CALL PINTI('t1',NT1,1,NT2) CALL PINTI('t2',NT2,NT1,NOBS) CALL PINTI('Variables including lags',NZ,1,LD) NT=NT2-NT1+1 CALL UM0SET(NK,NK,ZZ,LDZ) DO 50 IT=1,NT2 DO 15 L=NLAG,1,-1 15 CALL DCOPY(NVARS,A2(L,1),LD,A2(L+1,1),LD) READ(NFILE,*,END=110,ERR=120)(V1(IVARS),IVARS=1,NVARS) CALL DCOPY(NVARS,V1,1,A2,LD) IF(IT.LT.NT1)GO TO 50 IF(IT.GT.NT1)GO TO 30 DO 20 IK=1,NK A1(1,IK)=V1(IK) A1(2,IK)=V1(IK) A1(3,IK)=0.0D0 20 A1(4,IK)=V1(IK) 30 J=0 DO 40 IK=1,NK IF(V1(IK).LT.A1(1,IK))A1(1,IK)=V1(IK) IF(V1(IK).GT.A1(2,IK))A1(2,IK)=V1(IK) A1(3,IK)=A1(3,IK)+V1(IK) DO 40 L=0,NLAG J=J+1 40 V2(J)=A2(L+1,IK) DO 45 IZ=1,NZ DO 45 JZ=IZ,NZ 45 ZZ(IZ,JZ)=ZZ(IZ,JZ)+V2(IZ)*V2(JZ) 50 CONTINUE CALL DFULL(NZ,ZZ,LDZ) WRITE(6,60) 60 FORMAT(/,' Col',5X,'Minimum',5X,'Maximum',8X,'Mean',2X, 1 'Stan. dev.',7X,'First',8X,'Last') DO 70 IK=1,NK AMEAN=A1(3,IK)/NT AVAR=ZZ(IK,IK)/NT-AMEAN**2 SDEV=0.0D0 IF(AVAR.GT.0.0D0)SDEV=DSQRT(AVAR) 70 WRITE(6,80)IK,A1(1,IK),A1(2,IK),AMEAN,SDEV,A1(4,IK),V1(IK) 80 FORMAT(I4,6(1PE12.4)) RETURN 110 CALL PEND('UNANTICIPATED END-OF-FILE READING DATA') 120 CALL PEND('ERROR READING DATA') END SUBROUTINE DATFMC(NFILE,NOPT,NM,AN,LDM) C This routine reads a first-order Markov chain model data file, C and returns the data matrix. C Inputs: C NFILE File number to read C NOPT If >0 print data matrix, otherwise not C NM Number of categories C LDM Leading dimension of A C Output: C AN(NMxNM) Data matrix PARAMETER(LD=100) IMPLICIT REAL*8 (A-H,O-Z) LOGICAL LSTOP COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION AN(LDM,NM) CALL FOPEN('Data',NFILE,.FALSE.) READ(NFILE,*,END=210,ERR=220)NM1,NM2 IF((NM1.NE.NM).OR.(NM2.NE.NM))GO TO 230 READ(NFILE,*,END=210,ERR=220)((AN(I,J),J=1,NM),I=1,NM) IF(NOPT.LE.0)GO TO 60 DO 50 JBASE=1,NM,8 JTOP=JBASE+7 IF(JTOP.GT.NM)JTOP=NM WRITE(6,20)(J,J=JBASE,JTOP) 20 FORMAT(8I9) DO 30 I=1,NM 30 WRITE(6,40)I,(AN(I,J),J=JBASE,JTOP) 40 FORMAT(I3,8F9.1) 50 CONTINUE 60 CALL DSET(NM,0.0D0,V1,1) CALL DSET(NM,0.0D0,V2,1) LSTOP=.FALSE. DO 70 I=1,NM DO 70 J=1,NM V1(I)=V1(I)+AN(I,J) V2(J)=V2(J)+AN(I,J) IF(AN(I,J).GE.0.0D0)GO TO 70 WRITE(6,65)I,J 65 FORMAT(/,' NEGATIVE DATA ENTRY ROW',I3,' COLUMN',I3) LSTOP=.TRUE. 70 CONTINUE IF(LSTOP)CALL TERM WRITE(6,75)DSUM(NM,V1,1) 75 FORMAT(F10.1,' observations in sample') DO 90 I=1,NM IF(V1(I).EQ.0.0D0)WRITE(6,80)I 80 FORMAT(' +++ WARNING: ROW SUM',I3,' IS 0.0') 90 IF(V2(I).EQ.0.0D0)WRITE(6,100)I 100 FORMAT(' +++ WARNING: COLUMN SUM',I3,' IS 0.0') IF(NOPT.LT.0)RETURN WRITE(6,110)(I,V1(I),V2(I),I=1,NM) 110 FORMAT(/,' Row sums Col sums',/,(I3,2F9.1)) RETURN 210 CALL PEND('UNEXPECTED END OF FILE READING MARKOV CHAIN DATA') 220 CALL PEND('ERROR READING MARKOV CHAIN DATA') 230 CALL PEND('NUMBER CATEGORIES IN CONTROL AND DATA FILES CONFLICT') END SUBROUTINE DFULL(N,A,LDA) C This routine converts a symmetric matrix in IMSL symmetric storeage mode C (upper triangular) to one in full storeage mode. C Inputs: C N Order of symmetric matrix C A Matrix in symmetric storeage mode C LDA Leading dimension of A C Output: C A Matrix in full storeage mode IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N) N1=N-1 DO 10 I=1,N1 DO 10 J=I+1,N 10 A(J,I)=A(I,J) RETURN END SUBROUTINE DMXXTU(N,A,LDA,B,LDB) C This routine computes B=AA', where A is upper triangular. C Only the upper triangle of B is returned (IMSL symmetric storeage mode) C Inputs: C N Order of upper triangular matrix C A Upper triangular matrix C LDA Leading dimension of A C LDB Leading dimension of B C Output C B B=AA' (upper triangle only) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),B(LDB,N) DO 20 I=1,N DO 20 J=I,N AA=0.0D0 DO 10 L=J,N 10 AA=AA+A(I,L)*A(J,L) 20 B(I,J)=AA RETURN END SUBROUTINE DPDCK(A,LDA,N,INFO) C This routine uses IMSL to check whether the matrix supplied is positive C definite. Returns INFO=0 if not INFO=1 if it is pd IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION A(LDA,*),EVAL(1000) CALL DEVLSF(N,A,LDA,EVAL) DO 10 I=1,N IF(EVAL(I).LE.0) GO TO 100 10 CONTINUE INFO=0 RETURN 100 INFO=1 RETURN END REAL*8 FUNCTION DTLGPD(N,A,LDA) C This function returns the log determinant of a positive definite C matrix A. C Inputs: C N Order of matrix C A Matrix C LDA Leading dimension of matrix C Output: C DTLGPD Log determinant of matrix IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION A(LDA,N) CALL DPDCK(A,LDA,N,INFO) IF(INFO.NE.0)GO TO 110 call dlftds(n,a,lda,a1,ld) AA=0.0D0 DO 10 I=1,N 10 AA=AA+DLOG(A1(I,I)) DTLGPD=2.0D0*AA RETURN 110 WRITE(6,115) 115 FORMAT(/,' *** NON-P.D. MATRIX IN DTLGPD',/,' PROGRAM STOPS.') STOP END REAL*8 FUNCTION FARSD1(ALAM) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This function evaluates the product of the spectral density function C of a univariate autogregressive process, and the function C cos(lambda*tau), at the frequencies lambdaC Inputs: C ALAM Lambda C Inputs (/AUAR/): C B Coefficients of autoregression C H Precision of shock to autoregression C TAU tau C NP Number of coefficients in autoregression C Output: C FARSD1 Value of product described above IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LP=100) COMMON/AUAR/B(0:LP),H,TAU,NP COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII AL=ALAM AC=B(0) AS=0.0D0 DO 10 IP=1,NP AC=AC+B(IP)*DCOS(AL*IP) 10 AS=AS+B(IP)*DSIN(AL*IP) FARSD1=DCOS(AL*TAU)/(H*(AS**2+AC**2)) RETURN END SUBROUTINE FOPEN(HTITLE,N,LNEW) C Read a file name from input file 5 and open the file. C Inputs: C HTITLE Information to print about opened file C N File number to open C LNEW .TRUE.: File is new; .FALSE.: File is old CHARACTER*(*)HTITLE COMMON/FILES/HFILE(4) CHARACTER*40 HFILE LOGICAL LNEW READ(5,10,ERR=110,END=110)HFILE(N) 10 FORMAT(A40) IF(.NOT.LNEW)GO TO 30 OPEN(UNIT=N,STATUS='NEW',FILE=HFILE(N),ERR=120) WRITE(6,20)HTITLE,HFILE(N) 20 FORMAT(/,' >>> ',A,' written to file ',A40) RETURN 30 OPEN(UNIT=N,STATUS='OLD',FILE=HFILE(N),ERR=130) WRITE(6,40)HTITLE,HFILE(N) 40 FORMAT(/,' <<< ',A,' read from file ',A40) RETURN 110 WRITE(6,115) 115 FORMAT(/,' *** ERROR READING FILE NAME') CALL TERM 120 WRITE(6,125)HFILE(N) 125 FORMAT(/,' *** FAILURE TO OPEN NEW FILE ',A40) CALL TERM 130 WRITE(6,135)HFILE(N) 135 FORMAT(/,' *** FAILURE TO OPEN OLD FILE ',A40) CALL TERM END SUBROUTINE GAU(K,AMU,H,LDH,X1,X2) C This routine generates two antithetic pseudo-random vectors x1 and x2 C from a Gaussian distribution with mean vector mu and precision C matrix H. C Inputs: C K Dimension of x C AMU Mean vector mu (kx1) C H Precision matrix H (kxk) C LDH Leading dimension of H C Outputs: C X1,X2 Random draws IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION AMU(K),H(LDH,K),X1(K),X2(K) C Entry GAU: Factor precision CALL DLFTDS(K,H,LDH,A1,LD) CALL DLINRT(K,A1,LD,2,A1,LD) C Entry GAUA: Precision already factored; C *** Can be no intervening calls that change A1 !! *** ENTRY GAUA(K,AMU,H,LDH,X1,X2) 10 CALL URNMVN(K,A1,LD,V1) CALL UVADD(K,AMU,V1,X1) CALL UVSUB(K,AMU,V1,X2) RETURN C Entry GAUB: Coming in with R: RR'=Variance, R upper triangular ENTRY GAUB(K,AMU,R,LDR,X1,X2) CALL UCRGRG(K,K,R,LDR,A1,LD) GO TO 10 END REAL*8 FUNCTION GAUK(K,AMU,H,LDH,X) C This function evaluates the log density kernel of a multivariate C normal random vector x with mean mu and precision matrix H. The C kernel density is exp(-.5(x-mu)'H(x-mu)) C Inputs: C K Dimension of x C AMU Mean vector mu (kx1) C H Precision matrix h (kxk) C LDH Leading dimension of H C X Point of evaluation, x (kx1) C Output: C GAUK Evaluation of log kernel density IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION H(LDH,K),X(K) CALL UVSUB(K,X,AMU,V1) GAUK=-0.5D0*UQUAF(K,V1,H,LDH) RETURN END REAL*8 FUNCTION GAUN(K,H,LDH) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This routine evaluates the log normalizing constant for the log C density kernel of a random vector computed in GAUK. This C normalizing constant, when multiplied by the kernel in GAUK, C yields the p.d.f. of the random vector. C Inputs: C K Dimension of random vector C H Precision matrix H (kxk) C LDH Leading dimension of H C Output: C GAUN Log normalizing factor IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION H(LDH,K) CALL DLFTDS(K,H,LDH,A1,LD) AA=0.0D0 DO 10 I=1,K 10 AA=AA+DLOG(A1(I,I)) GAUN=-K*HLG2PI+AA RETURN END SUBROUTINE GAU0(NFILE,NK,TITLE,H0,LDH,B0,HB0,ANORM,THETA) C This routine reads in the mean vector mu and precision matrix H of a C kx1 parameter vector theta. It computes and returns the product of C H and mu, and prints information about the prior. C Inputs: C NFILE File number for reading prior C NK Dimension of parameter vector = k C TITLE Identifier for distribution (printed) C LDH Leading dimension of H0 C Ouptuts: C H0 Prior precision matrix H C B0 Prior mean vector mu C HB0 Product of H and mu C ANORM Log normalizing constant for prior density C THETA Draw from prior IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) TITLE PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION H0(LDH,NK),HB0(NK),B0(NK),THETA(NK) READ(NFILE,*,END=110,ERR=120)NCODE1,NCODE2,NDIM C NCODE1=1 Diagonal, sqrts only C NCODE1=2 Sqrts diagonal followed by normalized matrix C NCODE1=3 Full matrix C NCODE2=1 Variance matrix C NCODE2=2 Precision matrix IF(NDIM.NE.NK)CALL PEND('ORDER OF GAUSSIAN PRIOR INCORRECT') CALL PINTI('Gaussian prior code 1',NCODE1,1,3) CALL PINTI('Gaussian prior code 2',NCODE2,1,2) READ(NFILE,*,END=110,ERR=120)(B0(IK),IK=1,NK) CALL DSET(NK,1.0D0,V1,1) IF(NCODE1.NE.3)CALL PLISTR(NFILE,NK,V1) CALL UMISET(NK,H0,LDH) IF(NCODE1.NE.1)CALL PDSYMR(NFILE,NK,H0,LDH) DO 10 IK=1,NK DO 10 JK=1,IK H0(IK,JK)=H0(IK,JK)*V1(IK)*V1(JK) 10 H0(JK,IK)=H0(IK,JK) IF(NCODE2.EQ.1)CALL DLINDS(NK,H0,LDH,H0,LDH) C Form Hx(mu). CALL DMURRV(NK,NK,H0,LDH,NK,B0,1,NK,HB0) C Normalizing constant and initial draw ANORM=GAUN(NK,H0,LDH) CALL GAU(NK,B0,H0,LDH,THETA,V2) C Variances as convenience to user CALL DLINDS(NK,H0,LDH,A1,LD) WRITE(6,30)TITLE,(JK,B0(JK),DSQRT(A1(JK,JK)),H0(JK,JK), 1 THETA(JK),JK=1,NK) 30 FORMAT(/,' Gaussian prior for ',A,/, 1 ' Parameter',6X,'Prior mean',6X,'Prior s.d.', 2 5X,'Prior prec.',3X,'Initial value',/, 3 (5X,I5,4(1PE16.6))) CALL PRCOR(' Prior correlation matrix',NK,A1,LD) CALL PRCOR(' Normalized prior precision matrix',NK,H0,LDH) RETURN 110 CALL PEND('UNANTICIPATED END-OF-FILE, GAUSSIAN PRIOR') 120 CALL PEND('ERROR READING GAUSSIAN PRIOR') END REAL*8 FUNCTION GAU1MX(NMIX,AMU,SIG,P) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AMU(NMIX),SIG(NMIX),P(NMIX) J=NDISC(NMIX,P) GAU1MX=AMU(J)+SIG(J)*DRNNOF() RETURN END SUBROUTINE GAU2(K,H0,LDH0,H1,LDH1,HB0,HB1,THETA1,THETA2) C This routine generates two, antithetic, drawings from the posterior C distribution, given a Gaussian prior and Gaussian likelihood. C Inputs: C K Order of distribution C H0 Prior precision matrix C LDH0 Leading dimension of H0 C H1 Data precision matrix C LDH1 Leading dimension of H1 C HB0 Product of prior precision matrix and mean C HB1 Product of data precision matrix and mean C Outputs: C THETA1 Drawing from postrior C THETA2 Antithetic of THETA1 IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION H0(LDH0,K),H1(LDH1,K),HB0(K),HB1(K), 1 THETA1(K),THETA2(K) C Compute posterior precision CALL UMADD(K,K,H0,LDH0,H1,LDH1,A2,LD) C Compute posterior mean CALL UVADD(K,HB0,HB1,V1) CALL DLSLDS(K,A2,LD,V1,V2) C Draw realization CALL GAU(K,V2,A2,LD,THETA1,THETA2) RETURN END SUBROUTINE GAUI(NP,AMU,H,LDH,BETA) C This routine generates a random vector from a Gaussian distribution with C mean vector mu and precision matrix H, subject to the restriction that the C vector correspond to an invertible lag operator. C Inputs: C NP Dimension of vector C AMU Mean vector C H Precision matrix H C LDH Leading dimension of H C Output: C BETA Random vecctor IMPLICIT REAL*8 (A-H,O-Z) LOGICAL LFLIP,LINOP PARAMETER(LD=100) DIMENSION AMU(NP),H(LDH,NP),BETA(NP) DIMENSION V1(LD),V2(LD) DATA NREP1,NREP2/250,250/ C Make independent draws from normal with rejection. CALL GAU(NP,AMU,H,LDH,V1,V2) IF(LINOP(NP,V1))GO TO 20 DO 10 IREP=1,NREP1 CALL GAUA(NP,AMU,H,LDH,V1,V2) 10 IF(LINOP(NP,V1))GO TO 20 GO TO 100 20 CALL DCOPY(NP,V1,1,BETA,1) RETURN C Make draws using independence Metropolis with uniform p.a.c.f. source. C Separate entry can be used to force independence Metropolis algorithm C (mostly for testing purposes). ENTRY GAUI1(NP,AMU,H,LDH,BETA) 100 WTLOG0=DMACH(8) DO 120 IREP=1,NREP2 DO 110 IP=1,NP 110 V1(IP)=-1.0D0+2.0D0*DRNUNF() CALL PACFAR(NP,V1,V2) WTLOG=GAUK(NP,AMU,H,LDH,V2)+ARPALJ(NP,V1) IF(.NOT.LFLIP(WTLOG-WTLOG0))GO TO 120 CALL DCOPY(NP,V2,1,BETA,1) WTLOG0=WTLOG 120 CONTINUE RETURN C Auxilliary entry to provide a weighted sample (mostly for testing purposes) C Additional output: C WTL Log weight, for importance sampling ENTRY GAUI2(NP,AMU,H,LDH,BETA,WTL) DO 210 IP=1,NP 210 V1(IP)=-1.0D0+2.0D0*DRNUNF() CALL PACFAR(NP,V1,BETA) WTL=GAUK(NP,AMU,H,LDH,BETA)+ARPALJ(NP,V1) RETURN END C REAL*8 FUNCTION GAUIN(NP,AMU,H,LDH,SD,TIME) C This function evaluates the log normalizing constant for the log density C kernel of a random vector computed in GAUK, given an invertibility C restriction. This normalizing constant, when multiplied by the kernel C in GAUK, yields the p.d.f. of the random vector. C Inputs: C NP Dimension of vector C AMU Mean vector C H Precision matrix H C LDH Leading dimension of H C SD Routine stops when standard deviation of approximation error C of log kernel reaches this point. C TIME Routine stops when this much CPU time has been used. C Output: C GAUIN Log normalizing constanat (approximate C SD Standard error of approximation C TIME CPU time used IMPLICIT REAL*8 (A-H,O-Z) REAL*4 CPSEC LOGICAL LINOP PARAMETER(LD=100) DIMENSION AMU(NP),H(LDH,NP) DIMENSION V1(LD),V2(LD) DATA MREP1,MREP2,JUMP1,JUMP2/500,1000,500,1000/ C Measure efficiency of sampling from untruncated prior. CLOCK=CPSEC() AA1=0.0D0 SS1=0.0D0 DO 10 IREP=1,MREP1 CALL GAU(NP,AMU,H,LDH,V1,V2) BB=0.0D0 IF(LINOP(NP,V1))BB=BB+0.5D0 IF(LINOP(NP,V2))BB=BB+0.5D0 AA1=AA1+BB 10 SS1=SS1+BB**2 P=AA1/MREP1 IF(P.EQ.0.0D0)GO TO 200 IF(P.EQ.1.0D0)GO TO 300 VARLOG=(SS1-MREP1*P**2)/(MREP1*(MREP1-1)*P**2) IF((CPSEC()-CLOCK).EQ.0)GO TO 100 EFF1=1.0D0/(VARLOG*(CPSEC()-CLOCK)) C Measure efficiency of sampling uniformly from p.a.c.f. CLOCK2=CPSEC() AA2=0.0D0 SS2=0.0D0 DO 30 IREP=1,MREP2 DO 20 IP=1,NP 20 V1(IP)=2.0D0*DRNUNF()-1.0D0 CALL PACFAR(NP,V1,V2) BB=DEXP(ARPALJ(NP,V1)+GAUK(NP,AMU,H,LDH,V2)) AA2=AA2+BB 30 SS2=SS2+BB**2 C=(2.0D0**NP)*DEXP(GAUN(NP,H,LDH)) P=C*AA2/MREP2 VARLOG=(SS2*C**2-MREP2*P**2)/(MREP2*(MREP2-1)*P**2) EFF2=1.0D0/(VARLOG*(CPSEC()-CLOCK2)) IF(EFF2.GT.EFF1)GO TO 200 C Approximate by sampling from untruncated prior. 100 NREP1=MREP1 110 DO 120 JUMP=1,JUMP1 CALL GAU(NP,AMU,H,LDH,V1,V2) BB=0.0D0 IF(LINOP(NP,V1))BB=BB+0.5D0 IF(LINOP(NP,V2))BB=BB+0.5D0 AA1=AA1+BB 120 SS1=SS1+BB**2 NREP1=NREP1+JUMP1 P=AA1/NREP1 SDLOG=DSQRT((SS1-NREP1*P**2)/(NREP1*(NREP1-1)*P**2)) IF(SDLOG.LT.SD)GO TO 300 IF(CPSEC()-CLOCK.GT.TIME)GO TO 300 GO TO 110 C Approximate by sampling uniformly from p.a.c.f. 200 NREP2=MREP2 210 DO 230 JUMP=1,JUMP2 DO 220 IP=1,NP 220 V1(IP)=2.0D0*DRNUNF()-1.0D0 CALL PACFAR(NP,V1,V2) BB=DEXP(ARPALJ(NP,V1)+GAUK(NP,AMU,H,LDH,V2)) AA2=AA2+BB 230 SS2=SS2+BB**2 NREP2=NREP2+JUMP2 P=C*AA2/NREP2 SDLOG=DSQRT((SS2*C**2-NREP2*P**2)/(NREP2*(NREP2-1)*P**2)) IF(SDLOG.LT.SD)GO TO 300 IF((CPSEC()-CLOCK).GT.TIME)GO TO 300 GO TO 210 C Return parameters and value. 300 SD=SDLOG TIME=CPSEC()-CLOCK GAUIN=GAUN(NP,H,LDH)-DLOG(P) RETURN END SUBROUTINE GAUI0(NFILE,NK,TITLE,H0,LDH,B0,HB0,ANORM,THETA) C This routine reads in the mean vector mu and precision matrix H of a C kx1 parameter vector theta. It computes and returns the product of C H and mu, and prints information about the prior. C Inputs: C NFILE File number for reading prior C NK Dimension of parameter vector = k C TITLE Identifier for distribution (printed) C LDH Leading dimension of H0 C Ouptuts: C H0 Prior precision matrix H C B0 Prior mean vector mu C HB0 Product of H and mu C ANORM Log normalizing constant for prior density C THETA Draw from prior IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) TITLE PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION H0(LDH,NK),HB0(NK),B0(NK),THETA(NK) READ(NFILE,*,END=110,ERR=120)NCODE1,NCODE2,NDIM C NCODE1=1 Diagonal, sqrts only C NCODE1=2 Sqrts diagonal followed by normalized matrix C NCODE1=3 Full matrix C NCODE2=1 Variance matrix C NCODE2=2 Precision matrix IF(NDIM.NE.NK)CALL PEND('ORDER OF GAUSSIAN PRIOR INCORRECT') CALL PINTI('Gaussian prior code 1',NCODE1,1,3) CALL PINTI('Gaussian prior code 2',NCODE2,1,2) READ(NFILE,*,END=110,ERR=120)(B0(IK),IK=1,NK) CALL DSET(NK,1.0D0,V1,1) IF(NCODE1.NE.3)CALL PLISTR(NFILE,NK,V1) CALL UMISET(NK,H0,LDH) IF(NCODE1.NE.1)CALL PDSYMR(NFILE,NK,H0,LDH) DO 10 IK=1,NK DO 10 JK=1,IK H0(IK,JK)=H0(IK,JK)*V1(IK)*V1(JK) 10 H0(JK,IK)=H0(IK,JK) IF(NCODE2.EQ.1)CALL DLINDS(NK,H0,LDH,H0,LDH) C Form Hx(mu). CALL DMURRV(NK,NK,H0,LDH,NK,B0,1,NK,HB0) C Normalizing constant and nitial draw TIME=60.0D0 SD=0.01D0 ANORM=GAUIN(NK,B0,H0,LDH,SD,TIME) CALL GAUI(NK,B0,H0,LDH,THETA) C Variances as convenience to user CALL DLINDS(NK,H0,LDH,A1,LD) WRITE(6,20)TITLE,(JK,B0(JK),DSQRT(A1(JK,JK)),H0(JK,JK), 1 THETA(JK),JK=1,NK) 20 FORMAT(/,' Invertible Gaussian prior for ',A,/, 1 ' Parameter',6X,'Prior mean',6X,'Prior s.d.',5X, 2 'Prior prec.',3X,'Initial value',/,(5X,I5,4(1PE16.6))) WRITE(6,30)SD,TIME 30 FORMAT(/,' .. Means and standard deviations do not reflect', 1 ' invertibility constraint.',/,' Standard deviation of', 2 ' approximation of log prior kernel:',1PE15.5,/, 3 ' Computation time:',0PF8.2) RETURN 110 CALL PEND('UNANTICIPATED END-OF-FILE, INIVERTIBLE GAUSSIAN PRIOR') 120 CALL PEND('ERROR READING GAUSSIAN PRIOR') END SUBROUTINE GAUT(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER,X) C This routine draws from a truncated multivariate normal distribution C without rejection. The truncation is formed by linear restrictions C equal in number to the dimension of the distribution, +infinity and C -infinity allowed. The truncated distribution is C X ~ N(AMU, SIGMA), A < D*X < B C References: J. Geweke, 1991. "Efficient Simulation from the Multivariate C Normal and Student-t Distributions Subject to Linear C Constraints," in E.M. Keramidas (ed.), Computing Science and C Statistics: Proceedings of the 23rd Symposium on the C Interface, pp. 571-578. Fairfax, VA: Interface Foundation C of North America, Inc. C Inputs: C N Dimension of the multivariate normal distribution C AMU Mean of the distribution C H(NxN) Precision matrix of distribution C LDX Leading dimension of H, D, and DINV in calling program C D(NxN) Matrix of linear combinations of X for constraints C DINV(NxN) Inverse of D C A(N) Lower bounds for linear combinations C B(N) Upper bounds for linear combinations C LA(N) If .TRUE., no corresponding lower bound C LB(N) If .TRUE., no corresponding upper bound C LE(N) If .TRUE., equality constraint C NITER Number of Gibbs iterations to perform C Input/output: C X Random vector IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) LOGICAL LA(N),LB(N),LE(N) COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION AMU(N),H(LDX,N),DINV(LDD,N),A(N),B(N),X(N) C Form mean of Z =Nu in V1 CALL DMRRRR(N,N,D,LDD,N,1,AMU,N,N,1,V1,N) C Form precision of Z in A1 CALL UXTAXF(N,N,DINV,LDD,H,LDX,A1,LD) C Retrieve current Z vector in V2 CALL DMRRRR(N,N,D,LDD,N,1,X,N,N,1,V2,N) CALL UVSUB(N,V2,V1,V2) C Make drawings for Z=D*X. DO 20 ITER=1,NITER DO 20 I=1,N IF(LE(I))GO TO 20 AA=0.0D0 DO 10 J=1,N 10 IF(I.NE.J)AA=AA-(A1(I,J)/A1(I,I))*V2(J) V2(I)=XNOT(AA,A1(I,I),A(I)-V1(I),B(I)-V1(I),LA(I),LB(I)) 20 CONTINUE C Transform from Z to X (IMSL-MATH p. 982). CALL DMRRRR(N,N,DINV,LDD,N,1,V2,N,N,1,V1,N) CALL UVADD(N,V1,AMU,X) RETURN END REAL*8 FUNCTION GAUTK(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER,X) C This function returns the log kernel of a truncated normal distribution. C If the constraints are satisfied then the kernel is the same as for the C untruncated distribution. If they are violated, a small negative number C is returned. C Inputs: C N Dimension of the multivariate normal distribution C AMU Mean of the distribution C H(NxN) Precision matrix of distribution C LDX Leading dimension of H, D, and DINV in calling program C D(NxN) Matrix of linear combinations of X for constraints C DINV(NxN) Inverse of D C A(N) Lower bounds for linear combinations C B(N) Upper bounds for linear combinations C LA(N) If .TRUE., no corresponding lower bound C LB(N) If .TRUE., no corresponding upper bound C LE(N) If .TRUE., equality constraint C NITER Number of Monte Carlo draws to use C X(N) Point of kernel evaluation C Output: C GAUTK Log kernel IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION AMU(N),H(LDX,N),DINV(LDD,N),A(N),B(N) LOGICAL LA(N),LB(N),LE(N) CALL DMURRV(N,N,D,LDD,N,X,1,N,V1) DO 10 I=1,N IF((.NOT.LA(I)).AND.(V1(I).LT.A(I)))GO TO 20 10 IF((.NOT.LB(I)).AND.(V1(I).GT.B(I)))GO TO 20 GAUTK=GAUK(N,AMU,H,LDX,X) RETURN 20 GAUTK=-700.0D0 RETURN END REAL*8 FUNCTION GAUTN(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER) C This function returns the log normalizing factor of a truncated normal C distribution. The normalilzing factor, when multiplied by the kernel C evaluated in GAUTK, yields the properly normalized probability density C of the truncated normal distribution. The normalizing factor is the C one for the untruncated normal distribution, divided by the probability C that the untruncated normal would satisfy the inequalities in the C truncated normal distribution. This probability is approximated using C the GHK simulation algorithm. C Inputs: C N Dimension of the multivariate normal distribution C AMU Mean of the distribution C H(NxN) Precision matrix of distribution C LDX Leading dimension of H, D, and DINV in calling program C D(NxN) Matrix of linear combinations of X for constraints C DINV(NxN) Inverse of D C A(N) Lower bounds for linear combinations C B(N) Upper bounds for linear combinations C LA(N) If .TRUE., no corresponding lower bound C LB(N) If .TRUE., no corresponding upper bound C LE(N) If .TRUE., equality constraint C NITER Number of Monte Carlo draws to use in the GHK algorithm C Output: C GAUTN Log normalizing factor IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION AMU(N),H(LDX,N),DINV(LDD,N),A(N),B(N) GAUTN=GAUN(N,H,LDX) 1 -DLOG(GHK(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER,SE)) RETURN END REAL*8 FUNCTION GHK(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,SD,TIME) C This routine approximates the log probability that a random variable C with the distribution X~N(AMU,H^-1) satisfies the inequality constraints C AA(i+1) for some i IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(N) DO 10 I=2,N 10 IF(A(I).LT.A(I-1))GO TO 20 LORDER=.TRUE. RETURN 20 LORDER=.FALSE. RETURN END REAL*8 FUNCTION UMIN(A,B) IMPLICIT REAL*8 (A-H,O-Z) C=A IF(C.GT.B)C=B UMIN=C RETURN END REAL*8 FUNCTION UMAX(A,B) IMPLICIT REAL*8 (A-H,O-Z) C=A IF(C.LT.B)C=B UMAX=C RETURN END INTEGER FUNCTION NDISC(N,P) C This routine randomly selects one from among N categories, where the C probabilities are proportional to the vector P. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N) DATA CC/1.00000008D0/ U=DRNUNF()*DSUM(N,P,1) AA=0.0D0 DO 10 I=1,N AA=AA+P(I) IF(U.LE.AA)GO TO 20 10 CONTINUE IF(U.GT.(AA*CC))GO TO 30 NDISC=N RETURN 20 NDISC=I RETURN 30 CALL TERM('LIMITS EXCEEDED IN NDISC') END SUBROUTINE NMXD(N,NMIX,LDMIX,IDCODE,MXCODE,X, 1 H0,HA0,S0NU0,ANU0,R,ALPHA,H,P,LT) C This routine draws from the posterior distribution for a discrete C mixture of normals density. C Inputs: C N Sample size C NMIX Number of normals in mixture C LDMIX Leading dimension of H0 C IDCODE =1: Mean vector ALPHA in nondescending order C =2: Precision vector H in nondescending order C =3: Probability vector P in nondescending order C Otherwise, no identifying restrictions beyond prior C distributions applied C MXCODE = 1: Scale mixture of normals (means fixed at 0) C = 2: Mean mixture of normals (precisions all the same) C Otherwise, both means and precisions are free C X Nx1 Data vector C H0 NMIX x NMIX prior precision for Alpha (ignored if MXCODE=1) C HA0 Product of prior precision and prior mean for Alpha C (ignored if MXCODE=1) C S0NU0 Vector of constant parameters of chi-square priors for H C (only first element if MXCODE=2) C ANU0 Vector of degrees-of-freedom of chi-square priors for H C (only first element if MXCODE=2) C R Vector of parameters of multivariate beta prior for P C Input/output: C ALPHA NMIX x 1 vector of means. If MXCODE not 1, NMXD C exchanges values in Gibbs sampler. If MXCODE=1, C all elements of ALPHA1 must be set the same on input; C values then left unchanged by NMXD C H NMIX x 1 vector of precisons. If MXCODE not 2, NMXD C exchanges values in Gibbs sampler. If MXCODE=2, C H(1) is the common precision and only H(1) is exchanged C P NMIX x 1 vector of probabilities. NMXD exchanges values C in Gibbs sampler C LT N x 1 vector of integers between 1 and NMIX inclusive C indicating mixture status for each observation C (Need not be input, since LT is drawn first) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=10) COMMON/AMCMC/WT,PRI,PDAT,PAR(250),ITER,NPAR,LRECRD,LWRITE,LERR COMMON/NMXDBL/D(LD,LD),DINV(LD,LD),A(LD),B(LD),LA(LD),LB(LD), 1 LE(LD) LOGICAL LA,LB,LE DIMENSION X(N),H0(LDMIX,NMIX),HA0(NMIX),S0NU0(NMIX),ANU0(NMIX), 1 R(NMIX),ALPHA(NMIX),H(NMIX),P(NMIX),LT(N) DIMENSION A1(LD,LD),A2(LD,LD),V1(LD),V2(LD),V3(LD),NT(LD) C Draw latent classifications 100 DO 110 J=1,NMIX V1(J)=P(J)*DSQRT(H(J)) V2(J)=-.5D0*H(J) 110 NT(J)=0 DO 130 I=1,N DO 120 J=1,NMIX 120 V3(J)=V1(J)*DEXP(V2(J)*(X(I)-ALPHA(J))**2) L=NDISC(NMIX,V3) LT(I)=L 130 NT(L)=NT(L)+1 C Draw mean vector 200 IF(MXCODE.EQ.1)GO TO 300 CALL UM0SET(NMIX,NMIX,A1,LD) CALL DSET(NMIX,0.00D0,V1,1) DO 210 I=1,N L=LT(I) 210 V1(L)=V1(L)+H(L)*X(I) DO 220 J=1,NMIX 220 A1(J,J)=H(J)*NT(J) NTRY=0 230 NTRY=NTRY+1 IF(IDCODE.EQ.1)GO TO 240 CALL GAU2(NMIX,H0,LDMIX,A1,LD,HA0,V1,ALPHA,V2) GO TO 300 240 CALL UMADD(NMIX,NMIX,H0,LDMIX,A1,LD,A2,LD) CALL UVADD(NMIX,HA0,V1,V2) CALL DLSLDS(NMIX,A2,LD,V2,V3) CALL GAUT(NMIX,V3,A2,LD,D,DINV,LD,A,B,LA,LB,LE,5,ALPHA) C Draw precision vector 300 IF(MXCODE.EQ.2)GO TO 340 CALL DCOPY(NMIX,S0NU0,1,V1,1) DO 310 I=1,N L=LT(I) 310 V1(L)=V1(L)+(X(I)-ALPHA(L))**2 DO 315 J=1,NMIX 315 V2(J)=ANU0(J)+NT(J) IF(IDCODE.EQ.2)GO TO 330 DO 320 J=1,NMIX 320 H(J)=CHI(V1(J),V2(J)) GO TO 400 330 CALL CHIO(NMIX,V1,V2,1,H) GO TO 400 340 AA=S0NU0(1) DO 350 I=1,N 350 AA=AA+(X(I)-ALPHA(1))**2 H(1)=CHI(AA,ANU0(1)+N) CALL DSET(NMIX-1,H(1),H(2),1) C Draw probability vector 400 DO 410 J=1,NMIX 410 V1(J)=R(J)+NT(J) 420 CALL BETAM(NMIX,V1,P) IF(IDCODE.NE.3)RETURN DO 430 J=2,NMIX 430 IF(P(J).LT.P(J-1))GO TO 420 RETURN END SUBROUTINE NMXD0(NFILE,TITLE,NMIX,IDCODE,MXCODE, 1 AH0,LDMIX,AB0,AHB0,ALPHA,ANU0,S0NU0,H,R,P,PRI0) C This routine reads in parameters for priors for a discrete mixture of C normals distribution. C Inputs: C NFILE File number for read C TITLE Informative banner for log file C NMIX Number of components of mixture C IDCODE Orderings indicator; see NMXD C MXCODE Free parameter indicator; see NMXD C LDMIX Leading dimension of AH0 C Outputs: C AH0 Prior precision matrix for coefficient Gaussian prior C AB0 Prior mean vector for coefficient Gaussian prior C AHB0 Product of prior precision matrix and mean vector C ALPHA Initial draw of mean vector C ANU0 Degrees of freedom parameters for precision priors C S0NU0 Constant parameters for precision priors C H Initial draw of precision parameters C R Prior parameters for probability beta distribution C P Initial draw of probability parameters C PRI0 Log normalization constant for prior IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) TITLE DIMENSION AH0(LDMIX,NMIX),AB0(NMIX),AHB0(NMIX),ALPHA(NMIX), 1 ANU0(NMIX),S0NU0(NMIX),H(NMIX),R(NMIX),P(NMIX) PARAMETER(LD=10) COMMON/NMXDBL/D(LD,LD),DINV(LD,LD),A(LD),B(LD),LA(LD),LB(LD), 1 LE(LD) LOGICAL LA,LB,LE READ(1,*)NMIX,IDCODE,MXCODE CALL PINTI('Number of mixtures',NMIX,2,LDMIX) CALL PINTI('Identification code',IDCODE,0,3) CALL PINTI('Mixture code',MXCODE,0,2) WRITE(6,10)TITLE 10 FORMAT(//,' Discrete normal mixture model for ',A) CALL PLINE1(1,' Disturbance distribution is a mixture of',NMIX, 1 'normals') IF(IDCODE.EQ.0) 1 CALL PLINE(0,' No inequality restrictions on coefficients') IF(IDCODE.EQ.1) 1 CALL PLINE(0,' Mean vector in nondecreasing order') IF(IDCODE.EQ.2) 1 CALL PLINE(0,' Precision vector in nondecreasing order') IF(IDCODE.EQ.3) 1 CALL PLINE(0,' Probability vector in nondecreasing order') IF(MXCODE.EQ.0)CALL PLINE(0,' Full mixture of normals') IF(MXCODE.EQ.1)CALL PLINE(0,' Scale mixture of normals') IF(MXCODE.EQ.2)CALL PLINE(0,' Mean mixture of normals') C Obtain priors for normal mixture model PRI0A=0.0D0 IF(MXCODE.EQ.1)GO TO 40 CALL GAU0(NFILE, 1 NMIX,'mixture coefficients',AH0,LDMIX,AB0,AHB0,PRI0A,ALPHA) IF(IDCODE.NE.1)GO TO 40 CALL UMISET(NMIX,D,LD) DO 20 I=1,NMIX IF(I.LT.NMIX)D(I,I+1)=-1.0D0 A(I)=-1.0D0 B(I)=0.0D0 LA(I)=.TRUE. 20 LB(I)=.FALSE. LB(NMIX)=.TRUE. CALL DLINRG(NMIX,D,LD,DINV,LD) SD=0.01D0 TIME=60.0D0 PRI0A=PRI0A-GHK(NMIX,AB0,AH0,LDMIX,D,DINV,LD,A,B,LA,LB,LE,SD,TIME) CALL GAUT(NMIX,AB0,AH0,LDMIX,D,DINV,LD,A,B,LA,LB,LE,100,ALPHA) WRITE(6,30)SD,TIME,(ALPHA(I),I=1,NMIX) 30 FORMAT(/,' ...Means and standard deviations etc. do not reflect', 1 ' ordering constraint.',/,' Standard deviation of', 2 ' approximation of log prior kernel:',1PE15.5,/, 3 ' Computation time:',0PF8.2,/, 4 ' Initial draw from ordered distribution:',/, 5 (5(1PE16.7))) 40 PRI0H=0.0D0 IF(MXCODE.EQ.2)GO TO 60 CALL CHI0(NFILE,NMIX,'mixture precisions',IDCODE.EQ.2, 1 ANU0,S0NU0,PRI0H,H) IF(IDCODE.NE.2)GO TO 70 CALL CHIO(NMIX,S0NU0,ANU0,20,H) WRITE(6,50)(H(I),I=1,NMIX) 50 FORMAT(' Initial draw from ordered distribution:',/,(5(1PE16.7))) GO TO 70 60 CALL CHI0(NFILE,1,'common disturbance precision',.FALSE., 1 ANU0,S0NU0,PRI0H,H) CALL DSET(NMIX-1,ANU0(1),ANU0(2),1) CALL DSET(NMIX-1,S0NU0(1),S0NU0(2),1) CALL DSET(NMIX-1,H(1),H(2),1) 70 PRI0B=0.0D0 CALL BET0(NFILE,1,'mixture probabilities',NMIX,R,1,PRI0B,P,1) PRI0=PRI0A+PRI0H+PRI0B RETURN END SUBROUTINE NMXD0P(NFILE,TITLE,NMIX,IDCODE,MXCODE, 1 AH0,LDMIX,AB0,AHB0,ALPHA,ANU0,S0NU0,H,R,P,JFIX,PRI0) C This routine reads in parameters for priors for a discrete mixture of C normals distribution, for the case in which one precision is fixed. C Inputs: C NFILE File number for read C TITLE Informative banner for log file C NMIX Number of components of mixture C IDCODE Orderings indicator; see NMXD C MXCODE Free parameter indicator; see NMXD C LDMIX Leading dimension of AH0 C Outputs: C AH0 Prior precision matrix for coefficient Gaussian prior C AB0 Prior mean vector for coefficient Gaussian prior C AHB0 Product of prior precision matrix and mean vector C ALPHA Initial draw of mean vector C ANU0 Degrees of freedom parameters for precision priors C S0NU0 Constant parameters for precision priors C H Initial draw of precision parameters C R Prior parameters for probability beta distribution C P Initial draw of probability parameters C JFIX Precision whose value is fixed C PRI0 Log normalization constant for prior IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) TITLE DIMENSION AH0(LDMIX,NMIX),AB0(NMIX),AHB0(NMIX),ALPHA(NMIX), 1 ANU0(NMIX),S0NU0(NMIX),H(NMIX),R(NMIX),P(NMIX) PARAMETER(LD=10) COMMON/NMXDBL/D(LD,LD),DINV(LD,LD),A(LD),B(LD),LA(LD),LB(LD), 1 LE(LD) LOGICAL LA,LB,LE READ(1,*)NMIX,IDCODE,MXCODE,JFIX,HFIX CALL PINTI('Number of mixtures',NMIX,2,LDMIX) CALL PINTI('Identification code',IDCODE,0,3) CALL PINTI('Mixture code',MXCODE,0,2) CALL PINTI('Fixed precision',JFIX,1,NMIX) WRITE(6,10)TITLE,JFIX,HFIX 10 FORMAT(//,' Discrete normal mixture model for ',A,/, 1 ' Precision',I2,' fixed at',1PE13.6) CALL PLINE1(1,' Disturbance distribution is a mixture of',NMIX, 1 'normals') IF(IDCODE.EQ.0) 1 CALL PLINE(0,' No inequality restrictions on coefficients') IF(IDCODE.EQ.1) 1 CALL PLINE(0,' Mean vector in nondecreasing order') IF(IDCODE.EQ.2) 1 CALL PLINE(0,' Precision vector in nondecreasing order') IF(IDCODE.EQ.3) 1 CALL PLINE(0,' Probability vector in nondecreasing order') IF(MXCODE.EQ.0)CALL PLINE(0,' Full mixture of normals') IF(MXCODE.EQ.1)CALL PLINE(0,' Scale mixture of normals') IF(MXCODE.EQ.2)CALL PLINE(0,' Mean mixture of normals') C Obtain priors for normal mixture model PRI0A=0.0D0 IF(MXCODE.EQ.1)GO TO 40 CALL GAU0(NFILE, 1 NMIX,'mixture coefficients',AH0,LDMIX,AB0,AHB0,PRI0A,ALPHA) IF(IDCODE.NE.1)GO TO 40 CALL UMISET(NMIX,D,LD) DO 20 I=1,NMIX IF(I.LT.NMIX)D(I,I+1)=-1.0D0 A(I)=-1.0D0 B(I)=0.0D0 LA(I)=.TRUE. 20 LB(I)=.FALSE. LB(NMIX)=.TRUE. CALL DLINRG(NMIX,D,LD,DINV,LD) SD=0.01D0 TIME=60.0D0 PRI0A=PRI0A-GHK(NMIX,AB0,AH0,LDMIX,D,DINV,LD,A,B,LA,LB,LE,SD,TIME) CALL GAUT(NMIX,AB0,AH0,LDMIX,D,DINV,LD,A,B,LA,LB,LE,100,ALPHA) WRITE(6,30)SD,TIME,(ALPHA(I),I=1,NMIX) 30 FORMAT(/,' ...Means and standard deviations etc. do not refect', 1 ' ordering constraint.',/,' Standard deviation of', 2 ' approximation of log prior kernel:',1PE15.5,/, 3 ' Computation time:',0PF8.2,/, 4 ' Initial draw from ordered distribution:',/, 5 (5(1PE16.7))) 40 IF(MXCODE.EQ.2)GO TO 60 CALL CHI0P(NFILE,NMIX,'mixture precisions',IDCODE.EQ.2,JFIX,HFIX, 1 ANU0,S0NU0,PRI0H,H) IF(IDCODE.NE.2)GO TO 70 CALL CHIOP(NMIX,S0NU0,ANU0,20,JFIX,H) WRITE(6,50)(H(I),I=1,NMIX) 50 FORMAT(' Initial draw from ordered distribution:',/,(5(1PE16.7))) GO TO 70 60 CALL DSET(NMIX,HFIX,H,1) 70 CALL BET0(NFILE,1,'mixture probabilities',NMIX,R,1,PRI0B,P,1) PRI0=PRI0A+PRI0H+PRI0B RETURN END SUBROUTINE NMXDP(N,NMIX,LDMIX,IDCODE,MXCODE,X, 1 H0,HA0,S0NU0,ANU0,R,JFIX,ALPHA,H,P,LT) C This routine draws from the posterior distribution for a discrete C mixture of normals density, with one precision fixed at a prespecified C value. C Inputs: C N Sample size C NMIX Number of normals in mixture C LDMIX Leading dimension of H0 C IDCODE =1: Mean vector ALPHA in nondescending order C =2: Precision vector H in nondescending order C =3: Probability vector P in nondescending order C Otherwise, no identifying restrictions beyond prior C distributions applied C MXCODE = 1: Scale mixture of normals (means fixed at 0) C = 2: Mean mixture of normals (precisions all the same) C Otherwise, both means and precisions are free C X Nx1 Data vector C H0 NMIX x NMIX prior precision for Alpha (ignored if MXCODE=1) C HA0 Product of prior precision and prior mean for Alpha C (ignored if MXCODE=1) C S0NU0 Vector of constant parameters of chi-square priors for H C (only first element if MXCODE=2) C ANU0 Vector of degrees-of-freedom of chi-square priors for H C (only first element if MXCODE=2) C R Vector of parameters of multivariate beta prior for P C JFIX Precision whose value is fixed C Input/output: C ALPHA NMIX x 1 vector of means. If MXCODE not 1, NMXD C exchanges values in Gibbs sampler. If MXCODE=1, C all elements of ALPHA1 must be set the same on input; C values then left unchanged by NMXD C H NMIX x 1 vector of precisons. If MXCODE not 2, NMXD C exchanges values in Gibbs sampler. If MXCODE=2, C H(1) is the common precision and only H(1) is exchanged C P NMIX x 1 vector of probabilities. NMXD exchanges values C in Gibbs sampler C LT N x 1 vector of integers between 1 and NMIX inclusive C indicating mixture status for each observation C (Need not be input, since LT is drawn first) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=10) COMMON/AMCMC/WT,PRI,PDAT,PAR(250),ITER,NPAR,LRECRD,LWRITE,LERR COMMON/NMXDBL/D(LD,LD),DINV(LD,LD),A(LD),B(LD),LA(LD),LB(LD), 1 LE(LD) LOGICAL LA,LB,LE DIMENSION X(N),H0(LDMIX,NMIX),HA0(NMIX),S0NU0(NMIX),ANU0(NMIX), 1 R(NMIX),ALPHA(NMIX),H(NMIX),P(NMIX),LT(N) DIMENSION A1(LD,LD),A2(LD,LD),V1(LD),V2(LD),V3(LD),NT(LD) C Draw latent classifications 100 DO 110 J=1,NMIX V1(J)=P(J)*DSQRT(H(J)) V2(J)=-.5D0*H(J) 110 NT(J)=0 DO 130 I=1,N DO 120 J=1,NMIX 120 V3(J)=V1(J)*DEXP(V2(J)*(X(I)-ALPHA(J))**2) L=NDISC(NMIX,V3) LT(I)=L 130 NT(L)=NT(L)+1 C Draw mean vector 200 IF(MXCODE.EQ.1)GO TO 300 CALL UM0SET(NMIX,NMIX,A1,LD) CALL DSET(NMIX,0.00D0,V1,1) DO 210 I=1,N L=LT(I) 210 V1(L)=V1(L)+H(L)*X(I) DO 220 J=1,NMIX 220 A1(J,J)=H(J)*NT(J) NTRY=0 230 NTRY=NTRY+1 IF(IDCODE.EQ.1)GO TO 240 CALL GAU2(NMIX,H0,LDMIX,A1,LD,HA0,V1,ALPHA,V2) GO TO 300 240 CALL UMADD(NMIX,NMIX,H0,LDMIX,A1,LD,A2,LD) CALL UVADD(NMIX,HA0,V1,V2) CALL DLSLDS(NMIX,A2,LD,V2,V3) CALL GAUT(NMIX,V3,A2,LD,D,DINV,LD,A,B,LA,LB,LE,5,ALPHA) C Draw precision vector 300 IF(MXCODE.EQ.2)GO TO 400 CALL DCOPY(NMIX,S0NU0,1,V1,1) DO 310 I=1,N L=LT(I) 310 V1(L)=V1(L)+(X(I)-ALPHA(L))**2 DO 315 J=1,NMIX 315 V2(J)=ANU0(J)+NT(J) IF(IDCODE.EQ.2)GO TO 330 DO 320 J=1,NMIX 320 IF(J.NE.JFIX)H(J)=CHI(V1(J),V2(J)) GO TO 400 330 CALL CHIOP(NMIX,V1,V2,1,JFIX,H) C Draw probability vector 400 DO 410 J=1,NMIX 410 V1(J)=R(J)+NT(J) 420 CALL BETAM(NMIX,V1,P) IF(IDCODE.NE.3)RETURN DO 430 J=2,NMIX 430 IF(P(J).LT.P(J-1))GO TO 420 RETURN END SUBROUTINE NOWPR(NFILE) C This routine prints the current day of the week, date, and time C to the unit NFILE (which of course must have been opened). CHARACTER*10 HWEEK(7) DATA HWEEK/' Sunday',' Monday',' Tuesday',' Wednesday', 1 ' Thursday',' Friday',' Saturday'/ CALL TDATE(IDATE,MONTH,IYEAR) CALL TIMDY(IHOUR,MINUTE,ISEC) IW=IDYWK(IDATE,MONTH,IYEAR) C THis next line is for those who are using the IMSL independent version C of JGLIB. The programs called here use IMSL to get machine dependent C values. IF(IW.EQ.0) GO TO 100 WRITE(NFILE,10)HWEEK(IW),MONTH,IDATE,IYEAR,IHOUR,MINUTE,ISEC 10 FORMAT(A10,I3,'/',I2,'/',2I4,2(':',I2)) 100 RETURN END SUBROUTINE NRMK(NK,NT,X,LDX,Y,H,HPRI,LDHPRI,HBPRI,BETA) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This routine simulates one draw from the coefficient posterior distribution C arising from a regression model with a normal prior for the coefficients. C Inputs: C NK Dimension of coefficient vector (= number of regressors) C NT Number of observations C X NTxNK matrix of regressors (DESTROYED) C LDX Leading dimension of X C Y NTx1 vector of dependent variables C H Disturbance precision C HPRI NKxNK prior precision matrix (symmetric) C LDHPRI Leading dimension of HPRI C HBPRI [HPRI]x[Coefficient prior mean vector] C Output: C BETA NKx1 Draw from posterior distribution C PDFLG Log posterior p.d.f. at the drawn value. C ... When finished, the upper triangle of the posterior precision is in A1, C the posterior mean is in V2, C and the Choleski decomposition of the posterior variance is in A2. IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) C COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION X(LDX,NK),Y(NT),HPRI(LDHPRI,NK),HBPRI(NK),BETA(NK) C Form upper triangle of posterior precision in A1 and RHS vector in V1. DO 10 I=1,NK V1(I)=HBPRI(I)+H*DDOT(NT,X(1,I),1,Y,1) DO 10 J=I,NK 10 A1(I,J)=HPRI(I,J)+H*DDOT(NT,X(1,I),1,X(1,J),1) C Form upper triangular R: R'R = Precision in A2, and invert. CALL DLFTDS(NK,A1,LD,A2,LD) CALL DLINRT(NK,A2,LD,2,A2,LD) C Form mean vector in V2. DO 30 I=1,NK AA=0.0D0 DO 20 J=I,NK DO 20 L=1,J 20 AA=AA+A2(I,J)*A2(L,J)*V1(L) 30 V2(I)=AA C Form basic shocks in V1. CALL DRNNOA(NK,V1) C PDFLG=-NK*HLG2PI-.5D0*DNRM2(NK,V2,1)**2 C Form final vector as sum of mean and deviation. DO 50 I=1,NK C PDFLG=PDFLG-DLOG(A2(I,I)) AA=0.0D0 DO 40 J=I,NK 40 AA=AA+A2(I,J)*V1(J) 50 BETA(I)=V2(I)+AA RETURN END SUBROUTINE OUT(A,N,NAME) DIMENSION A(N) CHARACTER*(*) NAME L=0 WRITE(6,5)NAME 5 FORMAT(1X,A) DO 20 K=1,N,5 L1=L+1 LAST=L+5 IF(LAST.GT.N)LAST=N WRITE(6,10)L1,(A(I),I=L1,LAST) 10 FORMAT(I4,1X,5(1PE15.8)) 20 L=L+5 RETURN END SUBROUTINE DOUT(A,N,NAME) REAL*8 A(N) CHARACTER*(*) NAME L=0 WRITE(6,5)NAME 5 FORMAT(1X,A) DO 20 K=1,N,5 L1=L+1 LAST=L+5 IF(LAST.GT.N)LAST=N WRITE(6,10)L1,(A(I),I=L1,LAST) 10 FORMAT(I4,1X,5(1PE15.8)) 20 L=L+5 RETURN END SUBROUTINE IOUT(M,N,NAME) DIMENSION M(N) CHARACTER*(*) NAME L=0 WRITE(6,5)NAME 5 FORMAT(1X,A) DO 20 K=1,N,15 L1=L+1 LAST=L+15 IF(LAST.GT.N)LAST=N WRITE(6,10)L1,(M(I),I=L1,LAST) 10 FORMAT(16I5) 20 L=L+15 RETURN END SUBROUTINE DOUT2(A,N1,N2,NAME) REAL*8 A(N1,N2) CHARACTER*(*) NAME WRITE(6,5)NAME 5 FORMAT(1X,A) DO 40 J=1,N1 WRITE(6,10)J 10 FORMAT(' Row',I3) L=0 DO 30 K=1,N2,5 L1=L+1 LAST=L+5 IF(LAST.GT.N2)LAST=N2 WRITE(6,20)L1,(A(J,I),I=L1,LAST) 20 FORMAT(I4,1X,5(1PE15.8)) 30 L=L+5 40 CONTINUE RETURN END SUBROUTINE COUT(A,N,NAME) COMPLEX*16 A(N) CHARACTER*(*) NAME L=0 WRITE(6,5)NAME 5 FORMAT(1X,A) DO 20 K=1,N,2 L1=L+1 LAST=L+2 IF(LAST.GT.N)LAST=N WRITE(6,10)L1,(A(I),I=L1,LAST) 10 FORMAT(I4,1X,2(1PE15.8),4X,2(1PE15.8)) 20 L=L+2 RETURN END SUBROUTINE I2OUT(M,N,NAME) IMPLICIT INTEGER*2 (I-N) DIMENSION M(N) CHARACTER*(*) NAME L=0 WRITE(6,5)NAME 5 FORMAT(1X,A) DO 20 K=1,N,15 L1=L+1 LAST=L+15 IF(LAST.GT.N)LAST=N WRITE(6,10)L1,(M(I),I=L1,LAST) 10 FORMAT(16I5) 20 L=L+5 RETURN END SUBROUTINE DOUT22(A,NR1,NR2,NC1,NC2,IDA,NAME) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) NAME DIMENSION A(IDA,NC2) WRITE(6,10)NAME,NC1,NC2 10 FORMAT(/,' Matrix ',A,' Columns',I3,' -',I3) DO 20 IR=NR1,NR2 DO 20 IFC=NC1,NC2,5 ILC=IFC+4 IF(ILC.GT.NC2)ILC=NC2 IF(IFC.EQ.NC1)WRITE(6,30)IR,(A(IR,J),J=IFC,ILC) 20 IF(IFC.GT.NC1)WRITE(6,40)(A(IR,J),J=IFC,ILC) 30 FORMAT(I4,1X,5(1PE15.7)) 40 FORMAT(5X,5(1PE15.7)) RETURN END SUBROUTINE LOUT(LA,N,NAME) LOGICAL LA(N) CHARACTER*(*) NAME CHARACTER*7 LSTAT(10) L=0 WRITE(6,10)NAME 10 FORMAT(1X,A) DO 30 I=1,N,10 LAST=I+9 IF(LAST.GT.N)LAST=N K=0 DO 20 J=I,LAST K=K+1 LSTAT(K)='FALSE ' 20 IF(LA(J))LSTAT(K)='TRUE ' 30 WRITE(6,40)I,(LSTAT(J),J=1,K) 40 FORMAT(I5,2X,10A7) RETURN END SUBROUTINE PDSYMR(NFILE,N,A,LDA) C This routine reads an n x n matrix A from a file, and checks to see C that it is symmetrix and positive definite. C Inputs: C NFILE File number from which to read matrix in full row and column C form, one row per record C N Order of matrix C LDA Leading dimension of matrix C Output: C A Matrix IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON /SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION A(LDA,N) C Check for symmetry DO 10 I=1,N 10 READ(NFILE,*,END=110,ERR=120)(A(I,J),J=1,N) DO 20 I=1,N DO 20 J=1,I 20 IF(A(I,J).NE.A(J,I))GO TO 130 C Check for positive definite CALL DPDCK(A,LDA,N,INFO) IF (INFO.NE.0) GO TO 140 RETURN C Termination 110 WRITE(6,115)I 115 FORMAT(/,' *** UNEXPECTED END-OF-FILE, MATRIX ROW',I3) CALL TERM 120 WRITE(6,125)I 125 FORMAT(/,' *** READ ERROR, MATRIX ROW',I3) CALL TERM 130 WRITE(6,135)J,I 135 FORMAT(/,' *** MATRIX ASYMETRY, ELEMENTS',I3,' AND',I3) CALL TERM 140 WRITE(6,145) 145 FORMAT(/,' *** MATRIX IS NOT POSITIVE DEFINITE') CALL TERM END REAL*8 FUNCTION PNORDF(AMU,H,X) C This function evaluates the c.d.f. of a univariate normal. C Inputs: C AMU Mean of normal distribution C H Precision of normal distribution C X Point of evaluation C Output: C PNORDF Value of c.d.f. IMPLICIT REAL*8 (A-H,O-Z) Z=(X-AMU)*DSQRT(H) PNORDF=DNORDF(Z) RETURN END SUBROUTINE PLINE2(NSKIP,A1,I1,A2,I2,A3) CHARACTER*(*)A1,A2,A3 IF(NSKIP.EQ.0)GO TO 30 DO 10 I=1,NSKIP 10 WRITE(6,20) 20 FORMAT(/) 30 WRITE(6,40)A1,I1,A2,I2,A3 40 FORMAT(1X,A,I4,1X,A,I4,1X,A) END C SUBROUTINE PLINE1(NSKIP,A1,I1,A2) CHARACTER*(*)A1,A2 IF(NSKIP.EQ.0)GO TO 30 DO 10 I=1,NSKIP 10 WRITE(6,20) 20 FORMAT(/) 30 WRITE(6,40)A1,I1,A2 40 FORMAT(1X,A,I4,1X,A) END C SUBROUTINE PLINE(NSKIP,A1) CHARACTER*(*)A1 IF(NSKIP.EQ.0)GO TO 30 DO 10 I=1,NSKIP 10 WRITE(6,20) 20 FORMAT(/) 30 WRITE(6,40)A1 40 FORMAT(1X,A) END C SUBROUTINE PINTI(A1,K1,K2,K3) CHARACTER*(*)A1 IF((K1.GE.K2).AND.(K1.LE.K3))RETURN WRITE(6,10)A1,K1,K2,K3 10 FORMAT(/,' *** ',A,I4,' NOT IN RANGE',I4,' to ',I5) CALL TERM END C SUBROUTINE PINTR(A1,R1,R2,R3) REAL*8 R1,R2,R3 CHARACTER*(*)A1 IF((R1.GT.R2).AND.(R1.LT.R3))RETURN WRITE(6,10)A1,R1,R2,R3 10 FORMAT(/,' *** ',A,1PE12.5, 1 1X,'NOT IN RANGE',1PE12.5,' to ',1PE12.5) CALL TERM END C SUBROUTINE PSIMP(A1,N,P) REAL*8 P(N) CHARACTER*(*)A1 AA=0.0D0 DO 10 I=1,N IF((P(I).LT.0.0D0).OR.(P(I).GT.1.0D0))GO TO 20 10 AA=AA+P(I) IF((AA.GT..99999D0).AND.(AA.LT.1.00001D0))RETURN 20 WRITE(6,30)A1,(P(I),I=1,N) 30 FORMAT(/,' ***',A,' NOT IN UNIT SIMPLEX:',/,(5(1PE15.3))) CALL TERM END C SUBROUTINE PEND(A1) CHARACTER*(*)A1 WRITE(6,10)A1 10 FORMAT(/,' *** ',A) CALL TERM END C SUBROUTINE TERM C This routine prints a terminal message and halts execution. WRITE(6,10) 10 FORMAT(' PROGRAM STOPS.') STOP END SUBROUTINE PROVAR(N,P,LDP,SIGMA,LDSIG) C This routine constructs the variance matrix, from a multivariate projection C matrix (for description of latter see VARPRO). This routine is the inverse C of VARPRO. C Inputs: C N Order of multivariate distribution (= number of variables) C P Projection matrix C LDP Leading dimension of P C LDSIG Leading dimension of SIGMA C Outputs: C SIGMA Variance matrix IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION P(LDP,N),SIGMA(LDSIG,N) CALL UM0SET(N,N,A1,LD) DO 20 I=1,N AA=1.0D0/P(I,I) A1(I,I)=AA IF(I.EQ.1)GO TO 20 DO 10 J=1,I-1 10 A1(I,J)=-AA*P(I,J) 20 CONTINUE CALL DLINRT(N,A1,LD,1,A1,LD) DO 40 I=1,N DO 40 J=1,I AA=0.0D0 DO 30 K=1,J 30 AA=AA+A1(I,K)*A1(J,K) SIGMA(I,J)=AA 40 SIGMA(J,I)=AA RETURN END SUBROUTINE VARPRO(N,SIGMA,LDSIG,P,LDP) C This routine constructs a recursive projection matrix P, given a variance C matrix. The linear projection of x(i) on x(1), ... x(i-1) is C x(i)=p(i,1)x(1)+...+p(i,i-1)x(i-1)+u(i); s.d.[u(i)]=p(i,i). C Inputs: C N Order of multivariate distribution (= number of variables) C SIGMA Variance matrix C LDSIG Leading dimension of SIGMA C LDP Leading dimension of P C Outputs: C P Projection matrix IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION SIGMA(LDSIG,N),P(LDP,N) CALL DLFTDS(N,SIGMA,LDSIG,A1,LD) CALL DLINRT(N,A1,LD,2,A1,LD) CALL UM0SET(N,N,P,LDP) DO 20 I=1,N AA=1.0D0/A1(I,I) P(I,I)=AA IF(I.EQ.1)GO TO 20 DO 10 J=1,I-1 10 P(I,J)=-AA*A1(J,I) 20 CONTINUE RETURN END SUBROUTINE PRSYM(HTITLE,N,A,LDA) C This routine prints the lower triangle of a symmetric n x n matrix. C Inputs: C HTITLE Header to identify matrix C N Order of matrix, n C A Matrix C LDA Leading dimension of matrix IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) HTITLE DIMENSION A(LDA,N) WRITE(6,5)HTITLE 5 FORMAT(/,1X,A) DO 20 IC=1,N,5 NR1=IC NR2=N NC1=IC NC2=IC+4 IF(NC2.GT.N)NC2=N WRITE(6,10)(NC,NC=NC1,NC2) 10 FORMAT(/,5I15) DO 20 IR=NR1,NR2 NCA=NC2 IF(NC2.GT.IR)NCA=IR 20 WRITE(6,30)IR,(A(IR,J),J=NC1,NCA) 30 FORMAT(I4,5(1PE15.6)) RETURN END SUBROUTINE PRCOR(HTITLE,N,A,LDA) C This routine prints the normalized lower triangle of a symmetric n x n C matrix; e.g., it converts a variance matrix to a correlation matrix before C printing. C Inputs: C HTITLE Header to identify matrix C N Order of matrix, n C A Matrix C LDA Leading dimension of matrix IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) HTITLE DIMENSION A(LDA,N) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) WRITE(6,5)HTITLE 5 FORMAT(/,1X,A) DO 20 IC=1,N,10 NR1=IC NR2=N NC1=IC NC2=IC+9 IF(NC2.GT.N)NC2=N WRITE(6,10)(NC,NC=NC1,NC2) 10 FORMAT(/I10,9I7) DO 20 IR=NR1,NR2 V1(IR)=DSQRT(A(IR,IR)) NCA=NC2 IF(NC2.GT.IR)NCA=IR 20 WRITE(6,30)IR,(A(IR,J)/(V1(IR)*V1(J)),J=NC1,NCA) 30 FORMAT(I5,10F7.3) RETURN END REAL*8 FUNCTION PTDF(AMU,H,ANU,X) C This function evaluates the c.d.f. of a Student t distribution C Inputs: C AMU Location parameter of Student t C H Inverse of squared scale parameter of Student t C ANU Degrees of freedom parameter of Student t C X Point of evaluation C Output: C PTDF Value of c.d.f. IMPLICIT REAL*8 (A-H,O-Z) T=(X-AMU)*DSQRT(H) PTDF=DTDF(T,ANU) RETURN END SUBROUTINE QUANT(N,NQ,W,G,P,Q) C From the sampled distribution of a random variable, compute quantiles. C Inputs: C N Size of sample C NQ Number of quantiles to compute C W(Nx1) Vector of log weights C G(Nx1) Sampled random variables C P(NQx1) Probabilities corresponding to quantiles C Outputs: C Q(NQx1) Quantiles C W(Nx1) Actual weights, unless all log weights were 0 C G(Nx1) Ordered random variables IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=10000) DIMENSION W(N),G(N),P(NQ),Q(NQ) DIMENSION NV(LD) IF(DASUM(N,W,1).GT.0.0D0)GO TO 20 C Procedure for weights the same -- CALL DSVRGN(N,G,G) DO 10 IQ=1,NQ JQ=P(IQ)*N+1 10 Q(IQ)=G(JQ) RETURN C Procedure for different weights -- 20 WMAX=W(IDMAX(N,W,1)) WTOT=0.0D0 DO 30 I=1,N NV(I)=I AA=0.0D0 IF(W(I).GT.-700.0D0)AA=DEXP(W(I)-WMAX) W(I)=AA 30 WTOT=WTOT+AA CALL DSVRGP(N,G,G,NV) JQ=1 WCUT=P(1)*WTOT WRUN=0.0D0 DO 50 I=1,N WRUN=WRUN+W(NV(I)) 40 IF(WRUN.LT.WCUT)GO TO 50 Q(JQ)=G(I) JQ=JQ+1 IF(JQ.GT.NQ)RETURN WCUT=P(JQ)*WTOT GO TO 40 50 CONTINUE DO 60 IQ=JQ,NQ 60 Q(IQ)=G(N) RETURN END REAL*8 FUNCTION ROZTST(NQ,A) C This routine conducts Rozanov's spectral teset for the invertibility C of a univarite lag operator. If the operator is invertible then the C result should be zero, otherwise it is positive. A good definition C of algorithmically zero with this routine is 10^(-6). This can be made C a little smaller, about 10^(-7), with stricter tolerances in the C quadrature routine, but at the expense of a 40% increase in computing time. C Inputs: C A A(j) contains the term on the j'th lag. The contemnporaneous term C in A(0) must be 1.0. C NQ Length of lag operator C Output: C ROZTST Rozanov's test criterion IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL QFUNC DIMENSION A(0:NQ) DATA PI/3.1415926535D0/ COMMON/QQ1/AA(0:200),NNQ PI=DCONST('PI') NNQ=NQ CALL DCOPY(NQ+1,A,1,AA,1) CALL DQDAG(QFUNC,0.0D0,PI,1.0D-6,1.0D-6,2,HOLD,ERR) ROZTST=HOLD/PI RETURN END REAL*8 FUNCTION QFUNC(ALAM) IMPLICIT REAL*8 (A-H,O-Z) COMMON/QQ1/A(0:200),NQ AA=0.0D0 BB=0.0D0 DO 10 IQ=0,NQ ARG=ALAM*IQ AA=AA+A(IQ)*DCOS(ARG) 10 BB=BB+A(IQ)*DSIN(ARG) QFUNC=DLOG(AA**2+BB**2) RETURN END REAL*8 FUNCTION RTMIN(N,A) C This function returns the smallest modulus of the roots of a polynomial C of order N. C Inputs: C N Order of polynomial C A Coefficients of polynomial, beginning with order 0 C Outputs: C URTMIN Smallest modulus IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/ROOT(LD),A1(LD),A2(LD) COMPLEX*16 ROOT CALL DZPORC(N,A,ROOT) AMIN=CDABS(ROOT(1)) IF(N.EQ.1)GO TO 20 DO 10 I=2,N RT=CDABS(ROOT(I)) 10 IF(RT.LT.AMIN)AMIN=RT 20 RTMIN=AMIN RETURN END SUBROUTINE RDSIM(NFILE,ITER,WPP,G) C Entry RDSIM: Read an iteration record from the previously opened C simulation input file on unit NFILE. C Inputs: C NFILE Unit number for file C Outputs: C ITER Iteration number wread from first line C WPP(3) Log weight, log prior, log data p.d.f. read from 1st line C G(NPARS) Vector of simulated values read from succeeding lines C Entry RDSIM0: Open and read the first record from the C simulation output file on unit NFILE C Inputs: C NFILE Unit number for file C LD Maximum number of parameters allowed in vector G C Outputs: C NITER Number of iterations C NPARS Number of parameters IMPLICIT REAL*8 (A-H,O-Z) COMMON/ARDSIM/LDA,NPARSA DIMENSION WPP(3),G(LDA) READ(NFILE,*,END=110,ERR=120)ITER,(WPP(J),J=1,3),(G(J),J=1,NPARSA) RETURN ENTRY RDSIM0(NFILE,LD,NITER,NPARS) LDA=LD CALL FOPEN('Posterior simulation matrix',NFILE,.FALSE.) READ(NFILE,*,END=130,ERR=140)NITER,NPARS IF((NITER.LE.0).OR.(NPARS.LE.0))GO TO 150 IF(NPARS.GT.LD)GO TO 160 NPARSA=NPARS RETURN 110 WRITE(6,115)ITER 115 FORMAT(/,' *** UNEXPECTED END OF FILE AFTER ITERATION',I16,/, 1 ' INPUT SIMULATION FILE ') CALL TERM 120 WRITE(6,125)ITER 125 FORMAT(/,' *** READ ERROR AFTER ITERATION',I16,/, 1 ' INPUT SIMULATION FILE ') CALL TERM 130 WRITE(6,135) 135 FORMAT(/,' *** UNEXPECTED END OF FILE, FIRST RECORD OF',/, 1 ' INPUT SIMULATION FILE ') CALL TERM 140 WRITE(6,145) 145 FORMAT(/,' *** READ ERROR, FIRST RECORD OF',/, 1 ' INPUT SIMULATION FILE ') CALL TERM 150 WRITE(6,155)NITER,NPARS 155 FORMAT(/,' *** ILLEGAL FIRST RECORD OF INPUT SIMULATION FILE ', 1 /,2I16) CALL TERM 160 WRITE(6,165)NPARS,LD 165 FORMAT(/,' *** INPUT SIMULATION FILE ',/,I8, 1 ' PARAMETERS EXCEEDS LIMIT OF',I4) CALL TERM END SUBROUTINE WTSIM(NFILE,ITER,WPP,G) C Entry WTSIM: Write an iteration record on the previously opened C simulation output file on unit NFILE. C Inputs: C NFILE Unit number for file C ITER Iteration number written in first line C WPP(3) Log weight, log prior, log data p.d.f. written in 1st line C G(NPARS) Vector of simulated values written in succeeding lines C Entry WTSIM0: Write the first record on the previously opened C simulation output file on unit NFILE C Inputs: C NFILE Unit number for file C NITER Number of iterations C NPARS Number of parameters IMPLICIT REAL*8 (A-H,O-Z) COMMON/ATWSIM/NPARSA DIMENSION WPP(3),G(NPARSA) WRITE(NFILE,10)ITER,(WPP(J),J=1,3),(G(J),J=1,NPARSA) 10 FORMAT(I16,3(1PE16.7),/,(5(1PE16.7))) RETURN ENTRY WTSIM0(NFILE,NITER,NPARS) NPARSA=NPARS WRITE(NFILE,20)NITER,NPARS 20 FORMAT(2I16) RETURN END SUBROUTINE STT(K,AMU,H,LDH,ANU,X1,X2) C This routine generates two antithetic pseudo-random vectors x1 and x2 C from a Student-t distribution with mean vector mu and precision C matrix H. C Inputs: C K Dimension of x C AMU Mean vector mu (kx1) C H Precision matrix H (kxk) C LDH Leading dimension of H C ANU Degrees of freedom parameter C Outputs: C X1,X2 Random draws IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION AMU(K),H(LDH,K),X1(K),X2(K) C Entry STT: Factor precision CALL DLFTDS(K,H,LDH,A1,LD) CALL DLINRT(K,A1,LD,2,A1,LD) C Entry STTA: Precision already factored; C *** Can be no intervening calls that change A1 !! *** ENTRY STTA(K,AMU,H,LDH,ANU,X1,X2) 10 CALL URNMVN(K,A1,LD,V1) CALL DRNCHI(1,ANU,AA) AA=1.0D0/DSQRT(AA/ANU) CALL DSCAL(K,AA,V1,1) CALL UVADD(K,AMU,V1,X1) CALL UVSUB(K,AMU,V1,X2) RETURN C Entry GAUB: Coming in with R: RR'=Variance, R upper triangular ENTRY STTB(K,AMU,R,LDR,ANU,X1,X2) CALL UCRGRG(K,K,R,LDR,A1,LD) GO TO 10 END REAL*8 FUNCTION STTK(K,AMU,H,LDH,ANU,X) C This function evaluates the log density kernel of a multivariate C normal random vector x with mean mu and precision matrix H. The C kernel density is [1+.5(x-mu)'H(x-mu)/nu]^[-(nu+k)/2] C Inputs: C K Dimension of x C AMU Mean vector mu (kx1) C H Precision matrix h (kxk) C LDH Leading dimension of H C ANU Degrees of freedom parameter C X Point of evaluation, x (kx1) C Output: C STTK Evaluation of log kernel density IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION H(LDH,K),X(K) CALL UVSUB(K,X,AMU,V1) STTK=-.5D0*(ANU+K)*DLOG(1.0D0+UQUAF(K,V1,H,LDH)/ANU) RETURN END REAL*8 FUNCTION STTN(K,H,LDH,ANU) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This routine evaluates the log normalizing constant for the log C density kernel of a random vector computed in GAUK. This C normalizing constant, when multiplied by the kernel in GAUK, C yields the p.d.f. of the random vector. C Inputs: C K Dimension of random vector C H Precision matrix H (kxk) C LDH Leading dimension of H C ANU Degrees of freedom parameter C Output: C STTN Log normalizing factor IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION H(LDH,K) STTN=DLNGAM(.5D0*(ANU+K))-DLNGAM(.5D0*ANU) 1 +.5D0*(DTLGPD(K,H,LDH)-K*DLOG(PI*ANU)) RETURN END IMPLICIT REAL*8 (A-H,O-Z) LOGICAL DIFNAN DIMENSION A(10) DO 10 I=1,8 10 A(I)=DMACH(I) CALL DOUT(A,8,'A') READ(5,*)B A(1)=1.0D0/B A(2)=B**2 A(3)=DEXP(B) A(4)=DLOG(B) CALL DOUT(A,4,'A') DO 20 I=1,4 20 IF(DIFNAN(A(I)))WRITE(6,30)I 30 FORMAT(' Entry',I2,' detected NaN by IMSL') END REAL*8 FUNCTION TRAB(NRA,NCA,A,LDA,NRB,NCB,B,LDB) C This routine computes tr(AB), where A:kxm, B:mxk. C Inputs C NRA Rows in A C NCA Columns in A C A Matrix A C LDA Leading dimension of A C NRB Rows in B C NCB Columns in B C B Matrix B C LDB Leading dimension of B IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,NCA),B(LDB,NCB) IF((NCA.NE.NRB).OR.(NRA.NE.NCB))GO TO 110 AA=0.0D0 DO 10 I=1,NRA DO 10 J=1,NCA 10 AA=AA+A(I,J)*B(J,I) TRAB=AA RETURN 110 WRITE(6,115)NRA,NCA,NRB,NCB 115 FORMAT(' Order inconsistency in TRAB:',4I6) STOP END SUBROUTINE TRMARK(LFWD,N,VP,VPTR,TRJAC) C In the forward transformation -- C This routine performs a normalizing transformation for a Markov transition C matrix. The tranformation is from pij to qij=log(pij/pii) for all j.ne.i. C It also returns the log Jacobian of transformation. This is what you need C to add to the log density in the new parameterization. C Inputs: C LFWD Set .TRUE. C N Order of Markov matrix C VP Markov matrix, stored in vectorized form (by columns) C Outputs: C VPTR Transformed parameters, organized by **rows** C TRJAC Log Jacobian of forward transformation C C In the backward transformation -- C Inverse transformation is performed; TRJAC is negative of forward value C Inputs: C LFWD Set .FALSE. C N Order of Markov matrix C VPTR Transformed parameters (qij), stored in vectorized form (by rows) C Outputs: C VP Markov matrix, stored in vectorized form (by columns) C TRJAC Log Jacobian of backward transformation IMPLICIT REAL*8 (A-H,O-Z) LOGICAL LFWD PARAMETER (LD=100) COMMON/SCRATCH/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION VP(N*N),VPTR(N*(N-1)) IF(.NOT.LFWD)GO TO 100 CALL UNVEC(N,N,VP,A1,LD) AA=0.0D0 L=0 DO 20 I=1,N BB=DLOG(A1(I,I)) DO 10 J=1,N IF(J.EQ.I)GO TO 10 L=L+1 CC=DLOG(A1(I,J)) VPTR(L)=CC-BB AA=AA+CC 10 CONTINUE 20 AA=AA+BB TRJAC=AA RETURN 100 L=0 AA=0.0D0 DO 120 I=1,N BB=0.0D0 DO 110 J=1,N A1(I,J)=1.0D0 IF(I.EQ.J)GO TO 110 L=L+1 AA=AA-VPTR(L) A1(I,J)=DEXP(VPTR(L)) 110 BB=BB+A1(I,J) AA=AA+N*DLOG(BB) DO 120 J=1,N 120 A1(I,J)=A1(I,J)/BB CALL VEC(N,N,A1,LD,VP) TRJAC=AA RETURN END SUBROUTINE TRPD(LFWD,N,VA,VATR,TRJAC) C In the forward transformation -- C This routine performs a normalizing transformation on a symmetric positive C definite matrix, and returns the log Jacobian of transformation. C The transformation is based on the Choleski factorization of p.d. A, C A = LL', where L is lower triangular. The transformation is to the logs C of the diagonal elemnts of L, and the untransformed elements in the strict C lower triangle of A. The log Jacobian is what you add to the log density C of the original distribution, to get the right log density for the C transformed variables. C Inputs: C LFWD Set .TRUE. C N Order of p.d. matrix C VA p.d. matrix in lower triangular storage, by columns C Outputs: C VATR Transformed Choleski decomposition, same storage convention as VA C TRJAC Log Jacobian of forward transformation C C In the backward transformation -- C Inverse transformation is performed; TRJAC is negative of forward value C Inputs: C LFWD Set .FALSE. C N Order of p.d. matrix C VATR Transformed Choleski decomposition of matrix, lower tri. by columns C Outputs: C VA p.d. matrix in lower triangular storage, by columns C TRJAC Log Jacobian of backward transformation IMPLICIT REAL*8 (A-H,O-Z) LOGICAL LFWD DIMENSION VA(N*(N+1)/2),VATR(N*(N+1)/2) DATA DLOG2/.693147181D0/ PARAMETER(LD=100) COMMON/SCRATCH/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) IF(.NOT.LFWD)GO TO 100 L=1 DO 10 J=1,N DO 10 I=J,N A1(J,I)=VA(L) 10 L=L+1 CALL DLFTDS(N,A1,LD,A2,LD) AA=N*DLOG2 K=N+1 L=0 IF(N.EQ.1)GO TO 40 DO 30 J=1,N-1 L=L+1 BB=DLOG(A2(J,J)) VATR(L)=BB AA=AA+K*BB DO 20 I=J+1,N L=L+1 20 VATR(L)=A2(J,I) 30 K=K-1 40 BB=DLOG(A2(N,N)) VATR(L+1)=BB TRJAC=AA+K*BB RETURN 100 L=0 K=N+1 AA=-N*DLOG2 DO 120 J=1,N DO 110 I=J,N L=L+1 110 A1(I,J)=VATR(L) AA=AA-K*A1(J,J) K=K-1 120 A1(J,J)=DEXP(A1(J,J)) TRJAC=AA L=0 DO 140 J=1,N DO 140 I=J,N L=L+1 AA=0.0D0 DO 130 M=1,J 130 AA=AA+A1(I,M)*A1(J,M) 140 VA(L)=AA RETURN END SUBROUTINE UCRGRG(M,N,A,LDA,B,LDB) C Copy A(MxN) to B(MxN). IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),B(LDB,N) DO 10 IM=1,M DO 10 JN=1,N 10 B(IM,JN)=A(IM,JN) RETURN END SUBROUTINE UCROSS(NRA,NCA,A,LDA,NRB,NCB,B,LDB,NRC,NCC,C,LDC) C This routine computes the Kronecker product of A and B and places the result C in C. C Inputs: C NRA Number of rows in A C NCA Number of columns in A C A Matrix A C LDA Leading dimension of A C NRB Number of rows in B C NCB Number of columns in B C B Matrix B C LDB Leading dimension of B C NRC Number of rows in C (Must be NRA*NRB) C NCC Number of columns in C (Must be NCA*NCB) C LDC Leading dimension of C C Output: C C Matrix C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,NCA),B(LDB,NCB),C(LDC,NCC) IF((NRA*NRB).NE.NRC)GO TO 20 IF((NCA*NCB).NE.NCC)GO TO 30 IRC=0 DO 10 IRA=1,NRA DO 10 IRB=1,NRB IRC=IRC+1 ICC=0 DO 10 ICA=1,NCA DO 10 ICB=1,NCB ICC=ICC+1 10 C(IRC,ICC)=A(IRA,ICA)*B(IRB,ICB) RETURN 20 WRITE(6,25)NRA,NRB,NRC 25 FORMAT 1 (/,' Order error in UCROSS. NRA =',I3,' NRB =',I3,'NRC =',I4) STOP 30 WRITE(6,35)NCA,NCB,NCC 35 FORMAT 1 (/,' Order error in UCROSS. NCA =',I3,' NCB =',I3,'NCC =',I4) STOP END SUBROUTINE ULFTDS(N,A,LDA,R,LDR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),R(LDR,N) CALL DLFTDS(N,A,LDA,R,LDR) DO 10 I=2,N DO 10 J=1,I-1 10 R(I,J)=0.0D0 RETURN END SUBROUTINE UMADD(M,N,A,LDA,B,LDB,C,LDC) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),B(LDB,N),C(LDC,N) DO 10 I=1,M DO 10 J=1,N 10 C(I,J)=A(I,J)+B(I,J) RETURN END SUBROUTINE UMSUB(M,N,A,LDA,B,LDB,C,LDC) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),B(LDB,N),C(LDC,N) DO 20 I=1,M DO 20 J=1,N 20 C(I,J)=A(I,J)-B(I,J) RETURN END SUBROUTINE UMSADD(M,N,A,LDA,B,LDB,C,LDC,S) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),B(LDB,N),C(LDC,N) DO 10 I=1,M DO 10 J=1,N 10 C(I,J)=S*A(I,J)+B(I,J) RETURN END SUBROUTINE UMISET(N,A,LDA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N) DO 20 I=1,N DO 10 J=1,N 10 A(I,J)=0.0D0 20 A(I,I)=1.0D0 RETURN END SUBROUTINE UM0SET(M,N,A,LDA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N) DO 10 I=1,M DO 10 J=1,N 10 A(I,J)=0.0D0 RETURN END SUBROUTINE UNVEC(M,N,B,A,LDA) C This routine unvectorizes a vector B into an M x N matrix A. C The vector B is assumed to have the columns of A stored successively. C This routine is the inverse of the routine VEC. C Inputs: C M Row dimension of matrix C N Column dimension of matrix C B Vectorized matrix; vector is of length MN. C LDA Leading dimension of A in calling program C Outputs: C A Matrix that was vectorized in B. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),B(500) L=0 DO 10 JN=1,N DO 10 IM=1,M L=L+1 10 A(IM,JN)=B(L) RETURN END REAL*8 FUNCTION UQUAF(N,X,A,LDA) C This routine forms the quadratic form Q = x'Ax, A being symmetric. C Inputs: C N Dimension of x and order of A C X x vector C A A matrix C LDA Leading dimension of A in calling program C Outputs: C UQUAF x'Ax IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N),A(LDA,N) AA=0.0D0 DO 10 I=1,N DO 10 J=1,I BB=X(I)*A(J,I)*X(J) AA=AA+BB 10 IF(J.LT.I)AA=AA+BB UQUAF=AA RETURN END SUBROUTINE URNMVN(K,U,LDU,X) C This routine generates one pseudo-random vector from a multivariate C normal distribution with mean 0 and variance UU'. C References: This routine patches up some anomolies in the IMSL library. C Update history: Brought into ulib January 21, 1992. C Inputs: C K Length of normal random vector C U Upper triangular K x K matrix U: UU' = Variance of normal C LDU Leading dimension of U in calling routine C Output: C X Random draw from N(0, UU') IMPLICIT REAL*8 (A-H,O-Z) DIMENSION U(LDU,K),X(K) CALL DSET(K,0.0D0,X,1) DO 10 J=1,K EPS=DRNNOF() DO 10 I=1,J 10 X(I)=X(I)+U(I,J)*EPS RETURN END REAL*8 FUNCTION URTMIN(N,A) C This function returns the smallest modulus of the roots of a polynomial C of order N. C Inputs: C N Order of polynomial C A Coefficients of polynomial, beginning with order 0 C Outputs: C URTMIN Smallest modulus IMPLICIT REAL*8 (A-H,O-Z) COMPLEX*16 ROOT(100) CALL DZPORC(N,A,ROOT) RTMIN=CDABS(ROOT(1)) IF(N.EQ.1)GO TO 20 DO 10 I=2,N RT=CDABS(ROOT(I)) 10 IF(RT.LT.RTMIN)RTMIN=RT 20 URTMIN=RTMIN RETURN END SUBROUTINE USCAL2(NA,NB,X,A,LDA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,NB) XX=X DO 10 IA=1,NA DO 10 IB=1,NB 10 A(IA,IB)=XX*A(IA,IB) RETURN END SUBROUTINE UVADD(N,A,B,C) C This routine adds two vectors, C = A + B. C Inputs C N Order of the vectors C A A vector C B B vector C Outputs: C C C vector, C = A + B IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(N),B(N),C(N) DO 10 I=1,N 10 C(I)=A(I)+B(I) RETURN END SUBROUTINE UVSUB(N,A,B,C) C This routine subtracts one vector from another, C = A - B. C Inputs C N Order of the vectors C A A vector C B B vector C Outputs: C C C vector, C = A - B IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(N),B(N),C(N) DO 10 I=1,N 10 C(I)=A(I)-B(I) RETURN END SUBROUTINE UVSADD(N,A,S,B,C) C This routine adds a scalar times a vector and a vector, C = S*A + B. C Inputs C N Order of the vectors C A A vector C B B vector C Outputs: C C C vector, C = S*A + B IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(N),B(N),C(N) DO 10 I=1,N 10 C(I)=S*A(I)+B(I) RETURN END SUBROUTINE UXAXTF(N,M,X,LDX,A,LDA,B,LDB) C This routine forms the matrix product B = XAX', A being symmetric. C Inputs: C N Column dimension of X and order of A C M Row dimension of X and order of B C X X matrix C LDX Leading dimension of X in calling program C A A matrix (full storeage) C LDA Leading dimension of A in calling program C LDB Leading dimension of B in calling program C Outputs: C B B matrix, B = X'AX IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(LDX,N),A(LDA,N),B(LDB,M) DO 20 I=1,M DO 20 J=1,I AA=0.0D0 DO 10 K=1,N DO 10 L=1,N 10 AA=AA+X(I,K)*A(K,L)*X(J,L) B(I,J)=AA 20 B(J,I)=AA RETURN END SUBROUTINE UXTAXF(N,M,X,LDX,A,LDA,B,LDB) C This routine forms the matrix product B = X'AX, A being symmetric. C Inputs: C N Row dimension of X and order of A C M Column dimension of X and order of B C X X matrix C LDX Leading dimension of X in calling program C A A matrix (full storeage) C LDA Leading dimension of A in calling program C LDB Leading dimension of B in calling program C Outputs: C B B matrix, B = X'AX IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(LDX,M),A(LDA,N),B(LDB,M) DO 20 I=1,M DO 20 J=1,I AA=0.0D0 DO 10 K=1,N DO 10 L=1,N 10 AA=AA+X(K,I)*A(K,L)*X(L,J) B(I,J)=AA 20 B(J,I)=AA RETURN END SUBROUTINE VEC(M,N,A,LDA,B) C This routine vectorizes a matrix, B = vec(A): B is a vector C containing the columns of A stored successively with no gaps C between columns. C Inputs: C M Row dimension of matrix C N Column dimension of matrix C A Matrix to be vectorized C LDA Leading dimension of A in calling program C Outputs: C B Vectorized matrix; vector is of length MN. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,100),B(500) IF((M.EQ.0).OR.(N.EQ.0))RETURN L=0 DO 10 JN=1,N DO 10 IM=1,M L=L+1 10 B(L)=A(IM,JN) RETURN END SUBROUTINE VECLT(N,A,LDA,B) C This routine vectorizes the lower triangle of a square matrix of order N. C The first N elements of the vector contain the first column of A, C the next N-1 elements of the vector contain the last N-1 elements of the C second column of A, ... , the next N-j elements of the vector contain C the last N-j elements of the (j+1)'th column of A, ..., and the last C element of the vector contains the element in the last row and column of A. C Inputs: C N Order of input matrix C A Input matrix whose lower triangle is to be vectorized C LDA Leading dimension of A in calling program C Outputs: C B Vectorized lower triangle of A IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LDA,N),B(100) L=0 DO 10 JN=1,N DO 10 IN=JN,N L=L+1 10 B(L)=A(IN,JN) RETURN END C SUBROUTINE WISH(M,SIG,LDSIG,ANU,A,LDA) C This routine generates a pseudo-random, mxm matrix A from a Wishart C distribution with matrix parameter Sigma and degrees of freedom C parameter Nu. C Inputs: C M Order of matrix, m C SIG Matrix parameter, Sigma C LDSIG Leading dimension of SIG C ANU Degrees of freedom parameter, Nu C LDA Leading dimension of A C Output: C A Random matrix A IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION SIG(LDSIG,M),A(LDA,M),R(LDR,M) CALL ULFTDS(M,SIG,LDSIG,A1,LD) C Entry WISHA: Matrix parameter already factored C *** Can be no interveneing calls that change A1 !! *** ENTRY WISHA(M,SIG,LDSIG,ANU,A,LDA) 5 DF=ANU DO 20 I=1,M DO 10 J=1,I A(J,I)=0.0D0 10 A(I,J)=DRNNOF() CALL DRNCHI(1,DF,AA) A(I,I)=DSQRT(AA) 20 DF=DF-1.0D0 CALL DMXTYF(M,M,A1,LD,M,M,A,LDA,M,M,A2,LD) CALL DMXYTF(M,M,A2,LD,M,M,A2,LD,M,M,A,LDA) RETURN C Entry WISHB: Coming in with matrix parameter already factored C R: R'R=Matrix parameter, R upper triangular with zeros set below ENTRY WISHB(M,R,LDR,ANU,A,LDA) CALL UCRGRG(M,M,R,LDR,A1,LD) GO TO 5 END REAL*8 FUNCTION WISHK(M,P,LDP,ANU,A,LDA) C This function evaluates the log density kernel of a Wishart random C matrix A with matrix parameter Sigma and degrees of freedom parameter C Nu. The log density kernel is .5(Nu+1-m)logdet(A)-.5tr[(Sigma^-1)A]. C Inputs: C M Order of matrix, m C P Square matrix P: PP'=Sigma^-1 C LDP Leading dimension of P C ANU Degrees of freedom parameter, Nu C A Matrix point of evaluation of log density kernel C LDA Leading dimension of A C Output: C WISHK Evaluation of log kernel density IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) CALL DLFTDS(M,A,LDA,A1,LD) CALL UXTAXF(M,M,P,LDP,A,LDA,A2,LD) AA=0.0D0 BB=0.0D0 DO 10 I=1,M AA=AA+DLOG(A1(I,I)) 10 BB=BB+A2(I,I) WISHK=(ANU-1-M)*AA-0.5D0*BB RETURN END REAL*8 FUNCTION WISHN(M,SIG,LDSIG,ANU,P,LDP) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This routine evaluates the log normalizing constant for the log C density kernel of a random matrix computed in WISHK. This C normalizing constant, when multiplied by the kernel in WISHK, C yields the p.d.f. of the random matrix. C Inputs: C M Order of matrix, m C SIG Matrix parameter, Sigma C LDSIG Leading dimension of SIG C ANU Degrees of freedom parameter, Nu C LDP Leading dimension of P C Outputs: C WISHN Log normalizing factor C P Upper triangular mxm matrix: PP'=Sigma^-1 IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) CALL DLFTDS(M,SIG,LDSIG,A1,LD) AA=0.0D0 BB=0.0D0 ANU1=ANU+1.0D0 DO 10 I=1,M AA=AA+DLOG(A1(I,I)) 10 BB=BB+DLNGAM(0.5D0*(ANU1-I)) WISHN=-0.5D0*M*(ANU*DLOG2+0.5D0*(M-1)*DLGPI)-ANU*AA-BB CALL DLINRT(M,A1,LD,2,P,LDP) RETURN END SUBROUTINE WISH0(NFILE,M,TITLE,SIG0,LDSIG,ANU0,ANORM,P,LDP,A,LDA) C This routine reads in a matrix parameter S0 and degress of freedom C parameter ANU0, for the distribution H ~ W[S0^-1, ANU0]. C (For an unambiguous definition of the Wishart distribution, see WISHK.) C This routine reads in the matrix and degrees-of-freedom parameter of C an m x m matrix with a Wishart distribution. In computes and returns C the normalizing constant of the density, a related matrix for use in C drawing from the distribution, and one draw from the distribution. C Inputs: C NFILE File number for reading input information C M Order of matrix, m C LDSIG Leading dimension of SIG0 C LDP Leading dimension of P C LDA Leading dimension of A C Outputs: C SIG0 Matrix parameter, Sigma C ANU Degrees of freedom parameter, Nu C ANORM Log normalizing constant for distribution C P Upper triangr triangular mxm matrix: PP'=S0 C A Drawing from distribution IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) TITLE PARAMETER(LD=100) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION A3(LD,LD) DM=M READ(NFILE,*,END=110,ERR=120)ANU0 WRITE(6,20)TITLE,ANU0 20 FORMAT(/,' Wishart prior distribution for ',A,/,' with',1PE16.7, 1 ' degrees of freedom') IF(ANU0.LE.DM)GO TO 130 CALL PDSYMR(NFILE,M,SIG0,LDSIG) CALL DLINDS(M,SIG0,LDSIG,A3,LD) ANORM=WISHN(M,A3,LD,ANU0,P,LDP) CALL WISH(M,A3,LD,ANU0,A,LDA) CALL PRSYM('Matrix parameter S0: H~W[S0^-1, ANU0]',M,SIG0,LDSIG) CALL PRSYM('Inverse of matrix parameter S0',M,A3,LD) CALL PRSYM('Initial draw of H',M,A,LDA) CALL DLINDS(M,A,LDA,A1,LD) CALL PRSYM('Inverse of initial draw of H',M,A1,LD) RETURN 110 CALL PEND('UNANTICIPATED END-OF-FILE READING WISHART PRIOR') 120 CALL PEND('ERROR READING WISHART PRIOR DEGREES OF FREEDOM') 130 CALL PINTR('WISHART DEGREES OF FREEDOM',ANU0,DM,1.D1000) END REAL*8 FUNCTION XCHIT(S2,ANU,A,B,LB) C This function generates from a doubly truncated chi-square distribution, C (s2)h ~ chisq(nu), a 2. C Inputs: C S2 Distribution parameter s2 C ANU Degrees of freedom C A Lower bound C B Upper bound if LB=.FALSE. C LB If .TRUE., no upper bound C Output: C XCHIT Random variable h IMPLICIT REAL*8 (A-H,O-Z) LOGICAL DIFNAN COMMON/CHIT/C1 EXTERNAL CHILD,CHILD1,CHILD2 LOGICAL LB IF(ANU.LE.2.0D0)CALL PEND(' DF <=2 in xchit') A1=A*S2 IF(.NOT.LB)B1=B*S2 IF(LB)B1=DMAX1(ANU,A1)+DSQRT(2.0D0*ANU)*(6.0D0+100.0D0/ANU) C1=.5D0*(ANU-2.0D0) ITRY=1 10 XCHIT=YLC(CHILD,CHILD1,CHILD2,A1,B1,.TRUE.,ANU-2.0D0)/S2 IF(XCHIT.EQ.DMACH(7))WRITE(6,15)S2,ANU,A,B 15 FORMAT(' +Infinity in XCHIT:',4(1PE15.5)) IF(.NOT.DIFNAN(XCHIT))RETURN NB=0 IF(LB)NB=1 IF(ITRY.EQ.1)WRITE(6,20)S2,ANU,A,B,NB 20 FORMAT(' NaN in XCHIT:',4(1PE14.5),I2) ITRY=ITRY+1 IF(ITRY.LE.2)GO TO 10 CALL PEND(' 2 tries for XCHIT') RETURN END REAL*8 FUNCTION CHILD(X) IMPLICIT REAL*8 (A-H,O-Z) COMMON/CHIT/C1 CHILD=C1*DLOG(X)-.5D0*X RETURN END REAL*8 FUNCTION CHILD1(X) IMPLICIT REAL*8 (A-H,O-Z) COMMON/CHIT/C1 CHILD1=C1/X-.5D0 RETURN END REAL*8 FUNCTION CHILD2(X) IMPLICIT REAL*8 (A-H,O-Z) COMMON/CHIT/C1 CHILD2=-C1/X**2 RETURN END REAL*8 FUNCTION XNOT(AMU,H,A,B,LA,LB) C This function generates a univariate normal random variable subject C to the constraint that it be in an interval (a,b), where the C endpoints may be finite or infinite. C References: J. Geweke, 1991. "Efficient Simulation from the C Multivariate Normal and Student-t Distributions Subject C to Linear Constraints," in E.M. Keramidas (ed.), C Computing Science and Statistics: Proceedings of the C 23rd Symposium on the Interface, pp. 571-578. Fairfax, C VA: Interface Foundation of North America, Inc. C Inputs: C AMU Mean of parent univariate normal distribution C H Precision of parent univariate normal distribution C A, B Endpoints of interval; A < B if LA = LB = .FALSE. C LA .TRUE. if left endpoint is - infinity; in this C case A is ignored. C LB .TRUE. if right endpoint is + infinity; in this C case B is ignored. C Output: C XNOT Random variable IMPLICIT REAL*8 (A-H,O-Z) LOGICAL LA,LB,LFLIP DATA EPS,T1,T2,T3,T4/2.0D0,.375D0,2.18D0,.725D0,.45D0/ F(X)=DEXP(-.5D0*X**2) H1=DSQRT(H) IF(LA.AND.LB)GO TO 160 LFLIP=.FALSE. IF(LA.OR.LB)GO TO 100 IF(B.LE.A)GO TO 170 C ******* Finite interval C1=H1*(A-AMU) C2=H1*(B-AMU) IF((C1*C2).GT.0.0D0)GO TO 30 C ++++ (A,B) includes 0 IF((C1.GT.-T1).AND.(C2.LT.T1))GO TO 20 C -- F(A) or F(B) small: full normal with rejection 10 X=DRNNOF() IF(X.LT.C1)GO TO 10 IF(X.GT.C2)GO TO 10 GO TO 150 C -- F(A) and F(B) large: uniform importance sampling 20 CDEL=C2-C1 25 X=C1+CDEL*DRNUNF() IF(DRNUNF().GT.F(X))GO TO 25 GO TO 150 C ++++ (A,B) excludes 0 C -- Transform to both positive 30 IF(C1.GT.0.0D0)GO TO 40 C=C1 C1=-C2 C2=-C LFLIP=.TRUE. 40 F1=F(C1) F2=F(C2) IF(F2.LT.EPS)GO TO 60 IF((F1/F2).GT.T2)GO TO 60 C -- F(A)/F(B) not large: uniform importance sampling 50 CDEL=C2-C1 55 X=C1+CDEL*DRNUNF() IF(DRNUNF().GT.(F(X)/F1))GO TO 55 GO TO 140 60 IF(C1.GT.T3)GO TO 80 C -- P(X>A) and F(A)/F(B) large: half-normal with rejection 70 X=DABS(DRNNOF()) IF(X.LT.C1)GO TO 70 IF(X.GT.C2)GO TO 70 GO TO 140 C -- P(X>A) small, F(A)/F(B) large: exponential importance C sampling with rejection 80 C=C2-C1 90 Z=-DLOG(DRNUNF())/C1 IF(Z.GT.C)GO TO 90 IF(DRNUNF().GT.F(Z))GO TO 90 X=C1+Z GO TO 140 C ****** Half-line interval 100 C1=H1*(A-AMU) C -- Transform to bound from below if A = -infinity IF(LB)GO TO 110 C1=-H1*(B-AMU) LFLIP=.TRUE. 110 IF(C1.GT.T4)GO TO 130 C -- A not large: full normal with rejection 120 X=DRNNOF() IF(X.LT.C1)GO TO 120 GO TO 140 C -- A small: exponential importance sampling 130 Z=-DLOG(DRNUNF())/C1 IF(DRNUNF().GT.F(Z))GO TO 130 X=C1+Z 140 IF(LFLIP)X=-X 150 XNOT=AMU+X/H1 RETURN C ****** Whole interval 160 XNOT=AMU+DRNNOF()/H1 RETURN C ***** Singleton 170 XNOT=A RETURN END REAL*8 FUNCTION XNOTK(AMU,H,A,B,LA,LB,X) C This function returns the logarithm of the kernel of a univariate C normal distribution truncated to the interval (a,b). If x is in the C interval the kernel is -h*(x-mu)^2; otherwise it is zero. C Inputs: C AMU Mean mu of parent univariate normal distribution C H Precision h of parent univariate normal distribution C A, B Endpoints of interval; A < B if LA = LB = .FALSE. C LA .TRUE. if left endpoint is - infinity; in this C case A is ignored. C LB .TRUE. if right endpoint is + infinity; in this C case B is ignored. C X Point of evaluation of kernel C Output: C XNOTK Log kernel IMPLICIT REAL*8 (A-H,O-Z) LOGICAL LA,LB AA=-700.0D0 IF((.NOT.LA).AND.(X.LT.A))GO TO 10 IF((.NOT.LB).AND.(X.GT.B))GO TO 10 AA=-.5D0*H*(X-AMU)**2 10 XNOTK=AA RETURN END REAL*8 FUNCTION XNOTN(AMU,H,A,B,LA,LB) C*********************************************************************** C In the calling routine a call to INIT must be made prior to calling * C this routine for the first time. * C*********************************************************************** C This function returns the logarithm of the normalizing factor of a C univariate normal distribution truncated to the interval (a,b). The C normalizing factor corresponds to the log-kernel evaluation in C XNOTK. C Inputs: C AMU Mean mu of parent univariate normal distribution C H Precision h of parent univariate normal distribution C A, B Endpoints of interval; A < B if LA = LB = .FALSE. C LA .TRUE. if left endpoint is - infinity; in this C case A is ignored. C LB .TRUE. if right endpoint is + infinity; in this C case B is ignored. C Output: C XNOTN Log normalizing factor IMPLICIT REAL*8 (A-H,O-Z) LOGICAL LA,LB COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII H1=DSQRT(H) AA=0.0D0 BB=1.0D0 FAC=-HLG2PI+.5D0*DLOG(H) IF(.NOT.LA)AA=DNORDF(H1*(A-AMU)) IF(.NOT.LB)BB=DNORDF(H1*(B-AMU)) XNOTN=FAC-DLOG(BB-AA) RETURN END REAL*8 FUNCTION YLC(G,G1,G2,D1,D2,LXHAT,XHAT) C This function draws from a distribution whose log density is strictly C concave. Documentation is in notes on the ARMA(p,q) paper. C Inputs: C G Function G(X) that has input ordinate X and returns log of density C kernel, G C G1 Function G1(X) that has input ordinate X and returns first C derivative of log of density kernel, G1 C G2 Function G2(X) that has input ordinate X and returns second C derivative of log of density kernel, G2 C D1 Lower end of finite support C D2 Upper end of finite support C LXHAT .TRUE. if mode of log density kernel provided on input C Input/output: C XHAT Mode of log density kernel C Output: C YLC Draw from distribution IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LROW=500) DIMENSION T(LROW,5),V1(2),V2(2) C DIMENSION TEMP(5) EXTERNAL G,G1,G2 LOGICAL LXHAT C1=D1 C2=D2 IF(C2.LT.C1)GO TO 1050 IF(C2.GT.C1)GO TO 10 YLC=C1 RETURN 10 CALL DSET(5*LROW,0.0D0,T,1) C Step 1 G1C1=G1(C1) IF(G1C1.GE.0.0D0)GO TO 110 X1=UMIN(C1-1.0D0/G1C1,C1+.5D0*(C2-C1)) X2=UMIN(X1-1.0D0/G1(X1),C1+.9D0*(C2-C1)) GO TO 150 110 G1C2=G1(C2) IF(G1C2.LE.0.0D0)GO TO 120 X2=UMAX(C2-1.0D0/G1C2,C1+.5D0*(C2-C1)) X1=UMAX(X2-1.0D0/G1(X2),C1+.1D0*(C2-C1)) GO TO 150 120 IF(LXHAT)GO TO 130 MAXFN=1000 AA=C1 XHAT=C2 CALL DZBREN(G1,0.1D0,0.0D0,AA,XHAT,MAXFN) 130 S=1.0D0/DSQRT(-G2(XHAT)) X1=DMAX1(XHAT-S,C1+.5D0*(XHAT-C1),C1+.1D0*(C2-C1)) X2=DMIN1(XHAT+S,XHAT+.5D0*(C2-XHAT),C2+.1D0*(C1-C2)) C Robustify to scaling problems so that subsequent exponentiation will not C overflow or underflow. 150 XBAR=0.5D0*(X1+X2) SCLOG=-G(XBAR) SC=DSQRT(-G2(XBAR)) X1=SC*X1 X2=SC*X2 C1=SC*C1 C2=SC*C2 SC1=1.0D0/SC C Step 2 200 T(1,1)=X1 T(2,1)=X2 T(1,3)=G(SC1*X1)+SCLOG T(2,3)=G(SC1*X2)+SCLOG T(1,4)=G1(SC1*X1)*SC1 T(2,4)=G1(SC1*X2)*SC1 T(1,2)=C1 T(2,2)=(T(2,3)-T(1,3)+T(1,4)*T(1,1)-T(2,4)*T(2,1)) 1 /(T(1,4)-T(2,4)) T(3,2)=C2 T(1,5)=(DEXP(T(1,3)+T(1,4)*(T(2,2)-T(1,1))) 1 -DEXP(T(1,3)+T(1,4)*(T(1,2)-T(1,1))))/T(1,4) T(2,5)=T(1,5)+(DEXP(T(2,3)+T(2,4)*(T(3,2)-T(2,1))) 1 -DEXP(T(2,3)+T(2,4)*(T(2,2)-T(2,1))))/T(2,4) C Step 3 300 JCAP=2 C Step 4 400 TSTAR=T(JCAP,5)*DRNUNF() C Step 5 500 DO 510 I=1,JCAP IF(TSTAR.GT.T(I,5))GO TO 510 J=I GO TO 600 510 CONTINUE GO TO 1010 C Step 6 600 ZJ=0.0D0 IF(J.GT.1)ZJ=T(J-1,5) C V1(1)=T(J,4)*(TSTAR-ZJ) C V1(2)=1.0D0 C V2(1)=-T(J,3) C V2(2)=T(J,4)*(T(J,2)-T(J,1)) C X=T(J,1)+DLGLCE(2,V1,V2)/T(J,4) X=T(J,1) + (DLOG(T(J,4)*(TSTAR-ZJ)*DEXP(-T(J,3)) 1 +DEXP(T(J,4)*(T(J,2)-T(J,1))))) /T(J,4) C Step 7 700 S3=G(SC1*X)+SCLOG IF(DRNUNF().GE.DEXP(S3-T(J,3)-T(J,4)*(X-T(J,1))))GO TO 800 C CALL DOUT(SC1,1,'SC1') C CALL DOUT(X,1,'X') ZZZ=T(J,4)*(TSTAR-ZJ)*DEXP(-T(J,3)) 1 +DEXP(T(J,4)*(T(J,2)-T(J,1))) C CALL DOUT(ZZZ,1,'ZZZ') CALL DCOPY(5,T(J,1),LROW,TEMP,1) C CALL DOUT(TEMP,5,'T line') C CALL DOUT(TSTAR,1,'TSTAR') C CALL DWRRRN('T',JCAP,5,T,LROW,0) YLC=SC1*X RETURN C Step 8 800 JT=J IF(X.GT.T(J,1))JT=J+1 JTM=JT-1 JTP=JT+1 C Step 8(a) JCAP=JCAP+1 IF(JCAP.EQ.LROW)GO TO 1030 DO 810 J=JCAP,JT,-1 J1=J+1 DO 810 K=1,5 810 T(J1,K)=T(J,K) C Step 8(b) T(JT,1)=X T(JT,3)=S3 T(JT,4)=G1(SC1*X)*SC1 C Step 8(c) IF(JT.GT.1)T(JT,2) 1 =(T(JT,3)-T(JTM,3)+T(JTM,4)*T(JTM,1)-T(JT,4)*T(JT,1)) 2 /(T(JTM,4)-T(JT,4)) IF(JT.LT.JCAP)T(JTP,2) 1 =(T(JTP,3)-T(JT,3)+T(JT,4)*T(JT,1)-T(JTP,4)*T(JTP,1)) 2 /(T(JT,4)-T(JTP,4)) C Step 8(d) IF(JT.GT.1)YA=(DEXP(T(JTM,3)+T(JTM,4)*(T(JT,2)-T(JTM,1))) 2 -DEXP(T(JTM,3)+T(JTM,4)*(T(JTM,2)-T(JTM,1))))/T(JTM,4) YB=(DEXP(T(JT,3)+T(JT,4)*(T(JTP,2)-T(JT,1))) 2 -DEXP(T(JT,3)+T(JT,4)*(T(JT,2)-T(JT,1))))/T(JT,4) IF(JT.LT.JCAP)YC=(DEXP(T(JTP,3)+T(JTP,4)*(T(JT+2,2)-T(JTP,1))) 2 -DEXP(T(JTP,3)+T(JTP,4)*(T(JTP,2)-T(JTP,1))))/T(JTP,4) C Step 8(e) IF(JT.EQ.1)GO TO 820 A1=0.0D0 IF(JT.GE.3)A1=T(JT-2,5) T(JTM,5)=A1+YA C Step 8(f) 820 A1=0.0D0 IF(JT.GE.2)A1=T(JTM,5) T(JT,5)=A1+YB C Step 8(g) IF(JT.EQ.JCAP)GO TO 400 DEL=T(JT,5)+YC-T(JTP,5) DO 830 J=JTP,JCAP 830 T(J,5)=T(J,5)+DEL C Step 9 GO TO 400 C Errors 1010 WRITE(6,1020) 1020 FORMAT(/,' Adding-up violations in YLC.') CALL DWRRRN('T',JCAP+1,5,T,LROW,0) CALL DOUT(X1,1,'X1 ') CALL DOUT(X2,1,'X2 ') CALL TERM 1030 CALL PEND(' Storage capacity in YLC exhausted') 1050 WRITE(6,1060)C1,C2 1060 FORMAT(/,' Constraint interval [',1PE15.5,',',1PE15.5,'] in YLC') CALL TERM END REAL*8 FUNCTION YUM(F,C1,C2,XSTAR,D2LGF) C This subroutine draws from a distribution with a unimodal density on a C finite support. Documentation is in notes on the ARMA(p,q) paper. C Inputs: C F Function F(X) that has input ordinate X and returns density kernel F. C C1 Lower end of finite support C C2 Upper end of finite support C XSTAR Mode of density kernel C D2LGF Second derivative of log density kernel at mode C Output: C YUM Drawing from distribution IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL F DATA TSTAR,E1,E2,IR/10.0D0,0.0D0,1.0D-8,2/ C Step 1 100 CALL DQDAG(F,C1,XSTAR,E1,E2,IR,R1,E3) CALL DQDAG(F,XSTAR,C2,E1,E2,IR,R2,E3) R=R1+R2 A=1.0D0/DSQRT(-D2LGF) B1P=XSTAR B2M=XSTAR C Step 2 200 SSTAR=R*DRNUNF() S=0.0D0 C Step 3 300 IF(B1P.GE.C2)GO TO 400 310 Q=A F1P=F(B1P) QMAX=C2-B1P IF(Q.GT.QMAX)Q=QMAX 320 B2P=B1P+Q IF((TSTAR*F(B2P)).GT.F1P)GO TO 330 Q=.5D0*Q GO TO 320 330 CALL DQDAG(F,B1P,B2P,E1,E2,IR,SINC,E3) S=S+SINC IF(S.GE.SSTAR)GO TO 340 B1P=B2P GO TO 400 340 T=F1P B1=B1P B2=B2P GO TO 500 C Step 4 400 IF(B2M.GT.C1)GO TO 410 IF(B1P.GE.C2)GO TO 610 GO TO 310 410 Q=A F2M=F(B2M) QMAX=B2M-C1 IF(Q.GT.QMAX)Q=QMAX 420 B1M=B2M-Q IF((TSTAR*F(B1M)).GT.F2M)GO TO 430 Q=.5D0*Q GO TO 420 430 CALL DQDAG(F,B1M,B2M,E1,E2,IR,SINC,E3) S=S+SINC IF(S.GT.SSTAR)GO TO 440 B2M=B1M GO TO 300 440 T=F2M B1=B1M B2=B2M C Step 5 500 X=B1+(B2-B1)*DRNUNF() IF(T*DRNUNF().GT.F(X))GO TO 500 YUM=X RETURN C Errors 610 WRITE(6,620) 620 FORMAT(/,' Adding-up violations in YUM.') STOP END SUBROUTINE ZM1(NA,A,NB,B,NC,C) C Multiply two z-transforms, C(z)=A(z)*B(z) C Inputs: C NA Highest power of A C A Coefficients A(0) through A(NA) C NB Hightest power of B C B Coefficients B(0) through B(NB) C NC Highest power of C to compute C Outputs: C C Coefficients C(0) through C(NC) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(0:NA),B(0:NB),C(0:NC) MC=MIN0(NA+NB,NC) IF(MC.LT.NC)CALL DSET(NC-MC,0.0D0,C(MC+1),1) DO 20 IC=0,MC AA=0.0D0 K1=MAX0(0,IC-NB) K2=MIN0(IC,NA) DO 10 K=K1,K2 10 AA=AA+A(K)*B(IC-K) 20 C(IC)=AA RETURN END SUBROUTINE ZQ1(NA,A,NB,B,NC,C) C Quotient of two z-transforms, C(z)=A(z)/B(z) C Inputs: C NA Highest power of A C A Coefficients A(0) through A(NA) C NB Hightest power of B C B Coefficients B(0) through B(NB) C NC Highest power of C to compute C Outputs: C C Coefficients C(0) through C(NC) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(0:NA),B(0:NB),C(0:NC) C(0)=A(0)/B(0) DO 20 IC=1,NC AA=0.0D0 IF(IC.LE.NA)AA=A(IC) K2=MIN0(IC,NB) DO 10 K=1,K2 10 AA=AA-B(K)*C(IC-K) 20 C(IC)=AA/B(0) RETURN END SUBROUTINE ZI1(NA,A,NB,B) C Inverse of z-transform, B(z)=1/A(z) C Inputs: C NA Highest power of A C A Coefficients A(0) through A(NA) C NB Highest power of B to compute C Outputs: C B Coefficients B(0) through B(NB) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(0:NA),B(0:NB) B(0)=1.0D0/A(0) DO 20 IB=1,NB AA=0.0D0 K2=MIN0(IB,NA) DO 10 K=1,K2 10 AA=AA+A(K)*B(IB-K) 20 B(IB)=-B(0)*AA RETURN END SUBROUTINE ZM2(NA,NRA,NCA,A,LDRA,LDCA,NB,NRB,NCB,B,LDRB,LDCB, 1 NC,NRC,NCC,C,LDRC,LDCC) C Multiply two matrix z-transforms, C(z)=A(z)*B(z) C Inputs: C NA Highest power of A C NRA Number of rows in A C NCA Number of columns in A C A Coefficients A(0) through A(NA) C LDRA Leading row dimension of A C LDCA Leading column dimension of A C NB Hightest power of B C NRB Number of rows in B C NCB Number of columns in B C B Coefficients B(0) through B(NB) C LDRB Leading row dimension of B C LDCB Leading column dimension of B C NC Highest power of C to compute C NRC Number of rows in C C NCC Number of columns in C C LDRC Leading row dimension of C C LDCC leading column dimension of C C Outputs: C C Coefficients C(0) through C(NC) PARAMETER(LD=100) IMPLICIT REAL*8 (A-H,O-Z) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION A(LDRA,LDCA,0:NA),B(LDRB,LDCB,0:NB),C(LDRC,LDCC,0:NC) MC=MIN0(NA+NB,NC) DO 10 IC=0,NC 10 CALL UM0SET(NRC,NCC,C(1,1,IC),LDRC) DO 20 IC=0,MC K1=MAX0(0,IC-NB) K2=MIN0(IC,NA) DO 20 K=K1,K2 CALL DMRRRR(NRA,NRB,A(1,1,K),LDRA,NRB,NCB,B(1,1,IC-K),LDRB, 1 NRC,NCC,A1,LD) 20 CALL UMADD(NRC,NCC,C(1,1,IC),LDRC,A1,LD,C(1,1,IC),LDRC) RETURN END SUBROUTINE ZQ2(NA,NRA,NCA,A,LDRA,LDCA,NB,NRB,NCB,B,LDRB,LDCB, 1 NC,NRC,NCC,C,LDRC,LDCC) C Quotient of two matrix z-transforms, C(z)=[B(z)^-1]*A(z) C Inputs: C NA Highest power of A C NRA Number of rows in A C NCA Number of columns in A C A Coefficients A(0) through A(NA) C LDRA Leading row dimension of A C LDCA Leading column dimension of A C NB Hightest power of B C NRB Number of rows in B C NCB Number of columns in B C B Coefficients B(0) through B(NB) C LDRB Leading row dimension of B C LDCB Leading column dimension of B C NC Highest power of C to compute C NRC Number of rows in C C NCC Number of columns in C C LDRC Leading row dimension of C C LDCC leading column dimension of C C Outputs: C C Coefficients C(0) through C(NC) PARAMETER(LD=100) IMPLICIT REAL*8 (A-H,O-Z) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION A(LDRA,LDCA,0:NA),B(LDRB,LDCB,0:NB),C(LDRC,LDCC,0:NC) CALL DLINRG(NRB,B,LDRB,A2,LD) CALL DMRRRR(NRB,NCB,A2,LD,NRA,NCA,A,LDRA,NRC,NCC,C,LDRC) DO 20 IC=1,NC CALL UM0SET(NRA,NCA,A1,LD) IF(IC.LE.NA)CALL UCRGRG(NRA,NCA,A(1,1,IC),LDRA,A1,LD) K2=MIN0(IC,NB) DO 10 K=1,K2 CALL DMRRRR(NRB,NCB,B(1,1,K),LDRB,NRC,NCC,C(1,1,IC-K),LDRC, 1 NRC,NCC,C(1,1,IC),LDRC) 10 CALL UMSUB(NRC,NCC,A1,LD,C(1,1,IC),LDRC,A1,LD) 20 CALL DMRRRR(NRB,NCB,A2,LD,NRC,NCC,A1,LD,NRC,NCC,C(1,1,IC),LDRC) RETURN END SUBROUTINE ZI2(M,NA,A,LDRA,LDCA,NB,B,LDRB,LDCB) C Invert a matrix z-transform, B(z)=A(z)^-1 C Inputs: C NA Highest power of A C NRA Number of rows in A C NCA Number of columns in A C A Coefficients A(0) through A(NA) C LDRA Leading row dimension of A C LDCA Leading column dimension of A C NB Hightest power of B to compute C NRB Number of rows in B C NCB Number of columns in B C LDRB Leading row dimension of B C LDCB Leading column dimension of B C Outputs: C B Coefficients B(0) through B(NB) PARAMETER(LD=100) IMPLICIT REAL*8 (A-H,O-Z) COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD) DIMENSION A(LDRA,LDCA,0:NA),B(LDRB,LDCB,0:NB) CALL DLINRG(M,A,LDRA,B,LDRB) DO 20 IB=1,NB CALL UM0SET(M,M,A1,LD) K2=MIN0(IB,NA) DO 10 K=1,K2 CALL DMRRRR(M,M,A(1,1,K),LDRA,M,M,B(1,1,IB-K),LDRB, 1 M,M,B(1,1,IB),LDRB) 10 CALL UMSUB(M,M,A1,LD,B(1,1,IB),LDRB,A1,LD) 20 CALL DMRRRR(M,M,B,LDRB,M,M,A1,LD,M,M,B(1,1,IB),LDRB) RETURN END