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

Syntax

C#
public static void g02qg(
	int sorder,
	int ic1,
	int n,
	int m,
	double[,] dat,
	int[] isx,
	int ip,
	double[] y,
	double[] wt,
	int ntau,
	double[] tau,
	out double df,
	double[,] b,
	double[,] bl,
	double[,] bu,
	double[,,] ch,
	double[,] res,
	G02..::..g02qgOptions options,
	G05..::..G05State g05state,
	int[] info,
	out int ifail
)
Visual Basic
Public Shared Sub g02qg ( _
	sorder As Integer, _
	ic1 As Integer, _
	n As Integer, _
	m As Integer, _
	dat As Double(,), _
	isx As Integer(), _
	ip As Integer, _
	y As Double(), _
	wt As Double(), _
	ntau As Integer, _
	tau As Double(), _
	<OutAttribute> ByRef df As Double, _
	b As Double(,), _
	bl As Double(,), _
	bu As Double(,), _
	ch As Double(,,), _
	res As Double(,), _
	options As G02..::..g02qgOptions, _
	g05state As G05..::..G05State, _
	info As Integer(), _
	<OutAttribute> ByRef ifail As Integer _
)
Visual C++
public:
static void g02qg(
	int sorder, 
	int ic1, 
	int n, 
	int m, 
	array<double,2>^ dat, 
	array<int>^ isx, 
	int ip, 
	array<double>^ y, 
	array<double>^ wt, 
	int ntau, 
	array<double>^ tau, 
	[OutAttribute] double% df, 
	array<double,2>^ b, 
	array<double,2>^ bl, 
	array<double,2>^ bu, 
	array<double,3>^ ch, 
	array<double,2>^ res, 
	G02..::..g02qgOptions^ options, 
	G05..::..G05State^ g05state, 
	array<int>^ info, 
	[OutAttribute] int% ifail
)
F#
static member g02qg : 
        sorder : int * 
        ic1 : int * 
        n : int * 
        m : int * 
        dat : float[,] * 
        isx : int[] * 
        ip : int * 
        y : float[] * 
        wt : float[] * 
        ntau : int * 
        tau : float[] * 
        df : float byref * 
        b : float[,] * 
        bl : float[,] * 
        bu : float[,] * 
        ch : float[,,] * 
        res : float[,] * 
        options : G02..::..g02qgOptions * 
        g05state : G05..::..G05State * 
        info : int[] * 
        ifail : int byref -> unit 

Parameters

sorder
Type: System..::..Int32
On entry: determines the storage order of variates supplied in dat.
Constraint: sorder=1 or 2.
ic1
Type: System..::..Int32
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.
ic1=1
An intercept will be included in the model.
ic1=0
An intercept will not be included in the model.
Constraint: ic1=0 or 1.
n
Type: System..::..Int32
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.
m
Type: System..::..Int32
On entry: m, the total number of variates in the dataset.
Constraint: m0.
dat
Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [dim1, _sddat]
Note: dim1 must satisfy the constraint:
  • if sorder=1, dim1n;
  • otherwise dim1m.
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-1,j-1] if sorder=1, and
  • dat[j-1,i-1] if sorder=2.
The design matrix X is constructed from dat, isx and ic1.
isx
Type: array<System..::..Int32>[]()[][]
An array of size [m]
On entry: indicates which independent variables are to be included in the model.
isx[j-1]=0
The jth variate, supplied in dat, is not included in the regression model.
isx[j-1]=1
The jth variate, supplied in dat, is included in the regression model.
Constraints:
  • isx[j-1]=0 or 1, for j=1,2,,m;
  • if ic1=1, exactly ip-1 values of isx must be set to 1;
  • if ic1=0, exactly ip values of isx must be set to 1.
ip
Type: System..::..Int32
On entry: p, the number of independent variables in the model, including the intercept, see ic1, if present.
Constraints:
  • 1ip<n;
  • if ic1=1, 1ipm+1;
  • if ic1=0, 1ipm.
y
Type: array<System..::..Double>[]()[][]
An array of size [n]
On entry: y, observations on the dependent variable.
wt
Type: array<System..::..Double>[]()[][]
An array of size [_lwt]
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-1]=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-1]0.0, for i=1,2,,n;
  • The effective number of observations 2.
