    PROGRAM g05hmfe

!      G05HMF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g05hmf, g05kbf, g13fef, g13fff, 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 = 2000, num1 = 3000
       INTEGER, PARAMETER              :: ldcovr = nparmx
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: df, fac1, gamma, hp, lgf, mean, tol, &
                                          xterm
       INTEGER                         :: i, iflag, igen, ip, iq, 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,*) 'G05HMF Example Program Results'
       WRITE (nout,*)

       iseed(1) = 111
       igen = 0

       lwork = num1*3 + nparmx + nregmx*num1 + 1
       nreg = 0
       ldx = num1
       df = 5.1E0_nag_wp
       gamma = 0.1E0_nag_wp
       bx(1) = 1.5E0_nag_wp
       bx(2) = 2.5E0_nag_wp
       bx(3) = 3.0E0_nag_wp
       mean = 4.0E0_nag_wp

       DO i = 1, num
          fac1 = real(i,kind=nag_wp)*0.01E0_nag_wp
          x(i,2) = 0.01E0_nag_wp + 0.7E0_nag_wp*sin(fac1)
          x(i,1) = 0.5E0_nag_wp + fac1*0.1E0_nag_wp
          x(i,3) = 1.0E0_nag_wp
       END DO

       mn = 1
       nreg = 2
       gamma = 0.1E0_nag_wp
       ip = 1
       iq = 1
       npar = iq + ip + 1
       param(1) = 0.4E0_nag_wp
       param(2) = 0.1E0_nag_wp
       param(3) = 0.7E0_nag_wp
       fcall = .TRUE.

       dist = 'N'
       CALL g05kbf(igen,iseed)
       iflag = 1
       CALL g05hmf(dist,200,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,iflag)

       IF (iflag==0) THEN
          fcall = .FALSE.
          iflag = 1
          CALL g05hmf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen, &
             iseed,rwsav,iflag)
       END IF
       IF (iflag/=0) THEN
          WRITE (nout,99997) iflag
          GO TO 20
       END IF

       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

       iflag = -1

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

       maxit = 100
       tol = 1.0E-5_nag_wp

       DO i = 1, npar
          theta(i) = param(i)*0.5E0_nag_wp
       END DO
       theta(npar+1) = gamma*0.5E0_nag_wp
       IF (mn==1) THEN
          theta(npar+mn+1) = mean*0.5E0_nag_wp
       END IF
       DO i = 1, nreg
          theta(npar+mn+1+i) = bx(i)*0.5E0_nag_wp
       END DO

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

       IF (iflag<0) GO TO 20
       WRITE (nout,*)
       WRITE (nout,*) 'Normal distribution'
       WRITE (nout,*)
       WRITE (nout,*) '          Parameter       Standard        Correct'
       WRITE (nout,*) '          estimates       errors          values'

       DO i = 1, npar
          WRITE (nout,99999) theta(i), se(i), param(i)
       END DO

       WRITE (nout,99999) theta(npar+1), se(npar+1), gamma

       IF (mn==1) THEN
          WRITE (nout,99999) theta(npar+2), se(npar+2), mean
       END IF

       DO i = 1, nreg
          WRITE (nout,99999) theta(npar+mn+i+1), se(npar+mn+i+1), bx(i)
       END DO

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

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

       dist = 'T'
       mean = 3.0E0_nag_wp

       DO i = 1, num
          fac1 = real(i,kind=nag_wp)*0.01E0_nag_wp
          x(i,2) = 0.01E0_nag_wp + 0.7E0_nag_wp*sin(fac1)
          x(i,1) = 0.5E0_nag_wp + fac1*0.1E0_nag_wp
          x(i,3) = 1.0E0_nag_wp
       END DO

       mn = 1
       nreg = 2
       gamma = 0.09E0_nag_wp
       ip = 1
       iq = 1

       npar = iq + ip + 1
       param(1) = 0.05E0_nag_wp
       param(2) = 0.1E0_nag_wp
       param(3) = 0.8E0_nag_wp

       fcall = .TRUE.
       iseed(1) = 111
       CALL g05kbf(igen,iseed)

       iflag = 0
       CALL g05hmf(dist,200,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,iflag)
       fcall = .FALSE.
       iflag = 0
       CALL g05hmf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,iflag)
       iflag = 0
       CALL g05hmf(dist,num,ip,iq,param,gamma,df,ht,yt,fcall,rvec,igen,iseed, &
          rwsav,iflag)

       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

       iflag = -1

       maxit = 100
       tol = 1.0E-5_nag_wp

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

       theta(npar+1) = gamma*0.5E0_nag_wp

       theta(npar+2) = df*0.65E0_nag_wp

       IF (mn==1) THEN
          theta(npar+2+mn) = mean*0.5E0_nag_wp
       END IF
       DO i = 1, nreg
          theta(npar+2+mn+i) = bx(i)*0.5E0_nag_wp
       END DO

       copts(1) = .TRUE.
       copts(2) = .TRUE.
       hp = 0.5E0_nag_wp

       npar2 = 3 + ip + iq + mn + nreg

       CALL g13fef(dist,yt,x,ldx,num,ip,iq,nreg,mn,npar2,theta,se,sc,covr, &
          ldcovr,hp,etm,htm,lgf,copts,maxit,tol,work,lwork,iflag)

       WRITE (nout,*)
       WRITE (nout,*) 'Students t distribution'
       WRITE (nout,*)
       WRITE (nout,*) '          Parameter       Standard        Correct'
       WRITE (nout,*) '          estimates       errors          values'

       DO i = 1, npar
          WRITE (nout,99999) theta(i), se(i), param(i)
       END DO

       WRITE (nout,99999) theta(npar+1), se(npar+1), gamma

       WRITE (nout,99999) theta(npar+2), se(npar+2), df

       IF (mn==1) THEN
          WRITE (nout,99999) theta(npar+2+mn), se(npar+2+mn), mean
       END IF

       DO i = 1, nreg
          WRITE (nout,99999) theta(npar+2+mn+i), se(npar+2+mn+i), bx(i)
       END DO

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

       WRITE (nout,*)
       WRITE (nout,99998) 'Volatility forecast = ', cvar(nt)
       WRITE (nout,*)
20     CONTINUE

99999  FORMAT (1X,3F16.3)
99998  FORMAT (1X,A,F12.3)
99997  FORMAT (1X,' ** G05HMF returned with IFAIL = ',I5)
    END PROGRAM g05hmfe
