NAG Library Manual, Mark 28.4
```    Program g13bjfe

!     G13BJF Example Program Text

!     Mark 28.4 Release. NAG Copyright 2022.

!     .. Use Statements ..
Use nag_library, Only: g13bjf, nag_wp, x04caf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
Integer                          :: i, ifail, imwa, isttf, iwa, kfc,     &
kzef, ldparx, ldxxy, n, ncf, nev,    &
nfv, nis, npara, nparx, nser, nsttf
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: fsd(:), fva(:), para(:), parx(:,:),  &
rmsxy(:), sttf(:), wa(:), xxy(:,:)
Integer                          :: mr(7)
Integer, Allocatable             :: mrx(:,:), mt(:,:), mwa(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: max, sum
!     .. Executable Statements ..
Write (nout,*) 'G13BJF Example Program Results'
Write (nout,*)

!     Skip heading in data file

Read (nin,*) kfc, nev, nfv, nser, kzef

!     Number of input series
nis = nser - 1

Allocate (mt(4,nser))

!     Read in the orders for the output noise

Do i = 1, nis
End Do

!     Calculate NPARA
npara = 0
Do i = 1, nis
npara = npara + mt(2,i) + mt(3,i)
End Do
npara = npara + mr(1) + mr(3) + mr(4) + mr(6) + nser

!     Calculate array sizes
n = nev + nfv
ldxxy = n
ncf = 0
Do i = 1, nis
If (mt(4,i)>1) Then
ncf = sum(mt(1:3,i))
End If
End Do
isttf = mr(4)*mr(7) + mr(2) + mr(5)*mr(7) + mr(3) +                      &
max(mr(1),mr(6)*mr(7)) + ncf
iwa = 0
imwa = 0
Allocate (para(npara),xxy(ldxxy,nser),rmsxy(nser),mrx(7,nser),fva(nfv),  &
fsd(nfv),sttf(isttf),wa(iwa),mwa(imwa))

!     Read in multi-input model parameters

!     Read in the observed values for the input and output series

!     Read in the future values for the input series

If (nis>=1) Then
!       Read in residual variance of input series

!       Read in orders for input series ARIMA where available
!       (i.e. where residual variance is not zero)
ldparx = 0
Do i = 1, nis
If (rmsxy(i)/=0.0E0_nag_wp) Then
nparx = mrx(1,i) + mrx(3,i) + mrx(4,i) + mrx(6,i)
ldparx = max(ldparx,nparx)
Else
mrx(1:7,i) = 0
End If
End Do
Else
!       No input series
ldparx = 1
End If

Allocate (parx(ldparx,nser))

!     Read in parameters for each input series ARIMA
If (nis>0) Then
Do i = 1, nis
If (rmsxy(i)/=0.0E0_nag_wp) Then
nparx = mrx(1,i) + mrx(3,i) + mrx(4,i) + mrx(6,i)
If (nparx>0) Then
End If
End If
End Do
End If

ifail = 0
Call g13bjf(mr,nser,mt,para,npara,kfc,nev,nfv,xxy,ldxxy,kzef,rmsxy,mrx,  &
parx,ldparx,fva,fsd,sttf,isttf,nsttf,wa,iwa,mwa,imwa,ifail)

!     Display results
Write (nout,99999) 'After processing', nev, ' sets of observations'
Write (nout,99998) nsttf, ' values of the state set are derived'
Write (nout,*)
Write (nout,99997) sttf(1:nsttf)
Write (nout,*)
Write (nout,*) 'The residual mean square for the output'
Write (nout,99996) 'series is also derived and its value is',            &
rmsxy(nser)
Write (nout,*)
Write (nout,*) 'The forecast values and their standard errors are'
Write (nout,*)
Write (nout,*) '   I       FVA       FSD'
Write (nout,*)
Write (nout,99995)(i,fva(i),fsd(i),i=1,nfv)
Write (nout,*)
Flush (nout)
ifail = 0
Call x04caf('General',' ',n,nser,xxy,ldxxy,                              &
'The values of z(t) and n(t) are',ifail)
Write (nout,99994) 'The first ', nis,                                    &
' columns hold the z(t) and the last column the n(t)'

99999 Format (1X,A,I3,A)
99998 Format (1X,I3,A)
99997 Format (1X,6F10.4)
99996 Format (1X,A,F10.4)
99995 Format (1X,I4,F10.3,F10.4)
99994 Format (1X,A,I0,A)
End Program g13bjfe
```