C.hr spcode1.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     snpcode1.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 SETPRD
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     SETPRD     12/5/88
C
C     PURPOSE
C     COMPUTE THE INTEGER*4 ARRAYS IY, IDY, IX WHICH ALLOW EFFICIENT
C     COMPUTATION OF A MULTIVARIATE POLYNOMIAL AND ITS DERIVATIVES IN
C     LNHYX.   THEY ARE STORED IN THE COMMON BLOCK
C     COMMON /ZSETPR/ IY,IDY,IX
C     IY   - IY(2*(I-1)+1) POINTS TO A MONOMIAL WHICH WHEN MULTIPLIED BY
C            Y(IY(2*(I-1)+2)) WILL PRODUCE THE I-TH MONOMIAL IN Y.
C     IDY  - IDY(M*(I-1)+K) POINTS TO A MONOMIAL WHICH WHEN MULTIPLIED
C            BY THE K-TH ELEMENT OF THE I-TH MULTI-INDEX WILL PRODUCE
C            THE K-TH PARTIAL DERIVATIVE OF THE I-TH MONOMIAL IN Y.
C     IX   - IX(2*(I-1)+1) POINTS TO A MONOMIAL WHICH WHEN MULTIPLIED BY
C            X(IX(2*(I-1)+2)) WILL PRODUCE THE I-TH MONOMIAL IN X.
C
C     USAGE
C     CALL SETPRD(M,L,KY,KX,IDEGY,IDEGX)
C
C     WARNING
C     IF ZZLIM IS CHANGED HERE IT MUST BE CHANGED IN SUBROUTINE LNHYX
C     AND EHYX AS WELL.
C
C     SUBROUTINES CALLED
C     MULTIY, MULTIX, SETMYX
C
C     ARGUMENTS
C     M    - DIMENSION OF THE TIME SERIES, I.E. 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
C     REMARK
C     TYPICAL CODE TO COMPUTE THE POLYNOMIAL P AND ITS GRADIENT DP IS:
C     WY(1)=1
C     DO 10 I=2,KY
C10   WY(I)=WY(IY(2*(I-1)+1))*Y(IY(2*(I-1)+2))
C     WX(1)=1
C     DO 20 I=2,KX
C20   WX(I)=WX(IX(2*(I-1)+1))*X(IX(2*(I-1)+2))
C     P=0.D0
C     DO 30 K=1,M
C30   DP(K)=0.D0
C     DO 40 J=1,KY
C     CALL MULTIY(J,ALPHA)
C     DO 40 I=1,KX
C     P=P+A(I,J)*WX(I)*WY(J)
C     DO 40 K=1,M
C40   DP(K)=DP(K)+A(I,J)*WX(I)*ALPHA(K)*WY(IDY(M*(I-1)+K))
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 SETPRD(M,L,KY,KX,IDEGY,IDEGX)
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 ALPHA(100),BETA(100),GAMMA(100)
C-----------------------------------------------------------------------
C     DECLARATIONS SHARED BETWEEN SETPRD, PZX, EHYX, AND LNHYX
C-----------------------------------------------------------------------
      INTEGER*4  ZZLIM
      PARAMETER (ZZLIM=500)
      INTEGER*4 IY(ZZLIM),IDY(ZZLIM),IX(ZZLIM)
      COMMON /ZSETPR/ IY,IDY,IX
      SAVE
C-----------------------------------------------------------------------
C
C     LOGICAL:
C     INTEGER*4 IY(2,KY),IDY(M,KY),IX(2,KX)
C
C     ERROR TRAPS
C
      NOUT=3
      IF((2*KY.GT.ZZLIM).OR.(M*KY.GT.ZZLIM).OR.(2*KX.GT.ZZLIM)) THEN
      WRITE(NOUT,3001)
      STOP
      ENDIF
C
C     INITIALIZATION
C
      DO 10 IALPHA=1,KY
      IY(2*(IALPHA-1)+1)=0
      IY(2*(IALPHA-1)+2)=0
      DO 10 INDEX=1,M
10    IDY(M*(IALPHA-1)+INDEX)=1
      DO 20 IBETA=1,KX
      IX(2*(IBETA-1)+1)=0
20    IX(2*(IBETA-1)+2)=0

C
C     OUTER LOOP OVER ALL MULTI-INDEXES ALPHA
C
      DO 190 IALPHA=1,KY
      CALL MULTIY(IALPHA,ALPHA)
C
C     LOOK FOR A PREDECESSOR TO THE MULTI-INDEX ALPHA
C
      DO 180 ITEST=1,IALPHA
      CALL MULTIY(ITEST,GAMMA)
C
      DO 130 INDEX=1,M
      ALPHA(INDEX)=ALPHA(INDEX)-1
C
      DO 120 K=1,M
120   IF(ALPHA(K).NE.GAMMA(K)) GO TO 130

C     FALL THROUGH LOOP 120 MEANS THAT ITEST IS A PREDECESSOR
      IY(2*(IALPHA-1)+1)=ITEST
      IY(2*(IALPHA-1)+2)=INDEX
      IDY(M*(IALPHA-1)+INDEX)=ITEST

130   ALPHA(INDEX)=ALPHA(INDEX)+1
C
180   CONTINUE
C
190   CONTINUE
C
C
C     OUTER LOOP OVER ALL MULTI-INDEXES BETA
C
      DO 290 IBETA=1,KX
      CALL MULTIX(IBETA,BETA)
C
C     LOOK FOR A PREDECESSOR TO THE MULTI-INDEX BETA
C
      DO 280 ITEST=1,IBETA
      CALL MULTIX(ITEST,GAMMA)
C
      DO 230 INDEX=1,M*L
      BETA(INDEX)=BETA(INDEX)-1
C
      DO 220 K=1,M*L
220   IF(BETA(K).NE.GAMMA(K)) GO TO 230

C     FALL THROUGH LOOP 220 MEANS THAT ITEST IS A PREDECESSOR
      IX(2*(IBETA-1)+1)=ITEST
      IX(2*(IBETA-1)+2)=INDEX

230   BETA(INDEX)=BETA(INDEX)+1
C
280   CONTINUE
C
290   CONTINUE
C
      RETURN
3001  FORMAT
     &(' ERROR RETURN FROM SETPRD, 2*KY, M*KY, OR 2*KX EXCEEDS ZZLIM')
      END
