SUBROUTINE CPD
C  This routine implements a Gibbs/Metropoils algorithm for sampling from the
C  posterior distribution of a linear regression model with stationary AR(P
C  disturbances.  The prior for coefficients (Beta) is multivariate normal,
C  the prior for the disturbance (h) is chi-square, and the prior for the
C  AR coefficients (Phi) is multivariate normal truncated to the invertible
C  region.
      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL LFLIP,LFORE,LINOP
      PARAMETER(LDK=20,LDP=10,NCD=30,LDT=2000)
      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),E(LDT,LDP),EPS(LDT),V1(LDT),X(LDT,LDK),
     1          Z(LDT,LDK)
      DIMENSION BETA(LDK),B0(LDK),HB0(LDK),H0(LDK,LDK)
      DIMENSION NV1(LDK)
      DIMENSION APT(LDP,LDP),APTCAN(LDP,LDP),A1(LDP,LDP),
     1          H0PHI(LDP,LDP),VP(LDP,LDP)
      DIMENSION B0PHI(LDP),HB0PHI(LDP),PHI(LDP),PHICAN(LDP),V(0:LDP)
C  Entry point CPD
C  Draw Phi
  100 DO 110 IP=1,NP
  110 CALL DCOPY(NTMNP,EPS(NP1-IP),1,E(1,IP),1)
      CALL NRMK(NP,NTMNP,E,LDT,EPS(NP1),H,H0PHI,LDP,HB0PHI,PHICAN)
      IF(.NOT.LINOP(NP,PHICAN))GO TO 200
      CALL ARACF(NP,PHICAN,1.0D0,V)
      DO 120 IP=1,NP
      DO 120 JP=1,NP
  120 VP(IP,JP)=V(IABS(IP-JP))
      CALL DLFTDS(NP,VP,LDP,A1,LDP)
      CALL DLINRT(NP,A1,LDP,2,APTCAN,LDP)
      AACAN=0.0D0
      BBOLD=0.0D0
      BBCAN=0.0D0
      DO 130 IP=1,NP
      AACAN=AACAN+DLOG(APTCAN(IP,IP))
      BBOLD=BBOLD+DDOT(IP,APT(1,IP),1,EPS,1)**2
  130 BBCAN=BBCAN+DDOT(IP,APTCAN(1,IP),1,EPS,1)**2
      IF(.NOT.LFLIP(AAOLD-AACAN+.5D0*H*(BBOLD-BBCAN)))GO TO 200
      NEW=NEW+1
      AAOLD=AACAN
      CALL DCOPY(NP,PHICAN,1,PHI,1)
      CALL DCRGRG(NP,APTCAN,LDP,APT,LDP)
C  Draw Beta
  200 DO 220 IK=1,NK1
      DO 210 IT=1,NP
  210 X(IT,IK)=DDOT(IT,APT(1,IT),1,Z(1,IK),1)
      DO 220 IT=NP1,NT
  220 X(IT,IK)=Z(IT,IK)-DDOT(NP,PHI,1,Z(IT-NP,IK),-1)
      CALL NRMK(NK,NT,X,LDT,X(1,NK1),H,H0,LDK,HB0,BETA)
      CALL DMURRV(NT,NK,Z,LDT,NK,BETA,1,NT,EPS)
      CALL UVSUB(NT,Z(1,NK1),EPS,EPS)
C  Draw h
  300 CALL DMURRV(NT,NK,X,LDT,NK,BETA,1,NT,V1)
      CALL UVSUB(NT,X(1,NK1),V1,V1)
      S21=DNRM2(NT,V1,1)**2
      H=CHI(S20+S21,ANU2)
C  Record
  400 IF(.NOT.(LRECRD.AND.LWRITE))RETURN
      CALL DCOPY(NK,BETA,1,PAR,1)
      PAR(LH)=H
      CALL DCOPY(NP,PHI,1,PAR(LPHI),1)
      PDAT=PDAT0+HNT*DLOG(H)-AAOLD-0.5D0*S21*H
      PRI=PRI0+GAUK(NK,B0,H0,LDK,BETA)
     1        +GAUK(NP,B0PHI,H0PHI,LDP,PHI)+CHIK(S20,ANU0,H)
      IF(.NOT.LFORE)RETURN
      SD=1.0D0/DSQRT(H)
      DO 410 IF=1,NF
      JT=NT+IF
      EPS(JT)=DDOT(NP,PHI,1,EPS(JT-NP),-1)+SD*DRNNOF()
  410 PAR(LF+IF-1)=DDOT(NK,BETA,1,Z(IT,1),LDT)+EPS(JT)
      RETURN
