SUBROUTINE CPD
C  This routine implements sampling for probit1, the normal linear probit model
C  with a rpoper, multivariate normal coefficient prior.  The algorithm is
C  Gibbs sampling - data augmentation with T+1 blocks: the coefficient vector,
C  and one block for each probit.
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LDK=20,NCD=20,LDT=20021)
      LOGICAL LFORE,LY(LDT)
      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 D(LDT,NCD),Z(LDT,LDK),Y(LDT)
      DIMENSION H0(LDK,LDK),XX(LDK,LDK)
      DIMENSION B0(LDK),HB0(LDK),HB1(LDK),BETA(LDK),BETA2(LDK)
      DIMENSION NV1(LDK)
C  Entry point PD0
C  Draw probits.
  100 CALL DMURRV(NT,NK,Z,LDT,NK,BETA,1,NT,Y)
      DO 110 IT=1,NT
      AA=Y(IT)
  110 Y(IT)=XNOT(AA,1.0D0,0.0D0,0.0D0,.NOT.LY(IT),LY(IT))
C  Draw Beta.
  200 CALL DMURRV(NT,NK,Z,LDT,NT,Y,2,NK,HB1)
      CALL GAU2(NK,H0,LDK,XX,LDK,HB0,HB1,BETA,BETA2)
C  Record.
  300 IF(.NOT.(LRECRD.AND.LWRITE))RETURN
      CALL DCOPY(NK,BETA,1,PAR,1)
      IF(.NOT.LFORE)GO TO 320
      DO 310 IF=1,NF
      AA=DDOT(NK,Z(NT+IF,1),LDT,BETA,1)
  310 PAR(NK+IF)=1.0D0-PNORDF(AA,1.0D0,0.0D0)
  320 PRI=PRI0+GAUK(NK,B0,H0,LDK,BETA)
      PDAT=0.0D0
      DO 330 IT=1,NT
      AA=DDOT(NK,Z(IT,1),LDT,BETA,1)
      P=PNORDF(AA,1.0D0,0.0D0)
      IF(LY(IT))P=1.0D0-P
      PLOG=-1.0D100
      IF(P.GT.0.0D0)PLOG=DLOG(P)
      IF(P.EQ.0.0D0)WRITE(6,325)ITER,IT
  325 FORMAT(' Zero prob at iter',I6,'  t=',I4)
  330 PDAT=PDAT+PLOG
      RETURN
C  Entry point PD0
      ENTRY PD0
      WT=0.0D0
C  Obtain model specification.
      WRITE(6,1010)
 1010 FORMAT(10X,'Normal linear univariate probit model')
      CALL FOPEN('Model specification',1,.FALSE.)
      READ(1,*,END=2010,ERR=2020)NTA,NTB,NK,NF
      CALL PINTI('t1',NTA,1,NTB)
      CALL PINTI('t2',NTB,NTA,LDT)
      NTA1=NTA-1
      NT=NTB-NTA1
      CALL PINTI('k',NK,1,MIN0(LDK,LDPAR-1))
      NK1=NK+1
      CALL PINTI('f',NF,0,MIN0(LDT-NT,LDPAR-NK-1))
      NTF=NT+NF
      LFORE=.TRUE.
      IF(NF.EQ.0)LFORE=.FALSE.
C  Obtain data.
      CALL DAT1(2,1,NTB+NF,NVARS,D,LDT,NCD)
      CALL LISTR(1,NK1,NVARS,.TRUE.,NV1)
      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 explanatory variabiles (k):',I5,/,
     3         '            Dependent variable column:',I5,/,
     4         '         Explanatory variable columns:',9I5,/,
     5           (30X,9I5))
      DO 1030 IK=1,NK
 1030 CALL DCOPY(NTF,D(NTA,NV1(IK)),1,Z(1,IK),1)
      CALL DMXTXF(NT,NK,Z,LDT,NK,XX,LDK)
      J=NV1(NK1)
      DO 1040 IT=1,NT
      AA=D(NTA1+IT,J)
      LY(IT)=.FALSE.
      IF(AA.EQ.0.0D0)GO TO 1040
      IF(AA.NE.1.0D0)GO TO 2030
      LY(IT)=.TRUE. 
 1040 CONTINUE
C  Obtain prior
      CALL GAU0(1,NK,'covariate coefficients',H0,LDK,B0,HB0,PRI0,BETA)
      IF(.NOT.LWRITE)RETURN
C  Simulation vector organization
      NPAR=NK+NF
      CALL PINTI('Number of parameters',NPAR,1,LDPAR)
      WRITE(6,1050)1
 1050 FORMAT(/,' Organization of simulation vector',/,
     1         ' Covariate coefficient vector:',I5)
      IF(.NOT.LFORE)RETURN
      LF=NK+1
      WRITE(6,1060)NF,LF
 1060 FORMAT(I7,' probability forecasts:',I5)
      RETURN
 2010 CALL PEND('UNANTICIPATED END OF FILE READING MODEL')
 2020 CALL PEND('ERROR READING MODEL')
 2030 WRITE(6,2035)NTA1+IT,J,AA
 2035 FORMAT(/,'   VALUE OF OBSERVATION',I5,' OF DEP VAR',I3,' IS',
     1       1PE15.6)
      CALL TERM
C  Entry point PD1
      ENTRY PD1
      RETURN
      END