C.hr SETUPM
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
      SUBROUTINE SETUPM(IDEGY,M,KY)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C     NOTE:    ZZMNTS IS MAXIMUM LENGTH OF ARRAY OF GAUSSIAN MOMENTS.
C     WARNING: IF ZZMNTS IS CHANGED HERE IT MUST BE CHANGED IN
C              SUBROUTINES EHYX AND EYLX AS WELL
C
      INTEGER     ZZMNTS
      PARAMETER ( ZZMNTS = 3000 )
      INTEGER*4 ALPHA(100),GAMMA(100),OMEGA(100)
      REAL*8 MOMENT(34),MNTS(ZZMNTS)
      COMMON /ZMNTS/MNTS
      SAVE
      DATA MOMENT/
     &  1.00000000000000D+00,  0.00000000000000D+00,
     &  1.00000000000000D+00,  0.00000000000000D+00,
     &  3.00000000000000D+00,  0.00000000000000D+00,
     &  1.50000000000000D+01,  0.00000000000000D+00,
     &  1.05000000000000D+02,  0.00000000000000D+00,
     &  9.45000000000000D+02,  0.00000000000000D+00,
     &  1.03950000000000D+04,  0.00000000000000D+00,
     &  1.35135000000000D+05,  0.00000000000000D+00,
     &  2.02702500000000D+06,  0.00000000000000D+00,
     &  3.44594250000000D+07,  0.00000000000000D+00,
     &  6.54729075000000D+08,  0.00000000000000D+00,
     &  1.37493105750000D+10,  0.00000000000000D+00,
     &  3.16234143225000D+11,  0.00000000000000D+00,
     &  7.90585358062500D+12,  0.00000000000000D+00,
     &  2.13458046676875D+14,  0.00000000000000D+00,
     &  6.19028335362938D+15,  0.00000000000000D+00,
     &  1.91898783962511D+17,  0.00000000000000D+00/
C
C     ERROR TRAP
C
      IF ((M.GT.100).OR.(KY*KY.GT.ZZMNTS).OR.(IDEGY.GT.16)) THEN
      WRITE(3,3002)
      RETURN
      ENDIF
C
C     GENERATE MOMENTS AND STORE IN MATRIX MNTS
C
      DO 30 IGAMMA=1,KY
      CALL MULTIY(IGAMMA,GAMMA)
C
      DO 30 IALPHA=1,KY
      CALL MULTIY(IALPHA,ALPHA)
C
      DO 10 I=1,M
10    OMEGA(I)=ALPHA(I)+GAMMA(I)
C
      DINT=1.D0
      DO 20 I=1,M
20    DINT=DINT*MOMENT(OMEGA(I)+1)
C
      MNTS(KY*(IGAMMA-1)+IALPHA)=DINT
30    CONTINUE
C
      RETURN
3002  FORMAT(' ERROR RETURN FROM SETUPM:'/
     &' M EXCEEDS 100, KY*KY EXCEEDS ZZMNTS, OR IDEGY EXCEEDS 16')
      END
