Program g13bafe

!     G13BAF Example Program Text

!     Mark 28.4 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: g13ajf, g13baf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: cx, cy, rms
      Integer                          :: i, idd, ifail, ifv, ii, ij, ipar,    &
                                          iqxd, ist, iw, nb, nmr, npar, nparx, &
                                          nst, nwa, nx, ny, pp, qp, sy
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:), fsd(:), fva(:), par(:),        &
                                          parx(:), st(:), w(:), wa(:), x(:),   &
                                          y(:)
      Integer                          :: isf(4), mrx(7)
      Integer, Allocatable             :: mr(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, min, mod
!     .. Executable Statements ..
      Write (nout,*) 'G13BAF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size
      Read (nin,*) nx

!     Read univariate ARIMA for series
      Read (nin,*) mrx(1:7)
      Read (nin,*) cx

!     Calculate number of backforecasts required
      iqxd = mrx(3) + mrx(6)*mrx(7)
      If (iqxd/=0) Then
        nmr = 14
      Else
        nmr = 7
      End If

!     Back forecasts will be stored in first IQXD elements
!     of Y, the series will be stored in last NX elements of
!     Y, so calculate start point for the series
      sy = iqxd + 1

!     Calculate length of series with back forecasts
      ny = nx + iqxd

      Allocate (y(ny),mr(nmr))

!     Read in the series into the end of Y
      Read (nin,*) y(sy:ny)

!     Get back forecasts if required
      If (iqxd/=0) Then

!       Calculate number of parameters in ARIMA model
        nparx = mrx(1) + mrx(3) + mrx(4) + mrx(6)

        ist = mrx(4) + mrx(7) + mrx(2) + mrx(5) + mrx(3) +                     &
          max(mrx(1),mrx(6)*mrx(7))
        ifv = max(1,iqxd)
        qp = mrx(6)*mrx(7) + mrx(3)
        pp = mrx(4)*mrx(7) + mrx(1)
        iw = 6*nx + 5*nparx + qp*(qp+11) + 3*pp + 7
        Allocate (parx(nparx),x(nx),st(ist),fva(ifv),fsd(ifv),w(iw))

!       Read in initial values
        Read (nin,*) parx(1:nparx)

!       Reverse series
        x(nx:1:-1) = y(sy:ny)

!       Possible sign reversal for ARIMA constant
        idd = mrx(2) + mrx(5)
        If (mod(idd,2)/=0) Then
          cx = -cx
        End If

!       Calculate back forecasts
        ifail = 0
        Call g13ajf(mrx,parx,nparx,cx,1,x,nx,rms,st,ist,nst,iqxd,fva,fsd,ifv,  &
          isf,w,iw,ifail)

!       Move back forecasts into Y, in reverse order
        y(1:iqxd) = fva(iqxd:1:-1)

!       Reverse sign for ARIMA constant back again
        If (mod(idd,2)/=0) Then
          cx = -cx
        End If
      End If

!     Read model by which to filter series
      Read (nin,*) mr(1:7)

!     Calculate NPAR
      ipar = mr(1) + mr(3) + mr(4) + mr(6)
      npar = ipar + nparx

      Allocate (par(npar))

!     Read in initial parameter values
      Read (nin,*) par(1:ipar)

      If (iqxd/=0) Then
!       Move ARIMA series into MR
        mr(8:14) = mrx(1:7)

!       Move parameters of ARIMA for Y into PAR
        par((ipar+1):(ipar+nparx)) = parx(1:nparx)
      End If

!     Move constant
      cy = cx

!     Set parameters for call to filter routine G13BAF
      If (nmr==14) Then
        nwa = mr(3) + mr(6)*mr(7) + mr(8) + mr(9) + (mr(11)+mr(12))*mr(14)
        nwa = nwa*(nwa+2)
        nb = ny + max(mr(3)+mr(6)*mr(7),mr(1)+mr(2)+(mr(4)+mr(5))*mr(7))
      Else
        nwa = 1
        nb = ny
      End If
      Allocate (wa(nwa),b(nb))

!     Filter series by call to G13BAF
      ifail = 0
      Call g13baf(y,ny,mr,nmr,par,npar,cy,wa,nwa,b,nb,ifail)

!     Display results
      If (iqxd/=0) Then
        Write (nout,*) '                 Original        Filtered'
        Write (nout,*) 'Backforecasts    y-series         series'
        ij = -iqxd
        Do i = 1, iqxd
          Write (nout,99999) ij, y(i), b(i)
          ij = ij + 1
        End Do
        Write (nout,*)
      End If
      Write (nout,*)                                                           &
        '       Filtered        Filtered        Filtered        Filtered'
      Write (nout,*)                                                           &
        '        series          series          series          series'
      Do i = iqxd + 1, ny, 4
        Write (nout,99998)(ii-iqxd,b(ii),ii=i,min(ny,i+3))
      End Do

99999 Format (1X,I8,F17.4,F15.4)
99998 Format (1X,I5,F9.4,I7,F9.4,I7,F9.4,I7,F9.4)
    End Program g13bafe