SUBROUTINE ARPACF(NP,BETA,R)
C  This routine computes the p partial autocorrelation function coefficients,
C  given an autoregression of order p, y(t)=Beta(1)y(t-1)+...+Beta(p)y(t-p)
C  + Epsilon(t).
C  Inputs:
C  NP     Order of autoregression p
C  BETA   AR coefficients, BETA(j)=Beta(j)
C  Output:
C  R      Vector of partial autocorrelation coefficients, in order.
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION BETA(NP),R(NP)
      COMMON/SCRA/V1(10100),V2(10100)
      CALL DCOPY(NP,BETA,1,V2,1)
      IF(NP.EQ.1)GO TO 30
      DO 20 K=NP,2,-1
      RR=V2(K)
      AA=1.0D0/(1.0D0-RR**2)
      K1=K-1
      DO 10 I=1,K1
   10 V1(I)=(V2(I)+RR*V2(K-I))*AA
   20 CALL DCOPY(K1,V1,1,V2,1)
   30 CALL DCOPY(NP,V2,1,R,1)
      RETURN
      END 
C
      SUBROUTINE PACFAR(NP,R,BETA)
C  This routine computes the autoregression of order p corresponding to the
C  first p autocorrrelation function coefficients.  It is the inverse of
C  ARPACF.
C  Inputs:
C  NP     Order of autoregression p
C  R      Vector of partial autocorrelation coefficients, in order.
C  Output:
C  BETA   AR coefficients, in order.
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(NP),BETA(NP)
      COMMON/SCRA/V1(10100),V2(10100)
      V1(1)=R(1)
      IF(NP.EQ.1)GO TO 30
      DO 20 K=2,NP
      RR=R(K)
      DO 10 I=1,K-1
   10 V2(I)=V1(I)-RR*V1(K-I)
      V2(K)=RR
   20 CALL DCOPY(K,V2,1,V1,1)
   30 CALL DCOPY(NP,V1,1,BETA,1)
      RETURN
      END
C      
      REAL*8 FUNCTION ARPALJ(NP,R)
C  This function finds the log determinant Jacobian of transformation from p
C  partial autocorrelation function coefficients to the p coefficients in an
C  AR(p) -- i.e., log{det[partial(beta)/partial(r)]}.
C  Inputs:
C  NP     Order of autoregression p
C  R      Vector of partial autocorrelation coefficients, in order.
C  Output:
C  ARPALJ Log determinant Jacobian of transformation
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(NP)
      AA=0.0D0
      IF(NP.EQ.1)GO TO 30
      AA=DLOG(1.0D0-R(2))
      IF(NP.EQ.2)GO TO 30
      DO 10 J=3,NP,2
   10 AA=AA+((J-1)/2)*DLOG(1.0D0-R(J)**2)
C  10 AA=AA+DLOG(1.0D0-R(J)**(J-1))
      IF(NP.EQ.3)GO TO 30
      DO 20 J=4,NP,2
   20 AA=AA+((J-2)/2)*DLOG(1.0D0-R(J)**2)+DLOG(1.0D0-R(J))
C  20 AA=AA+DLOG((1.0D0-R(J))*(1.0D0-R(J)**(J-2)))
   30 ARPALJ=AA
      RETURN
      END