1
$\begingroup$

I am trying to compare the results of estimators from least squares and quantile regression from the perspective of the mean square error (MSE). The steps for performing this Monte Carlo approach are:

  1. Set parameters;
  2. Generate data (x and y to be more specific, for three different betas). xnod and tau*quantile will be used later to build true estimator;
  3. Run QuantileRegressionFit and extract parameters;
  4. Run LinearModelFit and extract parameters, including errors. These errors are needed to build the quantile for least squares results (just add the quantile of errors to the final equation).
  5. Rebuild the true y's at a fixed point (xnod) and compare with yqr (rebuilt using quantile regression parameters) and with yls at the same xnod.
  6. Compare calculating MSE.

My supervisor believes that my results are wrong because the MSEs are too large, except for the first, $\beta_{low}$. I have examined this code and fine-tuned it many times. Can someone please help me by pointing out where are my errors and suggestions to fix them?

AppendTo[$Path, FileNameJoin[{"D:", "Wolfram Mathematica"}]];
Needs["QuantileRegression`"]

Set n, m and τ

n = 1000;
τ = 0.95;
columns = 10000;
m = columns;

Generate data

SeedRandom[41];
xdata = Table[RandomVariate[NormalDistribution[], n], columns];
xnod = InverseCDF[NormalDistribution[], τ];
ϵdata = Table[RandomVariate[NormalDistribution[], n], columns];
τquantile = InverseCDF[NormalDistribution[], τ];
βdata = {1/3, 1, 3};
numβ = Length[βdata];
ydata = Table[xdata βdata[[k]] + ϵdata, {k, numβ}];
data = Table[Transpose[{ydata[[q, k]], xdata[[k]]}], {q, numβ}, {k, columns}];

Quantile Regression

funcs = {1, x};
qrFuncs = Table[QuantileRegressionFit[#, funcs, x, {τ}] & /@ data[[q]], {q, numβ}];
qrFuncsFlat = Flatten[qrFuncs];
parametersQRflat = 
  MapThread[Join, 
    {Map[Cases[#, _?NumberQ] &, qrFuncsFlat] /. {} -> {0}, 
     Outer[Coefficient, qrFuncsFlat, DeleteCases[funcs, _?NumberQ]]}];
parametersQR = Partition[parametersQRflat, columns];

Run Least Squares

lsFunc = Table[LinearModelFit[#, x, x] & /@ data[[q]], {q, numβ}];
parametersLs = Table[lsFunc[[q, j]]["BestFitParameters"], {q, numβ}, {j, columns}];
residualsLs = Table[lsFunc[[q, j]]["FitResiduals"], {q, numβ}, {j, columns}];

Estimated quantiles

qϵhat = Table[Quantile[#, τ] & /@ residualsLs[[q]], {q, numβ}];

Quantile regression estimator, least squares estimator and true estimator

yqr = 
  Table[
    parametersQR[[q, k, 1]] + parametersQR[[q, k, 2]] xnod, 
    {q, numβ}, {k, columns}];
yls = 
   Table[
     qϵhat[[q, k]] + parametersLs[[q, k, 2]] xnod + parametersLs[[q, k, 1]], 
     {q, numβ}, {k, columns}];
yactual = τquantile + xnod βdata;

MSE

mseQR = Table[Total[(yqr[[q]] - yactual[[q]])^2]/m, {q, numβ}];
mseLS = Table[Total[(yls[[q]] - yactual[[q]])^2]/m, {q, numβ}];

Bias^2

biasQR = Table[(Total[yqr[[q]] - yactual[[q]]]/m)^2, {q, numβ}];
biasLS = Table[(Total[yls[[q]] - yactual[[q]]]/m)^2, {q, numβ}];

Variance

varQR = mseQR - biasQR;
varLS = mseLS - biasLS;

Conjecture and table of results

TableResults

10000 columns took my computer about 6 - 7 hours to run. I have included the issues that occurred with QuantileRegressionFit during calculation:

qrFuncs = Table[QuantileRegressionFit[#, funcs, x, {τ}] & /@ data[[q]], {q, numβ}];

LinearProgramming::lpipcv: Warning: the interior point algorithm cannot converge to the tolerance of 1.4901161193847656*^-8. The best residual achieved is 0.00004395507383878355, and the solution at that residual has been returned. Setting the option Method -> RevisedSimplex should give a more definite answer, though large problems may take longer computing time.

LinearProgramming::lpipcv: Warning: the interior point algorithm cannot converge to the tolerance of 1.4901161193847656*^-8. The best residual achieved is 6.671193676426196*^-7, and the solution at that residual has been returned. Setting the option Method -> RevisedSimplex should give a more definite answer, though large problems may take longer computing time.

LinearProgramming::lpipcv: Warning: the interior point algorithm cannot converge to the tolerance of 1.4901161193847656*^-8. The best residual achieved is 2.2673711534277855*^-6, and the solution at that residual has been returned. Setting the option Method -> RevisedSimplex should give a more definite answer, though large problems may take longer computing time.

$\endgroup$

closed as off-topic by Kuba May 27 '18 at 7:53

  • The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    $\begingroup$ I figured it out!! When Transposing the data to prepare them for the regression functions, I switched the order. X must come first!! Cheers $\endgroup$ – Thadeu Freitas Filho May 27 '18 at 6:47
  • $\begingroup$ I'm voting to close this question as off-topic because OP found a mistake and the topic is too specific to help future visitors. $\endgroup$ – Kuba May 27 '18 at 7:53