Program g13bbfe

!     G13BBF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: g13ajf, g13bbf, 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, iwa, nb, nmr, npar,   &
                                          nparx, nst, 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,*) 'G13BBF Example Program Results'
      Write (nout,*)

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

!     Read in 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 = 10
      Else
        nmr = 3
      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 series
      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,0,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:3)

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

      Allocate (par(npar))

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

      If (iqxd/=0) Then
!       Move ARIMA for series into MR
        mr(4:10) = 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 G13BBF
      If (nmr==10) Then
        iwa = mr(3) + mr(4) + mr(5) + (mr(7)+mr(8))*mr(10)
        iwa = npar + iwa*(iwa+2)
        nb = ny + max(mr(1)+mr(2),mr(3))
      Else
        iwa = mr(1) + npar
        nb = ny
      End If

      Allocate (wa(iwa),b(nb))

!     Filter series by call to G13BBF
      ifail = 0
      Call g13bbf(y,ny,mr,nmr,par,npar,cy,wa,iwa,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.1,F16.1)
99998 Format (1X,I5,F10.1,I6,F10.1,I6,F10.1,I6,F10.1)
    End Program g13bbfe