NAG Library Manual, Mark 28.4
```    Program g05pjfe

!     G05PJF Example Program Text

!     Mark 28.4 Release. NAG Copyright 2022.

!     .. Use Statements ..
Use nag_library, Only: g05kff, g05pjf, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: lseed = 1, nin = 5, nout = 6
!     .. Local Scalars ..
Integer                          :: genid, i, ifail, ii, ip, iq, j, k,   &
k2, l, ldvar, ldx, lphi, lr, lstate, &
ltheta, mode, n, nreal, rn, subid
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: phi(:), r(:), theta(:), var(:,:),    &
x(:,:), xmean(:)
Integer                          :: seed(lseed)
Integer, Allocatable             :: state(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: max
!     .. Executable Statements ..
Write (nout,*) 'G05PJF Example Program Results'
Write (nout,*)

!     Skip heading in data file

!     Read in the base generator information and seed

!     Initial call to initializer to get size of STATE array
lstate = 0
Allocate (state(lstate))
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!     Reallocate STATE
Deallocate (state)
Allocate (state(lstate))

!     Initialize the generator to a repeatable sequence
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!     Read in the sample size and number of realizations

!     Read in the number of coefficients

k2 = k**2
rn = max(ip,iq)
l = k*(k+1)/2
If (ip>0) Then
l = l + (ip-1)*k2
End If
If (k>=6) Then
lr = (5*rn**2+1)*k2 + (4*rn+3) + 4
Else
lr = ((ip+iq)**2+1)*k2 + (4*(ip+iq)+3)*k + max(k*rn*(k*rn+2),k2*(ip+iq &
)**2+l*(l+3)+k2*(iq+1)) + 4
End If
lphi = ip*k*k
ltheta = iq*k*k
ldvar = k
ldx = k
Allocate (phi(lphi),theta(ltheta),var(ldvar,k),r(lr),x(ldx,n),xmean(k))

!     Read in the AR parameters
Do l = 1, ip
Do i = 1, k
ii = (l-1)*k*k + i
End Do
End Do

!     Read in the MA parameters
Do l = 1, iq
Do i = 1, k
ii = (l-1)*k*k + i
End Do
End Do

!     Read in the variance / covariance matrix

!     For the first realization we need to set up the reference vector
!     as well as generate the series
mode = 2

!     Generate NREAL realizations
d_lp: Do rn = 1, nreal

ifail = 0
Call g05pjf(mode,n,k,xmean,ip,phi,iq,theta,var,ldvar,r,lr,state,x,ldx, &
ifail)

!       Display the results
Write (nout,99999) ' Realization Number ', rn
Do i = 1, k
Write (nout,*)
Write (nout,99999) '  Series number ', i
Write (nout,*) '  -------------'
Write (nout,*)
Write (nout,99998) x(i,1:n)
End Do
Write (nout,*)

!       For subsequent realizations we use previous reference vector
mode = 3

End Do d_lp

99999 Format (1X,A,I0)
99998 Format (8(2X,F8.3))
End Program g05pjfe
```