PROGRAM g05pmfe

!      G05PMF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g01amf, g01faf, g05kff, g05pmf, g13amf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: lseed = 1, nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: ad, alpha, dv, tmp, var, z
       INTEGER                         :: en, genid, i, ifail, itype, k, le,   &
                                          linit, lparam, lr, lstate, mode, n,  &
                                          nf, nsim, p, smode, subid
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: blim(:,:), bsim(:,:), e(:), fse(:),  &
                                          fv(:), glim(:,:), gsim(:,:),         &
                                          init(:), param(:), r(:), res(:),     &
                                          tsim1(:), tsim2(:), y(:), yhat(:)
       REAL (KIND=nag_wp)              :: q(2)
       INTEGER                         :: seed(lseed)
       INTEGER, ALLOCATABLE            :: state(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'G05PMF Example Program Results'
       WRITE (nout,*)

!      Skip heading in data file
       READ (nin,*)

!      Read in the base generator information and seed
       READ (nin,*) genid, subid, seed(1)

!      Initial call to initialiser to get size of STATE array
       lstate = 0
       ALLOCATE (state(lstate))
       ifail = 0
       CALL g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!      Reallocate STATE
       DEALLOCATE (state)
       ALLOCATE (state(lstate))

!      Initialize the generator to a repeatable sequence
       ifail = 0
       CALL g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!      Read in the initial arguments and check array sizes
       READ (nin,*) mode, itype, n, nf, nsim, alpha

       SELECT CASE (itype)
       CASE (1)
          lparam = 1
          p = 0
          linit = 1
       CASE (2)
          lparam = 2
          p = 0
          linit = 2
       CASE (3)
          lparam = 3
          p = 0
          linit = 2
       CASE DEFAULT
          lparam = 4

!         Read in seasonal order
          READ (nin,*) p

          linit = p + 2
       END SELECT
       lr = 13 + p
!      Not using the E array for the bootstrap
       le = 0
       ALLOCATE (param(lparam),init(linit),r(lr),e(le),fv(nf),fse(nf),yhat(n), &
          res(n),blim(2,nf),glim(2,nf),tsim1(nf),tsim2(nf),gsim(nsim,nf), &
          bsim(nsim,nf),y(n))

!      Read in series to be smoothed
       READ (nin,*) y(1:n)

!      Read in parameters
       READ (nin,*) param(1:lparam)

!      Read in the MODE dependent arguments (skipping headings)
       SELECT CASE (mode)
       CASE (0)
!         User supplied initial values
          READ (nin,*) init(1:linit)
       CASE (1)
!         Continuing from a previously saved R
          READ (nin,*) r(1:(p+13))
       CASE (2)
!         Initial values calculated from first K observations
          READ (nin,*) k
       END SELECT

!      Fit a smoothing model (parameter R in G05PMF and STATE in G13AMF
!      are in the same format)
       ifail = 0
       CALL g13amf(mode,itype,p,param,n,y,k,init,nf,fv,fse,yhat,res,dv,ad,r, &
          ifail)

!      Simulate forecast values from the model, and don't update R
       smode = 2
       var = dv*dv

!      Simulate NSIM forecasts
       DO i = 1, nsim
!         Not using E array for gaussian errors
          en = 0

!         Simulations assuming gaussian errors
          ifail = 0
          CALL g05pmf(smode,nf,itype,p,param,init,var,r,state,e,en,tsim1, &
             ifail)

!         For bootstrapping error, we are using RES from call to G13AMF as the
!         errors, and length of RES is N
          en = n

!         Bootstrapping errors
          ifail = 0
          CALL g05pmf(smode,nf,itype,p,param,init,0.0E0_nag_wp,r,state,res,en, &
             tsim2,ifail)

!         Copy and transpose the simulated values
          gsim(i,1:nf) = tsim1(1:nf)
          bsim(i,1:nf) = tsim2(1:nf)
       END DO

!      Calculate CI based on the quantiles for each simulated forecast
       q(1) = alpha/2.0E0_nag_wp
       q(2) = 1.0E0_nag_wp - q(1)
       DO i = 1, nf
          ifail = 0
          CALL g01amf(nsim,gsim(1,i),2,q,glim(1,i),ifail)
          ifail = 0
          CALL g01amf(nsim,bsim(1,i),2,q,blim(1,i),ifail)
       END DO

!      Display the forecast values and associated prediction intervals
       WRITE (nout,*) 'Initial values used:'
       WRITE (nout,99998) init(1:linit)
       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,*)
       WRITE (nout,99997) (i,y(i),yhat(i),res(i),i=1,n)
       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 i = 1, nf
          tmp = z*fse(i)
          WRITE (nout,99996) n + i, fv(i), fv(i) - tmp, fv(i) + tmp, &
             glim(1,i), glim(2,i), blim(1,i), blim(2,i)
       END DO
       WRITE (nout,99995) 100.0E0_nag_wp*(1.0E0_nag_wp-alpha), &
          '% CIs were produced'

99999  FORMAT (A,E12.4)
99998  FORMAT (F12.3)
99997  FORMAT (I4,1X,F12.3,1X,F12.3,1X,F12.3)
99996  FORMAT (I3,7(1X,F10.3))
99995  FORMAT (1X,F5.1,A)
    END PROGRAM g05pmfe