    PROGRAM g05hkfe

!      G05HKF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g05hkf, g05kbf, g13faf, g13fbf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER   :: zero = 0.0E0_nag_wp
       INTEGER, PARAMETER              :: nout = 6, nparmx = 10, nregmx = 10,  &
                                          num = 1500, num1 = 3000
       INTEGER, PARAMETER              :: ldcovr = nparmx
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: df, fac1, gamma, hp, lgf, mean, tol, &
                                          xterm
       INTEGER                         :: i, ifail, igen, ip, iq, isym, k,     &
                                          ldx, lwork, maxit, mn, npar, npar2,  &
                                          nreg, nt
       LOGICAL                         :: fcall
       CHARACTER (1)                   :: dist
!      .. Local Arrays ..
       REAL (KIND=nag_wp)              :: bx(10), covr(ldcovr,nparmx),         &
                                          cvar(100), etm(num1), ht(num1+10),   &
                                          htm(num1), param(nparmx), rvec(40),  &
                                          rwsav(9), sc(nparmx), se(nparmx),    &
                                          theta(nparmx),                       &
                                      work(num1*3+nparmx+nregmx*num1+20*20+1), &
                                          x(num1,10), yt(num1+10)
       INTEGER                         :: iseed(4)
       LOGICAL                         :: copts(2)
!      .. Intrinsic Functions ..
       INTRINSIC                          real, sin
!      .. Executable Statements ..
       WRITE (nout,*) 'G05HKF Example Program Results'
       WRITE (nout,*)

!      Initialize the generator
       iseed(1) = 111
       igen = 0
       CALL g05kbf(igen,iseed)

!      Example using Normal errors
       nreg = 0
       ldx = num1
       bx(1) = 1.5E0_nag_wp
       bx(2) = 2.5E0_nag_wp
       bx(3) = 3.0E0_nag_wp
       mean = 3.0E0_nag_wp
       DO i = 1, num
          fac1 = real(i,kind=nag_wp)*0.01E0_nag_wp
          x(i,1) = 0.01E0_nag_wp + 0.7E0_nag_wp*sin(fac1)
          x(i,2) = 0.5E0_nag_wp + fac1*0.1E0_nag_wp
          x(i,3) = 1.0E0_nag_wp
       END DO

       isym = 1
       mn = 1
       gamma = -0.4E0_nag_wp
       ip = 0
       iq = 3
       param(1) = 0.8E0_nag_wp
       param(2) = 0.6E0_nag_wp
       param(3) = 0.2E0_nag_wp
       param(4) = 0.1E0_nag_wp
       npar = 1 + iq + ip
       lwork = nreg*num + 3*num + npar + isym + mn + nreg + 403
       dist = 'N'

!      First call to routine
       fcall = .TRUE.
       ifail = 0
       CALL g05hkf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,ifail)
!      Discard the first NUM values and call the routine again
       fcall = .FALSE.
       ifail = 0
       CALL g05hkf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,ifail)

       DO i = 1, num
          xterm = zero
          DO k = 1, nreg
             xterm = xterm + x(i,k)*bx(k)
          END DO
          IF (mn==1) THEN
             yt(i) = mean + xterm + yt(i)
          ELSE
             yt(i) = xterm + yt(i)
          END IF
       END DO

       DO i = 1, npar
          theta(i) = param(i)*0.5E0_nag_wp
       END DO

       IF (isym==1) THEN
          theta(npar+isym) = gamma*0.5E0_nag_wp
       END IF

!      Fit a GARCH model to the generated data
       copts(1) = .TRUE.
       copts(2) = .TRUE.
       maxit = 200
       tol = 1.0E-5_nag_wp
       npar2 = 1 + iq + ip + isym + mn + nreg
       ifail = 0
       CALL g13faf(dist,yt,x,ldx,num,ip,iq,nreg,mn,isym,npar2,theta,se,sc, &
          covr,ldcovr,hp,etm,htm,lgf,copts,maxit,tol,work,lwork,ifail)
