NAG FL Interface
g02qgf (quantile_​linreg)

Note: this routine uses optional parameters to define choices in the problem specification and in the details of the algorithm. If you wish to use default settings for all of the optional parameters, you need only read Sections 1 to 10 of this document. If, however, you wish to reset some or all of the settings please refer to Section 11 for a detailed description of the algorithm, to Section 12 for a detailed description of the specification of the optional parameters and to Section 13 for a detailed description of the monitoring information produced by the routine.
Settings help

FL Name Style:


FL Specification Language:


1 Purpose

g02qgf performs a multiple linear quantile regression. Parameter estimates and, if required, confidence limits, covariance matrices and residuals are calculated. g02qgf may be used to perform a weighted quantile regression. A simplified interface for g02qgf is provided by g02qff.

2 Specification

Fortran Interface
Subroutine g02qgf ( sorder, intcpt, weight, n, m, dat, lddat, isx, ip, y, wt, ntau, tau, df, b, bl, bu, ch, res, iopts, opts, state, info, ifail)
Integer, Intent (In) :: sorder, n, m, lddat, isx(m), ip, ntau, iopts(*)
Integer, Intent (Inout) :: state(*), ifail
Integer, Intent (Out) :: info(ntau)
Real (Kind=nag_wp), Intent (In) :: dat(lddat,*), y(n), wt(*), tau(ntau), opts(*)
Real (Kind=nag_wp), Intent (Inout) :: b(ip,ntau), bl(ip,*), bu(ip,*), ch(ip,ip,*), res(n,*)
Real (Kind=nag_wp), Intent (Out) :: df
Character (1), Intent (In) :: intcpt, weight
C Header Interface
#include <nag.h>
void  g02qgf_ (const Integer *sorder, const char *intcpt, const char *weight, const Integer *n, const Integer *m, const double dat[], const Integer *lddat, const Integer isx[], const Integer *ip, const double y[], const double wt[], const Integer *ntau, const double tau[], double *df, double b[], double bl[], double bu[], double ch[], double res[], const Integer iopts[], const double opts[], Integer state[], Integer info[], Integer *ifail, const Charlen length_intcpt, const Charlen length_weight)
The routine may be called by the names g02qgf or nagf_correg_quantile_linreg.

3 Description

Given a vector of n observed values, y = { y i :i=1,2,,n} , an n×p design matrix X, a column vector, x, of length p holding the ith row of X and a quantile τ (0,1) , g02qgf estimates the p-element vector β as the solution to
minimize β p i=1 n ρ τ ( y i -xiTβ) (1)
where ρ τ is the piecewise linear loss function ρ τ (z) = z (τ-I(z<0)) , and I (z<0) is an indicator function taking the value 1 if z<0 and 0 otherwise. Weights can be incorporated by replacing X and y with WX and Wy respectively, where W is an n×n diagonal matrix. Observations with zero weights can either be included or excluded from the analysis; this is in contrast to least squares regression where such observations do not contribute to the objective function and are, therefore, always dropped.
g02qgf uses the interior point algorithm of Portnoy and Koenker (1997), described briefly in Section 11, to obtain the parameter estimates β^, for a given value of τ.
Under the assumption of Normally distributed errors, Koenker (2005) shows that the limiting covariance matrix of β^-β has the form
Σ = τ (1-τ) n H n -1 J n H n -1  
where Jn = n-1 i=1 n x i xiT and Hn is a function of τ, as described below. Given an estimate of the covariance matrix, Σ^, lower (β^L) and upper (β^U) limits for an (100×α)% confidence interval can be calculated for each of the p parameters, via
β^ Li = β^ i - t n-p , (1+α) / 2 Σ^ ii , β^ Ui = β^ i + t n-p , (1+α) / 2 Σ^ ii  
where tn-p,0.975 is the 97.5 percentile of the Student's t distribution with n-k degrees of freedom, where k is the rank of the cross-product matrix XTX.
Four methods for estimating the covariance matrix, Σ, are available:
  1. (i)Independent, identically distributed (IID) errors
    Under an assumption of IID errors the asymptotic relationship for Σ simplifies to
    Σ = τ (1-τ) n (s(τ)) 2 (XTX)-1  
    where s is the sparsity function. g02qgf estimates s(τ) from the residuals, ri = yi - xiT β^ and a bandwidth hn.
  2. (ii)Powell Sandwich
    Powell (1991) suggested estimating the matrix Hn by a kernel estimator of the form
    H^ n = (ncn)-1 i=1 n K( r i c n ) xi xiT  
    where K is a kernel function and cn satisfies lim n cn 0 and lim n n cn . When the Powell method is chosen, g02qgf uses a Gaussian kernel (i.e., K=ϕ) and sets
    cn = min(σr, (qr3-qr1) / 1.34 ) × (Φ-1(τ+hn)-Φ-1(τ-hn))  
    where hn is a bandwidth, σr,qr1 and qr3 are, respectively, the standard deviation and the 25% and 75% quantiles for the residuals, ri.
  3. (iii)Hendricks–Koenker Sandwich
    Koenker (2005) suggested estimating the matrix Hn using
    H^ n = n-1 i=1 n [ 2 hn xiT (β^(τ+hn)-β^(τ-hn)) ] xi xiT  
    where hn is a bandwidth and β^(τ+hn) denotes the parameter estimates obtained from a quantile regression using the (τ+hn)th quantile. Similarly with β^(τ-hn).
  4. (iv)Bootstrap
    The last method uses bootstrapping to either estimate a covariance matrix or obtain confidence intervals for the parameter estimates directly. This method, therefore, does not assume Normally distributed errors. Samples of size n are taken from the paired data {yi,xi} (i.e., the independent and dependent variables are sampled together). A quantile regression is then fitted to each sample resulting in a series of bootstrap estimates for the model parameters, β. A covariance matrix can then be calculated directly from this series of values. Alternatively, confidence limits, β^L and β^U, can be obtained directly from the (1-α)/2 and (1+α)/2 sample quantiles of the bootstrap estimates.
