Program g05pjfe
! G05PJF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
! .. 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 (nin,*)
! Read in the base generator information and seed
Read (nin,*) genid, subid, seed(1)
! Initial call to initialiser 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 (nin,*) n, nreal
! Read in the number of coefficients
Read (nin,*) k, ip, iq
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
Read (nin,*)(phi(ii+k*(j-1)),j=1,k)
End Do
End Do
! Read in the MA parameters
Do l = 1, iq
Do i = 1, k
ii = (l-1)*k*k + i
Read (nin,*)(theta(ii+k*(j-1)),j=1,k)
End Do
End Do
! Read in the means
Read (nin,*) xmean(1:k)
! Read in the variance / covariance matrix
Read (nin,*)(var(i,1:i),i=1,k)
! 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