*     G05PMF Example Program Text
*     Mark 22 Revised. NAG Copyright 2008.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NMAX, NFMAX, PMAX, MSTATE, MSEED, NSMAX
      PARAMETER        (NMAX=200,NFMAX=12,PMAX=5,MSTATE=633,MSEED=1,
     +                 NSMAX=100)
      DOUBLE PRECISION ONE, TWO
      PARAMETER        (ONE=1.0D0,TWO=2.0D0)
*     .. Local Scalars ..
      DOUBLE PRECISION AD, ALPHA, DV, TMP, VAR, Z
      INTEGER          EN, GENID, I, IFAIL, ITYPE, IVAL, J, K, LSEED,
     +                 LSTATE, MODE, N, NF, NSIM, P, SMODE, SUBID
*     .. Local Arrays ..
      DOUBLE PRECISION BLIM(2,NFMAX), BSIM(NSMAX,NFMAX), E(1),
     +                 FSE(NFMAX), FV(NFMAX), GLIM(2,NFMAX),
     +                 GSIM(NSMAX,NFMAX), INIT(PMAX+2), PARAM(4), Q(2),
     +                 R(PMAX+13), RES(NMAX), TSIM1(NFMAX),
     +                 TSIM2(NFMAX), Y(NMAX), YHAT(NMAX)
      INTEGER          SEED(MSEED), STATE(MSTATE)
*     .. External Functions ..
      DOUBLE PRECISION G01FAF
      EXTERNAL         G01FAF
*     .. External Subroutines ..
      EXTERNAL         G01AMF, G05KFF, G05PMF, G13AMF
*     .. Executable Statements ..
      CONTINUE
      WRITE (NOUT,*) 'G05PMF Example Program Results'
*     Initialize the seed to a repeatable sequence
      SEED(1) = 1762543
*     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) THEN
         WRITE (NOUT,99995) IFAIL
         GO TO 140
      END IF
*     Skip headings in data file
      READ (NIN,*)
*     Read in the initial arguments and check array sizes
      READ (NIN,*) MODE, ITYPE, N, NF, NSIM, ALPHA
      IF (N.GT.NMAX .OR. NF.GT.NFMAX .OR. NSIM.GT.NSMAX) THEN
         WRITE (NOUT,99991)
         GO TO 140
      END IF
      READ (NIN,*) (Y(I),I=1,N)
*     Read in the ITYPE dependent arguments (skipping headings)
      IF (ITYPE.EQ.1) THEN
         READ (NIN,*) PARAM(1)
         P = 0
         IVAL = 1
      ELSE IF (ITYPE.EQ.2) THEN
         READ (NIN,*) PARAM(1), PARAM(2)
         P = 0
         IVAL = 2
      ELSE IF (ITYPE.EQ.3) THEN
         READ (NIN,*) PARAM(1), PARAM(2), PARAM(3)
         P = 0
         IVAL = 2
      ELSE
         READ (NIN,*) PARAM(1), PARAM(2), PARAM(3), PARAM(4), P
         IVAL = P + 2
         IF (P.GT.PMAX) THEN
            WRITE (NOUT,99991)
            GO TO 140
         END IF
      END IF
*     Read in the MODE dependent arguments (skipping headings)
      IF (MODE.EQ.0) THEN
*        User supplied initial values
         READ (NIN,*) (INIT(I),I=1,IVAL)
      ELSE IF (MODE.EQ.1) THEN
*        Continuing from a previously saved R
         READ (NIN,*) (R(I),I=1,P+13)
      ELSE IF (MODE.EQ.2) THEN
*        Initial values calculated from first K observations
         READ (NIN,*) K
      END IF
*
*     Fit a smoothing model (parameter R in G05PMF and STATE in G13AMF
*     are in the same format)
      IFAIL = 1
      CALL G13AMF(MODE,ITYPE,P,PARAM,N,Y,K,INIT,NF,FV,FSE,YHAT,RES,DV,
     +            AD,R,IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,99994) IFAIL
         GO TO 140
      END IF