Further details of the algorithms used to calculate the covariance matrices can be found in Section 11.
All three asymptotic estimates of the covariance matrix require a bandwidth, hn. Two alternative methods for determining this are provided:
  1. (i)Sheather–Hall
    hn = ( 1.5 (Φ-1(αb)ϕ(Φ-1(τ))) 2 n (2Φ-1(τ)+1) ) 13  
    for a user-supplied value αb,
  2. (ii)Bofinger
    hn = ( 4.5 (ϕ(Φ-1(τ))) 4 n (2Φ-1(τ)+1) 2 ) 15  
g02qgf allows optional parameters to be supplied via the iopts and opts arrays (see Section 12 for details of the available options). Prior to calling g02qgf the optional parameter arrays, iopts and opts must be initialized by calling g02zkf with optstr set to Initialize=g02qgf (see Section 12 for details on the available options). If bootstrap confidence limits are required (Interval Method=BOOTSTRAP XY) then one of the random number initialization routines g05kff (for a repeatable analysis) or g05kgf (for an unrepeatable analysis) must also have been previously called.

4 References

Koenker R (2005) Quantile Regression Econometric Society Monographs, Cambridge University Press, New York
Mehrotra S (1992) On the implementation of a primal-dual interior point method SIAM J. Optim. 2 575–601
Nocedal J and Wright S J (2006) Numerical Optimization (2nd Edition) Springer Series in Operations Research, Springer, New York
Portnoy S and Koenker R (1997) The Gaussian hare and the Laplacian tortoise: computability of squared-error versus absolute error estimators Statistical Science 4 279–300
Powell J L (1991) Estimation of monotonic regression models under quantile restrictions Nonparametric and Semiparametric Methods in Econometrics Cambridge University Press, Cambridge

5 Arguments

1: sorder Integer Input
On entry: determines the storage order of variates supplied in dat.
Constraint: sorder=1 or 2.
2: intcpt Character(1) Input
On entry: indicates whether an intercept will be included in the model. The intercept is included by adding a column of ones as the first column in the design matrix, X.
intcpt='Y'
An intercept will be included in the model.
intcpt='N'
An intercept will not be included in the model.
Constraint: intcpt='N' or 'Y'.
3: weight Character(1) Input
On entry: indicates if weights are to be used.
weight='W'
A weighted regression model is fitted to the data using weights supplied in array wt.
weight='U'
An unweighted regression model is fitted to the data and array wt is not referenced.
Constraint: weight='U' or 'W'.
4: n Integer Input
On entry: the total number of observations in the dataset. If no weights are supplied, or no zero weights are supplied or observations with zero weights are included in the model then n=n. Otherwise n=n+ the number of observations with zero weights.
Constraint: n2.
5: m Integer Input
On entry: m, the total number of variates in the dataset.
Constraint: m0.
6: dat(lddat,*) Real (Kind=nag_wp) array Input
Note: the second dimension of the array dat must be at least m if sorder=1, and at least n otherwise.
On entry: the ith value for the jth variate, for i=1,2,,n and j=1,2,,m, must be supplied in
  • dat(i,j) if sorder=1, and
  • dat(j,i) if sorder=2.
