(* :Title: QuantileRegression *) (* :Author: Johannes Ludsteck *) (* $Log: QuantileRegression.m,v $ Revision 2.0 2004/05/25 05:50:19 katz LinearProgramming is replaced by DualLinearProgramming making safe use of interior point methods possible. Method -> "InteriorPoint" in DualLinearProgramming is now default and cannot be changed by the user. QuantileRegression is now extremely fast. *) (* :Mathematica Version: 5.0.1 *) (* :References: e.g. Portnoy, Stephen & Koenker, Roger (1997): The Gaussian Hare and the Laplacian Tortoise: Computability of Squared-Errors versus Absolute Errors Estimators in: Statistical Science, Vol 12, No. 4, 279-300. *) (* :Limitations: [1] The program does not perform sophisticated input checks. We check only whether the regressor matrix has full rank. I.e. the user is responsible for sensible input data and additional diagnostic checks. *) BeginPackage["QuantileRegression`QuantileRegression`", "LinearAlgebra`MatrixManipulation`"]; Unprotect[{QuantileRegression, RegressionQuantile}]; QuantileRegression::usage = "QuantileRegression[x,y] computes the coefficients of the linear quantile regression of y on x. If the option Quantile-> q with 0 < q < 1 is given, the q-quantile is used. Otherwise the default value q == 0.5 is used, i.e. the median-regression is computed. y is the dependent variable and x the matrix of regressors with Dimensions[x] == {n,p} where n denotes the number of observations and p the number of regressors."; Options[QuantileRegression] = { RegressionQuantile -> 0.5, CovarianceMatrixEstimate -> "Asymptotic" }; Begin["`Private`"]; < "InteriorPoint"]; MakeOutput[x, y, -sol[[2]], q, opts]]; MakeOutput[x_, y_, beta_, q_, opts___]:= Module[{covopt}, covopt = CovarianceMatrixEstimate /. {opts} /. Options[QuantileRegression]; Which[covopt === None, beta, covopt === "Asymptotic", {beta, AsymptoticCovariance[x,y,beta,q]}, MatchQ[covopt,{"Bootstrap", r_Integer/;Positive[r]}], {beta, BSCovariance[x,y, Last[covopt], opts]}, True, Message[QuantileRegression::"covopt"]; beta]]; QuantileRegression::"bsfail" = "Estimation failed in `1` percent of \ bootstrap replications. A failure ratio above 5 percent may indicate \ severe estimation problems and may make the consistency \ of bootstrap results questionable."; QuantileRegression::"covopt" = "Options for covariance matrix estimation are incorrect. \ Only the coefficient vector estimate is returned. \ To obtain an estimate of the covariance matrix check your options and retry."; QuantileRegression::"collin" = "Regressor Matrix does not have full rank indicating collinearity of regressors. Check regressor matrix or the model specification and retry."; AsymptoticCovariance[x_, y_, b_, q_]:= Module[{u=y-x.b, n=Length[y], k}, k = Min[n, Round[3.0 * Sqrt[n]]]; q*(1-q) / (kNNDens[u, q, k]^2 * Transpose[x].x)] (* kNN-Estimate of the density of residual vector u in point F[q] *) kNNDens[x_, q_, k_] := k/(2 Length[x] Sort[Abs[x-Quantile[x,q]]][[k]]); BSCovariance[x_, y_, nReplic_, opts___] := Module[{n = Length[y], s, t, rep = 0, fail = 0, result, bsopts}, bsopts = DeleteCases[{opts}, HoldPattern[_->"Bootstrap"] | HoldPattern[_->{"Bootstrap", _ }]]; bsopts = Append[bsopts, CovarianceMatrixEstimate-> None]; t = Reap[While[rep < nReplic, s = Table[Random[Integer, {1, n}], {n}]; result = Check[QuantileRegression[x[[s]], y[[s]], Sequence@@bsopts], $Failed, LinearSolve::"luc", LinearSolve::"nosol"]; If[result === $Failed, fail++, Sow[result]; rep++]]][[2,1]]; Message[QuantileRegression::"bsfail", 100.0 fail/nReplic]; CovarianceMatrix[t]]; collinear[x_] := If[MatrixRank[x,Tolerance->0.00001] == Min[Dimensions[x]], False, Message[QuantileRegression::"collin"]; True]; SetAttributes[{QuantileRegression}, {Protected, ReadProtected}]; End[] EndPackage[]