*
*     Simulate forecast values from the model, and don't update R
      SMODE = 2
      VAR = DV*DV
      EN = N
*     Simulate NSIM forecasts
      DO 40 I = 1, NSIM
*        Simulations assuming gaussian errors
         IFAIL = 1
         CALL G05PMF(SMODE,NF,ITYPE,P,PARAM,INIT,VAR,R,STATE,E,0,TSIM1,
     +               IFAIL)
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99993) IFAIL
            GO TO 140
         END IF
*        Bootstrapping errors
         IFAIL = 1
         CALL G05PMF(SMODE,NF,ITYPE,P,PARAM,INIT,0.0D0,R,STATE,RES,EN,
     +               TSIM2,IFAIL)
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99993) IFAIL
            GO TO 140
         END IF
*        Copy and transpose the simulated values
         DO 20 J = 1, NF
            GSIM(I,J) = TSIM1(J)
            BSIM(I,J) = TSIM2(J)
   20    CONTINUE
   40 CONTINUE
*
*     Calculate CI based on the quantiles for each simulated forecast
      Q(1) = ALPHA/TWO
      Q(2) = ONE - Q(1)
      DO 60 I = 1, NF
         IFAIL = 1
         CALL G01AMF(NSIM,GSIM(1,I),2,Q,GLIM(1,I),IFAIL)
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99992) IFAIL
            GO TO 140
         END IF
         IFAIL = 1
         CALL G01AMF(NSIM,BSIM(1,I),2,Q,BLIM(1,I),IFAIL)
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99992) IFAIL
            GO TO 140
         END IF
   60 CONTINUE
*
*     Display the forecast values and associated prediction intervals
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Initial values used:'
      DO 80 I = 1, IVAL
         WRITE (NOUT,99998) I, INIT(I)
   80 CONTINUE
      WRITE (NOUT,*)
      WRITE (NOUT,99999) 'Mean Deviation     = ', DV
      WRITE (NOUT,99999) 'Absolute Deviation = ', AD
      WRITE (NOUT,*)
      WRITE (NOUT,*) '         Observed      1-Step'
      WRITE (NOUT,*) ' Period   Values      Forecast      Residual'
      WRITE (NOUT,*)
      DO 100 I = 1, N
         WRITE (NOUT,99998) I, Y(I), YHAT(I), RES(I)
  100 CONTINUE
      WRITE (NOUT,*)
      WRITE (NOUT,*) '                                     '//
     +  '      Simulated CI         Simulated CI'
      WRITE (NOUT,*) 'Obs.  Forecast      Estimated CI   '//
     +  '     (Gaussian Errors)    (Bootstrap Errors)'
      Z = G01FAF('L',Q(2),IFAIL)
      DO 120 I = 1, NF
         TMP = Z*FSE(I)
         WRITE (NOUT,99997) N + I, FV(I), FV(I) - TMP, FV(I) + TMP,
     +     GLIM(1,I), GLIM(2,I), BLIM(1,I), BLIM(2,I)
  120 CONTINUE
      WRITE (NOUT,99996) 100.0D0*(ONE-ALPHA), '% CIs were produced'
*
  140 CONTINUE
*
99999 FORMAT (A,E12.4)
99998 FORMAT (I4,3(1X,F12.3))
99997 FORMAT (I3,7(1X,F10.3))
99996 FORMAT (1X,F5.1,A)
99995 FORMAT (1X,' ** G05KFF returned with IFAIL = ',I5)
99994 FORMAT (1X,' ** G13AMF returned with IFAIL = ',I5)
99993 FORMAT (1X,' ** G05PMF returned with IFAIL = ',I5)
99992 FORMAT (1X,' ** G01AMF returned with IFAIL = ',I5)
99991 FORMAT (1X,' ** Problem size too large, increase array limits')
      END