The design matrix X is constructed from dat, isx and intcpt.
7: lddat Integer Input
On entry: the first dimension of the array dat as declared in the (sub)program from which g02qgf is called.
Constraints:
  • if sorder=1, lddatn;
  • otherwise lddatm.
8: isx(m) Integer array Input
On entry: indicates which independent variables are to be included in the model.
isx(j)=0
The jth variate, supplied in dat, is not included in the regression model.
isx(j)=1
The jth variate, supplied in dat, is included in the regression model.
Constraints:
  • isx(j)=0 or 1, for j=1,2,,m;
  • if intcpt='Y', exactly ip-1 values of isx must be set to 1;
  • if intcpt='N', exactly ip values of isx must be set to 1.
9: ip Integer Input
On entry: p, the number of independent variables in the model, including the intercept, see intcpt, if present.
Constraints:
  • 1ip<n;
  • if intcpt='Y', 1ipm+1;
  • if intcpt='N', 1ipm.
10: y(n) Real (Kind=nag_wp) array Input
On entry: y, the observations on the dependent variable.
11: wt(*) Real (Kind=nag_wp) array Input
Note: the dimension of the array wt must be at least n if weight='W'.
On entry: if weight='W', wt must contain the diagonal elements of the weight matrix W. Otherwise wt is not referenced.
When
Drop Zero Weights=YES
If wt(i)=0.0, the ith observation is not included in the model, in which case the effective number of observations, n, is the number of observations with nonzero weights. If Return Residuals=YES, the values of res will be set to zero for observations with zero weights.
Drop Zero Weights=NO
All observations are included in the model and the effective number of observations is n, i.e., n=n.
Constraints:
  • If weight='W', wt(i)0.0, for i=1,2,,n;
  • The effective number of observations 2.
12: ntau Integer Input
On entry: the number of quantiles of interest.
Constraint: ntau1.
13: tau(ntau) Real (Kind=nag_wp) array Input
On entry: the vector of quantiles of interest. A separate model is fitted to each quantile.
Constraint: ε<tau(j)<1-ε where ε is the machine precision returned by x02ajf, for j=1,2,,ntau.
14: df Real (Kind=nag_wp) Output
On exit: the degrees of freedom given by n-k, where n is the effective number of observations and k is the rank of the cross-product matrix XTX.
15: b(ip,ntau) Real (Kind=nag_wp) array Input/Output
On entry: if Calculate Initial Values=NO, b(i,l) must hold an initial estimates for β^i, for i=1,2,,ip and l=1,2,,ntau. If Calculate Initial Values=YES, b need not be set.
On exit: b(i,l), for i=1,2,,ip, contains the estimates of the parameters of the regression model, β^, estimated for τ=tau(l).
If intcpt='Y', b(1,l) will contain the estimate corresponding to the intercept and b(i+1,l) will contain the coefficient of the jth variate contained in dat, where isx(j) is the ith nonzero value in the array isx.
If intcpt='N', b(i,l) will contain the coefficient of the jth variate contained in dat, where isx(j) is the ith nonzero value in the array isx.
16: bl(ip,*) Real (Kind=nag_wp) array Output
Note: the second dimension of the array bl must be at least ntau if Interval MethodNONE.
On exit: if Interval MethodNONE, bl(i,l) contains the lower limit of an (100×α)% confidence interval for b(i,l), for i=1,2,,ip and l=1,2,,ntau.
If Interval Method=NONE, bl is not referenced.
The method used for calculating the interval is controlled by the optional parameters Interval Method and Bootstrap Interval Method. The size of the interval, α, is controlled by the optional parameter Significance Level.
17: bu(ip,*) Real (Kind=nag_wp) array Output
Note: the second dimension of the array bu must be at least ntau if Interval MethodNONE.
On exit: if Interval MethodNONE, bu(i,l) contains the upper limit of an (100×α)% confidence interval for b(i,l), for i=1,2,,ip and l=1,2,,ntau.
If Interval Method=NONE, bu is not referenced.
The method used for calculating the interval is controlled by the optional parameters Interval Method and Bootstrap Interval Method. The size of the interval, α is controlled by the optional parameter Significance Level.
18: ch(ip,ip,*) Real (Kind=nag_wp) array Output
Note: the last dimension of the array ch must be at least ntau if Interval MethodNONEandMatrix Returned=COVARIANCE and at least ntau+1 if Interval MethodNONE, IID or BOOTSTRAP XYandMatrix Returned=H INVERSE.
On exit: depending on the supplied optional parameters, ch will either not be referenced, hold an estimate of the upper triangular part of the covariance matrix, Σ, or an estimate of the upper triangular parts of nJn and n-1Hn-1.
If Interval Method=NONE or Matrix Returned=NONE, ch is not referenced.
If Interval Method=BOOTSTRAP XY or IID and Matrix Returned=H INVERSE, ch is not referenced.
Otherwise, for i,j=1,2,,ip,ji and l=1,2,,ntau:
  • If Matrix Returned=COVARIANCE, ch(i,j,l) holds an estimate of the covariance between b(i,l) and b(j,l).
  • If Matrix Returned=H INVERSE, ch(i,j,1) holds an estimate of the (i,j)th element of nJn and ch(i,j,l+1) holds an estimate of the (i,j)th element of n-1Hn-1, for τ=tau(l).
