Program g05pmfe

!     G05PMF Example Program Text

!     Mark 28.3 Release. NAG Copyright 2022.

!     .. 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 initializer 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