C.hr snpcode5.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     snpcode5.f
C
C     Copyright (C) 1990.
C
C     A. Ronald Gallant & George E. Tauchen 
C     P.O. Box 5513, Raleigh NC 27650-5513, USA   
C
C     Permission to use, copy, modify, and distribute this
C     software and its documentation for any purpose and
C     without fee is hereby granted, provided that the
C     above copyright notice appear in all copies and that
C     both that copyright notice and this permission notice
C     appear in supporting documentation.
C
C     This software is provided "as is" without any
C     expressed or implied warranty.
C
C....*...1.........2.........3.........4.........5.........6.........7.*
C.hr WHX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     WHX      1/6/90
C
C     PURPOSE
C
C     COMPUTE WEIGHTS V(J) SUCH THAT THE HERMITE CONDITIONAL DENSITY
C     CAN BE COMPUTED AS
C
C       H(Z|X,THETA)=[SUM_J: V(J)* PRD_K: Z(K)**L(K,J)]**2 * N(Z)
C
C     WHERE L(K,J) DENOTES THE ELEMENTS OF APPROPRIATE MULTI-INDEXES
C     COMPUTED USING MULTIY.
C
C     COMPUTE WEIGHTS W(I,J) SUCH THAT
C
C       G(Z|X,THETA) = SUM_IJ: W(I,J) * PRD_K: C[Z(K)| L(K,I)+L(K,J)+1]
C
C     IS A W(I,J) WEIGHTED MIXTURE OF INDEPENDENT CHI RANDOM VARIABLES
C     AS DEFINED IN MONAHAN:
C
C        C(U|A) = [2**(1-A/2)/GAMMA(A/2)] * U**(A-1) * EXP(-X**2/2)
C
C     NOTE
C     ONE CAN DRAW FROM G(Z|X,THETA), CHANGE THE SIGN OF EACH COMPONENT
C     OF Z WITH PROBABILITY 1/2, AND USE
C
C       B(Z|X,THETA)=[SUM_J: V(J)* PRD_K: |Z(K)|**L(K,J)]**2 * N(Z)
C
C     TO MAJORIZE H(Z|X,THETA) TO GET A SAMPLE FROM H(Z|X,THETA) BY THE
C     REJECTION METHOD.  SEE KENNEDY AND GENTLE.
C
C     USAGE
C     CALL WHX(V,W,X,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,LWORK)
C
C     REFERENCES
C     MONAHAN, JOHN F. (1987) "AN ALGORITHM FOR GENERATING CHI RANDOM 
C     VARIABLES," ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOL. 13, 
C     NO. 2, P. 168-172. ALGORITHM 651. (1988) CORRIGENDUM, VOL. 13, 
C     NO. 1  P. 111.
C     KENNEDY, WILLIAM J. JR. AND JAMES E. GENTLE (1980), STATISTICAL
C     COMPUTING, NEW YORK: MARCEL DEKKER, INC., P. 186.
C
C     SUBROUTINES CALLED
C     MULTIY,MULTIX,SETMYX,EHYX,GAMMLN
C
C     ARGUMENTS
C     V    - KY-VECTOR OF COMPUTED WEIGHTS
C            OUTPUT, REAL*8
C     W    - KY*KY-VECTOR OF COMPUTED WEIGHTS
C            OUTPUT, REAL*8
C     X    - M*L VECTOR OF LAGGED VALUES OF Y, STORED END TO END
C            X=(Y'_T-L, ..., Y'_T-1), M*L MUST BE LESS THAN 100
C            INPUT, REAL*8
C     A    - COEFFICIENTS OF THE POLYNOMIAL, A KX BY KY MATRIX, STORED
C            COLUMNWISE
C            INPUT, REAL*8
C     B0   - M-VECTOR OF COEFFICIENTS, THE CONSTANT TERM OF THE VAR
C            INPUT, REAL*8
C     B    - COEFFICIENTS FOR THE REGRESSION ON X, AN M BY M*L MATRIX,
C            STORED COLUMNWISE
C            INPUT, REAL*8
C     R    - AN M*(M+1)/2-VECTOR CONTAINING THE UPPER TRIANGLE OF THE
C            FACTOR OF THE VARIANCE-COVARIANCE MATRIX OF THE VAR
C            STORED COLUMNWISE, EG. R=(R11,R12,R22,R13,R23,R33,...,RMM)
C            SIGMA=R*R', THIS IS NOT THE STANDARD CHOLESKY FACTORIZATION
C            AS IT IS UPPER TRIANGULAR, NOT LOWER, USE SUBROUTINE FACTOR
C            TO COMPUTE IT FROM SIGMA, GIVEN SIGMA.
C            INPUT, REAL*8
C     M    - DIMENSION OF THE TIME SERIES, IE. LENGTH OF Y_T, MUST BE
C            LESS THAN 100
C            INPUT, INTEGER*4
C     L    - NUMBER OF LAGS IN THE VAR, M*L MUST BE LESS THAN 100
C            INPUT, INTEGER*4
C     KY   - NUMBER OF TERMS IN THE POLYNOMIAL IN Y, COMPUTED BY SETMYX,
C            SEE REMARK BELOW.
C            INPUT, INTEGER*4
C     KX   - NUMBER OF TERMS IN THE POLYNOMIAL IN X, COMPUTED BY SETMYX,
C            SEE REMARK BELOW.
C            INPUT, INTEGER*4
C     IDEGY- DEGREE OF THE POLYNOMIAL IN Y, SEE REMARK BELOW
C            INPUT, INTEGER*4
C     IDEGX- DEGREE OF THE POLYNOMIAL IN X, SEE REMARK BELOW
C            INPUT, INTEGER*4
C     WORK - WORK VECTOR OF LENGTH M+KX*KY+M+M*M*L+M*(M+1)/2+2*KY.
C            REAL*8
C     LWORK- DECLARED DIMENSION OF WORK IN CALLING PROGRAM.
C            INPUT, INTEGER*4
C
C     REMARK
C     THE CALL
C     CALL SETMYX(M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX,KY,KX,IDUMMY,IWORK)
C     MUST PRECEDE USAGE OF THIS ROUTINE
C
      SUBROUTINE WHX(V,W,X,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,LWORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 V(KY),W(KY,KY),X(M*L)
      REAL*8 A(KX*KY),B0(M),B(M,M*L),R(M*(M+1)/2)
      REAL*8 WORK(LWORK)
      INTEGER*4 ALPHA(100),BETA(100),LI(100),LJ(100)
      SAVE
C
C     ERROR TRAP FOR: ALPHA(M), BETA(M*L), LI(M), LJ(M)
C
      IF ((M.GT.100).OR.(M*L.GT.100)) THEN
        WRITE(3,3002)
        STOP
      ENDIF
C
C     COMPUTE LENGTH OF STACKED PARAMETER VECTOR
C
      IBA  = 0
      IBB0 = IBA +KX*KY
      IBB  = IBB0+M
      IBR  = IBB +M*M*L
      LTHETA=IBR +M*(M+1)/2
C
C     COMPUTE POINTERS TO WORK VECTOR
C
      IBY    = 0
      IBDEL  = IBY   + M
      IBW    = IBDEL + LTHETA
      LENGTH = IBW   + 2*KY
C
      IAY    = IBY   + 1
      IADEL  = IBDEL + 1
      IAW    = IBW   + 1
C
      IF (LENGTH.GT.LWORK) THEN
        WRITE(3,3003)
        STOP
      ENDIF
C
C     INITIALIZE WEIGHTS
C
      DO 20 J=1,KY
        V(J)=0.D0
      DO 20 I=1,KY
        W(I,J)=0.D0
20    CONTINUE
C
C     GET NORMALIZING CONSTANT EH
C
      CALL EHYX(WORK(IAY),X,EH,WORK(IADEL),A,B0,B,R,M,L,
     &  KY,KX,IDEGY,IDEGX,WORK(IAW))
      ROOTEH=DSQRT(EH)
C
C     COMPUTE V
C
      DO 120 IBETA=1,KX
C       GET THE EXPONENT FOR THE MONOMIAL XBETA=X**BETA
        CALL MULTIX(IBETA,BETA)
C       COMPUTE XBETA
        XBETA=1.D0
        DO 100 I=1,M*L
          IF (BETA(I).GT.0) XBETA=XBETA*(X(I)**BETA(I))
100     CONTINUE
      DO 120 IALPHA=1,KY
C       ACCUMULATE THE MONOMIALS TO GET V(J)
        V(IALPHA)=V(IALPHA)+A(KX*(IALPHA-1)+IBETA)*XBETA
120   CONTINUE
C
      DO 125 IALPHA=1,KY
        V(IALPHA)=V(IALPHA)/ROOTEH
125   CONTINUE
C
C     COMPUTE PRODUCT OF GAMMA DENSITY NORMALIZATION FACTORS
C     USE W(I,J) AS TEMPORARY STORAGE
C
      ADJ=0.D0
      DLNB=DLOG(2.D0)
      DO 165 J=1,KY
      IF (V(J).EQ.0.D0) GO TO 165
C       GET LAMBDA_J
        CALL MULTIY(J,LJ)
      DO 160 I=1,J
      IF (V(I).EQ.0.D0) GO TO 160
C       GET LAMBDA_I
        CALL MULTIY(I,LI)
C       COMPUTE SUM OF LI AND LJ
        DO 130 II=1,M
          ALPHA(II)=LI(II)+LJ(II)
130     CONTINUE
        DLG=0.D0
        DO 140 II=1,M
          IA=LI(II)+LJ(II)+1
          DI=IA/2.D0
          DLG=DLG+GAMMLN(DI)-DLNB*(1.D0-DI)
140     CONTINUE
        IF (DLG.GT.ADJ) ADJ=DLG
        W(I,J)=DLG
        W(J,I)=DLG
160   CONTINUE
165   CONTINUE
C
C     COMPUTE W(I,J)
C
      SUM=0.D0
      DO 185 J=1,KY
      IF (V(J).EQ.0.D0) GO TO 185
      DO 180 I=1,KY
      IF (V(I).EQ.0.D0) GO TO 180
        W(I,J)=DABS(V(J))*DABS(V(I))*DEXP(W(I,J)-ADJ)
        SUM=SUM+W(I,J)
180   CONTINUE
185   CONTINUE
      DO 190 J=1,KY
      DO 190 I=1,KY
        W(I,J)=W(I,J)/SUM
190   CONTINUE
C
      RETURN
3002  FORMAT(' ERROR RETURN FROM WHX, M OR M*L EXCEEDS 100')
3003  FORMAT(' ERROR RETURN FROM WHX, LWORK NOT BIG ENOUGH')
      END
C.hr SAMPY
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     SAMPY    1/17/90
C
C     PURPOSE
C     SIMULATE A Y SCALED TO ORIGINAL UNITS GIVEN AN X SCALED IN 
C     ORIGINAL UNITS
C
C     USAGE
C
C     DO 20 I=1,NREPS
C     code to fetch an X in original units
C     CALL SAMPY (0,ISEED,Y,MUSTAT,CSTAT,X,NLOG,C,
C    &            A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,LWORK)
C     DO 10 J=1,NSTEP
C     code to compute R
C     CALL SAMPY (1,ISEED,Y,MUSTAT,CSTAT,X,NLOG,C,
C    &            A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,LWORK)
C     other statements
C10   CONTINUE  
C20   CONTINUE  
C
C     SUBROUTINES CALLED
C     WHX,MULTIY,MULTIX,SETMYX,EHYX,GAMMLN,DRINV,BZ,HZ,SAMPZ,GCHIRV,RAN
C
C     ARGUMENTS
C     ISW -   SWITCH, ISW=0: INITIALIZE X.
C                     ISW=1: COMPUTE Y, NEW X
C             INPUT, INTEGER*4
C     ISEED - SEED FOR THE RANDOM NUMBER GENERATOR RAN
C             INPUT, INTEGER*4
C     Y     - A SIMULATED VALUE OF Y IN ORIGINAL UNITS.
C             OUTPUT, REAL*8
C     MUSTAT- M-VECTOR, SAMPLE MEAN OF THE DATA IN ORIGINAL UNITS.
C             INPUT, REAL*8
C     CSTAT - M BY M MATRIX CONTAINING THE FACTORED SAMPLE VARIANCE OF
C             THE DATA IN ORIGINAL UNITS, VAR = (CSTAT)*(CSTAT)'. CSTAT
C             IS ASSUMED TO HAVE ALL ZEROES BELOW THE DIAGONAL.
C             INPUT, REAL*8
C     X     - M*L VECTOR OF LAGGED VALUES OF Y, STORED END TO END IN
C             TIME SERIAL ORDER. X=(Y'_T-L, ..., Y'_T-1), M*L MUST BE
C             LESS THAT 100
C             INPUT & OUTPUT, REAL*8
C     NLOG  - SWITCH, IF NLOG.EQ.1 THE LOGISTIC TRANSFORM IS NOT USED 
C             INPUT, INTEGER*4                                        
C     C     - LOGISTIC TRANSFORM SCALE FACTOR.
C             INPUT, REAL*8
C     A    -  COEFFICIENTS OF THE POLYNOMIAL, A KX BY KY MATRIX, STORED
C             COLUMNWISE
C             INPUT, REAL*8
C     B0   -  M-VECTOR OF COEFFICIENTS, THE CONSTANT TERM OF THE VAR
C             INPUT, REAL*8
C     B    -  COEFFICIENTS FOR THE REGRESSION ON X, AN M BY M*L MATRIX,
C             STORED COLUMNWISE
C             INPUT, REAL*8
C     R    -  AN M*(M+1)/2-VECTOR CONTAINING THE UPPER TRIANGLE OF THE
C             FACTOR OF THE VARIANCE-COVARIANCE MATRIX OF THE VAR PART
C             STORED COLUMNWISE, EG. R=(R11,R12,R22,R13,R23,R33,...,RMM)
C             SIGMA=R*R', THIS IS NOT THE STANDARD CHOLESKY FACTORIZATION
C             AS IT IS UPPER TRIANGULAR, NOT LOWER, USE SUBROUTINE FACTOR
C             TO COMPUTE IT FROM SIGMA, GIVEN SIGMA.
C             INPUT, REAL*8
C     M    -  DIMENSION OF THE TIME SERIES, IE. LENGTH OF Y_T, MUST BE
C             LESS THAN 100
C             INPUT, INTEGER*4
C     L    -  NUMBER OF LAGS IN THE VAR, M*L MUST BE LESS THAN 100
C             INPUT, INTEGER*4
C     KY   -  NUMBER OF TERMS IN THE POLYNOMIAL IN Y, COMPUTED BY SETMYX,
C             SEE REMARK BELOW.
C             INPUT, INTEGER*4
C     KX   -  NUMBER OF TERMS IN THE POLYNOMIAL IN X, COMPUTED BY SETMYX,
C             SEE REMARK BELOW.
C             INPUT, INTEGER*4
C     IDEGY-  DEGREE OF THE POLYNOMIAL IN Y, SEE REMARK BELOW
C             INPUT, INTEGER*4
C     IDEGX-  DEGREE OF THE POLYNOMIAL IN X, SEE REMARK BELOW
C             INPUT, INTEGER*4
C     WORK -  WORK VECTOR OF LENGTH 
C             M + KX*KY + M + M*M*L +M*(M+1)/2 + 2*KY + KY + KY*KY
C     LWORK-  DECLARED DIMENSION OF WORK IN CALLING PROGRAM.
C             INPUT, INTEGER*4
C
      SUBROUTINE SAMPY(ISW,ISEED,Y,MUSTAT,CSTAT,X,NLOG,C,
     &                  A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,LWORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Y(M),MUSTAT(M),CSTAT(M*M),X(M*L)
      REAL*8 A(KX*KY),B0(M),B(M,M*L),R(M*(M+1)/2)
      REAL*8 WORK(LWORK)
      SAVE
C
C     COMPUTE POINTERS TO WORKSPACE
C
      IBV    = 0
      IBW    = IBV    + KY
      IBWORK = IBW    + KY*KY
      LENGTH = IBWORK + M + KX*KY + M + M*M*L +M*(M+1)/2 + 2*KY
C
      IAV    = IBV    + 1
      IAW    = IBW    + 1
      IAWORK = IBWORK + 1
C
      IF(ISW.EQ.0) GO TO 10
      IF(ISW.GT.0) GO TO 100
C
C     INITIALIZATION
C
10    CONTINUE
C
C     ERROR CHECKS
C
      IF (LENGTH.GT.LWORK) THEN
        WRITE(3,3003)
        STOP
      ENDIF
      IF (L.LE.0) THEN
        WRITE(3,3004) 
        STOP
      ENDIF
C
C     POINTERS FOR THIS BLOCK
C
      IBP = 0
      IBC = IBP + M*M
C
      IAP = IBP+1
      IAC = IBC+1
C
C     LOCATION SCALE SHIFT OF X OFF OF CSTAT AND MUSTAT
C
      DO 35 J=1,M
      DO 35 I=1,J
        WORK(IBC+J*(J-1)/2+I)=CSTAT(M*(J-1)+I)
35    CONTINUE
      CALL DRINV(WORK(IAC),M)
      DO 40 J=1,M
      DO 40 I=1,M
        WORK(IBP+M*(J-1)+I)=0.D0
        IF(I.LE.J) WORK(IBP+M*(J-1)+I)=WORK(IBC+J*(J-1)/2+I)
40    CONTINUE
      DO 50 J=1,L
      DO 50 I=1,M
        WORK(IBC+M*(J-1)+I)=X(M*(J-1)+I)-MUSTAT(I)
50    CONTINUE
      DO 60 J=1,L
        CALL DGMPRD(WORK(IAP),WORK(IBC+M*(J-1)+1),X(M*(J-1)+1),M,M,1)
60    CONTINUE
C
C     LOGISTIC TRANSFORM, IF REQUIRED
C
      IF (NLOG.NE.1) THEN                                               
        DO 70 I=1,M*L                                                  
          ECX=DEXP(C*X(I))                                     
          X(I)=(4.D0/C)*ECX/(1.D0+ECX)-2.D0/C                  
70      CONTINUE                                                       
      ENDIF
      RETURN
C
C     COMPUTE Y
C
100   CONTINUE
C
C     COMPUTE WEIGHTS V & W
C
      CALL WHX(WORK(IAV),WORK(IAW),X,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,
     &         WORK(IAWORK),LWORK-IBWORK)
C
C     COMPUTE Y
C
      CALL SAMPZ(ISEED,WORK(IAWORK),WORK(IAV),WORK(IAW),M,KY,KX)
C
      DO 120 I=1,M
        Y(I)=0.D0
120   CONTINUE
      DO 130 J=1,M
      DO 130 I=1,J
        Y(I)=Y(I)+R(J*(J-1)/2+I)*WORK(IBWORK+J)
130   CONTINUE
      CALL DGMPRD(B,X,WORK(IAWORK),M,M*L,1)
      DO 140 I=1,M
        Y(I)=Y(I)+B0(I)+WORK(IBWORK+I)
140   CONTINUE
C
C     PUSH ELEMENTS OF X BACK ONE LAG
C
      IF (L.GE.2) THEN
        DO 150 I=1,M*(L-1)
          X(I) = X(M+I)
150     CONTINUE
      ENDIF
C
C     ADD NEW, LAGGED, TRANSFORMED Y TO THE END OF X
C
      IF (NLOG.NE.1) THEN                                               
        DO 170 I=1,M                                                  
          ECX=DEXP(C*Y(I))                                     
          X(M*(L-1)+I)=(4.D0/C)*ECX/(1.D0+ECX)-2.D0/C                  
170     CONTINUE                                                       
      ELSE
        DO 180 I=1,M                                                  
          X(M*(L-1)+I)=Y(I)
180     CONTINUE                                                       
      ENDIF
C
C     CONVERT Y TO ORIGINAL UNITS
C
      CALL DGMPRD(CSTAT,Y,WORK(IAWORK),M,M,1)
      DO 190 I=1,M
        Y(I)=WORK(IBWORK+I)+MUSTAT(I)
190   CONTINUE
      RETURN
3003  FORMAT(' ERROR RETURN FROM SAMPY, LWORK NOT BIG ENOUGH')
3004  FORMAT(' ERROR RETURN FROM SAMPY, L .LE. 0')
      END
C.hr SAMPZ
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     SAMPZ      1/9/90
C
C     PURPOSE
C     SAMPLE FROM THE HERMITE CONDITIONAL DENSITY USING OUTPUT FROM WHX.  
C     SEE WHX FOR DOCUMENTATION.
C
C     USAGE
C     CALL SAMPZ(ISEED,Z,V,W,M,KY,KX)
C
C     NOTE
C     THIS ROUTINE RELIES ON WHX FOR ERROR CHECKING.
C
C     SUBROUTINES CALLED
C     GCHIVR, RAN, HZ, BZ
C
      SUBROUTINE SAMPZ(ISEED,Z,V,W,M,KY,KX)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Z(M),V(KY),W(KY,KY)
      INTEGER*4 ALPHA(100),LI(100),LJ(100)
      REAL*4 GCHIRV,RAN,SA
      SAVE
C
C     TOP OF "DO ... UNTIL ... " LOOP
C
100   CONTINUE
C
      U=RAN(ISEED)
      SUM=0.D0
      DO 10 J=1,KY
      DO 10 I=1,KY
        II=I
        JJ=J
        SUM=SUM+W(I,J)
        IF (SUM.GT.U) GO TO 20
10    CONTINUE
20    CONTINUE
C
      CALL MULTIY(II,LI)
      CALL MULTIY(JJ,LJ)
C
C       COMPUTE SUM OF LI AND LJ
        DO 50 I=1,M
          ALPHA(I)=LI(I)+LJ(I)
50      CONTINUE
C
C       GENERATE M CHI'S AND RANDOMLY CHANGE SIGN
        DO 60 I=1,M
          SA=ALPHA(I)+1
          Z(I)=GCHIRV(SA,ISEED)
          IF (RAN(ISEED).GT.0.5E0) Z(I)=-Z(I)
60      CONTINUE
C       
C     UNTIL CONDITION
C
C       ACCEPT OR REJECT
        U=RAN(ISEED)
        IF ( U*BZ(Z,V,M,KY,KX) .GT. HZ(Z,V,M,KY,KX) ) GO TO 100
C
C     BOTTOM OF "DO ... UNTIL ... " LOOP
C
      RETURN
      END
C.hr BZ
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     BZ       1/9/90
C
C     PURPOSE
C     MAJORIZE HZ USING OUTPUT FROM WHX.  SEE WHX FOR DOCUMENTATION.
C     REAL*8 FUNCTION.
C
C     USAGE
C     A=BZ(Z,V,M,KY,KX)
C
C     THIS ROUTINE RELIES ON WHX FOR ERROR CHECKING.
C
      REAL*8 FUNCTION BZ(Z,V,M,KY,KX)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Z(M),V(KY)
      INTEGER*4 ALPHA(100)
      SAVE
C
      PI=3.1415926535897932D0
C
      SUM=0.D0
      DO 20 I=1,KY
        CALL MULTIY(I,ALPHA)
        ZALPHA=1.D0
        DO 10 J=1,M
          IF(ALPHA(J).GT.0) ZALPHA=ZALPHA*DABS(Z(J))**ALPHA(J)
10      CONTINUE
        SUM=SUM+DABS(V(I))*ZALPHA
20    CONTINUE
C
      SSZ=0.D0
      DO 30 J=1,M
        SSZ=SSZ+Z(J)**2
30    CONTINUE
C
      BZ=(SUM**2)*DEXP(-0.5D0*SSZ)/DSQRT(2.D0*PI)**M
C
      RETURN
      END
C.hr HZ
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     HZ       1/9/90
C
C     PURPOSE
C     COMPUTE THE HERMITE CONDITIONAL DENSITY USING OUTPUT FROM WHX.  
C     SEE WHX FOR DOCUMENTATION.  REAL*8 FUNCTION.
C
C     USAGE
C     A=HZ(Z,V,M,KY,KX)
C
C     THIS ROUTINE RELIES ON WHX FOR ERROR CHECKING.
C
      REAL*8 FUNCTION HZ(Z,V,M,KY,KX)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Z(M),V(KY)
      INTEGER*4 ALPHA(100)
      SAVE
C
      PI=3.1415926535897932D0
C
      SUM=0.D0
      DO 20 I=1,KY
        CALL MULTIY(I,ALPHA)
        ZALPHA=1.D0
        DO 10 J=1,M
          IF(ALPHA(J).GT.0) ZALPHA=ZALPHA*Z(J)**ALPHA(J)
10      CONTINUE
        SUM=SUM+V(I)*ZALPHA
20    CONTINUE
C
      SSZ=0.D0
      DO 30 J=1,M
        SSZ=SSZ+Z(J)**2
30    CONTINUE
C
      HZ=(SUM**2)*DEXP(-0.5D0*SSZ)/DSQRT(2.D0*PI)**M
C
      RETURN
      END
C.hr GCHIRV
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
      REAL FUNCTION GCHIRV(A,ISEED)
C
C  ALGORITHM TO GENERATE RANDOM VARIABLES WITH THE CHI DISTRIBUTION
C  WITH A DEGREES OF FREEDOM,  FOR A GREATER THAN OR EQUAL TO ONE
C
C  J F MONAHAN (MAY, 1986) DEPT OF STATISTICS, NCSU, RALEIGH, NC USA
C  MODIFIED 1/24/90 BY A R GALLANT TO SET SEED DIFFERENTLY
C
      REAL A
      REAL ALPHA,ALPHM1,BETA,U,V,VMIN,VDIF,Z,ZZ,RNUM,W,S,EMHLF,VMAXU,
     *     RSQRT2,EMHLF4,EQTRT2,C
      SAVE
C
      DATA EMHLF/ .6065307 /, VMAXU/ .8577639 /, RSQRT2/ .7071068 /
      DATA EMHLF4/ .1516327 /, EQTRT2/ 2.568051 /, C/ 1.036961 /
C
      DATA ALPHA/ 0. /
C                     IS THIS ALPHA THE SAME AS THE LAST ONE?
      IF( A .EQ. ALPHA ) GO TO 1
C                     DO A LITTLE SETUP
      ALPHA = A
      ALPHM1 = ALPHA - 1.
      BETA = SQRT( ALPHM1 )
C                     GET DIMENSIONS OF BOX
      VMAXU = EMHLF * ( RSQRT2 + BETA )/( .5 + BETA )
      VMIN =  -BETA
      IF( BETA .GT. 0.483643 ) VMIN = EMHLF4/ALPHA - EMHLF
      VDIF = VMAXU - VMIN
C                     START ( AND RESTART ) ALGORITHM HERE
  1   CONTINUE
      U = RAN(ISEED)
      V = VDIF*RAN(ISEED) + VMIN
      Z = V / U
C                     DO SOME QUICK REJECT CHECKS FIRST
      IF( Z .LE. - BETA ) GO TO 1
      ZZ = Z * Z
      RNUM = 2.5 - ZZ
      IF( Z .LT. 0. ) RNUM = RNUM + ZZ * Z / ( 3. * (Z + BETA ) )
C                     DO QUICK INNER CHECK
      IF( U .LT. RNUM/EQTRT2 ) GO TO 9
      IF( ZZ .GE. C / U + 1.4 ) GO TO 1
C                     ABOVE WAS KNUTH'S NORMAL OUTER CHECK
      W = 2. * ALOG( U )
C                     NOW THE REAL CHECK
      S = - ( ZZ / 2. + Z * BETA )
      IF( BETA .GT. 0. ) S = ALPHM1*ALOG(1.+Z/BETA) + S
      IF( W .GT. S ) GO TO 1
C                     SUCCESS HERE -- TRANSFORM AND GO
  9   GCHIRV = Z + BETA
      RETURN
      END
C.hr GAMMLN
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C 
C     THIS IS A MODIFICATION OF THE ROUTINE OF THE SAME NAME IN PRESS, 
C     WILLIAM H., BRIAN P. FLANNERY, SAUL A. TEUKOLSKY, AND WILLIAM T.
C     VETTERLING (1986), NUMERICAL RECIPES, CAMBRIDGE UNIVERSITY PRESS,
C     CAMBRIDGE, U.K., P. 157.
C
      REAL*8 FUNCTION GAMMLN(XX)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER
      SAVE
      DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0,
     *    -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
      DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/
      X=XX-ONE
      TMP=X+FPF
      TMP=(X+HALF)*DLOG(TMP)-TMP
      SER=ONE
      DO 11 J=1,6
        X=X+ONE
        SER=SER+COF(J)/X
11    CONTINUE
      GAMMLN=TMP+DLOG(STP*SER)
      RETURN
      END
C.hr RAN
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
      REAL FUNCTION RAN(IX)
C  UNIFORM PSEUDORANDOM NUMBER GENERATOR
C  FORTRAN VERSION OF LEWIS, GOODMAN, MILLER
C  SCHRAGE,  ACM TOMS V.5 (1979) P132
C  MODIFIED 1/24/90 BY A R GALLANT TO SET SEED DIFFERENTLY
      INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K
      SAVE
      DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/
      XHI=IX/B16
      XALO=(IX-XHI*B16)*A
      LEFTLO=XALO/B16
      FHI=XHI*A+LEFTLO
      K=FHI/B15
      IX=(((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16)+K
      IF(IX.LT.0) IX=IX+P
      RAN=FLOAT(IX)*4.656612875E-10
      RETURN
      END