The method used for calculating Σ and Hn-1 is controlled by the optional parameter Interval Method.
19: res(n,*) Real (Kind=nag_wp) array Output
Note: the second dimension of the array res must be at least ntau if Return Residuals=YES.
On exit: if Return Residuals=YES, res(i,l) holds the (weighted) residuals, ri, for τ=tau(l), for i=1,2,,n and l=1,2,,ntau.
If weight='W' and Drop Zero Weights=YES, the value of res will be set to zero for observations with zero weights.
If Return Residuals=NO, res is not referenced.
20: iopts(*) Integer array Communication Array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument iopts in the previous call to g02zkf.
On entry: optional parameter array, as initialized by a call to g02zkf.
21: opts(*) Real (Kind=nag_wp) array Communication Array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument opts in the previous call to g02zkf.
On entry: optional parameter array, as initialized by a call to g02zkf.
22: state(*) Integer array Communication Array
Note: the actual argument supplied must be the array state supplied to the initialization routines g05kff or g05kgf.
If Interval Method=BOOTSTRAP XY, state contains information about the selected random number generator. Otherwise state is not referenced.
23: info(ntau) Integer array Output
On exit: info(i) holds additional information concerning the model fitting and confidence limit calculations when τ=tau(i).
Code Warning
0 Model fitted and confidence limits (if requested) calculated successfully
1 The routine did not converge. The returned values are based on the estimate at the last iteration. Try increasing Iteration Limit whilst calculating the parameter estimates or relaxing the definition of convergence by increasing Tolerance.
2 A singular matrix was encountered during the optimization. The model was not fitted for this value of τ.
4 Some truncation occurred whilst calculating the confidence limits for this value of τ. See Section 11 for details. The returned upper and lower limits may be narrower than specified.
8 The routine did not converge whilst calculating the confidence limits. The returned limits are based on the estimate at the last iteration. Try increasing Iteration Limit.
16 Confidence limits for this value of τ could not be calculated. The returned upper and lower limits are set to a large positive and large negative value respectively as defined by the optional parameter Big.
It is possible for multiple warnings to be applicable to a single model. In these cases the value returned in info is the sum of the corresponding individual nonzero warning codes.
24: ifail Integer Input/Output
On entry: ifail must be set to 0, −1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of −1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value −1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or −1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=11
On entry, sorder=value.
Constraint: sorder=1 or 2.
ifail=21
On entry, intcpt=value was an illegal value.
ifail=31
On entry, weight had an illegal value.
Constraint: weight='U' or 'W'.
ifail=41
On entry, n=value.
Constraint: n2.
ifail=51
On entry, m=value.
Constraint: m0.
ifail=71
On entry, lddat=value and n=value.
Constraint: lddatn.
ifail=72
On entry, lddat=value and m=value.
Constraint: lddatm.
ifail=81
On entry, isx(value)=value.
Constraint: isx(i)=0 or 1, for all i.
ifail=91
On entry, ip=value and n=value.
Constraint: 1ip<n.
ifail=92
On entry, ip is not consistent with isx or intcpt: ip=value, expected value=value.
ifail=111
On entry, wt(value)=value.
Constraint: wt(i)0.0, for all i.
ifail=112
On entry, effective number of observations=value.
Constraint: effective number of observationsvalue.
ifail=121
On entry, ntau=value.
Constraint: ntau1.
ifail=131
On entry, tau(value)=value.
Constraint: ε<tau(l)<1-ε where ε is the machine precision returned by x02ajf, for all ntau.
ifail=201
On entry, either the option arrays have not been initialized or they have been corrupted.
ifail=221
On entry, state vector has been corrupted or not initialized.
ifail=231
A potential problem occurred whilst fitting the model(s).
Additional information has been returned in info.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

