*     G13FEF Example Program Text
*     Mark 22 Revised. NAG Copyright 2006.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          NPARMX, LDCOVR, NUM, MSTATE, MSEED, NUM1, NREGMX,
     +                 LRVEC
      PARAMETER        (NPARMX=10,LDCOVR=NPARMX,NUM=2000,MSTATE=200,
     +                 MSEED=1,NUM1=3000,NREGMX=10,LRVEC=40)
      DOUBLE PRECISION ZERO
      PARAMETER        (ZERO=0.0D0)
*     .. Local Scalars ..
      DOUBLE PRECISION FAC1, GAMMA, HP, LGF, MEAN, TOL, XTERM
      INTEGER          D, DF, GENID, I, IFAIL, IP, IQ, J, K, LDX, LSEED,
     +                 LSTATE, LWORK, MAXIT, MN, NPAR, NPAR2, NREG, NT,
     +                 SUBID
      LOGICAL          FCALL
      CHARACTER        DIST
*     .. Local Arrays ..
      DOUBLE PRECISION BX(10), COVR(LDCOVR,NPARMX), CVAR(100),
     +                 ETM(NUM1), HT(NUM1+10), HTM(NUM1), PARAM(NPARMX),
     +                 RVEC(LRVEC), SC(NPARMX), SE(NPARMX),
     +                 THETA(NPARMX), WORK(NUM1*3+NPARMX+NREGMX*NUM1+20*
     +                 20+1), X(NUM1,10), YT(NUM1+10)
      INTEGER          SEED(MSEED), STATE(MSTATE)
      LOGICAL          COPTS(2)
*     .. External Subroutines ..
      EXTERNAL         G05KFF, G05PFF, G13FEF, G13FFF
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE, SIN
*     .. Executable Statements ..
      WRITE (NOUT,*) 'G13FEF Example Program Results'
      WRITE (NOUT,*)
*     Initialize the seed
      SEED(1) = 111
*     GENID and SUBID identify the base generator
      GENID = 1
      SUBID = 1
*     Initialize the generator to a repeatable sequence
      LSTATE = MSTATE
      LSEED = MSEED
      IFAIL = -1
      CALL G05KFF(GENID,SUBID,SEED,LSEED,STATE,LSTATE,IFAIL)
*
      IF (IFAIL.NE.0) GO TO 200
*
*     Initialize the time dependent, exogenous, variables
      LDX = NUM1
      GAMMA = 0.1D0
      BX(1) = 1.5D0
      BX(2) = 2.5D0
      BX(3) = 3.0D0
      DO 20 I = 1, NUM
         FAC1 = DBLE(I)*0.01D0
         X(I,2) = 0.01D0 + 0.7D0*SIN(FAC1)
         X(I,1) = 0.5D0 + FAC1*0.1D0
         X(I,3) = 1.0D0
   20 CONTINUE
*
      DO 180 D = 1, 2
         IF (D.EQ.1) THEN
*           Errors distributed as a normal distribution
            WRITE (NOUT,*)
            WRITE (NOUT,*) 'Normal distribution'
            WRITE (NOUT,*)
            DIST = 'N'
*           Initialize the series parameters required
*           for generation of some test data
            MEAN = 4.0D0
            MN = 1
            NREG = 2
            GAMMA = 0.1D0
            IP = 1
            IQ = 1
            PARAM(1) = 0.4D0
            PARAM(2) = 0.1D0
            PARAM(3) = 0.7D0
         ELSE
*           Errors distributed as a students t distribution
            WRITE (NOUT,*)
            WRITE (NOUT,*) 'Students T distribution'
            WRITE (NOUT,*)
            DIST = 'T'
*           Initialize the series parameters required
*           for generation of some test data
            MEAN = 3.0D0
            DF = 5
            MN = 1
            NREG = 2
            GAMMA = 0.09D0
            IP = 1
            IQ = 1
            PARAM(1) = 0.05D0
            PARAM(2) = 0.1D0
            PARAM(3) = 0.8D0
         END IF
*
*        Start of data generation
*        Generate the errors
         FCALL = .TRUE.
         CALL G05PFF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,FCALL,RVEC,
     +               LRVEC,STATE,IFAIL)
