Program g05pjfe ! G05PJF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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