Not applicable.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
g02qgf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g02qgf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

g02qgf allocates internally approximately the following elements of real storage: 13n+ np+ 3p2+ 6p+ 3(p+1)×ntau. If Interval Method=BOOTSTRAP XY then a further np elements are required, and this increases by p×ntau×Bootstrap Iterations if Bootstrap Interval Method=QUANTILE. Where possible, any user-supplied output arrays are used as workspace and so the amount actually allocated may be less. If sorder=2, weight='U', intcpt='N' and ip=m an internal copy of the input data is avoided and the amount of locally allocated memory is reduced by np.

10 Example

A quantile regression model is fitted to Engels 1857 study of household expenditure on food. The model regresses the dependent variable, household food expenditure, against two explanatory variables, a column of ones and household income. The model is fit for five different values of τ and the covariance matrix is estimated assuming Normal IID errors. Both the covariance matrix and the residuals are returned.

10.1 Program Text

Program Text (g02qgfe.f90)

10.2 Program Data

Program Data (g02qgfe.d)

10.3 Program Results

Program Results (g02qgfe.r)
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 0 500 1000 1500 2000 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Household Food Expenditure Household Income Example Program Quantile Regression Engels 1857 Study of Household Expenditure on Food gnuplot_plot_1 gnuplot_plot_2 τ = 0.10 gnuplot_plot_3 τ = 0.25 gnuplot_plot_4 τ = 0.50 gnuplot_plot_5 τ = 0.75 gnuplot_plot_6 τ = 0.90

11 Algorithmic Details

By the addition of slack variables the minimization (1) can be reformulated into the linear programming problem
minimize (u,v,β) + n × + n × p τ e T u + (1-τ) e T v ​   subject to   y = Xβ + u - v (2)
and its associated dual
maximize d y T d ​   subject to   XTd=0, d [τ-1,τ] n (3)
where e is a vector of n 1s. Setting a=d+(1-τ)e gives the equivalent formulation
maximize a y T a ​   subject to   XT a = (1-τ) X T e , a [0,1] n . (4)
The algorithm introduced by Portnoy and Koenker (1997) and used by g02qgf, uses the primal-dual formulation expressed in equations (2) and (4) along with a logarithmic barrier function to obtain estimates for β. The algorithm is based on the predictor-corrector algorithm of Mehrotra (1992) and further details can be obtained from Portnoy and Koenker (1997) and Koenker (2005). A good description of linear programming, interior point algorithms, barrier functions and Mehrotra's predictor-corrector algorithm can be found in Nocedal and Wright (2006).

11.1 Interior Point Algorithm

In this section a brief description of the interior point algorithm used to estimate the model parameters is presented. It should be noted that there are some differences in the equations given here – particularly (7) and (9) – compared to those given in Koenker (2005) and Portnoy and Koenker (1997).

11.1.1 Central path

Rather than optimize (4) directly, an additional slack variable s is added and the constraint a [0,1] n is replaced with a+s=e , a i 0 , s i 0 , for i= 1,2,,n .
The positivity constraint on a and s is handled using the logarithmic barrier function
B (a,s,μ) = yT a + μ i=1 n (logai+logsi) .  
The primal-dual form of the problem is used giving the Lagrangian
L (a,s,β,u,μ) = B (a,s,μ) - βT (XTa-(1-τ)XTe) - uT (a+s-e)  
whose central path is described by the following first order conditions
XTa = (1-τ) XTe a+s = e Xβ+u-v = y SUe = μe AVe = μe (5)
where A denotes the diagonal matrix with diagonal elements given by a, similarly with S,U and V. By enforcing the inequalities on s and a strictly, i.e., ai>0 and si>0 for all i we ensure that A and S are positive definite diagonal matrices and hence A-1 and S-1 exist.
Rather than applying Newton's method to the system of equations given in (5) to obtain the step directions δβ,δa,δs,δu and δv , Mehrotra substituted the steps directly into (5) giving the augmented system of equations
XT (a+δa) = (1-τ)XTe (a+δa) + (s+δs) = e X (β+δβ) + (u+δu) - (v+δv) = y (S+Δs) (U+Δu) e = μe (A+Δa) (V+Δv) e = μe (6)
where Δa,Δs,Δu and Δv denote the diagonal matrices with diagonal elements given by δa,δs,δu and δv respectively.

