```    Program g13befe

!     G13BEF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2017.

!     .. Use Statements ..
Use nag_library, Only: g13bef, nag_wp, x04abf, x04caf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: iset = 1, nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: d, s
Integer                          :: dp, i, ifail, imwa, inc, isttf, itc, &
iwa, kef, kfc, kpriv, kzef, kzsp,    &
ldcm, ldxxy, mx, nadv, ncd, nce,     &
ncf, ncg, ndf, ndv, nis, nit, npara, &
nser, nsttf, nxxy, qp, qx, smx
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: cm(:,:), para(:), res(:), sd(:),     &
sttf(:), wa(:), xxy(:,:)
Real (Kind=nag_wp)               :: zsp(4)
Integer                          :: mr(7)
Integer, Allocatable             :: mt(:,:), mwa(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: max, sum
!     .. Executable Statements ..
Write (nout,*) 'G13BEF Example Program Results'
Write (nout,*)

!     Skip heading in data file

Read (nin,*) kzef, kfc, nxxy, nser, kef, nit, kzsp, kpriv
If (kzsp/=0) Then
End If

!     Number of input series
nis = nser - 1

!     Set the advisory channel to NOUT for monitoring information
If (kpriv/=0) Then
End If

Allocate (mt(4,nser))

Do i = 1, nis
End Do

!     Calculate NPARA and various other quantities required
!     for calculate array sizes
npara = 0
ncg = 0
qx = 0
smx = 0
ncf = nser
inc = 1
Do i = 1, nis
npara = npara + mt(2,i) + mt(3,i)
If (mt(4,i)>1) Then
ncg = ncg + sum(mt(1:3,i))
If (mt(4,i)==3) Then
mx = max(mt(1,i)+mt(2,i),mt(3,i))
qx = max(qx,mx)
smx = smx + mx
End If
Else If (mt(4,i)==1 .And. kef==3) Then
If (mt(3,i)>0) Then
ncf = ncf + 1
End If
inc = inc + 1
End If
End Do
npara = npara + mr(1) + mr(3) + mr(4) + mr(6) + nser

!     Calculate size of arrays
isttf = mr(4)*mr(7) + mr(2) + mr(5)*mr(7) + mr(3) +                      &
max(mr(1),mr(6)*mr(7)) + ncg
ldxxy = nxxy
ldcm = npara
qp = mr(3) + mr(6)*mr(7)
dp = mr(2) + mr(5)*mr(7)
If (mr(3)>0 .And. kef>1) Then
inc = inc + 1
End If
If (kfc>0 .And. kef==3) Then
inc = inc + 1
End If
qx = qp
ncd = npara + kfc + smx
If (mr(1)>0) Then
ncf = ncf + inc
End If
If (mr(3)>0) Then
ncf = ncf + inc
End If
If (mr(4)>0) Then
ncf = ncf + inc
End If
If (mr(6)>0) Then
ncf = ncf + inc
End If
If (qx>0) Then
ncf = ncf + 1
End If
If (kfc>0) Then
ncf = ncf + 1
End If
ncd = ncd + qx
nce = nxxy + dp + 6*qx
iwa = 2*ncd**2 + nce*(ncf+4)
iwa = 2*iwa
imwa = 16*nser + 7*ncd + 3*npara + 3*kfc + 27

Allocate (xxy(ldxxy,nser),para(npara),sd(npara),cm(ldcm,npara),          &
res(nxxy),sttf(isttf),wa(iwa),mwa(imwa))

!     Read in rest of data

ifail = -1
Call g13bef(mr,nser,mt,para,npara,kfc,nxxy,xxy,ldxxy,kef,nit,kzsp,zsp,   &
itc,sd,cm,ldcm,s,d,ndf,kzef,res,sttf,isttf,nsttf,wa,iwa,mwa,imwa,      &
kpriv,ifail)
If (ifail/=0) Then
If (ifail/=8 .And. ifail/=9) Then
Go To 100
End If
End If

!     Display results
Write (nout,99999) 'The number of iterations carried out is', itc
Write (nout,*)
If (ifail/=8 .And. ifail/=9) Then
Write (nout,*)                                                         &
'Final values of the parameters and their standard deviations'
Write (nout,*)
Write (nout,*) '   I            PARA(I)                 SD'
Write (nout,*)
Write (nout,99998)(i,para(i),sd(i),i=1,npara)
Write (nout,*)
Flush (nout)
ifail = 0
Call x04caf('General',' ',npara,npara,cm,ldcm,                         &
'The correlation matrix is',ifail)
End If
Write (nout,*)
Write (nout,*) 'The residuals and the z and n values are'
Write (nout,*)
Write (nout,*) '   I         RES(I)          z(t)           n(t)'
Write (nout,*)
ndv = nxxy - mr(2) - mr(5)*mr(7)
Do i = 1, nxxy
If (i<=ndv) Then
Write (nout,99997) i, res(i), xxy(i,1:nser)
Else
Write (nout,99996) i, xxy(i,1:nser)
End If
End Do
If (mr(2)/=0 .Or. mr(5)/=0) Then
Write (nout,*)
Write (nout,*)                                                         &
'** Note that the residuals relate to differenced values **'
End If
Write (nout,*)
Write (nout,99995) 'The state set consists of', nsttf, ' values'
Write (nout,*)
Write (nout,99994) sttf(1:nsttf)
Write (nout,*)
Write (nout,99999) 'The number of degrees of freedom is', ndf

100   Continue

99999 Format (1X,A,I4)
99998 Format (1X,I4,2F20.6)
99997 Format (1X,I4,3F15.3)
99996 Format (1X,I4,F30.3,F15.3)
99995 Format (1X,A,I4,A)
99994 Format (1X,6F10.4)
End Program g13befe
```