SUBROUTINE CPD
C  This routine implements sampling for uvr1, the normal linear model
C  with a proper, multivariate normal coefficient prior and a chi square
C  precision prior.  The algorithm is Gibbs sampling, with two blocks:
C  coefficient vector, and precision parameter.  This version works directly
C  with sums of squares and cross products, and has no forecasting option.
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100,LDK=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
      DIMENSION NV1(LD)
      DIMENSION H0(LDK,LDK),H1(LDK,LDK)
      DIMENSION B0(LDK),HB0(LDK),HB1(LDK),BETA1(LDK),BETA2(LDK)
      DIMENSION ZZ(LDK,LDK)
C  Entry point CDP
C  Draw Beta
      DO 10 IK=1,NK
      HB1(IK)=H*ZZ(IK,NK1)
      DO 10 JK=1,IK
      H1(IK,JK)=H*ZZ(IK,JK)
   10 H1(JK,IK)=H1(IK,JK)
      CALL GAU2(NK,H0,LDK,H1,LDK,HB0,HB1,BETA1,BETA2)
C  Draw h
      SSR=UQUAF(NK1,BETA1,ZZ,LDK)
      H=CHI(S20+SSR,ANU2)
C  Record
      IF(.NOT.(LRECRD.AND.LWRITE))RETURN
      CALL DCOPY(NK,BETA1,1,PAR,1)
      CALL DCOPY(NK,BETA2,1,PAR(NK1),1)
      PAR(NPAR)=H
      PRI=PRI0+GAUK(NK,B0,H0,LDK,BETA1)+CHIK(S20,ANU0,H)
      PDAT=PDAT0+HNT*DLOG(H)-.5D0*H*SSR
      RETURN
C  Entry point PD0
      ENTRY PD0
      WRITE(6,1010)
 1010 FORMAT(10X,'Linear model with normal disturbances')
C  Obtain basic model specification parameters
      CALL FOPEN('Model specification',1,.FALSE.)
      READ(1,*,END=2010,ERR=2020)NTA,NTB,NK
      CALL PINTI('t1',NTA,1,NTB)
      CALL PINTI('k',NK,1,MIN0(LD,LDK-1))
      NK1=NK+1
C  Obtain data (No limit on NT)
      CALL DAT0(2,NTA,NTB,NVARS,ZZ,LDK)
      NT=NTB-NTA+1
      CALL LISTR(1,NK1,NVARS,.TRUE.,NV1)
      CALL CPROD(NK1,NK1,NV1,NV1,ZZ,LDK,ZZ,LDK)
      WRITE(6,1020)NTA,NTB,NK,NV1(NK1),(NV1(IK),IK=1,NK)
 1020 FORMAT(/,'       First observation (t1):',I5,/,
     1         '        Last observation (t2):',I5,/,
     2         '     Number of regressors (k):',I5,/,
     3         '    Dependent variable column:',I5,/,
     4         ' Explanatory variable columns:',9I5,/,(30X,9I5))
      CALL GAU0(1,NK,'coefficients',H0,LDK,B0,HB0,AA,BETA1)
      CALL CHI0(1,1,'precision',.FALSE.,ANU0,S20,BB,H)
      BETA1(NK1)=-1.0D0
      ANU2=NT+ANU0
      IF(.NOT.LWRITE)RETURN
C  Simulation vector organization
      NPAR=2*NK+1
      CALL PINTI('Number of parameters',NPAR,1,LDPAR)
      WT=0.0D0
      PRI0=AA+BB
      HNT=0.5D0*NT
      PDAT0=-NT*HLG2PI
      WRITE(6,1030)1,NK1,NPAR
 1030 FORMAT(/,' Organiation of simulation vector',/,
     1         ' Covariate coefficient vector (Beta):',I5,/,
     2         '    Antithetic of coefficient vector:',I5,/,
     3         '           Disturbance precision (h):',I5)
      RETURN
C  Entry point PD1
      ENTRY PD1
      RETURN
C  Error processing
 2010 CALL PEND('UNANTICIPATED END OF FILE READING MODEL')
 2020 CALL PEND('ERROR READING MODEL')
      END