11.1.2 Affine scaling step

The affine scaling step is constructed by setting μ=0 in (5) and applying Newton's method to obtain an intermediate set of step directions
(XTWX) δ β = XT W (y-Xβ) + (τ-1) XT e + XT a δ a = W (y-Xβ-X δ β ) δ s = - δ a δ u = S-1 U δ a - U e δ v = A-1 V δ s - V e (7)
where W = (S-1U+A-1V)-1 .
Initial step sizes for the primal (γ^P) and dual (γ^D) parameters are constructed as
γ^P = σ × min{ min i , δ a i < 0 {ai/ δ a i }, min i , δ s i < 0 {si/ δ s i }} γ^D = σ × min{ min i , δ u i < 0 {ui/ δ u i }, min i , δ v i < 0 {vi/ δ v i }} (8)
where σ is a user-supplied scaling factor. If γ^P × γ^D 1 then the nonlinearity adjustment, described in Section 11.1.3, is not made and the model parameters are updated using the current step size and directions.

11.1.3 Nonlinearity Adjustment

In the nonlinearity adjustment step a new estimate of μ is obtained by letting
g^(γ^P,γ^D) = (s+γ^Pδs) T (u+γ^Dδu) + (a+γ^Pδa) T (v+γ^Dδv)  
and estimating μ as
μ = ( g^(γ^P,γ^D) g^(0,0) ) 3 g^(0,0) 2n .  
This estimate, along with the nonlinear terms (Δu, Δs, Δa and Δv) from (6) are calculated using the values of δa,δs,δu and δv obtained from the affine scaling step.
Given an updated estimate for μ and the nonlinear terms the system of equations
(XTWX) δ β = XT W (y-Xβ+μ(S-1-A-1)e+S-1 Δ s Δ u e-A-1 Δ a Δ v e) + (τ-1) XT e + XT a δ a = W (y-Xβ-X δ β +μ(S-1-A-1)) δ s = - δ a δ u = μ S-1 e + S-1 U δ a - U e - S-1 Δs Δu e δ v = μ A-1 e + A-1 V δ s - V e - A-1 Δa Δv e (9)
are solved and updated values for δβ,δa,δs,δu,δv,γ^P and γ^D calculated.

11.1.4 Update and convergence

At each iteration the model parameters (β,a,s,u,v) are updated using step directions, (δβ,δa,δs,δu,δv) and step lengths (γ^P,γ^D) .
Convergence is assessed using the duality gap, that is, the differences between the objective function in the primal and dual formulations. For any feasible point (u,v,s,a) the duality gap can be calculated from equations (2) and (3) as
τ eT u + (1-τ) eT v - dT y = τ eT u + (1-τ) eT v - (a-(1-τ)e) T y = sT u + aT v = eT u - aT y + (1-τ) eT X β  
and the optimization terminates if the duality gap is smaller than the tolerance supplied in the optional parameter Tolerance.

11.1.5 Additional information

Initial values are required for the parameters a,s,u,v and β. If you have not supplied them, initial values for β are calculated from a least squares regression of y on X. This regression is carried out by first constructing the cross-product matrix XTX and then using a pivoted QR decomposition as performed by f08bff. In addition, if the cross-product matrix is not of full rank, a rank reduction is carried out and, rather than using the full design matrix, X, a matrix formed from the first p-rank columns of XP is used instead, where P is the pivot matrix used during the QR decomposition. Parameter estimates, confidence intervals and the rows and columns of the matrices returned in the argument ch (if any) are set to zero for variables dropped during the rank-reduction. The rank reduction step is performed irrespective of whether initial values are supplied by the user.
Once initial values have been obtained for β, the initial values for u and v are calculated from the residuals. If |ri|<εu then a value of ±εu is used instead, where εu is supplied in the optional parameter Epsilon. The initial values for the a and s are always set to 1-τ and τ respectively.
The solution for δβ in both (7) and (9) is obtained using a Bunch–Kaufman decomposition, as implemented in f07mdf.

11.2 Calculation of Covariance Matrix

