Program g05pmfe ! G05PMF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: g01amf, g01faf, g05kff, g05pmf, g13amf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: lseed = 1, nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: ad, alpha, dv, tmp, var, z Integer :: en, genid, i, ifail, itype, k, le, & linit, lparam, lr, lstate, mode, n, & nf, nsim, p, smode, subid ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: blim(:,:), bsim(:,:), e(:), fse(:), & fv(:), glim(:,:), gsim(:,:), & init(:), param(:), r(:), res(:), & tsim1(:), tsim2(:), y(:), yhat(:) Real (Kind=nag_wp) :: q(2) Integer :: seed(lseed) Integer, Allocatable :: state(:) ! .. Executable Statements .. Write (nout,*) 'G05PMF 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 initial arguments and check array sizes Read (nin,*) mode, itype, n, nf, nsim, alpha Select Case (itype) Case (1) lparam = 1 p = 0 linit = 1 Case (2) lparam = 2 p = 0 linit = 2 Case (3) lparam = 3 p = 0 linit = 2 Case Default lparam = 4 ! Read in seasonal order Read (nin,*) p linit = p + 2 End Select lr = 13 + p ! Not using the E array for the bootstrap le = 0 Allocate (param(lparam),init(linit),r(lr),e(le),fv(nf),fse(nf),yhat(n), & res(n),blim(2,nf),glim(2,nf),tsim1(nf),tsim2(nf),gsim(nsim,nf), & bsim(nsim,nf),y(n)) ! Read in series to be smoothed Read (nin,*) y(1:n) ! Read in parameters Read (nin,*) param(1:lparam) ! Read in the MODE dependent arguments (skipping headings) Select Case (mode) Case (0) ! User supplied initial values Read (nin,*) init(1:linit) Case (1) ! Continuing from a previously saved R Read (nin,*) r(1:(p+13)) Case (2) ! Initial values calculated from first K observations Read (nin,*) k End Select ! Fit a smoothing model (parameter R in G05PMF and STATE in G13AMF ! are in the same format) ifail = 0 Call g13amf(mode,itype,p,param,n,y,k,init,nf,fv,fse,yhat,res,dv,ad,r, & ifail) ! Simulate forecast values from the model, and don't update R smode = 2 var = dv*dv ! Simulate NSIM forecasts Do i = 1, nsim ! Not using E array for gaussian errors en = 0 ! Simulations assuming gaussian errors ifail = 0 Call g05pmf(smode,nf,itype,p,param,init,var,r,state,e,en,tsim1,ifail) ! For bootstrapping error, we are using RES from call to G13AMF as the ! errors, and length of RES is N en = n ! Bootstrapping errors ifail = 0 Call g05pmf(smode,nf,itype,p,param,init,0.0E0_nag_wp,r,state,res,en, & tsim2,ifail) ! Copy and transpose the simulated values gsim(i,1:nf) = tsim1(1:nf) bsim(i,1:nf) = tsim2(1:nf) End Do ! Calculate CI based on the quantiles for each simulated forecast q(1) = alpha/2.0E0_nag_wp q(2) = 1.0E0_nag_wp - q(1) Do i = 1, nf ifail = 0 Call g01amf(nsim,gsim(1,i),2,q,glim(1,i),ifail) ifail = 0 Call g01amf(nsim,bsim(1,i),2,q,blim(1,i),ifail) End Do ! Display the forecast values and associated prediction intervals Write (nout,*) 'Initial values used:' Write (nout,99998) init(1:linit) Write (nout,*) Write (nout,99999) 'Mean Deviation = ', dv Write (nout,99999) 'Absolute Deviation = ', ad Write (nout,*) Write (nout,*) ' Observed 1-Step' Write (nout,*) ' Period Values Forecast Residual' Write (nout,*) Write (nout,99997)(i,y(i),yhat(i),res(i),i=1,n) Write (nout,*) Write (nout,*) ' ' // & ' Simulated CI Simulated CI' Write (nout,*) 'Obs. Forecast Estimated CI ' // & ' (Gaussian Errors) (Bootstrap Errors)' z = g01faf('L',q(2),ifail) Do i = 1, nf tmp = z*fse(i) Write (nout,99996) n + i, fv(i), fv(i) - tmp, fv(i) + tmp, glim(1,i), & glim(2,i), blim(1,i), blim(2,i) End Do Write (nout,99995) 100.0E0_nag_wp*(1.0E0_nag_wp-alpha), & '% CIs were produced' 99999 Format (A,E12.4) 99998 Format (F12.3) 99997 Format (I4,1X,F12.3,1X,F12.3,1X,F12.3) 99996 Format (I3,7(1X,F10.3)) 99995 Format (1X,F5.1,A) End Program g05pmfe