!      Display the output
       WRITE (nout,*)
       WRITE (nout,*) 'Normal distribution'
       WRITE (nout,*)
       WRITE (nout,*) '             Correct    Parameter   Standard'
       WRITE (nout,*) '             values     estimates    errors'
       DO i = 1, npar
          WRITE (nout,99998) 'PARAM(', i, ')', param(i), theta(i), se(i)
       END DO
       IF (isym==1) THEN
          WRITE (nout,99999) 'Gamma   ', gamma, theta(npar+1), se(npar+1)
       END IF
       IF (mn==1) THEN
          WRITE (nout,99999) 'Mean    ', mean, theta(npar+isym+1), &
             se(npar+isym+1)
       END IF
       DO i = 1, nreg
          WRITE (nout,99998) 'BX(', i, ')   ', bx(i), theta(npar+isym+mn+i), &
             se(npar+isym+mn+i)
       END DO

       nt = 4
       ifail = 0
       CALL g13fbf(num,nt,ip,iq,theta,gamma,cvar,htm,etm,ifail)
       WRITE (nout,*)
       WRITE (nout,99997) 'Volatility forecast = ', cvar(nt)
       WRITE (nout,*)

!      Example using Students T errors
       dist = 'T'
       nreg = 2
       mn = 1
       df = 4.1E0_nag_wp
       ip = 1
       iq = 2
       isym = 1
       gamma = -0.2E0_nag_wp

       npar = iq + ip + 1
       lwork = nreg*num + 3*num + npar + isym + mn + nreg + 404

       param(1) = 0.1E0_nag_wp
       param(2) = 0.2E0_nag_wp
       param(3) = 0.3E0_nag_wp
       param(4) = 0.4E0_nag_wp

       igen = 0
       iseed(1) = 111
       CALL g05kbf(igen,iseed)

!      First call to the routine
       ifail = 0
       fcall = .TRUE.
       CALL g05hkf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,ifail)
!      Second call to the routine
       fcall = .FALSE.
       CALL g05hkf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,ifail)
!      Discard the previously generated values, and call the routine again
       fcall = .FALSE.
       CALL g05hkf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,ifail)

       DO i = 1, num
          xterm = zero
          DO k = 1, nreg
             xterm = xterm + x(i,k)*bx(k)
          END DO
          IF (mn==1) THEN
             yt(i) = mean + xterm + yt(i)
          ELSE
             yt(i) = xterm + yt(i)
          END IF
       END DO

       copts(1) = .TRUE.
       copts(2) = .TRUE.

       maxit = 200
       tol = 1.0E-5_nag_wp

       DO i = 1, npar
          theta(i) = param(i)*0.5E0_nag_wp
       END DO

       theta(npar+isym) = gamma*0.5E0_nag_wp
       theta(npar+isym+1) = df*0.65E0_nag_wp

       npar2 = 2 + ip + iq + isym + mn + nreg
       ifail = 0
       CALL g13faf(dist,yt,x,ldx,num,ip,iq,nreg,mn,isym,npar2,theta,se,sc, &
          covr,ldcovr,hp,etm,htm,lgf,copts,maxit,tol,work,lwork,ifail)

       WRITE (nout,*)
       WRITE (nout,*) 'Students t distribution'
       WRITE (nout,*)
       WRITE (nout,*) '             Correct    Parameter   Standard'
       WRITE (nout,*) '             values     estimates    errors'
       DO i = 1, npar
          WRITE (nout,99998) 'PARAM(', i, ')', param(i), theta(i), se(i)
       END DO
       IF (isym==1) THEN
          WRITE (nout,99999) 'Gamma   ', gamma, theta(npar+isym), &
             se(npar+isym)
       END IF
       WRITE (nout,99999) 'DF      ', df, theta(npar+isym+1), se(npar+isym+1)
       IF (mn==1) THEN
          WRITE (nout,99999) 'Mean    ', mean, theta(npar+isym+2), &
             se(npar+isym+2)
       END IF
       DO i = 1, nreg
          WRITE (nout,99998) 'BX(', i, ')   ', bx(i), theta(npar+isym+1+mn+i), &
             se(npar+isym+1+mn+i)
       END DO

       nt = 4
       ifail = 0
       CALL g13fbf(num,nt,ip,iq,theta,gamma,cvar,htm,etm,ifail)

       WRITE (nout,*)
       WRITE (nout,99997) 'Volatility forecast = ', cvar(nt)

99999  FORMAT (1X,A,3(F10.2,2X))
99998  FORMAT (1X,A,I1,A,3(F10.2,2X))
99997  FORMAT (1X,A,F12.2)
    END PROGRAM g05hkfe
