Example description
!   G13DKF Example Program Text
!   Mark 27.1 Release. NAG Copyright 2020.

    Module g13dkfe_mod

!     G13DKF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fprint
!     .. Parameters ..
      Integer, Parameter, Public       :: iset = 1, nin = 5, nout = 6
    Contains
      Subroutine fprint(k,nm,lmax,predz,sefz,ldsefz,nout)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: k, ldsefz, lmax, nm, nout
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: predz(ldsefz,lmax),                 &
                                          sefz(ldsefz,lmax)
!       .. Local Scalars ..
        Integer                        :: i, i2, j, l, l2, loop
!       .. Intrinsic Procedures ..
        Intrinsic                      :: min, mod
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*) ' FORECAST SUMMARY TABLE'
        Write (nout,*) ' ----------------------'
        Write (nout,*)
        Write (nout,99999) ' Forecast origin is set at t = ', nm
        Write (nout,*)
        loop = lmax/5
        If (mod(lmax,5)/=0) Then
          loop = loop + 1
        End If
        Do j = 1, loop
          i2 = (j-1)*5
          l2 = min(i2+5,lmax)
          Write (nout,99998) 'Lead Time ', (i,i=i2+1,l2)
          Write (nout,*)
          i = 1
          Write (nout,99997) 'Series ', i, ' : Forecast    ',                  &
            (predz(1,l),l=i2+1,l2)
          Write (nout,99996) ' : Standard Error ', (sefz(1,l),l=i2+1,l2)
          Do i = 2, k
            Write (nout,99997) 'Series ', i, ' : Forecast    ',                &
              (predz(i,l),l=i2+1,l2)
            Write (nout,99996) ' : Standard Error ', (sefz(i,l),l=i2+1,l2)
          End Do
          Write (nout,*)
        End Do

        Return

99999   Format (1X,A,I4)
99998   Format (1X,A,12X,5I10)
99997   Format (1X,A,I2,A,5F10.2)
99996   Format (10X,A,4(F7.2,3X),F7.2)
      End Subroutine fprint
    End Module g13dkfe_mod
    Program g13dkfe

!     G13DKF Example Main Program

!     .. Use Statements ..
      Use g13dkfe_mod, Only: fprint, iset, nin, nout
      Use nag_library, Only: g13ddf, g13djf, g13dkf, g13dlf, nag_wp, x04abf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: cgetol, rlogl
      Integer                          :: d, i, ifail, ip, iprint, iq, ishow,  &
                                          k, kmax, ldcm, liwork, lmax, lpar,   &
                                          lref, lwork, m, maxcal, mlast, n,    &
                                          nadv, nd, niter, r, tddelta, tdv
      Logical                          :: exact, meanl
      Character (1)                    :: mean
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: cm(:,:), delta(:,:), g(:), par(:),   &
                                          predz(:,:), qq(:,:), ref(:),         &
                                          sefz(:,:), v(:,:), w(:,:), work(:),  &
                                          workl(:), z(:,:)
      Integer, Allocatable             :: id(:), iwork(:)
      Logical, Allocatable             :: parhld(:)
      Character (1), Allocatable       :: tr(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, maxval
!     .. Executable Statements ..
      Write (nout,*) 'G13DKF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size
      Read (nin,*) k, n

      Allocate (id(k))

!     Read in differencing
      Read (nin,*) id(1:k)

      d = maxval(id(1:k))
      tddelta = max(d,1)
      nd = n - d
      kmax = k
      Allocate (z(kmax,n),tr(k),delta(kmax,tddelta),w(kmax,nd),workl(k*n))

!     Read in series and the transformation flag
      Read (nin,*)(z(i,1:n),i=1,k)
      Read (nin,*) tr(1:k)

!     If required, read in delta
      If (d>0) Then
        Read (nin,*)(delta(i,1:id(i)),i=1,k)
      End If

!     Difference and / or transform series
      ifail = 0
      Call g13dlf(k,n,z,kmax,tr,id,delta,w,nd,workl,ifail)

!     Read in information on the VARMA
      Read (nin,*) ip, iq, mean, lmax

!     Calculate number of parameters for the VARMA
      lpar = (ip+iq)*k*k
      If (mean=='M' .Or. mean=='m') Then
        lpar = lpar + k
        meanl = .True.
      Else
        meanl = .False.
      End If

!     Read in control parameters
      Read (nin,*) iprint, cgetol, maxcal, ishow

!     Read in exact likelihood flag
      Read (nin,*) exact

      ldcm = lpar
      tdv = nd
      Allocate (par(lpar),parhld(lpar),qq(kmax,k),v(kmax,tdv),g(lpar),         &
        cm(ldcm,lpar))

!     Read in initial parameter estimates and free parameter flags
      Read (nin,*) par(1:lpar)
      Read (nin,*) parhld(1:lpar)

!     Read in initial values for covariance matrix Q
      Read (nin,*)(qq(i,1:i),i=1,k)

!     Set the advisory channel to NOUT for monitoring information
      If (iprint>=0 .Or. ishow/=0) Then
        nadv = nout
        Call x04abf(iset,nadv)
      End If

!     Fit a VARMA model
      ifail = -1
      Call g13ddf(k,nd,ip,iq,meanl,par,lpar,qq,kmax,w,parhld,exact,iprint,     &
        cgetol,maxcal,ishow,niter,rlogl,v,g,cm,ldcm,ifail)
      If (ifail/=0) Then
        If (ifail<4) Then
          Go To 100
        End If
      End If

      lref = (lmax-1)*k*k + 2*k*lmax + k
      r = max(ip,iq)
      lwork = max(k*r*(k*r+2),(ip+d+2)*k**2+(n+lmax)*k)
      liwork = k*max(ip,iq)
      Allocate (predz(kmax,lmax),sefz(kmax,lmax),ref(lref),work(lwork),        &
        iwork(liwork))

!     Forecast from VARMA
      ifail = 0
      Call g13djf(k,n,z,kmax,tr,id,delta,ip,iq,mean,par,lpar,qq,v,lmax,predz,  &
        sefz,ref,lref,work,lwork,iwork,liwork,ifail)

!     Display results
      Call fprint(k,n,lmax,predz,sefz,kmax,nout)

!     Update forecasts
      mlast = 0
d_lp: Do
        Read (nin,*,Iostat=ifail) m
        If (ifail/=0) Then
          Exit d_lp
        End If
        Read (nin,*,Iostat=ifail)(z(1:k,i),i=1,m)
        If (ifail/=0) Then
          Exit d_lp
        End If

!       Reallocate V if required
        If (tdv<m) Then
          Deallocate (v)
          Allocate (v(kmax,m))
        End If

!       Reallocate WORK if required
        If (lwork<k*m) Then
          Deallocate (work)
          Allocate (work(lwork))
        End If

!       Update forecast
        ifail = 0
        Call g13dkf(k,lmax,m,mlast,z,kmax,ref,lref,v,predz,sefz,work,ifail)

        Call fprint(k,n+mlast,lmax,predz,sefz,kmax,nout)
      End Do d_lp

100   Continue

    End Program g13dkfe