C.hr snpcode4.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     snpcode4.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 HYX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     HYX      10/28/86
C
C     PURPOSE
C     COMPUTE AN UNNORMALIZED HERMITE CONDITIONAL DENSITY,
C
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
C     AND ITS DERIVATIVE,
C
C     DELH(Y|X,THETA) = (D/D_THETA)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 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 HYX(Y,X,H,DELH,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK)
C
C     SUBROUTINES CALLED
C     DRINV,MULTIY,MULTIX,SETMYX
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 DENSITY
C            OUTPUT, REAL*8
C     DELH - COMPUTED GRADIENT OF H 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 POYNOMIAL IN Y, SEE 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 1+M+M*(M+1)/2+M+M. ON RETURN CONTAINS
C            ( N(Y|B0+B*X,SIGMA) | (Y-B0-B*X)'INV(SIGMA) | INV(R) |
C            (Y-B0-B*X)'INV(R') | (D/DZ')P )
C            REAL*8
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 HYX(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),BETA(100)
      SAVE
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),BETA(M*L)
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)) 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     LOOP 160 ACCUMULATES MONOMIALS TO GET P AND (D/DZ)P;
C     IT ALSO PUTS THEM IN DELH TO GET (D/DA)P
      DO 160 IALPHA=1,KY
C        GET EXPONANT FOR MONOMIAL YALPHA
         CALL MULTIY(IALPHA,ALPHA)
C        INITIALIZE ACCUMULATORS FOR YALPHA AND (D/DZ)YALPHA
         Y1=1.D0
         Y2=1.D0
         DO 105 I=1,M
105      WORK(IBSZ+I)=ALPHA(I)
C        LOOP 115 ACCUMULATES DEGREE ONE FACTORS, IF ANY
         DO 115 I=1,M
         IF (ALPHA(I).GT.0) THEN
         Y1=Y1*WORK(IBZ+I)
         DO 110 J=1,M
110      IF (I.NE.J) WORK(IBSZ+J)=WORK(IBSZ+J)*WORK(IBZ+I)
         ENDIF
115      CONTINUE
C        LOOPS 120 & 130 ACCUMULATE HIGHER ORDER FACTORS, IF ANY
         DO 120 I=1,M
         NTIMES=ALPHA(I)-1
120      IF(NTIMES.GT.0) Y2=Y2*(WORK(IBZ+I)**NTIMES)
         YALPHA=Y1*Y2
         DO 130 I=1,M
130      WORK(IBSZ+I)=WORK(IBSZ+I)*Y2
      DO 160 IBETA=1,KX
C        GET EXPONANT FOR MONOMIAL XBETA
         CALL MULTIX(IBETA,BETA)
C        INITIALIZE ACCUMULATOR
         XBETA=1.D0
C        ACCUMULATE
         DO 140 I=1,ML
140      IF (BETA(I).GT.0) XBETA=XBETA*(X(I)**BETA(I))
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     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
      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)=EXPQ*DELH(IBB0+I)+H*EXPQ*WORK(IBSZ+I)
      DO 390 J=1,ML
      DO 390 I=1,M
      IJ=M*(J-1)+I
390   DELH(IBB+IJ)=EXPQ*DELH(IBB+IJ)+H*EXPQ*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)=EXPQ*DELH(IBR+IJ)+H*EXPQ*WORK(IBSZ+I)*WORK(IBZ+J)
C
C     COMPUTE P**2*DABS[1/DET(R)]*EXP(Q) AND ADJUST DERIVATIVES
C
      H=H*EXPQ
      DO 420 I=1,LA
420   DELH(IBA+I)=EXPQ*DELH(IBA+I)
C
C     COMPUTE H=P**2*DEXP[1/DET(R)]*EXP(Q)*{DSQRT[1/(2*PI)]}**M
C     AND ADJUST DERIVATIVES
C
      TWOPI=8.D0*DATAN(1.D0)
      SCALE=(1.D0/DSQRT(TWOPI))**M
      H=SCALE*H
      DO 430 I=1,LTHETA
430   DELH(IBA+I)=SCALE*DELH(IBA+I)
C
      WORK(IAN)=SCALE*SR*EXPQ
C
      RETURN
3002  FORMAT(' ERROR RETURN FROM HYX, M OR M*L EXCEEDS 100')
      END