g02qgf supplies four methods to calculate the covariance matrices associated with the parameter estimates for β. This section gives some additional detail on three of the algorithms, the fourth, (which uses bootstrapping), is described in Section 3.
  1. (i)Independent, identically distributed (IID) errors
    When assuming IID errors, the covariance matrices depend on the sparsity, s(τ), which g02qgf estimates as follows:
    1. (a)Let ri denote the residuals from the original quantile regression, that is ri = yi - xiT β^ .
    2. (b)Drop any residual where |ri| is less than εu, supplied in the optional parameter Epsilon.
    3. (c)Sort and relabel the remaining residuals in ascending order, by absolute value, so that εu < |r1| < |r2| < .
    4. (d)Select the first l values where l=hnn, for some bandwidth hn.
    5. (e)Sort and relabel these l residuals again, so that r1 < r2 < < rl and regress them against a design matrix with two columns (p=2) and rows given by xi = {1, i/(n-p) } using quantile regression with τ=0.5.
    6. (f)Use the resulting estimate of the slope as an estimate of the sparsity.
  2. (ii)Powell Sandwich
    When using the Powell Sandwich to estimate the matrix Hn, the quantity
    cn = min(σr, (qr3-qr1) / 1.34 ) × (Φ-1(τ+hn)-Φ-1(τ-hn))  
    is calculated. Dependent on the value of τ and the method used to calculate the bandwidth (hn), it is possible for the quantities τ±hn to be too large or small, compared to machine precision (ε). More specifically, when τ-hnε, or τ+hn1-ε, a warning flag is raised in info, the value is truncated to ε or 1-ε respectively and the covariance matrix calculated as usual.
  3. (iii)Hendricks–Koenker Sandwich
    The Hendricks–Koenker Sandwich requires the calculation of the quantity di= xiT (β^(τ+hn)-β^(τ-hn)) . As with the Powell Sandwich, in cases where τ-hnε, or τ+hn1-ε, a warning flag is raised in info, the value truncated to ε or 1-ε respectively and the covariance matrix calculated as usual.
    In addition, it is required that di>0, in this method. Hence, instead of using 2 hn / di in the calculation of Hn, max( 2 hn / (di+εu) ,0) is used instead, where εu is supplied in the optional parameter Epsilon.

12 Optional Parameters

Several optional parameters in g02qgf control aspects of the optimization algorithm, methodology used, logic or output. Their values are contained in the arrays iopts and opts; these must be initialized before calling g02qgf by first calling g02zkf with optstr set to Initialize=g02qgf.
Each optional parameter has an associated default value; to set any of them to a non-default value, use g02zkf. The current value of an optional parameter can be queried using g02zlf.
The remainder of this section can be skipped if you wish to use the default values for all optional parameters.
The following is a list of the optional parameters available. A full description of each optional parameter is provided in Section 12.1.

12.1 Description of the Optional Parameters

