SUBROUTINE CPD
C  This routine implements sampling for mvr1, the normal multivariate
C  regression model with a proper, multivariate normal coefficient
C  prior and a Wishart precision matrix prior.  The algorithm is Gibbs
C  sampling, with two blocks: coefficients, and precision matrix.  This
C  version works directly with sums of squares and cross products, and has
C  no forecasting option.
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100,LDK=100,LDM=10,LDKM=100,LDZ=100)
      PARAMETER(LDPAR=250)
      COMMON/AMCMC/WT,PRI,PDAT,PAR(LDPAR),ITER,NPAR,LRECRD,LWRITE,LERR
      LOGICAL LRECRD,LWRITE,LERR
      COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII
      COMMON/ISCRA/NV1(LD),NV2(LD),NA1(LD,LD),NA2(LD,LD)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION H(LDM,LDM),P(LDM,LDM),
     1  S0(LDM,LDM),S1(LDM,LDM),S2(LDM,LDM)
      DIMENSION B0(LDKM),HB0(LDKM),HB1(LDKM),BETA1(LDKM),BETA2(LDKM)
      DIMENSION H0(LDKM,LDKM),H1(LDKM,LDKM)
      DIMENSION ZZ(LDZ,LDZ),B(LDZ,LDM)
C  Entry point CPD
C  Draw Beta
      CALL UCROSS(NM,NM,H,LDM,NK,NK,ZZ,LDZ,NKM,NKM,H1,LDKM)
      CALL DMRRRR(NK,NM,ZZ(1,NK1),LDZ,NM,NM,H,LDM,NK,NM,A1,LD)
      CALL VEC(NK,NM,A1,LD,HB1)
      CALL GAU2(NKM,H0,LDKM,H1,LDKM,HB0,HB1,BETA1,BETA2)
C  Draw H
      CALL UNVEC(NK,NM,BETA1,B,LDZ)
      CALL UXTAXF(NZ,NM,B,LDZ,ZZ,LDZ,S1,LDM)
      CALL UMADD(NM,NM,S0,LDM,S1,LDM,S2,LDM)
      CALL DLINDS(NM,S2,LDM,S2,LDM)
      CALL WISH(NM,S2,LDM,ANU2,H,LDM)
C  Record
      IF(.NOT.(LRECRD.AND.LWRITE))RETURN
      PRI=PRI0+GAUK(NKM,B0,H0,LDKM,BETA1)+WISHK(NM,P,LDM,ANU0,H,LDM)
      PDAT=PDAT0+HNT*DTLGPD(NM,H,LDM)
     1      -.5D0*TRAB(NM,NM,S1,LDM,NM,NM,H,LDM)
      CALL DCOPY(NKM,BETA1,1,PAR,1)
      CALL DCOPY(NKM,BETA2,1,PAR(LB2),1)
      CALL VECLT(NM,H,LDM,PAR(LH))
      RETURN
C  Entry point PD0
C  Get prior and data
      ENTRY PD0
      WRITE(6,1010)
 1010 FORMAT(10X,'Multivariate linear model with normal disturbances')
      CALL FOPEN('Model specification',1,.FALSE.)
      READ(1,*,END=2010,ERR=2020)NTA,NTB,NK,NM
      CALL PINTI('t1',NTA,1,NTB)
      CALL PINTI('k',NK,1,MIN0(LDK-1,LD))
      NK1=NK+1
      CALL PINTI('m',NM,1,MIN0(LDM,LD))
      NZ=NK+NM
      CALL PINTI('k+m',NZ,1,LDKM)
      NKM=NK*NM
      CALL PINTI('k*m',NKM,1,LDKM)
      CALL DAT0(2,NTA,NTB,NVARS,ZZ,LDZ)
      NT=NTB-NTA+1
C  Get regressors.  
      CALL LISTR(1,NZ,NVARS,.TRUE.,NV1)
      CALL CPROD(NZ,NZ,NV1,NV1,ZZ,LDZ,ZZ,LDZ)
      WRITE(6,1020)NTA,NTB,NK,(NV1(NK+IM),IM=1,NM)
 1020 FORMAT(/,'       First observation (t1):',I5,/,
     1         '        Last observation (t2):',I5,/,
     2         ' Number of dep. variables (m):',I5,/,
     3         '     Number of regressors (k):',I5,/,
     4         '   Dependent variable columns:',9I5,/,(30X,9I5))   
      WRITE(6,1030)(NV1(IK),IK=1,NK)
 1030 FORMAT(' Explanatory variable columns:',9I5,/,(30X,9I5))
C  Get priors.
      CALL GAU0(1,NKM,'coefficients',H0,LDKM,B0,HB0,AA,BETA1)
      CALL WISH0(1,NM,'precision matrix',S0,LDM,ANU0,BB,P,LDM,H,LDM)
      CALL UM0SET(NM,NM,B(NK1,1),LDZ)
      DO 1040 IM=1,NM
 1040 B(NK+IM,IM)=-1.0D0
      ANU2=NT+ANU0
C  Simulation vector organization
      IF(.NOT.LWRITE)RETURN
      LB2=1+NKM
      LH=LB2+NKM
      NPAR=LH+NM*(NM+1)/2-1
      CALL PINTI('Number of parameters',NPAR,1,LDPAR)
      WT=0.0D0
      PRI0=AA+BB
      HNT=0.5D0*NT
      PDAT0=-NT*NM*HLG2PI
      WRITE(6,1050)1,LB2,LH
 1050 FORMAT(/,' Organization of simulation vector',/,
     1  '   Coefficient vector Beta (B matrix by columns):',I5,/,
     2  '                Antithetic of coefficient vector:',I5,/,
     3  ' Disturbance precision matrix H (lower triangle):',I5)
      RETURN
C  Entry PD1
      ENTRY PD1
      RETURN
C  Error processing
 2010 CALL PEND('UNANTICIPATED END OF FILE READING PRIOR')
 2020 CALL PEND('ERROR READING PRIOR')
      END