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 )
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 _ )
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 )
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
- Type: System..::..Int32On entry: determines the storage order of variates supplied in dat.Constraint: or .
- Type: System..::..Int32On 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, .
Constraint: or .
- An intercept will be included in the model.
- An intercept will not be included in the model.
- Type: System..::..Int32On 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 . Otherwise the number of observations with zero weights.Constraint: .
- Type: System..::..Int32On entry: , the total number of variates in the dataset.Constraint: .
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [dim1, _sddat]Note: dim1 must satisfy the constraint:
- if , ;
- otherwise .
- Type: array<System..::..Int32>()An array of size [m]On entry: indicates which independent variables are to be included in the model.
- Type: System..::..Int32On entry: , the number of independent variables in the model, including the intercept, see ic1, if present.Constraints:
- if , ;
- if , .
- Type: array<System..::..Double>()An array of size [n]On entry: , observations on the dependent variable.
- Type: array<System..::..Double>()An array of size [_lwt]On entry: if , wt must contain the diagonal elements of the weight matrix . Otherwise wt is not referenced.When
- If , the th observation is not included in the model, in which case the effective number of observations, , is the number of observations with nonzero weights. If , the values of res will be set to zero for observations with zero weights.
- All observations are included in the model and the effective number of observations is n, i.e., .
- If , , for ;
- The effective number of observations .
- Type: System..::..Int32On entry: the number of quantiles of interest.Constraint: .
- 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: where is the machine precision returned by x02aj, for .
- Type: System..::..Double%On exit: the degrees of freedom given by , where is the effective number of observations and is the rank of the cross-product matrix .
- Type: array<System..::..Double,2>[,](,)[,][,]On entry: if , must hold an initial estimates for , for and . If , b need not be set.On exit: , for , contains the estimates of the parameters of the regression model, , estimated for .If , will contain the estimate corresponding to the intercept and will contain the coefficient of the th variate contained in dat, where is the th nonzero value in the array isx.
- 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 if .On exit: if , contains the lower limit of an confidence interval for , for and .If , bl is not referenced.
- 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 if .On exit: if , contains the upper limit of an confidence interval for , for and .If , bu is not referenced.
- 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 if and and at least if , or and .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 and .If or , ch is not referenced.If or and , ch is not referenced.Otherwise, for and :
The method used for calculating and is controlled by the optional parameter Interval Method.
- If , holds an estimate of the covariance between and .
- If , holds an estimate of the th element of and holds an estimate of the th element of , for .
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [n, dim2]Note: dim2 must satisfy the constraint:
- Type: NagLibrary..::..G02..::..g02qgOptionsAn Object of type G02.g02qgOptions. Used to configure optional parameters to this method.
- Type: array<System..::..Int32>()An array of size On exit: holds additional information concerning the model fitting and confidence limit calculations when .
Code Warning Model fitted and confidence limits (if requested) calculated successfully The 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. A singular matrix was encountered during the optimization. The model was not fitted for this value of . Some 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. The 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. 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.
Given a vector of observed values, , an design matrix , a column vector, , of length holding the th row of and a quantile , g02qg estimates the -element vector as the solution to
where is the piecewise linear loss function , and is an indicator function taking the value if and otherwise. Weights can be incorporated by replacing and with and respectively, where is an 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
where and is a function of , as described below. Given an estimate of the covariance matrix, , lower () and upper () limits for an confidence interval can be calculated for each of the parameters, via
where is the percentile of the Student's distribution with degrees of freedom, where is the rank of the cross-product matrix .
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
where is the sparsity function. g02qg estimates from the residuals, and a bandwidth .
Powell (1991) suggested estimating the matrix by a kernel estimator of the form
where is a kernel function and satisfies and . When the Powell method is chosen, g02qg uses a Gaussian kernel (i.e., ) and sets
where is a bandwidth, and are, respectively, the standard deviation and the and quantiles for the residuals, .
Koenker (2005) suggested estimating the matrix using
where is a bandwidth and denotes the parameter estimates obtained from a quantile regression using the th quantile. Similarly with .
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 are taken from the paired data (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, and , can be obtained directly from the and 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, . Two alternative methods for determining this are provided:
for a user-supplied value ,
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 (see [Optional Parameters] for details on the available options). If bootstrap confidence limits are required () 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.
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
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.
- On entry, or .
- On entry, or .
- On entry, or .
- On entry, .
- On entry, .
- On entry, , .
- On entry, , .
- On entry, or .
- On entry, or .
- On entry, and for at least one .
- On entry, the effective number of observations is less than two.
- On entry, .
- On entry, tau is invalid.
- On entry, and state was not initialized or has been corrupted.
- On exit, problems were encountered whilst fitting at least one model. Additional information has been returned in info.
g02qg allocates internally approximately the following elements of real storage: . If then a further elements are required, and this increases by if . Where possible, any user-supplied output arrays are used as workspace and so the amount actually allocated may be less. If , , and an internal copy of the input data is avoided and the amount of locally allocated memory is reduced by .
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.
By the addition of slack variables the minimization (1) can be reformulated into the linear programming problem
and its associated dual
where is a vector of s. Setting gives the equivalent formulation
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).
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).
Rather than optimize (4) directly, an additional slack variable is added and the constraint is replaced with , for .
The positivity constraint on and is handled using the logarithmic barrier function
The primal-dual form of the problem is used giving the Lagrangian
whose central path is described by the following first order conditions
where denotes the diagonal matrix with diagonal elements given by , similarly with and . By enforcing the inequalities on and strictly, i.e., and for all we ensure that and are positive definite diagonal matrices and hence and exist.
Rather than applying Newton's method to the system of equations given in (5) to obtain the step directions and , Mehrotra substituted the steps directly into (5) giving the augmented system of equations
where and denote the diagonal matrices with diagonal elements given by and respectively.
The affine scaling step is constructed by setting in (5) and applying Newton's method to obtain an intermediate set of step directions
Initial step sizes for the primal () and dual () parameters are constructed as
where is a user-supplied scaling factor. If then the nonlinearity adjustment, described in [Nonlinearity Adjustment], is not made and the model parameters are updated using the current step size and directions.
In the nonlinearity adjustment step a new estimate of is obtained by letting
and estimating as
This estimate, along with the nonlinear terms (, , and ) from (6) are calculated using the values of and obtained from the affine scaling step.
Given an updated estimate for and the nonlinear terms the system of equations
are solved and updated values for and calculated.
At each iteration the model parameters are updated using step directions, and step lengths .
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 the duality gap can be calculated from equations (2) and (3) as
and the optimization terminates if the duality gap is smaller than the tolerance supplied in the optional parameter Tolerance.
Initial values are required for the parameters and . If not supplied by the user, initial values for are calculated from a least squares regression of on . This regression is carried out by first constructing the cross-product matrix and then using a pivoted 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, , a matrix formed from the first -rank columns of is used instead, where is the pivot matrix used during the 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 and are calculated from the residuals. If then a value of is used instead, where is supplied in the optional parameter Epsilon. The initial values for the and are always set to 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).
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, , which g02qg estimates as follows:
When using the Powell Sandwich to estimate the matrix , the quantity
is calculated. Dependent on the value of and the method used to calculate the bandwidth (), it is possible for the quantities to be too large or small, compared to machine precision (). More specifically, when , or , a warning flag is raised in info, the value is truncated to or respectively and the covariance matrix calculated as usual.
The Hendricks–Koenker Sandwich requires the calculation of the quantity . As with the Powell Sandwich, in cases where , or , a warning flag is raised in info, the value truncated to or respectively and the covariance matrix calculated as usual.
In addition, it is required that , in this method. Hence, instead of using in the calculation of , is used instead, where is supplied in the optional parameter Epsilon.
See the description of the optional argument Monitoring.