C.hr RESID
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     RESID    12/17/87
C
C     PURPOSE
C     COMPUTE A NORMALIZED RESIDUAL FROM A HERMITE CONDITIONAL DENSITY
C
C     RES = INV{HALF[VAR(Y|X)]}[Y-E(Y|X)]
C         = INV(P)[INV(R)(Y-B0-B*X)-E(Z|X)]
C
C     AND ITS DERIVATIVE
C
C     DER = (D/D_THETA)RES
C
C     WHERE
C     THE DENSITIES DEFINING EXPECTATIONS ARE:
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)/INTEGRAL_OF_NUMERATOR
C     H(Z|X,THETA)= [SUM_J  A(X,J)*Z**ALPHA(I)]**2
C                     * N(Z | 0, I)/INTEGRAL_OF_NUMERATOR
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     VAR(Z|X) = E[Z-E(Z|X)][Z-E(Z|X)]' = P*P'
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 RESID(Y,X,RES,DER,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,IWORK)
C
C     SUBROUTINES CALLED
C     MULTIY,MULTIX,SETMYX,SETUPM,DRINV,FACTOR,EHYX,EYLX
C
C     ARGUMENTS
C     Y     - M-VECTOR CONTAINING OBSERVED VALUE OF Y
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     RES   - COMPUTED RESIDUAL, AN M-VECTOR
C             OUTPUT, REAL*8
C     DER   - COMPUTED JACOBIAN OF RES WITH RESPECT TO THETA, AN M BY
C             LTHETA MATRIX WHERE LTHETA=KX*KY+M+M*M*L+M*(M+1)/2,
C             STORED COLUMNWISE.
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, 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
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 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
C             (M+M*(M+1)/2+1)*LTHETA+2*M+MAX0(2*KY,(M*(M+1)/2))
C             +(M*(M+1)/2+1)*(M*(M+1)/2)
C             WORKSPACE, REAL*8
C     IWORK - WORK VECTOR OF LENGTH M
C             WORKSPACE, INTEGER*4
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
C
      SUBROUTINE
     &   RESID(Y,X,RES,DER,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,IWORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Y(1),X(1),DER(1),RES(1)
      REAL*8 A(1),B0(1),B(1),R(1)
      REAL*8 WORK(1)
      INTEGER*4 IWORK(1)
      SAVE
C
C     LOGICAL:
C     REAL*8 Y(M),X(M*L),DER(M,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     OR:
C     REAL*8 Y(M),X(M,L),A(KX,KY)
C
      IOUT=3
C
C     ERROR TRAP
C
      IF ((M.GT.100).OR.(M*L.GT.100)) THEN
      WRITE(IOUT,3002)
      RETURN
      ENDIF
C
C     COMPUTE POINTERS TO STACKED PARAMETER VECTOR, USED TO POINT TO
C     ELEMENTS OF DER=(D/D_THETA)RES
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
      IBDMU =0
      IBDP  =IBDMU +M*LTHETA
      IBBOT =IBDP  +(M*(M+1)/2)*LTHETA
      IBZ   =IBBOT +LTHETA
      IBMUZ =IBZ   +M
      IBW   =IBMUZ +M
      IBPINV=IBW   +MAX0(2*KY,(M*(M+1)/2))
      IBS   =IBPINV+M*(M+1)/2
      LWORK =IBS   +(M*(M+1)/2)*(M*(M+1)/2)
C
      IADMU =IBDMU +1
      IADP  =IBDP  +1
      IABOT =IBBOT +1
      IAZ   =IBZ   +1
      IAMUZ =IBMUZ +1
      IAW   =IBW   +1
      IAPINV=IBPINV+1
      IAS   =IBS   +1
C
C     EQUIVALENCED POINTERS
C
      IBRINV=IBPINV
      IBPZ  =IBS
C
      IARINV=IBRINV+1
      IAPZ  =IBPZ  +1
C
C     COMPUTE LENGTHS OF VARIOUS ARRAYS
C
      LA=KY*KX
      LR=M*(M+1)/2
      ML=M*L
C
C     INITIALIZE DER TO ZERO
C
      DO 10 J=1,LTHETA
      DO 10 I=1,M
10    DER(M*(J-1)+I)=0.D0
C
C     COMPUTE Z=INV(R)(Y-B0-B*X)
C
      DO 20 I=1,M
20    RES(I)=Y(I)-B0(I)
      DO 30 J=1,ML
      DO 30 I=1,M
30    RES(I)=RES(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)*RES(J)
C
C     COMPUTE:
C
C     (D/DB0')Z  = (-1)*INV(R)
C     (D/DB')Z   = (-1)*INV(R)*(X'@I)
C     (D/DR_IJ)Z = (-1)*INV(R)*T(IJ)*INV(R)*(Y-B0-B*X)
C
C     WHERE
C     @ DENOTES KRONECKER PRODUCT AND
C     T(IJ) HAS A ONE IN THE IJ_TH ELEMENT AND ZEROES ELSEWHERE
C
      DO 70 I=1,M
      DO 70 K=1,I
70    DER(M*(IBB0+I-1)+K)=-WORK(IBRINV+I*(I-1)/2+K)
C
      DO 80 J=1,ML
      DO 80 I=1,M
      IJ=M*(J-1)+I
      DO 80 K=1,I
80    DER(M*(IBB+IJ-1)+K)=-X(J)*WORK(IBRINV+I*(I-1)/2+K)
C
      DO 90 J=1,M
      DO 90 I=1,J
      IJ=J*(J-1)/2+I
      DO 90 K=1,I
90    DER(M*(IBR+IJ-1)+K)=-WORK(IBRINV+I*(I-1)/2+K)*WORK(IBZ+J)
C
C     COMPUTE E(Z|X) AND PUT IT IN MUZ
C     COMPUTE THE TRANSPOSE OF (D/D_THETA)E(Z|X) AND PUT IT IN DMU
C     SUBTRACT E(Z|X) FROM Z AND (D/D_THETA)E(Z|X) FROM DER
C
      CALL EHYX(RES,X,EH,WORK(IABOT),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      DO 140 I=1,M
      DO 110 K=1,M
110   IWORK(K)=0
      IWORK(I)=1
      CALL EYLX(IWORK,RES,X,EYL,WORK(IADMU+LTHETA*(I-1)),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      RATIO=EYL/EH
      WORK(IBMUZ+I)=RATIO
      DO 115 J=1,LTHETA
115   WORK(IBDMU+LTHETA*(I-1)+J)
     &   =(WORK(IBDMU+LTHETA*(I-1)+J)/EH-(RATIO/EH)*WORK(IBBOT+J))
      WORK(IBZ+I)=WORK(IBZ+I)-RATIO
      DO 120 J=1,LTHETA
120   DER(M*(J-1)+I)=DER(M*(J-1)+I)-WORK(IBDMU+LTHETA*(I-1)+J)
140   CONTINUE
C
C     COMPUTE E(ZZ'|X)-E(Z|X)E(Z'|X), PUT IT IN PINV
C
C     COMPUTE THE TRANSPOSE OF (D/D_THETA)[E(ZZ'|X)-E(Z|X)E(Z'|X)]
C     AND PUT IT IN DP
C
      DO 160 J=1,M
      DO 160 I=1,J
      IJ=J*(J-1)/2+I
      DO 150 K=1,M
150   IWORK(K)=0
      IWORK(I)=1
      IWORK(J)=IWORK(J)+1
      CALL EYLX(IWORK,RES,X,EYL,WORK(IADP+LTHETA*(IJ-1)),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      WORK(IBS+IJ)=EYL/EH
      WORK(IBPINV+IJ)=EYL/EH-WORK(IBMUZ+I)*WORK(IBMUZ+J)
160   CONTINUE
C
      DO 165 J=1,M
      DO 165 I=1,J
      IJ=J*(J-1)/2+I
      DO 165 K=1,LTHETA
      WORK(IBDP+LTHETA*(IJ-1)+K)
     &   =WORK(IBDP+LTHETA*(IJ-1)+K)/EH
     &      -(WORK(IBS+IJ)/EH)*WORK(IBBOT+K)
     &      -WORK(IBMUZ+I)*WORK(IBDMU+LTHETA*(J-1)+K)
     &      -WORK(IBMUZ+J)*WORK(IBDMU+LTHETA*(I-1)+K)
165   CONTINUE
C
C     FACTOR E(ZZ'|X)-E(Z|X)E(Z'|X)=PP', INVERT P, PUT RESULT IN PINV
C     COMPUTE (D/D_THETA)P, PUT RESULT IN DP
C
      CALL FACTOR(WORK(IAPINV),WORK(IAS),M,1.D-13,IER)
      IF (IER.LT.0) THEN
      WRITE(IOUT,3003) IER
      RETURN
      ENDIF
      CALL DRINV(WORK(IAPINV),M)
C
      DO 174 J=1,LTHETA
      DO 171 I=1,LR
171   WORK(IBW+I)=0.D0
      DO 172 I=1,LR
      DO 172 K=1,LR
172   WORK(IBW+I)=WORK(IBW+I)
     &            +WORK(IBS+LR*(K-1)+I)*WORK(IBDP+LTHETA*(K-1)+J)
      DO 173 I=1,LR
173   WORK(IBDP+LTHETA*(I-1)+J)=WORK(IBW+I)
174   CONTINUE
C
C     COMPUTE (PINV)Z AND PUT IT IN PZ
C
      DO 175 I=1,M
175   WORK(IBPZ+I)=0.D0
      DO 177 J=1,M
      DO 177 I=1,J
177   WORK(IBPZ+I)=WORK(IBPZ+I)+WORK(IBPINV+J*(J-1)/2+I)*WORK(IBZ+J)
C
C     PUT DER TO (PINV)DER
C
      DO 184 K=1,LTHETA
      DO 181 I=1,M
181   RES(I)=0.D0
      DO 182 J=1,M
      DO 182 I=1,J
182   RES(I)=RES(I)+WORK(IBPINV+J*(J-1)/2+I)*DER(M*(K-1)+J)
      DO 183 I=1,M
183   DER(M*(K-1)+I)=RES(I)
184   CONTINUE
C
C     ADD [(D/D_VECP)(INVP*Z)]*(D/D_THETA)VECP TO DER
C
      DO 196 JJ=1,LTHETA
      DO 190 K=1,M
190   RES(K)=0.D0
      DO 192 J=1,M
      DO 192 I=1,J
      IJ=J*(J-1)/2+I
      DO 192 K=1,I
192   RES(K)=RES(K)
     &   -WORK(IBPINV+I*(I-1)/2+K)*WORK(IBPZ+J)
     &      *WORK(IBDP+LTHETA*(IJ-1)+JJ)
      DO 194 K=1,M
194   DER(M*(JJ-1)+K)=DER(M*(JJ-1)+K)+RES(K)
196   CONTINUE
C
C     PUT PZ IN RES
C
      DO 198 I=1,M
198   RES(I)=WORK(IBPZ+I)
C
      RETURN
3002  FORMAT(' ERROR RETURN FROM RESID, M OR M*L EXCEEDS 100')
3003  FORMAT(' ERROR RETURN FROM RESID, FACTOR RETURNED IER =',I5)
      END
C.hr MEAN
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     MEAN     12/17/87
C
C     PURPOSE
C     COMPUTE THE MEAN OF A HERMITE CONDITIONAL DENSITY
C
C     EY = E(Y|X)
C        = B0 + B*X + R*E(Z|X)
C
C     AND ITS DERIVATIVE
C
C     DER = (D/D_THETA)MEAN
C
C     WHERE
C     THE DENSITIES DEFINING EXPECTATIONS ARE:
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)/INTEGRAL_OF_NUMERATOR
C     H(Z|X,THETA)= [SUM_J  A(X,J)*Z**ALPHA(I)]**2
C                     * N(Z | 0, I)/INTEGRAL_OF_NUMERATOR
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 MEAN(Y,X,EY,DER,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,IWORK)
C
C     SUBROUTINES CALLED
C     MULTIY,MULTIX,SETMYX,SETUPM,EHYX,EYLX
C
C     ARGUMENTS
C     Y     - M-VECTOR CONTAINING OBSERVED VALUE OF Y
C             WORKSPACE, 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     EY    - COMPUTED MEAN, AN M-VECTOR
C             OUTPUT, REAL*8
C     DER   - COMPUTED JACOBIAN OF EY WITH RESPECT TO THETA, AN M BY
C             LTHETA MATRIX WHERE LTHETA=KX*KY+M+M*M*L+M*(M+1)/2,
C             STORED COLUMNWISE.
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, 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
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 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
C             (M+1)*LTHETA+M+MAX0(2*KY,(M*(M+1)/2))
C             WORKSPACE, REAL*8
C     IWORK - WORK VECTOR OF LENGTH M
C             WORKSPACE, INTEGER*4
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
C
      SUBROUTINE
     &   MEAN(Y,X,EY,DER,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,IWORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Y(1),X(1),DER(1),EY(1)
      REAL*8 A(1),B0(1),B(1),R(1)
      REAL*8 WORK(1)
      INTEGER*4 IWORK(1)
      SAVE
C
C     LOGICAL:
C     REAL*8 Y(M),X(M*L),DER(M,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     OR:
C     REAL*8 Y(M),X(M,L),A(KX,KY)
C
      IOUT=3
C
C     ERROR TRAP
C
      IF ((M.GT.100).OR.(M*L.GT.100)) THEN
      WRITE(IOUT,3002)
      RETURN
      ENDIF
C
C     COMPUTE POINTERS TO STACKED PARAMETER VECTOR, USED TO POINT TO
C     ELEMENTS OF DER=(D/D_THETA)EY
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
      IBDMU =0
      IBBOT =IBDMU +M*LTHETA
      IBMUZ =IBBOT +LTHETA
      IBW   =IBMUZ +M
      LWORK =IBW   +MAX0(2*KY,(M*(M+1)/2))
C
      IADMU =IBDMU +1
      IABOT =IBBOT +1
      IAMUZ =IBMUZ +1
      IAW   =IBW   +1
C
C     COMPUTE LENGTHS OF VARIOUS ARRAYS
C
      ML=M*L
C
C     INITIALIZE DER TO ZERO
C
      DO 10 J=1,LTHETA
      DO 10 I=1,M
10    DER(M*(J-1)+I)=0.D0
C
C     PUT EY = B0+B*X
C
      DO 20 I=1,M
20    EY(I)=B0(I)
      DO 30 J=1,ML
      DO 30 I=1,M
30    EY(I)=EY(I)+B(M*(J-1)+I)*X(J)
C
C     PUT DER = [(D/DB0')(B0+B*X)|(D/DB')(B0+B*X)]
C             = [ I | (X'@I) ]
C
      DO 70 I=1,M
70    DER(M*(IBB0+I-1)+I)=1.D0
C
      DO 80 J=1,ML
      IPTR=M*IBB+M*M*(J-1)
      DO 80 I=1,M
80    DER(IPTR+M*(I-1)+I)=X(J)
C
C     COMPUTE E(Z|X) AND PUT IT IN MUZ
C     COMPUTE THE TRANSPOSE OF (D/D_THETA)E(Z|X) AND PUT IT IN DMU
C
      CALL EHYX(Y,X,EH,WORK(IABOT),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      DO 140 I=1,M
      DO 110 K=1,M
110   IWORK(K)=0
      IWORK(I)=1
      CALL EYLX(IWORK,Y,X,EYL,WORK(IADMU+LTHETA*(I-1)),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      RATIO=EYL/EH
      WORK(IBMUZ+I)=RATIO
      DO 115 J=1,LTHETA
115   WORK(IBDMU+LTHETA*(I-1)+J)
     &   =(WORK(IBDMU+LTHETA*(I-1)+J)/EH-(RATIO/EH)*WORK(IBBOT+J))
140   CONTINUE
C
C     ADD R*E(Z|X) TO EY
C
      DO 150 J=1,M
      DO 150 I=1,J
      IJ=J*(J-1)/2+I
150   EY(I)=EY(I)+R(IJ)*WORK(IBMUZ+J)
C
C     ADD (D/D_R)(R*E) EVALUATED AT E = E(Z|X) TO DER
C
      DO 160 J=1,M
      IPTR=M*IBR+M*(J*(J-1)/2)
      DO 160 I=1,J
160   DER(IPTR+M*(I-1)+I)=DER(IPTR+M*(I-1)+I)+WORK(IBMUZ+J)
C
C     ADD R[(D/D_THETA)E(Z|X)] TO DER
C
      DO 170 K=1,LTHETA
      DO 170 J=1,M
      DO 170 I=1,J
      IJ=J*(J-1)/2+I
170   DER(M*(K-1)+I)
     &   =DER(M*(K-1)+I)+R(IJ)*WORK(IBDMU+LTHETA*(J-1)+K)
C
      RETURN
3002  FORMAT(' ERROR RETURN FROM MEAN, M OR M*L EXCEEDS 100')
      END
C.hr VARCOV
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     VARCOV   12/17/87
C
C     PURPOSE
C     COMPUTE THE VARIANCE-COVARIANCE MATRIX OF A HERMITE CONDITIONAL
C     DENSITY
C
C     VAR = VAR(Y|X)
C         = R*VAR(Z|X)*R'
C
C     AND ITS DERIVATIVE
C
C     DER = (D/D_THETA)VARCOV
C
C     WHERE
C     THE DENSITIES DEFINING EXPECTATIONS ARE:
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)/INTEGRAL_OF_NUMERATOR
C     H(Z|X,THETA)= [SUM_J  A(X,J)*Z**ALPHA(I)]**2
C                     * N(Z | 0, I)/INTEGRAL_OF_NUMERATOR
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     VAR(Z|X) = E[Z-E(Z|X)][Z-E(Z|X)]' = P*P'
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 VARCOV(Y,X,VAR,DER,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,IWORK)
C
C     SUBROUTINES CALLED
C     MULTIY,MULTIX,SETMYX,SETUPM,EHYX,EYLX
C
C     ARGUMENTS
C     Y     - M-VECTOR CONTAINING OBSERVED VALUE OF Y
C             WORKSPACE, 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     VAR   - COMPUTED UPPER TRIANGLE OF THE VARIANCE-COVARIANCE MATRIX,
C             STORED AS AN LR-VECTOR WHERE LR=M*(M+1)/2. THE ORDERING IS
C             VAR=(V11,V12,V22,V13,V23,V33,...,VMM).
C             OUTPUT, REAL*8
C     DER   - COMPUTED JACOBIAN OF VAR WITH RESPECT TO THETA, AN LR BY
C             LTHETA MATRIX WHERE LTHETA=KX*KY+M+M*M*L+M*(M+1)/2,
C             STORED COLUMNWISE.
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 LR-VECTOR CONTAINING THE UPPER TRIANGLE OF THE CHOLESKY
C             FACTOR OF THE VARIANCE-COVARIANCE MATRIX SIGMA=R*R' OF THE
C             VAR PART OF THE DENSITY STORED COLUMNWISE,
C             R=(R11,R12,R22,R13,R23,R33,...,RMM)
C             INPUT, REAL*8
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
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 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
C             M*M*(LTHETA+1)+LTHETA+MAX0(2*KY,(M*(M+1)/2))+M*(M+1)/2
C             WORKSPACE, REAL*8
C     IWORK - WORK VECTOR OF LENGTH M
C             WORKSPACE, INTEGER*4
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
C
      SUBROUTINE
     &   VARCOV(Y,X,VAR,DER,A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK,IWORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Y(1),X(1),DER(1),VAR(1)
      REAL*8 A(1),B0(1),B(1),R(1)
      REAL*8 WORK(1)
      INTEGER*4 ROW,COL,IWORK(1)
      SAVE
C
C     LOGICAL:
C     REAL*8 Y(M),X(M*L),DER(LR,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     OR:
C     REAL*8 Y(M),X(M,L),A(KX,KY)
C
      IOUT=3
C
C     ERROR TRAP
C
      IF ((M.GT.100).OR.(M*L.GT.100)) THEN
      WRITE(IOUT,3002)
      RETURN
      ENDIF
C
C     COMPUTE POINTERS TO STACKED PARAMETER VECTOR, USED TO POINT TO
C     ELEMENTS OF DER=(D/D_THETA)VAR
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
      IBGV  =0
      IBGD  =IBGV  +M*M
      IBBOT =IBGD  +M*M*LTHETA
      IBW   =IBBOT +LTHETA
      IBS   =IBW   +MAX0(2*KY,(M*(M+1)/2))
      LWORK =IBS   +M*(M+1)/2
C
      IAGV  =IBGV  +1
      IAGD  =IBGD  +1
      IABOT =IBBOT +1
      IAW   =IBW   +1
      IAS   =IBS   +1
C
C     EQUIVALENCED POINTERS
C
      IBMUZ =IBGV
      IBDMU =IBGD
C
      IAMUZ =IBMUZ +1
      IADMU =IBDMU +1
C
C     COMPUTE LENGTHS OF VARIOUS ARRAYS
C
      LA=KY*KX
      LR=M*(M+1)/2
      ML=M*L
      LGV=M*M
C
C     COMPUTE E(Z|X) AND PUT IT IN MUZ
C     COMPUTE THE TRANSPOSE OF (D/D_THETA)E(Z|X) AND PUT IT IN DMU
C
      CALL EHYX(Y,X,EH,WORK(IABOT),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      DO 140 I=1,M
      DO 110 K=1,M
110   IWORK(K)=0
      IWORK(I)=1
      CALL EYLX(IWORK,Y,X,EYL,WORK(IADMU+LTHETA*(I-1)),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      RATIO=EYL/EH
      WORK(IBMUZ+I)=RATIO
      DO 115 J=1,LTHETA
115   WORK(IBDMU+LTHETA*(I-1)+J)
     &   =(WORK(IBDMU+LTHETA*(I-1)+J)/EH-(RATIO/EH)*WORK(IBBOT+J))
140   CONTINUE
C
C     COMPUTE E(ZZ'|X)-E(Z|X)E(Z'|X), PUT IT IN VAR
C
C     COMPUTE THE TRANSPOSE OF (D/D_THETA)[E(ZZ'|X)-E(Z|X)E(Z'|X)]
C     AND PUT IT IN DER
C
      DO 160 J=1,M
      DO 160 I=1,J
      IJ=J*(J-1)/2+I
      DO 150 K=1,M
150   IWORK(K)=0
      IWORK(I)=1
      IWORK(J)=IWORK(J)+1
      CALL EYLX(IWORK,Y,X,EYL,DER(LTHETA*(IJ-1)+1),
     &   A,B0,B,R,M,L,KY,KX,IDEGY,IDEGX,WORK(IAW))
      WORK(IBS+IJ)=EYL/EH
      VAR(IJ)=EYL/EH-WORK(IBMUZ+I)*WORK(IBMUZ+J)
160   CONTINUE
C
      DO 165 J=1,M
      DO 165 I=1,J
      IJ=J*(J-1)/2+I
      DO 165 K=1,LTHETA
      DER(LTHETA*(IJ-1)+K)
     &   =DER(LTHETA*(IJ-1)+K)/EH
     &      -(WORK(IBS+IJ)/EH)*WORK(IBBOT+K)
     &      -WORK(IBMUZ+I)*WORK(IBDMU+LTHETA*(J-1)+K)
     &      -WORK(IBMUZ+J)*WORK(IBDMU+LTHETA*(I-1)+K)
165   CONTINUE
C
C     EXPAND CONTENTS OF VAR AND DER TO THE GENERAL MATRICES GV AND GD
C
      DO 175 J=1,M
      DO 175 I=1,J
      IJ1=J*(J-1)/2+I
      IJ0=M*(J-1)+I
      JI0=M*(I-1)+J
      WORK(IBGV+IJ0)=VAR(IJ1)
      WORK(IBGV+JI0)=VAR(IJ1)
      DO 170 K=1,LTHETA
      WORK(IBGD+LGV*(K-1)+IJ0)=DER(LTHETA*(IJ1-1)+K)
      WORK(IBGD+LGV*(K-1)+JI0)=DER(LTHETA*(IJ1-1)+K)
170   CONTINUE
175   CONTINUE
C
C     SET DER TO ZERO
C
      DO 180 J=1,LTHETA
      DO 180 I=1,LR
180   DER(LR*(J-1)+I)=0.D0
C
C     PUT VAR TO R*GV*R' AND PUT THE DIFFERENTIAL WRT R IN DER
C
      DO 185 COL=1,M
      DO 185 ROW=1,COL
      IJ=COL*(COL-1)/2+ROW
      VAR(IJ)=0.D0
      DO 185 I=ROW,M
      DO 185 J=COL,M
      VAR(IJ)=VAR(IJ)
     &   +R(I*(I-1)/2+ROW)*WORK(IBGV+M*(J-1)+I)*R(J*(J-1)/2+COL)
      DER(LR*(IBR+I*(I-1)/2+ROW-1)+IJ)=DER(LR*(IBR+I*(I-1)/2+ROW-1)+IJ)
     &   +WORK(IBGV+M*(J-1)+I)*R(J*(J-1)/2+COL)
      DER(LR*(IBR+J*(J-1)/2+COL-1)+IJ)=DER(LR*(IBR+J*(J-1)/2+COL-1)+IJ)
     &   +WORK(IBGV+M*(J-1)+I)*R(I*(I-1)/2+ROW)
185   CONTINUE
C
C     ADD R*GD(K)*R' TO DER
C
      DO 195 K=1,LTHETA
      DO 195 COL=1,M
      DO 195 ROW=1,COL
      IJ=COL*(COL-1)/2+ROW
      DO 190 I=ROW,M
      DO 190 J=COL,M
190   DER(LR*(K-1)+IJ)=DER(LR*(K-1)+IJ)
     & +R(I*(I-1)/2+ROW)*WORK(IBGD+LGV*(K-1)+M*(J-1)+I)*R(J*(J-1)/2+COL)
195   CONTINUE
C
      RETURN
3002  FORMAT(' ERROR RETURN FROM VARCOV, M OR M*L EXCEEDS 100')
      END
C.hr PCORD
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     PCORD    4/14/88
C
C     PURPOSE
C     GENERATE Y'S FOR PLOTTING PURPOSES.
C
C     NOTE
C     IN ORIGINAL UNITS, THE GRID GENERATED IS:
C     MU(I)-NGRD*SCALE*SE(I), ..., MU(I)+NGRD*SCALE*SE(I)
C     FOR I = 1, ..., M  WHERE MU(I) AND SE(I) ARE THE SAMPLE MEAN AND
C     VARIANCE OF THE ORIGINAL DATA. (2*NGRD+1)**M DISTINCT POINTS CAN
C     BE GENERATED IN SUCCESIVE CALLS TO PCORD WITH ISW=1.
C
C     USAGE
C     CALL PCORD(0,YPLT,Y,X,M,L,NGRD,S,MUSTAT,CSTAT,PSTAT,IY,YINC,WORK)
C     DO 10 I=1,(2*NGRD+1)**M
C     CALL PCORD(1,YPLT,Y,X,M,L,NGRD,S,MUSTAT,CSTAT,PSTAT,IY,YINC,WORK)
C     OTHER STATEMENTS
C10   CONTINUE
C
C     SUBROUTINES CALLED
C     DRINV, DGMPNT
C
C     ARGUMENTS
C     ISW   - SWITCH, ISW = 0: INITIALIZE X, PSTAT, IY, YINC;
C                     ISW = 1: COMPUTE THE NEXT YPLT AND Y VALUES.
C             INPUT, INTEGER*4
C     YPLT  - AN M-VECTOR OF PLOT POINTS IN ORIGINAL UNITS.  THE COMPUTE
C             GRID IS EQUALLY SPACED AND RECTANGULAR IN ORIGINAL UNITS.
C             OUTPUT, REAL*8
C     Y     - AN M-VECTOR OF PLOT POINTS IN TRANSFORMED UNITS AS
C             APPROPRIATE FOR CALLS TO HYX.
C             OUTPUT, 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     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     S     - SCALE FACTOR, 3.D0 IS A GOOD PLACE TO START.
C             INPUT, 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     PSTAT - M BY M INVERSE OF CSTAT.
C             OUTPUT, REAL*8
C     IY    - M-VECTOR USED TO REMEMBER WHICH YPLT'S HAVE ALREADY BEEN
C             COMPUTED.
C             OUTPUT, INTEGER*4
C     YINC  - M-VECTOR OF PLOT INCREMENTS.
C             OUTPUT, INTEGER*4
C     WORK  - WORKSPACE OF LENGTH M*M
C             REAL*8
C
      SUBROUTINE PCORD
     &(ISW,YPLT,Y,X,M,L,NGRD,S,MUSTAT,CSTAT,PSTAT,IY,YINC,WORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 YPLT(M),Y(M),X(M*L)
      REAL*8 MUSTAT(M),CSTAT(M*M),PSTAT(M*M),YINC(M)
      REAL*8 WORK(1)
      INTEGER*4 IY(M)
      SAVE
C
      IF(ISW.EQ.0) GO TO 10
      IF(ISW.EQ.1) GO TO 100
C
C     INITIALIZATION
C
10    CONTINUE
      DO 20 I=1,M
      IY(I)=-NGRD
20    YINC(I)=0.D0
      IY(1)=IY(1)-1
      DO 25 J=1,M
      DO 25 I=1,M
25    YINC(I)=YINC(I)+CSTAT(M*(J-1)+I)**2
      FNGRD=NGRD
      DO 30 I=1,M
30    YINC(I)=S*DSQRT(YINC(I))/FNGRD
      DO 35 J=1,M
      DO 35 I=1,J
35    WORK(J*(J-1)/2+I)=CSTAT(M*(J-1)+I)
      CALL DRINV(WORK,M)
      DO 40 J=1,M
      DO 40 I=1,M
      PSTAT(M*(J-1)+I)=0.D0
40    IF(I.LE.J) PSTAT(M*(J-1)+I)=WORK(J*(J-1)/2+I)
      DO 50 J=1,L
      DO 50 I=1,M
50    WORK(M*(J-1)+I)=X(M*(J-1)+I)-MUSTAT(I)
      DO 60 J=1,L
60    CALL DGMPRD(PSTAT,WORK(M*(J-1)+1),X(M*(J-1)+1),M,M,1)
      RETURN
C
C     COMPUTE YPLT AND Y
C
100   CONTINUE
      DO 110 I=1,M
      IF(IY(I).LT.NGRD) THEN
         IY(I)=IY(I)+1
         GO TO 120
      ELSE
         IY(I)=-NGRD
      ENDIF
110   CONTINUE
120   CONTINUE
      DO 130 I=1,M
      FIY=IY(I)
      WORK(I)=FIY*YINC(I)
130   YPLT(I)=MUSTAT(I)+WORK(I)
      CALL DGMPRD(PSTAT,WORK,Y,M,M,1)
      RETURN
      END