For each option, we give a summary line, a description of the optional parameter and details of constraints.
The summary line contains:
Keywords and character values are case and white space insensitive.
Band Width AlpharDefault =1.0
A multiplier used to construct the parameter αb used when calculating the Sheather–Hall bandwidth (see Section 3), with αb=(1-α)×Band Width Alpha. Here, α is the Significance Level.
Constraint: Band Width Alpha>0.0.
Band Width MethodaDefault =SHEATHER HALL
The method used to calculate the bandwidth used in the calculation of the asymptotic covariance matrix Σ and H-1 if Interval Method=HKS, KERNEL or IID (see Section 3).
Constraint: Band Width Method=SHEATHER HALL or BOFINGER.
BigrDefault =10.020
This parameter should be set to something larger than the biggest value supplied in dat and y.
Constraint: Big>0.0.
Bootstrap Interval MethodaDefault =QUANTILE
If Interval Method=BOOTSTRAP XY, Bootstrap Interval Method controls how the confidence intervals are calculated from the bootstrap estimates.
Bootstrap Interval Method=T
t intervals are calculated. That is, the covariance matrix, Σ = { σ ij :i,j=1,2,,p} is calculated from the bootstrap estimates and the limits calculated as β i ± t ( n - p , (1+α) / 2 ) σ ii where t ( n - p , (1+α) / 2 ) is the (1+α) / 2 percentage point from a Student's t distribution on n - p degrees of freedom, n is the effective number of observations and α is given by the optional parameter Significance Level.
Bootstrap Interval Method=QUANTILE
Quantile intervals are calculated. That is, the upper and lower limits are taken as the (1+α) / 2 and (1-α) / 2 quantiles of the bootstrap estimates, as calculated using g01amf.
Constraint: Bootstrap Interval Method=T or QUANTILE.
Bootstrap IterationsiDefault =100
The number of bootstrap samples used to calculate the confidence limits and covariance matrix (if requested) when Interval Method=BOOTSTRAP XY.
Constraint: Bootstrap Iterations>1.
Bootstrap MonitoringaDefault =NO
If Bootstrap Monitoring=YES and Interval Method=BOOTSTRAP XY, the parameter estimates for each of the bootstrap samples are displayed. This information is sent to the unit number specified by Unit Number.
Constraint: Bootstrap Monitoring=YES or NO.
Calculate Initial ValuesaDefault =YES
If Calculate Initial Values=YES then the initial values for the regression parameters, β, are calculated from the data. Otherwise they must be supplied in b.
Constraint: Calculate Initial Values=YES or NO.
Defaults
This special keyword is used to reset all optional parameters to their default values.
Drop Zero WeightsaDefault =YES
If a weighted regression is being performed and Drop Zero Weights=YES then observations with zero weight are dropped from the analysis. Otherwise such observations are included.
Constraint: Drop Zero Weights=YES or NO.
EpsilonrDefault =ε
εu, the tolerance used when calculating the covariance matrix and the initial values for u and v. For additional details see Section 11.2 and Section 11.1.5 respectively.
Constraint: Epsilon0.0.
Interval MethodaDefault =IID
The value of Interval Method controls whether confidence limits are returned in bl and bu and how these limits are calculated. This parameter also controls how the matrices returned in ch are calculated.
Interval Method=NONE
No limits are calculated and bl, bu and ch are not referenced.
Interval Method=KERNEL
The Powell Sandwich method with a Gaussian kernel is used.
Interval Method=HKS
The Hendricks–Koenker Sandwich is used.
Interval Method=IID
The errors are assumed to be identical, and independently distributed.
Interval Method=BOOTSTRAP XY
A bootstrap method is used, where sampling is done on the pair (yi,xi). The number of bootstrap samples is controlled by the parameter Bootstrap Iterations and the type of interval constructed from the bootstrap samples is controlled by Bootstrap Interval Method.
Constraint: Interval Method=NONE, KERNEL, HKS, IID or BOOTSTRAP XY.
Iteration LimitiDefault =100
The maximum number of iterations to be performed by the interior point optimization algorithm.
Constraint: Iteration Limit>0.
Matrix ReturnedaDefault =NONE
The value of Matrix Returned controls the type of matrices returned in ch. If Interval Method=NONE, this parameter is ignored and ch is not referenced. Otherwise:
Matrix Returned=NONE
No matrices are returned and ch is not referenced.
Matrix Returned=COVARIANCE
The covariance matrices are returned.
Matrix Returned=H INVERSE
If Interval Method=KERNEL or HKS, the matrices J and H-1 are returned. Otherwise no matrices are returned and ch is not referenced.
The matrices returned are calculated as described in Section 3, with the algorithm used specified by Interval Method. In the case of Interval Method=BOOTSTRAP XY the covariance matrix is calculated directly from the bootstrap estimates.
Constraint: Matrix Returned=NONE, COVARIANCE or H INVERSE.
MonitoringaDefault =NO
If Monitoring=YES then the duality gap is displayed at each iteration of the interior point optimization algorithm. In addition, the final estimates for β are also displayed.
The monitoring information is sent to the unit number specified by Unit Number.
Constraint: Monitoring=YES or NO.
QR TolerancerDefault =ε0.9
The tolerance used to calculate the rank, k, of the p×p cross-product matrix, XTX. Letting Q be the orthogonal matrix obtained from a QR decomposition of XTX, then the rank is calculated by comparing Qii with Q11×QR Tolerance.
If the cross-product matrix is rank deficient, the parameter estimates for the p-k columns with the smallest values of Qii are set to zero, along with the corresponding entries in bl, bu and ch, if returned. This is equivalent to dropping these variables from the model. Details on the QR decomposition used can be found in f08bff.
Constraint: QR Tolerance>0.0.
Return ResidualsaDefault =NO
If Return Residuals=YES, the residuals are returned in res. Otherwise res is not referenced.
Constraint: Return Residuals=YES or NO.
SigmarDefault =0.99995
The scaling factor used when calculating the affine scaling step size (see equation (8)).
Constraint: 0.0<Sigma<1.0.
Significance LevelrDefault =0.95
α, the size of the confidence interval whose limits are returned in bl and bu.
Constraint: 0.0<Significance Level<1.0.
TolerancerDefault =ε
Convergence tolerance. The optimization is deemed to have converged if the duality gap is less than Tolerance (see Section 11.1.4).
Constraint: Tolerance>0.0.
Unit Numberi Default =advisory message unit number
The unit number to which any monitoring information is sent.
Constraint: Unit Number>1.

13 Description of Monitoring Information

See the description of the optional parameter Monitoring.