*
*        Calculate the observed data
         DO 60 I = 1, NUM
            XTERM = ZERO
            DO 40 K = 1, NREG
               XTERM = XTERM + X(I,K)*BX(K)
   40       CONTINUE
            IF (MN.EQ.1) THEN
               YT(I) = MEAN + XTERM + YT(I)
            ELSE
               YT(I) = XTERM + YT(I)
            END IF
   60    CONTINUE
*        End of data generation
*
*        Calculate the number of parameters
         NPAR = IQ + IP + 1
         NPAR2 = NPAR + MN + NREG + 1
         IF (DIST.EQ.'T') NPAR2 = NPAR2 + 1
         LWORK = NREG*NUM + 3*NUM + NPAR2 + 403
*
*        Set the model fitting options
         COPTS(1) = .TRUE.
         COPTS(2) = .TRUE.
         MAXIT = 50
         TOL = 1.0D-5
*
*        Generate initial values for the alpha_i and beta_j parameters
         DO 80 I = 1, NPAR
            THETA(I) = PARAM(I)*0.5D0
   80    CONTINUE
*        Generate initial values for the asymmetry parameter, gamma
         THETA(NPAR+1) = GAMMA*0.5D0
*        Generate initial values for DF estimate
         J = NPAR + 2
         IF (DIST.EQ.'T') THEN
            THETA(J) = DF*0.65D0
            J = J + 1
         END IF
*        Generate the initial values for b_0
         IF (MN.EQ.1) THEN
            THETA(J) = MEAN*0.5D0
            J = J + 1
         END IF
*        Generate the initial values for the linear regression
*        coefficients b_i, if required
         IF ( .NOT. COPTS(2)) THEN
            DO 100 I = 1, NREG
               THETA(I+J-1) = BX(I)*0.5D0
  100       CONTINUE
         END IF
*
*        Call the analysis routine
         IFAIL = -1
         CALL G13FEF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,NPAR2,THETA,SE,SC,
     +               COVR,LDCOVR,HP,ETM,HTM,LGF,COPTS,MAXIT,TOL,WORK,
     +               LWORK,IFAIL)
*
*        Output the results
         WRITE (NOUT,*) '          Parameter       Standard',
     +     '        Correct'
         WRITE (NOUT,*) '          estimates       errors  ',
     +     '        values'
*        Output the coefficient alpha_0
         WRITE (NOUT,99999) THETA(1), SE(1), PARAM(1)
*        Output the coefficients alpha_i
         DO 120 I = 2, IQ
            WRITE (NOUT,99999) THETA(I), SE(I), PARAM(I)
  120    CONTINUE
*        Output the coefficients beta_j
         DO 140 I = IQ + 1, NPAR
            WRITE (NOUT,99999) THETA(I), SE(I), PARAM(I)
  140    CONTINUE
*        Output the estimated asymmetry parameter, gamma
         WRITE (NOUT,99999) THETA(NPAR+1), SE(NPAR+1), GAMMA
         J = NPAR + 2
*        Output the estimated degrees of freedom, df
         IF (DIST.EQ.'T') THEN
            WRITE (NOUT,99999) THETA(J), SE(J), DBLE(DF)
            J = J + 1
         END IF
*        Output the estimated mean term, b_0
         IF (MN.EQ.1) THEN
            WRITE (NOUT,99999) THETA(J), SE(J), MEAN
            J = J + 1
         END IF
*        Output the estimated linear regression coefficients, b_i
         DO 160 I = 1, NREG
            WRITE (NOUT,99999) THETA(I+J-1), SE(I+J-1), BX(I)
  160    CONTINUE
*        Calculate the volatility forecast
         NT = 4
         CALL G13FFF(NUM,NT,IP,IQ,THETA,GAMMA,CVAR,HTM,ETM,IFAIL)
*        Display the volatility forecast
         WRITE (NOUT,*)
         WRITE (NOUT,99998) 'Volatility forecast = ', CVAR(NT)
         WRITE (NOUT,*)
  180 CONTINUE
*
  200 CONTINUE
*
99999 FORMAT (1X,3F16.4)
99998 FORMAT (1X,A,F12.4)
      END