C  Entry point PD0
      ENTRY PD0
      WRITE(6,1010)
 1010 FORMAT(10X,'Linear model with normal AR(p) disturbances')
C  Obtain model specification
      CALL FOPEN('Model specification',1,.FALSE.)
      READ(1,*,END=2010,ERR=2020)NTA,NTB,NK,NP,NF
      CALL PINTI('t1',NTA,1,NTB)
      CALL PINTI('t2',NTB,NTA,LDT)
      NT=NTB-NTA+1
      CALL PINTI('k',NK,1,MIN0(LDK-1,LDPAR-2))
      NK1=NK+1
      CALL PINTI('p',NP,1,MIN0(NT,LDP,LDPAR-NK-1))
      NP1=NP+1
      NTMNP=NT-NP
      CALL PINTI('f',NF,0,MIN0(LDT-NT,LDPAR-NK-NP-1))
      NTF=NT+NF
      LFORE=.TRUE.
      IF(NF.EQ.0)LFORE=.FALSE.
C  Get data
      CALL DAT1(2,1,NTB+NF,NVARS,D,LDT,NCD)
C  Get regressors
      CALL LISTR(1,NK1,NVARS,.TRUE.,NV1)
      DO 1020 IK=1,NK1
 1020 CALL DCOPY(NTF,D(NTA,NV1(IK)),1,Z(1,IK),1)
      WRITE(6,1030)NTA,NTB,NK,NP,NV1(NK1),(NV1(IK),IK=1,NK)
 1030 FORMAT(/,'       First observation (t1):',I5,/,  
     1         '        Last observation (t2):',I5,/,
     2         '     Number of regressors (k):',I5,/,
     3         '  Order of autoregression (p):',I5,/,
     4         '    Dependent variable column:',I5,/,
     5         ' Explanatory variable columns:',9I5,/,(30X,9I5))
C  Get priors
      CALL GAU0(1,NK,'covariate coefficients',H0,LDK,B0,HB0,PRI0B,BETA)
      CALL DMURRV(NT,NK,Z,LDT,NK,BETA,1,NT,EPS)
      CALL UVSUB(NT,Z(1,NK1),EPS,EPS)
      CALL GAUI0(1,NP,'AR(p) coefficients',H0PHI,LDP,B0PHI,HB0PHI,
     1           PRI0PH,PHI)
      NEW=0
      AAOLD=DMACH(7) 
      CALL CHI0(1,1,'precision',.FALSE.,ANU0,S20,PRI0H,H)
      ANU2=ANU0+NT
C  Simulation vector organization
      IF(.NOT.LWRITE)RETURN
      LH=NK+1
      LPHI=LH+1
      NPAR=LPHI+NP-1
      CALL PINTI('Number of parameters',NPAR,1,LDPAR)
      WT=0.0D0
      HNT=0.5D0*NT
      PRI0=PRI0B+PRI0PH+PRI0H
      PDAT0=NT*HLG2PI
      WRITE(6,1040)1,LH,LPHI
 1040 FORMAT(/,' Organization of simulation vector',/,
     1         ' Covariate coefficient vector:',I5,/,
     2         '                    Precision:',I5,/,
     3         '     AR(p) coefficient vector:',I5)
      IF(.NOT.LFORE)RETURN
      CALL PINTI('Number of forecasts',NF,1,LDPAR-NPAR)
      LF=LPHI+NP
      NPAR=LF+NF-1
      WRITE(6,1050)NF,LF
 1050 FORMAT(I19,' forecasts:',I5)
      RETURN
 2010 CALL PEND('UNANTICIPATED END OF FILE READING LINE 1 OF MODEL')
 2020 CALL PEND('ERROR READING LINE 1 OF MODEL')
      ENTRY PD1
      WRITE(6,1510)NEW
 1510 FORMAT(/,I10,' accepted Metropolis candidates') 
      RETURN
      END