C.hr EHYX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C
C     EHYX     12/7/88
C
C     PURPOSE
C     COMPUTE THE INTEGRAL OF AN UNNORMALIZED HERMITE CONDITIONAL
C     DENSITY,
C
C     EH = I..I {SUM_J  A(X,J)*[INV(R)(Y-B0-B*X)]**ALPHA(J)}**2
C                * N(Y | B0 + B*X, SIGMA) DY_1...DY_M
C
C     AND ITS DERIVATIVE,
C
C     DELEH = (D/D_THETA)EH
C
C     WHERE
C     A(X,J) = SUM_I A(I,J) X**BETA(I), I = 1, ..., KX, J = 1, ..., KY
C     THETA = (VEC'(A)|B0'|VEC'(B)|R11,R12,R22,R13,R23,R33,...,RMM)'
C     SIGMA = R*R'
C
C     NOTE
C     R IS UPPER TRIANGULAR, STORED COLUMNWISE IN COMPRESSED FORM,
C     THE FACTOR N(Y | B0 + B*X, SIGMA) IS A VECTOR AUTOREGRESSION
C     (VAR), THE DEPENDENCE OF THE POLYNOMIAL IN Y ON X ACCOMMODATES
C     CONDITIONAL HETEROSKEDASTICITY AND NONLINEARITY
C
C     USAGE
C     CALL EHYX(Y,X,EH,DELEH,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
C
C     WARNING
C     IF ZZLIM IS CHANGED HERE IT MUST BE CHANGED IN SUBROUTINE SETPRD
C     AS WELL.
C
C     SUBROUTINES CALLED
C     SETMYX,SETUPM,SETPRD
C
C     ARGUMENTS
C     Y    - M-VECTOR USED AS WORKSPACE, M MUST BE LESS THAN 100
C            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     EH   - COMPUTED VALUE OF THE INTEGRAL OF THE DENSITY
C            OUTPUT, REAL*8
C     DELEH- COMPUTED GRADIENT OF EH WITH RESPECT TO THETA, A VECTOR OF
C            LENGTH KX*KY+M+M*M*L+M*(M+1)/2
C            OUTPUT, 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, MUST BE LESS THAN 16, SEE
C            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 2*KY
C            REAL*8
C
C     REMARK
C     THE CALLS:
C     CALL SETMYX(M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX,KY,KX,IDUMMY,IWORK)
C     CALL SETUPM(IDEGY,M,KY)
C     CALL SETPRD(M,L,KY,KX,IDEGY,IDEGX)
C     MUST PRECEDE USAGE OF THIS ROUTINE, IN THAT ORDER
C
      SUBROUTINE EHYX(Y,X,EH,DELEH,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C     NOTE:    ZZMNTS IS MAXIMUM LENGTH OF ARRAY OF GAUSSIAN MOMENTS.
C     WARNING: IF ZZMNTS IS CHANGED HERE IT MUST BE CHANGED IN
C              SUBROUTINES SETUPM AND EYLX AS WELL
C
      INTEGER     ZZMNTS
      PARAMETER ( ZZMNTS = 3000 )
      REAL*8 Y(1),X(1),EH,DELEH(1)
      REAL*8 A(1),B0(1),B(1),R(1)
      REAL*8 WORK(1)
      REAL*8 MNTS(ZZMNTS)
      COMMON /ZMNTS/MNTS
C
C     LOGICAL:
C     REAL*8 Y(M),X(M*L),DELEH(KX*KY+M+M*M*L+M*(M+1)/2)
C     REAL*8 A(KX*KY),B0(M),B(M,M*L),R(M*(M+1)/2)
C     OR:
C     REAL*8 Y(M),X(M,L),A(KX,KY)
C
C-----------------------------------------------------------------------
C     DECLARATIONS SHARED BETWEEN SETPRD, PZX, EHYX, AND LNHYX 
C-----------------------------------------------------------------------
      INTEGER*4  ZZLIM
      PARAMETER (ZZLIM=500)
      INTEGER*4 IY(ZZLIM),IDY(ZZLIM),IX(ZZLIM)
      COMMON /ZSETPR/ IY,IDY,IX
C-----------------------------------------------------------------------
C
C     LOGICAL:
C     INTEGER*4 IY(2,KY),IDY(M,KY),IX(2,KX)
C-----------------------------------------------------------------------
      SAVE
C
C     ERROR TRAP
C
      IF ((M.GT.100).OR.(M*L.GT.100).OR.(KY*KY.GT.ZZMNTS)
     &.OR.(IDEGY.GT.16)) THEN
      WRITE(3,3002)
      RETURN
      ENDIF
C
C     COMPUTE POINTERS TO STACKED PARAMETER VECTOR, USED TO POINT TO
C     ELEMENTS OF DELEH=(D/DTHETA)EH(Y|X,THETA)
C
      IBA  = 0
      IBB0 = IBA   +KX*KY
      IBB  = IBB0  + M
      IBR  = IBB   + M*M*L
      LTHETA=IBR+M*(M+1)/2
C
C     IAA  = IBA +1
C     IAB0 = IBB0+1
C     IAB  = IBB +1
C     IAR  = IBR +1
C
C     COMPUTE POINTERS TO WORK VECTOR
C
      IBAX = 0
      IBAXM= IBAX +KY
C     LWORK= IBAXM+KY
C
C     IAAX = IBAX +1
C     IAAXM= IBAXM+1
C
C     COMPUTE LENGTHS OF VARIOUS ARRAYS
C
      ML=M*L
C
C     INITIALIZE DELEH
C
      DO 10 I=1,LTHETA
 10   DELEH(IBA+I)=0.D0
C
C     STORE THE MONOMIALS FOR X IN FIRST COLUMN OF DELEH
C
      DELEH(IBA+1)=1.D0
      IF(KX.LE.1) GO TO 15
      DO 14 IBETA=2,KX
      IX1=IX(2*(IBETA-1)+1)
      IX2=IX(2*(IBETA-1)+2)
14    DELEH(IBA+IBETA)=DELEH(IBA+IX1)*X(IX2)
15    CONTINUE
C
C     COMPUTE A(X)
C
      DO 20 IALPHA=1,KY
20    WORK(IBAX+IALPHA)=0.D0
      DO 50 IBETA =1,KX
      XBETA=DELEH(IBA+IBETA)
      DO 40 IALPHA=1,KY
40    WORK(IBAX+IALPHA)=WORK(IBAX+IALPHA)+A(KX*(IALPHA-1)+IBETA)*XBETA
50    CONTINUE
C
C     COMPUTE A'(X)M
C
      DO 60 IGAMMA=1,KY
60    WORK(IBAXM+IGAMMA)=0.D0
      DO 70 IGAMMA=1,KY
      DO 70 IALPHA=1,KY
70    WORK(IBAXM+IGAMMA)=WORK(IBAXM+IGAMMA)
     &   +WORK(IBAX+IALPHA)*MNTS(KY*(IGAMMA-1)+IALPHA)
C
C     COMPUTE DELEH
C
      DO 90 IALPHA=1,KY
      JALPHA=KY-(IALPHA-1)
      DO 90 IBETA=1,KX
      XBETA=DELEH(IBA+IBETA)
90    DELEH(IBA+KX*(JALPHA-1)+IBETA)=2.D0*XBETA*WORK(IBAXM+JALPHA)
C
C     COMPUTE EH
C
      SUM=0.D0
      DO 110 IALPHA=1,KY
110   SUM=SUM+WORK(IBAXM+IALPHA)*WORK(IBAX+IALPHA)
C
      EH=SUM
C
      RETURN
 3002 FORMAT(' ERROR RETURN FROM EHYX:'/
     &' M OR M*L EXCEEDS 100, KY*KY EXCEEDS ZZMNTS, OR IDEGY'/
     &' EXCEEDS 16')
      END
C.hr EYLX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     EYLX     11/17/87
C
C     PURPOSE
C     COMPUTE THE LAMBDA'TH MOMENT OF  INV(R)(Y-B0-B*X)  OF AN
C     UNNORMALIZED HERMITE CONDITIONAL DENSITY
C
C     E{INV(R)(Y-B0-B*X)[**LAMBDA|X,THETA}
C        = I...I [INV(R)(Y-B0-B*X)]**LAMBDA * H(Y|X,THETA)
C
C     AND ITS DERIVATIVE
C
C     (D/D_THETA)E{INV(R)(Y-B0-B*X)[**LAMBDA|X,THETA}
C
C     WHERE
C     H(Y|X,THETA)= [SUM_J  A(X,J)*[INV(R)(Y-B0-B*X)]**ALPHA(I)]**2
C                     * N(Y | B0 + B*X, SIGMA)
C     A(X,J) = SUM_I A(I,J) X**BETA(I), I = 1, ..., KX, J = 1, ..., KY
C     THETA = (VEC'(A)|B0|VEC'(B)|R11,R12,R22,R13,R23,R33,...,RMM)'
C     SIGMA = R*R'
C
C     NOTE
C     R IS UPPER TRIANGULAR, STORED COLUMNWISE IN COMPRESSED FORM,
C     THE FACTOR N(Y | B0 + B*X, SIGMA) IS A VECTOR AUTOREGRESSION
C     (VAR), THE DEPENDENCE OF THE POLYNOMIAL IN Y ON X ACCOMMODATES
C     CONDITIONAL HETEROSKEDASTICITY AND NONLINEARITY
C
C     USAGE
C     CALL EYLX(LAMBDA,Y,X,EYL,DEL,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
C
C     SUBROUTINES CALLED
C     MULTIY,MULTIX,SETMYX,SETUPM,EHYX
C
C     REMARK:
C     THE FOLLOWING USAGE WILL COMPUTE THE MOMENTS OF THE NORMALIZED
C     DENSITY:
C     CALL EYLX(LAMBDA,Y,X,EYL,DEL,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
C     CALL EHYX(Y,X,EH,DELH,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
C     EYL=EYL/EH
C     DO 10 I=1,KX*KY+M+M*M*L+M*(M+1)/2
C10   DEL(I)=DEL(I)/EH-(EYL/EH)*DELH(I)
C
C     ARGUMENTS
C     LAMBDA- A MULTI-INDEX THAT DEFINES THE MOMENT.  I.E. LAMBDA IS
C             AN M-VECTOR OF NONNEGATIVE INTEGERS, THE I'TH ELEMENT
C             OF (Y-B0-B*X) IS TO BE RAISED TO THE POWER OF THE I'TH
C             ELEMENT OF LAMBDA.  M MUST BE LESS THAN 100.
C             INPUT, INTEGER*4
C     Y     - M-VECTOR USED AS WORKSPACE, M MUST BE LESS THAN 100
C             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     EYL   - COMPUTED LAMBDA'TH MOMENT ABOUT B0+B*X
C             OUTPUT, REAL*8
C     DEL   - COMPUTED GRADIENT OF EYL WITH RESPECT TO THETA, A VECTOR
C             OF LENGTH KX*KY+M+M*M*L+M*(M+1)/2
C             OUTPUT, 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             CHOLESKY FACTOR OF THE VARIANCE-COVARIANCE MATRIX OF THE
C             VAR STORED COLUMNWISE, R=(R11,R12,R22,R13,R23,R33,...,RMM)
C             SIGMA=R*R'
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
C             SETMYX, SEE REMARK BELOW.
C             INPUT, INTEGER*4
C     KX    - NUMBER OF TERMS IN THE POLYNOMIAL IN X, COMPUTED BY
C             SETMYX, SEE REMARK BELOW.
C             INPUT, INTEGER*4
C     IDEGY - DEGREE OF THE POYNOMIAL IN Y, MUST BE LESS THAN 16, SEE
C             REMARK BELOW
C             INPUT, INTEGER*4
C     IDEGX - DEGREE OF THE POYNOMIAL IN X, SEE REMARK BELOW
C             INPUT, INTEGER*4
C     WORK  - WORK VECTOR OF LENGTH 2*KY
C             REAL*8
C
C     REMARK
C     THE CALLS:
C     CALL SETMYX(M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX,KY,KX,IDUMMY,IWORK)
C     CALL SETUPM(IDEGY,M,KY)
C     MUST PRECEDE USAGE OF THIS ROUTINE, IN THAT ORDER
C
      SUBROUTINE
     &   EYLX(LAMBDA,Y,X,EYL,DEL,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
C     NOTE:    ZZMNTS IS MAXIMUM LENGTH OF ARRAY OF GAUSSIAN MOMENTS.
C     WARNING: IF ZZMNTS IS CHANGED HERE IT MUST BE CHANGED IN
C              SUBROUTINE SETUPM AND EHYX AS WELL
C
      INTEGER     ZZMNTS
      PARAMETER ( ZZMNTS = 3000 )
      REAL*8 Y(1),X(1),EYL,DEL(1)
      REAL*8 A(1),B0(1),B(1),R(1)
      REAL*8 WORK(1)
      INTEGER*4 ALPHA(100),BETA(100),GAMMA(100),OMEGA(100)
      INTEGER*4 LAMBDA(100),COPY(100),ENTRY
      REAL*8 MOMENT(34),LMNTS(ZZMNTS)
      SAVE
      DATA ENTRY/0/
      DATA MOMENT/
     &  1.00000000000000D+00,  0.00000000000000D+00,
     &  1.00000000000000D+00,  0.00000000000000D+00,
     &  3.00000000000000D+00,  0.00000000000000D+00,
     &  1.50000000000000D+01,  0.00000000000000D+00,
     &  1.05000000000000D+02,  0.00000000000000D+00,
     &  9.45000000000000D+02,  0.00000000000000D+00,
     &  1.03950000000000D+04,  0.00000000000000D+00,
     &  1.35135000000000D+05,  0.00000000000000D+00,
     &  2.02702500000000D+06,  0.00000000000000D+00,
     &  3.44594250000000D+07,  0.00000000000000D+00,
     &  6.54729075000000D+08,  0.00000000000000D+00,
     &  1.37493105750000D+10,  0.00000000000000D+00,
     &  3.16234143225000D+11,  0.00000000000000D+00,
     &  7.90585358062500D+12,  0.00000000000000D+00,
     &  2.13458046676875D+14,  0.00000000000000D+00,
     &  6.19028335362938D+15,  0.00000000000000D+00,
     &  1.91898783962511D+17,  0.00000000000000D+00/
C
C     LOGICAL:
C     REAL*8 Y(M),X(M*L),DEL(KX*KY+M+M*M*L+M*(M+1)/2)
C     REAL*8 A(KX*KY),B0(M),B(M,M*L),R(M*(M+1)/2)
C     INTEGER*4 ALPHA(M),BETA(M*L),GAMMA(M)
C     OR:
C     REAL*8 Y(M),X(M,L),A(KX,KY)
C
C
C     ERROR TRAP
C
      IF ((M.GT.100).OR.(M*L.GT.100).OR.(KY*KY.GT.ZZMNTS)
     &.OR.(IDEGY.GT.16)) THEN
      WRITE(3,3002)
      RETURN
      ENDIF
C
C     COMPUTE POINTERS TO STACKED PARAMETER VECTOR, USED TO POINT TO
C     ELEMENTS OF DEL
C
      DO 5 I=1,M
 5    IF (COPY(I).NE.LAMBDA(I)) ENTRY=0
C
      IF (ENTRY.EQ.0) THEN
C
C     GENERATE MOMENTS AND STORE IN MATRIX LMNTS
C
      DO 30 IGAMMA=1,KY
      CALL MULTIY(IGAMMA,GAMMA)
C
      DO 30 IALPHA=1,KY
      CALL MULTIY(IALPHA,ALPHA)
C
      DO 10 I=1,M
10    OMEGA(I)=ALPHA(I)+GAMMA(I)+LAMBDA(I)
C
      DINT=1.D0
      DO 20 I=1,M
20    DINT=DINT*MOMENT(OMEGA(I)+1)
C
      LMNTS(KY*(IGAMMA-1)+IALPHA)=DINT
30    CONTINUE
C
      ENDIF
C
      IBA  = 0
      IBB0 = IBA   +KX*KY
      IBB  = IBB0  + M
      IBR  = IBB   + M*M*L
      LTHETA=IBR+M*(M+1)/2
C
C     IAA  = IBA +1
C     IAB0 = IBB0+1
C     IAB  = IBB +1
C     IAR  = IBR +1
C
C     COMPUTE POINTERS TO WORK VECTOR
C
      IBAX = 0
      IBAXM= IBAX +KY
C     LWORK= IBAXM+KY
C
C     IAAX = IBAX +1
C     IAAXM= IBAXM+1
C
C     COMPUTE LENGTHS OF VARIOUS ARRAYS
C
      ML=M*L
C
C     INITIALIZE DEL TO ZERO
C
      DO 110 I=1,LTHETA
 110  DEL(IBA+I)=0.D0
C
C     COMPUTE A(X)
C
      DO 120 IALPHA=1,KY
 120  WORK(IBAX+IALPHA)=0.D0
      DO 150 IBETA =1,KX
      CALL MULTIX(IBETA,BETA)
      XBETA=1.D0
      DO 130 I=1,ML
 130  IF (BETA(I).GT.0) XBETA=XBETA*(X(I)**BETA(I))
      DO 140 IALPHA=1,KY
 140  WORK(IBAX+IALPHA)=WORK(IBAX+IALPHA)+A(KX*(IALPHA-1)+IBETA)*XBETA
 150  CONTINUE
C
C     COMPUTE A'(X)M
C
      DO 160 IGAMMA=1,KY
 160  WORK(IBAXM+IGAMMA)=0.D0
      DO 170 IGAMMA=1,KY
      DO 170 IALPHA=1,KY
 170  WORK(IBAXM+IGAMMA)=WORK(IBAXM+IGAMMA)
     &   +WORK(IBAX+IALPHA)*LMNTS(KY*(IGAMMA-1)+IALPHA)
C
C     COMPUTE DEL
C
      SUM=0.D0
      DO 200 IBETA =1,KX
      CALL MULTIX(IBETA,BETA)
      XBETA=1.D0
      DO 180 I=1,ML
 180  IF (BETA(I).GT.0) XBETA=XBETA*(X(I)**BETA(I))
C
      DO 190 IALPHA=1,KY
C
 190  DEL(IBA+KX*(IALPHA-1)+IBETA)=2.D0*XBETA*WORK(IBAXM+IALPHA)
C
 200  CONTINUE
C
C     COMPUTE EYL
C
      SUM=0.D0
      DO 210 IALPHA=1,KY
 210  SUM=SUM+WORK(IBAXM+IALPHA)*WORK(IBAX+IALPHA)
C
      EYL=SUM
C
      ENTRY=1
      DO 220 I=1,M
 220  COPY(I)=LAMBDA(I)
C
      RETURN
 3002 FORMAT(' ERROR RETURN FROM EYLX:'/
     &' M OR M*L EXCEEDS 100, KY*KY EXCEEDS ZZMNTS, OR IDEGY'/
     &' EXCEEDS 16')
      END
C.hr LNHYX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     LNHYX      12/5/88
C
C     PURPOSE
C     COMPUTE AN UNNORMALIZED LOG HERMITE CONDITIONAL DENSITY,
C
C     H(Y|X,THETA)
C        = LN {SUM_J  A(X,J)*[INV(R)*(Y-B0-B*X)]**ALPHA(J)}**2
C                  * N(Y | B0 + B*X, SIGMA)
C
C     AND ITS DERIVATIVE,
C
C     DELH(Y|X,THETA) = (D/D_THETA) LN H(Y|X,THETA)
C
C     WHERE
C     A(X,J) = SUM_I A(I,J) X**BETA(I), I = 1, ..., KX, J = 1, ..., KY
C     THETA = (VEC'(A)|B0'|VEC'(B)|R11,R12,R22,R13,R23,R33,...,RMM)'
C     SIGMA = R*R'
C
C     NOTE
C     R IS UPPER TRIANGULAR, STORED COLUMNWISE IN COMPRESSED
C     FORM,THE FACTOR N(Y | B0 + B*X, SIGMA) IS A VECTOR AUTOREGRESSION
C     (VAR), THE DEPENDENCE OF THE POLYNOMIAL IN Y ON X ACCOMMODATES
C     CONDITIONAL HETEROSKEDASTICITY AND NONLINEARITY
C
C     USAGE
C     CALL LNYX(Y,X,H,DELH,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
C
C     WARNING
C     IF ZZLIM IS CHANGED HERE IT MUST BE CHANGED IN SUBROUTINE SETPRD
C     AND EHYX AS WELL.
C     IF ZZMID IS CHANGED HERE IT MUST BE CHANGED IN SETMYX, MULTIY,
C     MULTIX, AND EHYX AS WELL.
C
C     SUBROUTINES CALLED
C     DRINV,MULTIY,MULTIX,SETMYX,SETPRD
C
C     ARGUMENTS
C     Y    - M-VECTOR OF CONTEMPORANEOUS VALUES Y_T OF THE TIME SERIES
C            M MUST BE LESS THAN 100
C            INPUT, 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     H    - COMPUTED VALUE OF THE LOG DENSITY
C            OUTPUT, REAL*8
C     DELH - COMPUTED GRADIENT OF H WITH RESPECT TO THETA, A VECTOR
C            OF LENGTH KX*KY+M+M*M*L+M*(M+1)/2
C            OUTPUT, 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 1+M+M*(M+1)/2+M+M.
C            REAL*8
C
C     REMARK
C     THE CALLS
C     CALL SETMYX(M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX,KY,KX,IDUMMY,IWORK)
C     CALL SETPRD(M,L,KY,KX,IDEGY,IDEGX)
C     MUST PRECEDE USAGE OF THIS ROUTINE, IN THAT ORDER.
C
      SUBROUTINE LNHYX(Y,X,H,DELH,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Y(1),X(1),DELH(1)
      REAL*8 A(1),B0(1),B(1),R(1)
      REAL*8 WORK(1)
      INTEGER*4 ALPHA(100)
C
C----Declarations for direct access of multi-indexes ------------------
C
      INTEGER ZZMID
      PARAMETER ( ZZMID = 5000 )
      INTEGER MIDYX(ZZMID)
      COMMON /ZSVMYX/ MIDYX
C
C     NOTE:
C        MIDYX(1..M*KY)            = Y MULTINDEXES STORED AS MxKY ARRAY
C        MIDYX(M*KY+1..M*KY+ML*KX) = X MULTINDEXES STORED AS MLxKX ARRAY
C
C     WARNING:
C        The value of the parameter ZZMID must agree with the value
C        set in the routines SETMYX, MULTIY, MULTIX, and also with
C        the value set in other routines containing these declarations.
C
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
C     DECLARATIONS SHARED BETWEEN SETPRD, PZX, EHYX, AND LNHYX
C-----------------------------------------------------------------------
      INTEGER*4  ZZLIM
      PARAMETER (ZZLIM=500)
      INTEGER*4 IY(ZZLIM),IDY(ZZLIM),IX(ZZLIM)
      COMMON /ZSETPR/ IY,IDY,IX
C-----------------------------------------------------------------------
C
C     LOGICAL:
C     INTEGER*4 IY(2,KY),IDY(M,KY),IX(2,KX)
C
C     11/9/86
C     THIS ROUTINE WAS CODED BY USING THE CODE FROM HYX, REMOVING ALL
C     MULTIPLICATIONS BY EXPQ, AND ADJUSTING THE LAST THREE STATEMENTS
C     BEFORE THE RETURN.  TO MAKE SENSE OF THE COMMENTS BELOW, READ
C     1.D0 FOR EXPQ.
C
C     12/5/88
C     LOOP 160 WAS MODIFIED TO IMPROVE EFFICIENCY.  LOGICALLY IT IS THE
C     SAME AS LOOP 160 IN HYX.
C
C     LOGICAL:
C     REAL*8 Y(M),X(M*L),DELH(KX*KY+M+M*M*L+M*(M+1)/2+1)
C     REAL*8 A(KX*KY),B0(M),B(M,M*L),R(M*(M+1)/2)
C     INTEGER*4 ALPHA(M)
C     OR:
C     REAL*8 Y(M),X(M,L),A(KX,KY)
C
C-----------------------------------------------------------------------
      SAVE
C
C     ERROR TRAP
C
      IF ((M.GT.100).OR.(M*L.GT.100)) THEN
      WRITE(3,3002)
      RETURN
      ENDIF
C
C     COMPUTE POINTERS TO STACKED PARAMETER VECTOR, USED TO POINT TO
C     ELEMENTS OF DELH=(D/DTHETA)H(Y|X,THETA)
C
      IBA  = 0
      IBB0 = IBA +KX*KY
      IBB  = IBB0+M
      IBR  = IBB +M*M*L
      LTHETA=IBR +M*(M+1)/2
C
      IAA  = IBA +1
      IAB0 = IBB0+1
      IAB  = IBB +1
      IAR  = IBR +1
C
C     COMPUTE POINTERS TO WORK VECTOR
C
      IBN   =0
      IBSZ  =IBN   +1
      IBRINV=IBSZ  +M
      IBZ   =IBRINV+M*(M+1)/2
      IBDZP =IBZ   +M
      LWORK =IBDZP +M
C
      IAN   =IBN   +1
      IASZ  =IBSZ  +1
      IARINV=IBRINV+1
      IAZ   =IBZ   +1
      IADZP =IBDZP +1
C
C     COMPUTE LENGTHS OF VARIOUS ARRAYS
C
      LA=KY*KX
      LR=M*(M+1)/2
      ML=M*L
C
C     INITIALIZE DELH TO ZERO
C
      DO 10 I=1,LTHETA
10    DELH(IBA+I)=0.D0
C
C     COMPUTE Z=INV(R)(Y-B0-B*X)
C
      DO 20 I=1,M
20    WORK(IBSZ+I)=Y(I)-B0(I)
      DO 30 J=1,ML
      DO 30 I=1,M
30    WORK(IBSZ+I)=WORK(IBSZ+I)-B(M*(J-1)+I)*X(J)
      DO 40 I=1,LR
40    WORK(IBRINV+I)=R(I)
      CALL DRINV(WORK(IARINV),M)
      DO 50 I=1,M
50    WORK(IBZ+I)=0.D0
      DO 60 J=1,M
      DO 60 I=1,J
60    WORK(IBZ+I)=WORK(IBZ+I)+WORK(IBRINV+J*(J-1)/2+I)*WORK(IBSZ+J)
C
C     COMPUTE P, (D/DA)P, AND (D/DZ)P
C
C     INITIALIZE ACCUMULATORS FOR P AND (D/DZ)P
      P=0.D0
      DO 100 I=1,M
100   WORK(IBDZP+I)=0.D0
C
C     ------------------------------------------------------------------
C     12/5/88 MODIFICATION OF LOOP 160
C     ------------------------------------------------------------------
C
C     INITIALIZE WY & WX, STORED AS LEADING ROW AND COLUMN OF (D/DA)P
C     WHICH IS THE LOGICAL KX*KY LEADING SUBMATRIX OF DELH.
C
      DELH(IBA+1)=1.D0
      IF(KY.LE.1) GO TO 103
      DO 102 IALPHA=2,KY
      IY1=IY(2*(IALPHA-1)+1)
      IY2=IY(2*(IALPHA-1)+2)
102   DELH(IBA+KX*(IALPHA-1)+1)=DELH(IBA+KX*(IY1-1)+1)*WORK(IBZ+IY2)
103   CONTINUE
      IF(KX.LE.1) GO TO 105
      DO 104 IBETA=2,KX
      IX1=IX(2*(IBETA-1)+1)
      IX2=IX(2*(IBETA-1)+2)
104   DELH(IBA+IBETA)=DELH(IBA+IX1)*X(IX2)
105   CONTINUE
C
C     COMPUTE P, (D/DZ)P, AND (D/DA)P
C
      DO 160 IALPHA=1,KY
C       DIRECT ACCESS OF MULTI-INDEXES, REPLACES
C       CALL MULTIY(IALPHA,ALPHA)
        DO 101 I=1,M
101     ALPHA(I)=MIDYX(M*(IALPHA-1)+I)
        YALPHA=DELH(IBA+KX*(IALPHA-1)+1)
        DO 130 I=1,M
        IDY0=IDY(M*(IALPHA-1)+I)
130     WORK(IBSZ+I)=ALPHA(I)*DELH(IBA+KX*(IDY0-1)+1)
      DO 160 IBETA=1,KX
        XBETA=DELH(IBA+IBETA)
C     ACCUMULATE THE MONOMIALS TO GET THE POLYNOMIAL P & (D/DZ)P
      P=P+A(KX*(IALPHA-1)+IBETA)*XBETA*YALPHA
      DO 150 I=1,M
150   WORK(IBDZP+I)=WORK(IBDZP+I)
     &             +A(KX*(IALPHA-1)+IBETA)*XBETA*WORK(IBSZ+I)
C     PUT THE MONOMIALS IN PROPER LOCATIONS TO GET (D/DA)P
      DELH(IBA+KX*(IALPHA-1)+IBETA)=XBETA*YALPHA
160   CONTINUE
C
C     ------------------------------------------------------------------
C
C     COMPUTE:
C     (D/DB0')P = [(D/DZ')P]*(-1)*INV(R)
C     (D/DB')P = [(D/DZ')P]*(-1)*INV(R)*(X'@I)
C        WHERE @ DENOTES KRONEKER PRODUCT
C     (D/DR_IJ)P
C        =[(D/DZ')P]*(-1)*INV(R)*T(IJ)*INV(R)*(Y-B0-B*X)
C        WHERE T(IJ) HAS A ONE IN THE IJ_TH ELEMENT AND ZEROES ELSEWHERE
C
      DO 165 I=1,M
165   WORK(IBSZ+I)=0.D0
      DO 170 J=1,M
      DO 170 I=1,J
170   WORK(IBSZ+J)=WORK(IBSZ+J)+WORK(IBRINV+J*(J-1)/2+I)*WORK(IBDZP+I)
      DO 180 I=1,M
180   DELH(IBB0+I)=-WORK(IBSZ+I)
      DO 190 J=1,ML
      DO 190 I=1,M
      IJ=M*(J-1)+I
190   DELH(IBB+IJ)=-X(J)*WORK(IBSZ+I)
      DO 195 J=1,M
      DO 195 I=1,J
      IJ=J*(J-1)/2+I
195   DELH(IBR+IJ)=-WORK(IBSZ+I)*WORK(IBZ+J)
C
C     COMPUTE P**2 AND (D/DTHETA)(P**2)
C
      DO 250 I=1,LTHETA
250   DELH(IBA+I)=2.D0*P*DELH(IBA+I)
      P2=P*P
C
C     COMPUTE DABS[1/DET(R)] AND (D/DTHETA)DABS[1/DET(R)]
C
      DETR=1.D0
      DO 260 I=1,M
260   DETR=DETR*R(I*(I-1)/2+I)
      DETR=DABS(DETR)
      SR=1.D0/DETR
      DELSR=-SR*SR
      DO 270 I=1,M
      IJ=I*(I-1)/2+I
      DELH(IBR+IJ)=DELH(IBR+IJ)+(P2/SR)*DELSR*DETR/R(IJ)
270   CONTINUE
C
C     COMPUTE P**2*DABS[1/DET(R)] AND ADJUST DERIVATIVES
C
      H=P2*SR
      DO 280 I=1,LTHETA
280   DELH(IBA+I)=DELH(IBA+I)*SR
C
C     COMPUTE EXP(Q) = EXP[(-1/2)(Y-B0-B*X)'INV(R')INV(R)(Y-B0-B*X)]
C        AND (D/DTHETA)EXP(Q)
C     (D/DB0')EXP(Q)=EXP(Q)(Y-B0-B*X)'*INV'(R)*INV(R)
C     (D/DB')EXP(Q) =EXP(Q)(Y-B0-B*X)'*INV'(R)*INV(R)*(X'@I)
C        WHERE @ DENOTES KRONEKER PRODUCT
C     (D/DR_IJ)EXP(Q)
C        =EXP(Q)[INV(R)*(Y-B0-B*X)]'*INV(R)*T(IJ)*[INV(R)*(Y-B0-B*X)]
C        WHERE T(IJ) HAS A ONE IN THE IJ_TH ELEMENT AND ZEROES ELSEWHERE
C
      Q=0.D0
      DO 350 I=1,M
350   Q=Q+WORK(IBZ+I)**2
      Q=-0.5D0*Q
C     EXPQ=DEXP(Q)
      DO 360 I=1,M
360   WORK(IBSZ+I)=0.D0
      DO 370 J=1,M
      DO 370 I=1,J
370   WORK(IBSZ+J)=WORK(IBSZ+J)+WORK(IBRINV+J*(J-1)/2+I)*WORK(IBZ+I)
      DO 380 I=1,M
380   DELH(IBB0+I)=DELH(IBB0+I)+H*WORK(IBSZ+I)
      DO 390 J=1,ML
      DO 390 I=1,M
      IJ=M*(J-1)+I
390   DELH(IBB+IJ)=DELH(IBB+IJ)+H*X(J)*WORK(IBSZ+I)
      DO 395 J=1,M
      DO 395 I=1,J
      IJ=J*(J-1)/2+I
395   DELH(IBR+IJ)=DELH(IBR+IJ)+H*WORK(IBSZ+I)*WORK(IBZ+J)
C
C     COMPUTE H=P**2*DEXP[1/DET(R)]*EXP(Q)*{DSQRT[1/(2*PI)]}**M
C     AND ADJUST DERIVATIVES
C
      SCALE=0.39894228040143267794D0**M
      DO 430 I=1,LTHETA
430   DELH(IBA+I)=DELH(IBA+I)/H
      H=DLOG(SCALE*H)+Q
C
C     WORK(IAN)=SCALE*SR*EXPQ
C
      RETURN
3002  FORMAT(' ERROR RETURN FROM LNHYX, M OR M*L EXCEEDS 100')
      END
C.hr PZX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     PZX      9/1/91
C
C     PURPOSE
C     COMPUTE THE POLYNOMIAL
C
C     P(Z|X,A)=[SUM_J  A(X,J)*Z**ALPHA(J)]**2
C
C     AND ITS DERIVATIVES
C
C     DELAP = (D/D_A)P(Z|X,A)
C     DELZP = (D/D_Z)P(Z|X,A)
C
C     WHERE
C     A(X,J) = SUM_I A(I,J)*X**BETA(I), I = 1, ..., KX, J = 1, ..., KY
C
C     USAGE
C     CALL PYZ(Z,X,P,DELAP,DELZP,A,M,L,KY,KX,WORK)
C
C     SUBROUTINES CALLED
C     SETMYX,SETPRD
C
C     ARGUMENTS
C     Z     - M-VECTOR OF NORMALIZED VALUES Z=INV(R)*(Y_T-B0-B*X) OF 
C             THE TIME SERIES Y. M MUST BE LESS THAN 100.
C             INPUT, 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     P     - COMPUTED VALUE OF THE POLYNOMIAL.
C             OUTPUT, REAL*8
C     DELAP - COMPUTED GRADIENT OF P WITH RESPECT TO A, A VECTOR OF
C             LENGTH KX*KY.
C             OUTPUT, REAL*8
C     DELZP - COMPUTED GRADIENT OF H WITH RESPECT TO THETA, A VECTOR OF
C             LENGTH M.
C             OUTPUT, REAL*8
C     A     - COEFFICIENTS OF THE POLYNOMIAL, A KX BY KY MATRIX, STORED
C             COLUMNWISE.
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     WORK  - WORK VECTOR OF LENGTH M.
C             REAL*8
C
C     REMARK
C     THE CALLS
C     CALL SETMYX(M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX,KY,KX,IDUMMY,IWORK)
C     CALL SETPRD(M,L,KY,KX,IDEGY,IDEGX)
C     MUST PRECEDE USE OF THIS ROUTINE
C
      SUBROUTINE PZX(Z,X,P,DELAP,DELZP,A,M,L,KY,KX,WORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      SAVE

      REAL*8 Z(M),X(M*L),DELAP(KX,KY),DELZP(M)
      REAL*8 A(KX,KY)
      REAL*8 WORK(M)
      REAL*8 P,TWOP,YALPHA,XBETA

      INTEGER*4 M,L,KY,KX
      INTEGER*4 IALPHA,IBETA,IY1,IY2,IDY0,IX1,IX2,I

C     MODIFIED 3/27/92 TO USE SETPRD (HORNER'S RULE) TO COMPUTE P.

C-----------------------------------------------------------------------
C     DECLARATIONS FOR DIRECT ACCESS OF MULTI-INDEXES 
C-----------------------------------------------------------------------
      INTEGER   ZZMID
      PARAMETER ( ZZMID = 5000 )
      INTEGER   MIDYX(ZZMID)
      COMMON /ZSVMYX/ MIDYX
C
C     NOTE:
C        MIDYX(1..M*KY)            = Y MULTINDEXES STORED AS MxKY ARRAY
C        MIDYX(M*KY+1..M*KY+ML*KX) = X MULTINDEXES STORED AS MLxKX ARRAY
C
C     WARNING:
C        The value of the parameter ZZMID must agree with the value
C        set in the routines SETMYX, MULTIY, MULTIX, and also with
C        the value set in other routines containing these declarations.
C
C-----------------------------------------------------------------------

C-----------------------------------------------------------------------
C     DECLARATIONS SHARED BETWEEN SETPRD, PZX, EHYX, AND LNHYX
C-----------------------------------------------------------------------
      INTEGER*4 ZZLIM
      PARAMETER (ZZLIM=500)
      INTEGER*4 IY(ZZLIM),IDY(ZZLIM),IX(ZZLIM)
      COMMON /ZSETPR/ IY,IDY,IX
C-----------------------------------------------------------------------

C     ZERO OUT ACCUMULATORS

      P=0.D0

      DO I=1,M
        DELZP(I)=0.D0
      END DO

C     INITIALIZE WY, STORED AS LEADING ROW OF (D/DA)P

      DELAP(1,1)=1.D0

      IF (KY.GT.1) THEN
        DO IALPHA=2,KY
          IY1=IY(2*(IALPHA-1)+1)
          IY2=IY(2*(IALPHA-1)+2)
          DELAP(1,IALPHA)=DELAP(1,IY1)*Z(IY2)
        END DO
      END IF

C     COMPUTE P, (D/DZ)P, AND (D/DA)P

      IF (KX.GT.1) THEN

C       INITIALIZE WX, STORED AS LEADING COLUMN OF (D/DA)P

        DO IBETA=2,KX
          IX1=IX(2*(IBETA-1)+1)
          IX2=IX(2*(IBETA-1)+2)
          DELAP(IBETA,1)=DELAP(IX1,1)*X(IX2)
        END DO

        DO IALPHA=1,KY

          YALPHA=DELAP(1,IALPHA)

          DO I=1,M
            IDY0=IDY(M*(IALPHA-1)+I)
            WORK(I)=DFLOAT(MIDYX(M*(IALPHA-1)+I))*DELAP(1,IDY0)
          END DO

        DO IBETA=1,KX

          XBETA=DELAP(IBETA,1)

          P=P+A(IBETA,IALPHA)*XBETA*YALPHA

          DO I=1,M
            DELZP(I)=DELZP(I)+A(IBETA,IALPHA)*XBETA*WORK(I)
          END DO

          DELAP(IBETA,IALPHA)=XBETA*YALPHA

        END DO
        END DO

      ELSE

        DO IALPHA=1,KY

          YALPHA=DELAP(1,IALPHA)

          DO I=1,M
            IDY0=IDY(M*(IALPHA-1)+I)
            DELZP(I)=DELZP(I)
     &        +DFLOAT(MIDYX(M*(IALPHA-1)+I))*(A(1,IALPHA)*DELAP(1,IDY0))
          END DO

          P=P+A(1,IALPHA)*YALPHA

        END DO

      END IF

C     COMPUTE P**2 AND (D/DA)(P**2)

      TWOP=2.D0*P

      DO IALPHA=1,KY
      DO IBETA=1,KX
        DELAP(IBETA,IALPHA)=TWOP*DELAP(IBETA,IALPHA)
      END DO
      END DO

      DO I=1,M
        DELZP(I)=TWOP*DELZP(I)
      END DO

      P = P*P

      RETURN
      END

