Program g02qgfe
! G02QGF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
! .. Use Statements ..
Use nag_library, Only: g02qgf, g02zkf, g02zlf, g05kff, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: lseed = 1, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp) :: df, rvalue
Integer :: genid, i, ifail, ip, ivalue, j, l, &
ldbl, lddat, ldres, liopts, lopts, &
lstate, lwt, m, n, ntau, optype, &
sorder, subid, tdch
Character (1) :: c1, weight
Character (30) :: cvalue, semeth
Character (100) :: optstr
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: b(:,:), bl(:,:), bu(:,:), ch(:,:,:), &
dat(:,:), opts(:), res(:,:), tau(:), &
wt(:), y(:)
Integer, Allocatable :: info(:), iopts(:), isx(:), state(:)
Integer :: seed(lseed)
! .. Intrinsic Procedures ..
Intrinsic :: count, len_trim, min
! .. Executable Statements ..
Write (nout,*) 'G02QGF Example Program Results'
Write (nout,*)
Flush (nout)
! Skip heading in data file
Read (nin,*)
! Read in the problem size
Read (nin,*) sorder, c1, weight, n, m, ntau
! Read in the data
If (weight=='W' .Or. weight=='w') Then
lwt = n
Else
lwt = 0
End If
Allocate (wt(lwt),isx(m),y(n),tau(ntau))
If (sorder==1) Then
! DAT(N,M)
lddat = n
Allocate (dat(lddat,m))
If (lwt==0) Then
Read (nin,*)(dat(i,1:m),y(i),i=1,n)
Else
Read (nin,*)(dat(i,1:m),y(i),wt(i),i=1,n)
End If
Else
! DAT(M,N)
lddat = m
Allocate (dat(lddat,n))
If (lwt==0) Then
Read (nin,*)(dat(1:m,i),y(i),i=1,n)
Else
Read (nin,*)(dat(1:m,i),y(i),wt(i),i=1,n)
End If
End If
! Read in variable inclusion flags
Read (nin,*) isx(1:m)
! Calculate IP
ip = count(isx(1:m)==1)
If (c1=='Y' .Or. c1=='y') Then
ip = ip + 1
End If
! Read in the quantiles required
Read (nin,*) tau(1:ntau)
liopts = 100
lopts = 100
Allocate (iopts(liopts),opts(lopts))
! Initialize the optional argument array
ifail = 0
Call g02zkf('INITIALIZE = G02QGF',iopts,liopts,opts,lopts,ifail)
c_lp: Do
! Read in any optional arguments. Reads in to the end of
! the input data, or until a blank line is reached
ifail = 1
Read (nin,99994,Iostat=ifail) optstr
If (ifail/=0) Then
Exit c_lp
Else If (len_trim(optstr)==0) Then
Exit c_lp
End If
! Set the supplied option
ifail = 0
Call g02zkf(optstr,iopts,liopts,opts,lopts,ifail)
End Do c_lp
! Assume that no intervals or output matrices are required
! unless the optional argument state differently
ldbl = 0
tdch = 0
ldres = 0
lstate = 0
! Query the optional arguments to see what output is required
ifail = 0
Call g02zlf('INTERVAL METHOD',ivalue,rvalue,cvalue,optype,iopts,opts, &
ifail)
semeth = cvalue
If (semeth/='NONE') Then
! Require the intervals to be output
ldbl = ip
If (semeth=='BOOTSTRAP XY') Then
! Need to find the length of the state array for the random
! number generator
! Read in the generator ID and a seed
Read (nin,*) genid, subid, seed(1)
! Query the length of the state array
Allocate (state(lstate))
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
! Deallocate STATE so that it can reallocated later
Deallocate (state)
End If
ifail = 0
Call g02zlf('MATRIX RETURNED',ivalue,rvalue,cvalue,optype,iopts,opts, &
ifail)
If (cvalue=='COVARIANCE') Then
tdch = ntau
Else If (cvalue=='H INVERSE') Then
If (semeth=='BOOTSTRAP XY' .Or. semeth=='IID') Then
! NB: If we are using bootstrap or IID errors then any request for
! H INVERSE is ignored
tdch = 0
Else
tdch = ntau + 1
End If
End If
ifail = 0
Call g02zlf('RETURN RESIDUALS',ivalue,rvalue,cvalue,optype,iopts,opts, &
ifail)
If (cvalue=='YES') Then
ldres = n
End If
End If
! Allocate memory for output arrays
Allocate (b(ip,ntau),info(ntau),bl(ldbl,ntau),bu(ldbl,ntau), &
ch(ldbl,ldbl,tdch),state(lstate),res(ldres,ntau))
If (lstate>0) Then
! Doing bootstrap, so initialise the RNG
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
End If
! Call the model fitting routine
ifail = -1
Call g02qgf(sorder,c1,weight,n,m,dat,lddat,isx,ip,y,wt,ntau,tau,df,b,bl, &
bu,ch,res,iopts,opts,state,info,ifail)
If (ifail/=0) Then
If (ifail==231) Then
Write (nout,*) 'Additional error information (INFO): ', info(1:ntau)
Else
Go To 100
End If
End If
! Display the parameter estimates
Do l = 1, ntau
Write (nout,99999) 'Quantile: ', tau(l)
Write (nout,*)
If (ldbl>0) Then
Write (nout,*) ' Lower Parameter Upper'
Write (nout,*) ' Limit Estimate Limit'
Else
Write (nout,*) ' Parameter'
Write (nout,*) ' Estimate'
End If
Do j = 1, ip
If (ldbl>0) Then
Write (nout,99998) j, bl(j,l), b(j,l), bu(j,l)
Else
Write (nout,99998) j, b(j,l)
End If
End Do
Write (nout,*)
If (tdch==ntau) Then
Write (nout,*) 'Covariance matrix'
Do i = 1, ip
Write (nout,99997) ch(1:i,i,l)
End Do
Write (nout,*)
Else If (tdch==ntau+1) Then
Write (nout,*) 'J'
Do i = 1, ip
Write (nout,99997) ch(1:i,i,1)
End Do
Write (nout,*)
Write (nout,*) 'H inverse'
Do i = 1, ip
Write (nout,99997) ch(1:i,i,l+1)
End Do
Write (nout,*)
End If
Write (nout,*)
End Do
If (ldres>0) Then
Write (nout,*) 'First 10 Residuals'
Write (nout,*) ' Quantile'
Write (nout,99995) 'Obs.', tau(1:ntau)
Do i = 1, min(n,10)
Write (nout,99996) i, res(i,1:ntau)
End Do
Else
Write (nout,*) 'Residuals not returned'
End If
Write (nout,*)
100 Continue
99999 Format (1X,A,F6.3)
99998 Format (1X,I3,3(3X,F7.3))
99997 Format (1X,10(E10.3,3X))
99996 Format (2X,I3,10(1X,F10.5))
99995 Format (1X,A,10(3X,F6.3,2X))
99994 Format (A100)
End Program g02qgfe