ntau
Type: System..::..Int32
On entry: the number of quantiles of interest.
Constraint: ntau1.
tau
Type: array<System..::..Double>[]()[][]
An array of size [ntau]
On entry: the vector of quantiles of interest. A separate model is fitted to each quantile.
Constraint: ε<tau[j-1]<1-ε where ε is the machine precision returned by x02aj, for j=1,2,,ntau.
df
Type: System..::..Double%
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.
b
Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [ip, ntau]
On entry: if Calculate Initial Values=NO, b[i-1,l-1] 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-1,l-1], for i=1,2,,ip, contains the estimates of the parameters of the regression model, β^, estimated for τ=tau[l-1].
If ic1=1, b[0,l-1] will contain the estimate corresponding to the intercept and b[i,l-1] will contain the coefficient of the jth variate contained in dat, where isx[j-1] is the ith nonzero value in the array isx.
If ic1=0, b[i-1,l-1] will contain the coefficient of the jth variate contained in dat, where isx[j-1] is the ith nonzero value in the array isx.
bl
Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [dim1, ntau]
Note: dim1 must satisfy the constraint:
Note: the second dimension of the array bl must be at least ntau if Interval MethodNONE.
On exit: if Interval MethodNONE, bl[i-1,l-1] contains the lower limit of an 100×α% confidence interval for b[i-1,l-1], 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.
bu
Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [dim1, ntau]
Note: dim1 must satisfy the constraint:
Note: the second dimension of the array bu must be at least ntau if Interval MethodNONE.
On exit: if Interval MethodNONE, bu[i-1,l-1] contains the upper limit of an 100×α% confidence interval for b[i-1,l-1], 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.
ch
Type: array<System..::..Double,3>[,](,)[,][,]
An array of size [dim1, dim2, dim3]
Note: dim1 must satisfy the constraint:
Note: dim2 must satisfy the constraint:
Note: dim3 must satisfy the constraint:
Note: the last dimension of the array ch must be at least ntau if Interval MethodNONE and Matrix Returned=COVARIANCE and at least ntau+1 if Interval MethodNONE, IID or BOOTSTRAP XY and Matrix 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-1,j-1,l-1] holds an estimate of the covariance between b[i-1,l-1] and b[j-1,l-1].
  • If Matrix Returned=H INVERSE, ch[i-1,j-1,0] holds an estimate of the i,jth element of nJn and ch[i-1,j-1,l+1-1] holds an estimate of the i,jth element of n-1Hn-1, for τ=tau[l-1].
The method used for calculating Σ and Hn-1 is controlled by the optional parameter Interval Method.
res
Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [n, dim2]
Note: dim2 must satisfy the constraint:
On exit: if Return Residuals=YES, res[i-1,l-1] holds the (weighted) residuals, ri, for τ=tau[l-1], for i=1,2,,n and l=1,2,,ntau.
If wtis notNULL 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.
options
Type: NagLibrary..::..G02..::..g02qgOptions
An Object of type G02.g02qgOptions. Used to configure optional parameters to this method.
g05state
Type: NagLibrary..::..G05..::..G05State
An Object of type G05.G05State.
info
Type: array<System..::..Int32>[]()[][]
An array of size [ntau]
On exit: info[i] holds additional information concerning the model fitting and confidence limit calculations when τ=tau[i].
CodeWarning
0Model fitted and confidence limits (if requested) calculated successfully
1The method 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.
2A singular matrix was encountered during the optimization. The model was not fitted for this value of τ.
4Some truncation occurred whilst calculating the confidence limits for this value of τ. See [Algorithmic Details] for details. The returned upper and lower limits may be narrower than specified.
8The method did not converge whilst calculating the confidence limits. The returned limits are based on the estimate at the last iteration. Try increasing Iteration Limit.
16Confidence 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.
ifail
Type: System..::..Int32%
On exit: ifail=0 unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

Description

