SUBROUTINE CPD
C  This routine implements sampling for sur1, the normal seemingly
C  unrelated regressions model with a proper, multivariate normal
C  coefficient prior and a Wishart precision matrix prior.  The 
C  algorithm is Gibbs sampling, with two blocks: coefficients, and
C  precision matrix.  This version works directly with sums of squares and
C  cross products, and has no forecasting option.
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100,LDK=100,LDM=10,LDKK=100,LDD=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(LDKK),HB0(LDKK),HB1(LDKK),BETA1(LDKK),
     1  BETA2(LDKK)
      DIMENSION H0(LDKK,LDKK),H1(LDKK,LDKK)
      DIMENSION DD(LDD,LDD),ZZ(LDK,LDK,LDM,LDM),B(LDK,LDM)
C  Entry point CPD
C  Draw Beta
      IB=0
      DO 30 IM=1,NM
      JB=IB
      DO 20 JM=IM,NM
      HH=H(IM,JM)
      DO 10 IK=1,NV1(IM)
      DO 10 JK=1,NV1(JM)
   10 H1(IB+IK,JB+JK)=HH*ZZ(IK,JK,IM,JM)
   20 JB=JB+NV1(JM)
   30 IB=IB+NV1(IM)
      CALL DFULL(KK,H1,LDKK)
      IB=0
      DO 50 IM=1,NM
      DO 50 IK=1,NV1(IM)
      AA=0.0D0
      DO 40 JM=1,NM
   40 AA=AA+H(IM,JM)*ZZ(IK,NV2(JM),IM,JM)
      IB=IB+1
   50 HB1(IB)=AA
      CALL GAU2(KK,H0,LDKK,H1,LDKK,HB0,HB1,BETA1,BETA2)
C  Draw H
      IB=1
      DO 110 IM=1,NM
      CALL DCOPY(NV1(IM),BETA1(IB),1,B(1,IM),1)
  110 IB=IB+NV1(IM)
      DO 120 IM=1,NM
      DO 120 JM=IM,NM
  120 S1(IM,JM)
     1  =DBLINF(NV2(IM),NV2(JM),ZZ(1,1,IM,JM),LDK,B(1,IM),B(1,JM))
      CALL DFULL(NM,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(KK,B0,H0,LDKK,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)
      IB=1
      DO 210 IM=1,NM
      NB=NV1(IM)
      CALL DCOPY(NB,BETA1(IB),1,PAR(IB),1)
      CALL DCOPY(NB,BETA2(IB),1,PAR(KK+IB),1)
  210 IB=IB+NB
      CALL VECLT(NM,H,LDM,PAR(LH))
      RETURN
C  Entry PD0
C  Get prior and data
      ENTRY PD0
      WRITE(6,1010)
 1010 FORMAT(10X,'Seemingly unrelated regressions model with normal',
     1  'disturbances')
      CALL FOPEN('Model specification',1,.FALSE.)
      READ(1,*,END=2010,ERR=2020)NTA,NTB,NM
      CALL PINTI('t1',NTA,1,NTB)
      CALL PINTI('m',NM,1,MIN0(LDM,LD))
      CALL DAT0(2,NTA,NTB,NVARS,DD,LDD)
      NT=NTB-NTA+1
      CALL LISTR(1,NM,LDK,.FALSE.,NV1)
      KK=0
      CALL IOUT(NM,1,'NM')
      DO 1020 IM=1,NM
      KK=KK+NV1(IM)
      NV2(IM)=NV1(IM)+1
 1020 CALL LISTR(1,NV2(IM),NVARS,.TRUE.,NA1(1,IM))
      CALL IOUT(NV2,NM,'NV2')
      DO 1030 IM=1,NM
      DO 1030 JM=1,NM
      CALL CPROD(NV2(IM),NV2(JM),NA1(1,IM),NA1(1,JM),DD,LDD,
     1  ZZ(1,1,IM,JM),LDK)
 1030 CONTINUE
      CALL PINTI('Total regressors',KK,1,LDKK)
      WRITE(6,1040)NTA,NTB
 1040 FORMAT(/,'       First observation (t1):',I5,/,
     1         '        Last observation (t2):',I5,//,
     2         ' Equation    Dep var    Regressors')
      DO 1050 IM=1,NM
 1050 WRITE(6,1060)IM,NA1(NV2(IM),IM),(NA1(J,IM),J=1,NV1(IM))
 1060 FORMAT(I5,I12,5X,15I3,/,(22X,15I3))
      CALL IOUT(KK,1,'KK')
      CALL GAU0(1,KK,'coefficients',H0,LDKK,B0,HB0,AA,BETA1)
      CALL WISH0(1,NM,'precision matrix',S0,LDM,ANU0,BB,P,LDM,H,LDM)
      DO 1070 IM=1,NM
 1070 B(NV2(IM),IM)=-1.0D0
      ANU2=NT+ANU0
C  Simulation vector organization
      IF(.NOT.LWRITE)RETURN
      LB2=1+KK
      LH=LB2+KK
      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,1080)1,LB2,LH
 1080 FORMAT(/,' Organization of simulation vector',/,
     1         '             Covariate coefficient vector (Beta):',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('UNANTICIP[ATED END OF FILE READING PRIOR')
 2020 CALL PEND('ERROR READING PRIOR')
      END