* G05PMF Example Program Text * Mark 22 Revised. NAG Copyright 2008. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX, NFMAX, PMAX, MSTATE, MSEED, NSMAX PARAMETER (NMAX=200,NFMAX=12,PMAX=5,MSTATE=633,MSEED=1, + NSMAX=100) DOUBLE PRECISION ONE, TWO PARAMETER (ONE=1.0D0,TWO=2.0D0) * .. Local Scalars .. DOUBLE PRECISION AD, ALPHA, DV, TMP, VAR, Z INTEGER EN, GENID, I, IFAIL, ITYPE, IVAL, J, K, LSEED, + LSTATE, MODE, N, NF, NSIM, P, SMODE, SUBID * .. Local Arrays .. DOUBLE PRECISION BLIM(2,NFMAX), BSIM(NSMAX,NFMAX), E(1), + FSE(NFMAX), FV(NFMAX), GLIM(2,NFMAX), + GSIM(NSMAX,NFMAX), INIT(PMAX+2), PARAM(4), Q(2), + R(PMAX+13), RES(NMAX), TSIM1(NFMAX), + TSIM2(NFMAX), Y(NMAX), YHAT(NMAX) INTEGER SEED(MSEED), STATE(MSTATE) * .. External Functions .. DOUBLE PRECISION G01FAF EXTERNAL G01FAF * .. External Subroutines .. EXTERNAL G01AMF, G05KFF, G05PMF, G13AMF * .. Executable Statements .. CONTINUE WRITE (NOUT,*) 'G05PMF Example Program Results' * Initialize the seed to a repeatable sequence SEED(1) = 1762543 * GENID and SUBID identify the base generator. GENID = 1 SUBID = 1 * Initialize the generator to a repeatable sequence LSTATE = MSTATE LSEED = MSEED IFAIL = 1 CALL G05KFF(GENID,SUBID,SEED,LSEED,STATE,LSTATE,IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99995) IFAIL GO TO 140 END IF * Skip headings in data file READ (NIN,*) * Read in the initial arguments and check array sizes READ (NIN,*) MODE, ITYPE, N, NF, NSIM, ALPHA IF (N.GT.NMAX .OR. NF.GT.NFMAX .OR. NSIM.GT.NSMAX) THEN WRITE (NOUT,99991) GO TO 140 END IF READ (NIN,*) (Y(I),I=1,N) * Read in the ITYPE dependent arguments (skipping headings) IF (ITYPE.EQ.1) THEN READ (NIN,*) PARAM(1) P = 0 IVAL = 1 ELSE IF (ITYPE.EQ.2) THEN READ (NIN,*) PARAM(1), PARAM(2) P = 0 IVAL = 2 ELSE IF (ITYPE.EQ.3) THEN READ (NIN,*) PARAM(1), PARAM(2), PARAM(3) P = 0 IVAL = 2 ELSE READ (NIN,*) PARAM(1), PARAM(2), PARAM(3), PARAM(4), P IVAL = P + 2 IF (P.GT.PMAX) THEN WRITE (NOUT,99991) GO TO 140 END IF END IF * Read in the MODE dependent arguments (skipping headings) IF (MODE.EQ.0) THEN * User supplied initial values READ (NIN,*) (INIT(I),I=1,IVAL) ELSE IF (MODE.EQ.1) THEN * Continuing from a previously saved R READ (NIN,*) (R(I),I=1,P+13) ELSE IF (MODE.EQ.2) THEN * Initial values calculated from first K observations READ (NIN,*) K END IF * * Fit a smoothing model (parameter R in G05PMF and STATE in G13AMF * are in the same format) IFAIL = 1 CALL G13AMF(MODE,ITYPE,P,PARAM,N,Y,K,INIT,NF,FV,FSE,YHAT,RES,DV, + AD,R,IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99994) IFAIL GO TO 140 END IF * * Simulate forecast values from the model, and don't update R SMODE = 2 VAR = DV*DV EN = N * Simulate NSIM forecasts DO 40 I = 1, NSIM * Simulations assuming gaussian errors IFAIL = 1 CALL G05PMF(SMODE,NF,ITYPE,P,PARAM,INIT,VAR,R,STATE,E,0,TSIM1, + IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99993) IFAIL GO TO 140 END IF * Bootstrapping errors IFAIL = 1 CALL G05PMF(SMODE,NF,ITYPE,P,PARAM,INIT,0.0D0,R,STATE,RES,EN, + TSIM2,IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99993) IFAIL GO TO 140 END IF * Copy and transpose the simulated values DO 20 J = 1, NF GSIM(I,J) = TSIM1(J) BSIM(I,J) = TSIM2(J) 20 CONTINUE 40 CONTINUE * * Calculate CI based on the quantiles for each simulated forecast Q(1) = ALPHA/TWO Q(2) = ONE - Q(1) DO 60 I = 1, NF IFAIL = 1 CALL G01AMF(NSIM,GSIM(1,I),2,Q,GLIM(1,I),IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99992) IFAIL GO TO 140 END IF IFAIL = 1 CALL G01AMF(NSIM,BSIM(1,I),2,Q,BLIM(1,I),IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99992) IFAIL GO TO 140 END IF 60 CONTINUE * * Display the forecast values and associated prediction intervals WRITE (NOUT,*) WRITE (NOUT,*) 'Initial values used:' DO 80 I = 1, IVAL WRITE (NOUT,99998) I, INIT(I) 80 CONTINUE 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,*) DO 100 I = 1, N WRITE (NOUT,99998) I, Y(I), YHAT(I), RES(I) 100 CONTINUE 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 120 I = 1, NF TMP = Z*FSE(I) WRITE (NOUT,99997) N + I, FV(I), FV(I) - TMP, FV(I) + TMP, + GLIM(1,I), GLIM(2,I), BLIM(1,I), BLIM(2,I) 120 CONTINUE WRITE (NOUT,99996) 100.0D0*(ONE-ALPHA), '% CIs were produced' * 140 CONTINUE * 99999 FORMAT (A,E12.4) 99998 FORMAT (I4,3(1X,F12.3)) 99997 FORMAT (I3,7(1X,F10.3)) 99996 FORMAT (1X,F5.1,A) 99995 FORMAT (1X,' ** G05KFF returned with IFAIL = ',I5) 99994 FORMAT (1X,' ** G13AMF returned with IFAIL = ',I5) 99993 FORMAT (1X,' ** G05PMF returned with IFAIL = ',I5) 99992 FORMAT (1X,' ** G01AMF returned with IFAIL = ',I5) 99991 FORMAT (1X,' ** Problem size too large, increase array limits') END