Given a vector of n observed values, y=yi: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, g02qg estimates the p-element vector β as the solution to
minimizeβp i=1 n ρτyi-xiTβ (1)
where ρτ is the piecewise linear loss function ρτz=zτ-Iz<0, and Iz<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.
g02qg uses the interior point algorithm of Portnoy and Koenker (1997), described briefly in [Algorithmic Details], 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-τnHn-1JnHn-1
where Jn=n-1 i=1 n xixiT 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-tn-p,1+α/2Σ^ii,β^Ui=β^i+tn-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:
(i) Independent, identically distributed (IID) errors
Under an assumption of IID errors the asymptotic relationship for Σ simplifies to
Σ=τ1-τnsτ2XTX-1
where s is the sparsity function. g02qg estimates sτ from the residuals, ri=yi-xiTβ^ and a bandwidth hn.
(ii) Powell Sandwich
Powell (1991) suggested estimating the matrix Hn by a kernel estimator of the form
H^n=ncn-1 i=1 n KricnxixiT
where K is a kernel function and cn satisfies limncn0 and limnncn. When the Powell method is chosen, g02qg 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.
(iii) Hendricks–Koenker Sandwich
Koenker (2005) suggested estimating the matrix Hn using
H^n=n-1 i=1 n 2hnxiTβ^τ+hn-β^τ-hnxixiT
where hn is a bandwidth and β^τ+hn denotes the parameter estimates obtained from a quantile regression using the τ+hnth quantile. Similarly with β^τ-hn.
(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 [Algorithmic Details].
All three asymptotic estimates of the covariance matrix require a bandwidth, hn. Two alternative methods for determining this are provided:
(i) Sheather–Hall
hn=1.5Φ-1αbϕΦ-1τ2n2Φ-1τ+113
for a user-supplied value αb,
(ii) Bofinger
hn=4.5ϕΦ-1τ4n2Φ-1τ+1215
g02qg allows optional arguments to be supplied via the iopts and opts arrays (see [Optional Parameters] for details of the available options). Prior to calling g02qg the optional parameter arrays, iopts and opts must be initialized by calling (G02ZKF not in this release) with optstr set to Initialize=g02qg (see [Optional Parameters] for details on the available options). If bootstrap confidence limits are required (Interval Method=BOOTSTRAP XY) then one of the random number initialization methods (G05KFF not in this release) (for a repeatable analysis) or (G05KGF not in this release) (for an unrepeatable analysis) must also have been previously called.

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 (1999) Numerical Optimization 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

Error Indicators and Warnings

Errors or warnings detected by the method:
Some error messages may refer to parameters that are dropped from this interface (LDDAT, RIP, TDCH, SDRES, LIOPTS, LOPTS, LSTATE) In these cases, an error in another parameter has usually caused an incorrect value to be inferred.
ifail=11
On entry, sorder1 or 2.
ifail=21
On entry, ic11 or 0.
ifail=31
On entry, _weight"U" or "W".
ifail=41
On entry, n<2.
ifail=51
On entry, m<0.
ifail=71
On entry, sorder=1, lddat<n.
ifail=72
On entry, sorder=2, lddat<m.
ifail=81
On entry, isx[j-1]0 or 1.
ifail=91
On entry, ip<1 or ipn.
ifail=92
On entry, ip is not consistent with isx and ic1.
ifail=111
On entry, _weight="W" and wt[i-1]<0.0 for at least one i.
ifail=112
On entry, the effective number of observations is less than two.
ifail=121
On entry, ntau<1.
ifail=131
On entry, tau is invalid.
ifail=201
On entry, one or more of the optional parameter arrays iopts and opts have not been initialized or have been corrupted.
ifail=221
On entry, Interval Method=BOOTSTRAP XY and state was not initialized or has been corrupted.
ifail=231
On exit, problems were encountered whilst fitting at least one model. Additional information has been returned in info.
ifail=-4000
Invalid dimension for array value
ifail=-8000
Negative dimension for array value
ifail=-6000
Invalid Parameters value

Accuracy

Not applicable.

Parallelism and Performance

None.

Further Comments

g02qg allocates internally approximately the following elements of real storage: 13n+np+3p2+6p+3p+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, wtisNULL, ic1=0 and ip=m an internal copy of the input data is avoided and the amount of locally allocated memory is reduced by np.

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.

Example program (C#): g02qge.cs

Example program data: g02qge.d

Example program results: g02qge.r

Algorithmic Details

By the addition of slack variables the minimization (1) can be reformulated into the linear programming problem
minimizeu,v,β+n×+n×pτeTu+1-τeTv​   subject to  y=Xβ+u-v (2)
and its associated dual
maximizedyTd​   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
maximizeayTa​   subject to  XTa=1-τXTe,a0,1n. (4)
The algorithm introduced by Portnoy and Koenker (1997) and used by g02qg, 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 (1999).

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).

Central path

Rather than optimize (4) directly, an additional slack variable s is added and the constraint a0,1n is replaced with a+s=e,ai0,si0, for i=1,2,,n.
The positivity constraint on a and s is handled using the logarithmic barrier function
Ba,s,μ=yTa+μ i=1 n logai+logsi.
The primal-dual form of the problem is used giving the Lagrangian
La,s,β,u,μ=Ba,s,μ-βTXTa-1-τXTe-uTa+s-e
whose central path is described by the following first order conditions
XTa=1-τXTea+s=eXβ+u-v=ySUe=μeAVe=μ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
XTa+δa=1-τXTea+δa+s+δs=eXβ+δβ+u+δu-v+δv=yS+ΔsU+Δue=μeA+ΔaV+Δve=μe (6)
where Δa,Δs,Δu and Δv denote the diagonal matrices with diagonal elements given by δa,δs,δu and δv respectively.

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δβ=XTWy-Xβ+τ-1XTe+XTaδa=Wy-Xβ-Xδβδs=-δaδu=S-1Uδa-Ueδv=A-1Vδs-Ve (7)
where W=S-1U+A-1V-1.
Initial step sizes for the primal (γ^P) and dual (γ^D) parameters are constructed as
γ^P=σ×minmini,δai<0ai/δai,mini,δsi<0si/δsiγ^D=σ×minmini,δui<0ui/δui,mini,δvi<0vi/δvi (8)
where σ is a user-supplied scaling factor. If γ^P×γ^D1 then the nonlinearity adjustment, described in [Nonlinearity Adjustment], is not made and the model parameters are updated using the current step size and directions.

Nonlinearity Adjustment

In the nonlinearity adjustment step a new estimate of μ is obtained by letting
g^γ^P,γ^D=s+γ^PδsTu+γ^Dδu+a+γ^PδaTv+γ^Dδv
and estimating μ as
μ=g^γ^P,γ^Dg^0,03g^0,02n.
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δβ=XTWy-Xβ+μS-1-A-1e+S-1ΔsΔue-A-1ΔaΔve+τ-1XTe+XTaδa=Wy-Xβ-Xδβ+μS-1-A-1δs=-δaδu=μS-1e+S-1Uδa-Ue-S-1ΔsΔueδv=μA-1e+A-1Vδs-Ve-A-1ΔaΔve (9)
are solved and updated values for δβ,δa,δs,δu,δv,γ^P and γ^D calculated.

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
τeTu+1-τeTv-dTy=τeTu+1-τeTv-a-1-τeTy=sTu+aTv=eTu-aTy+1-τeTXβ
and the optimization terminates if the duality gap is smaller than the tolerance supplied in the optional parameter Tolerance.

Additional information

Initial values are required for the parameters a,s,u,v and β. If not supplied by the user, 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 f08bf. 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 parameter 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 not in this release).

Calculation of Covariance Matrix

g02qg 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 [Description].
(i) Independent, identically distributed (IID) errors
When assuming IID errors, the covariance matrices depend on the sparsity, sτ, which g02qg estimates as follows:
(a) Let ri denote the residuals from the original quantile regression, that is ri=yi-xiTβ^.
(b) Drop any residual where ri is less than εu, supplied in the optional parameter Epsilon.
(c) Sort and relabel the remaining residuals in ascending order, by absolute value, so that εu<r1<r2<.
(d) Select the first l values where l=hnn, for some bandwidth hn.
(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.
(f) Use the resulting estimate of the slope as an estimate of the sparsity.
(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.
(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 2hn/di in the calculation of Hn, max2hn/di+εu,0 is used instead, where εu is supplied in the optional parameter Epsilon.

Description of Monitoring Information

See the description of the optional argument Monitoring.

See Also