PROGRAM g13bbfe

!      G13BBF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

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