"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/plugin/quantreg.c" (29 Aug 2019, 38231 Bytes) of package /linux/misc/gretl-2020e.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "quantreg.c" see the Fossies "Dox" file reference documentation.

    1 /*
    2  *  gretl -- Gnu Regression, Econometrics and Time-series Library
    3  *  Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
    4  *
    5  *  This program is free software: you can redistribute it and/or modify
    6  *  it under the terms of the GNU General Public License as published by
    7  *  the Free Software Foundation, either version 3 of the License, or
    8  *  (at your option) any later version.
    9  *
   10  *  This program is distributed in the hope that it will be useful,
   11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   13  *  GNU General Public License for more details.
   14  *
   15  *  You should have received a copy of the GNU General Public License
   16  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
   17  *
   18  */
   19 
   20 /* Please see the additional license, rq/uiuc-ncsa.txt: this license
   21    governs the Fortran code, authored by Roger Koenker, which is
   22    called by this module.
   23 */
   24 
   25 #include "libgretl.h"
   26 #include "version.h"
   27 #include "gretl_f2c.h"
   28 #include "usermat.h"
   29 #include "matrix_extra.h"
   30 #include "libset.h"
   31 
   32 #include <errno.h>
   33 
   34 #define QDEBUG 0
   35 
   36 /* Frisch-Newton algorithm: we use this if we're not computing
   37    rank-inversion confidence intervals.
   38 */
   39 
   40 extern int rqfnb_ (integer *n,    /* number of observations */
   41            integer *p,    /* number of parameters */
   42            double *a,     /* transposed matrix of regressors */
   43            double *y,     /* dependent variable vector */
   44            double *rhs,
   45            double *d,
   46            double *u,
   47            double *beta,
   48            double *eps,    /* tolerance */
   49            double *wn,     /* work array, length n */
   50            double *wp,     /* work array, length p */
   51            integer *nit,   /* iteration counts */
   52            integer *info, /* exit status */
   53            void (*callback)(void));
   54 
   55 /* Modified simplex, a la Barrodale-Roberts: this variant lets us get
   56    confidence intervals using the rank-inversion method.
   57 */
   58 
   59 extern int rqbr_ (int n,         /* number of observations */
   60           int p,         /* number of parameters */
   61           double *x,     /* matrix of regressors */
   62           double *y,     /* dependent var vector */
   63           double tau,    /* the desired quantile */
   64           double tol,    /* machine precision to the 2/3 */
   65           double *coeff, /* p-vector of parameter estimates */
   66           double *resid, /* n-vector of residuals */
   67           int *s,        /* int work array, length n */
   68           double *wa,    /* real work array, size (n+5) * (p+4) */
   69           double *wb,    /* and another, size n */
   70           double *sol,   /* primal solution array, size (p + 3) * 2 */
   71           double *dsol,  /* dual solution array, size n * 2 */
   72           int *h,        /* matrix of observation indices, size p * nsol */
   73           double *qn,    /* vector of residual variances from projection of
   74                     each column of x on the remaining columns */
   75           double cutoff, /* critical point for interval */
   76           double *ci,    /* matrix of confidence intervals, size 4 * p */
   77           double *tnmat, /* matrix of JGPK rank-test statistics */
   78           double big,    /* "large positive floating-point number" */
   79           int rmax,      /* max iterations (added for gretl) */
   80           int ci1,       /* doing confidence intervals? */
   81           void (*callback)(void));
   82 
   83 /* machine precision to the 2/3 */
   84 #define calc_eps23 (pow(2.22045e-16, 2/3.0))
   85 
   86 /* wrapper struct for use with Barrodale-Roberts */
   87 
   88 struct br_info {
   89     int warning;
   90     int rmax;
   91     int n, p;
   92     int n5, p3, p4;
   93     int nsol, ndsol;
   94     double tau;
   95     double tol;
   96     double big;
   97     double cut;
   98     double *rspace;
   99     double *coeff, *resid;
  100     double *wa, *wb;
  101     double *qn;
  102     double *sol, *dsol;
  103     int *ispace;
  104     int *s, *h;
  105     gretl_matrix *ci;
  106     gretl_matrix *tnmat;
  107     void (*callback)();
  108 };
  109 
  110 /* wrapper struct for use with Frisch-Newton method */
  111 
  112 struct fn_info {
  113     integer n, p;
  114     double tau;
  115     double beta;
  116     double eps;
  117     double *rspace;
  118     double *rhs;
  119     double *d;
  120     double *u;
  121     double *resid;
  122     double *coeff;
  123     integer nit[3];
  124     integer info;
  125     void (*callback)();
  126 };
  127 
  128 /* destructors for the above wrappers */
  129 
  130 static void fn_info_free (struct fn_info *rq)
  131 {
  132     free(rq->rspace);
  133 }
  134 
  135 static void br_info_free (struct br_info *rq)
  136 {
  137     free(rq->rspace);
  138     free(rq->ispace);
  139     gretl_matrix_free(rq->ci);
  140     gretl_matrix_free(rq->tnmat);
  141 }
  142 
  143 static int real_br_calc (gretl_matrix *y, gretl_matrix *X,
  144              double tau, struct br_info *rq,
  145              int calc_ci);
  146 
  147 /* Bandwidth selection for sparsity estimation, as per Hall, P. and
  148    Sheather, S. J. (1988), "On the distribution of a studentized
  149    quantile", Journal of the Royal Statistical Society, Series B, 50,
  150    381–391: rate = O(n^{-1/3})
  151 */
  152 
  153 static double hs_bandwidth (double tau, int n, int *err)
  154 {
  155     double alpha = 0.05;
  156     double x0 = normal_cdf_inverse(tau);
  157     double f0 = normal_pdf(x0);
  158     double b1 = pow(n, -1 / 3.0);
  159     double b2 = pow(normal_cdf_inverse(1 - alpha/2), 2 / 3.0);
  160     double b3 = (1.5 * f0 * f0) / (2 * x0 * x0 + 1);
  161     double h = b1 * b2 * pow(b3, 1 / 3.0);
  162 
  163 #if QDEBUG
  164     fprintf(stderr, "tau = %g, hs bandwidth = %g\n", tau, h);
  165 #endif
  166 
  167     if (err != NULL) {
  168     if (tau + h > 1) {
  169         gretl_errmsg_set("Hall-Sheather bandwidth is out of bounds");
  170         fprintf(stderr, "hs_bandwidth: tau + h > 1\n");
  171         *err = E_DATA;
  172     } else if (tau - h < 0) {
  173         gretl_errmsg_set("Hall-Sheather bandwidth is out of bounds");
  174         fprintf(stderr, "hs_bandwidth: tau - h < 0\n");
  175         *err = E_DATA;
  176     }
  177     }
  178 
  179     return h;
  180 }
  181 
  182 /* Compute the loglikelihood for a quantile model */
  183 
  184 static double rq_loglik (MODEL *pmod, double tau)
  185 {
  186     double R = 0.0;
  187     int n = pmod->nobs;
  188     int t;
  189 
  190     for (t=pmod->t1; t<=pmod->t2; t++) {
  191     if (!na(pmod->uhat[t])) {
  192         R += pmod->uhat[t] * (tau - (pmod->uhat[t] < 0));
  193     }
  194     }
  195 
  196     return n * (log(tau * (1-tau)) - 1 - log(R/n));
  197 }
  198 
  199 enum {
  200     RQ_STAGE_1,
  201     RQ_STAGE_2,
  202     RQ_LAD
  203 };
  204 
  205 /* Transcribe quantreg results into model struct.  Note: we have to be
  206    careful to replace or invalidate any values generated by the
  207    initial OLS run that are not appropriate in the quantreg context,
  208    so we don't serve up misleading values when the user calls a
  209    model-data accessor.
  210 */
  211 
  212 static void rq_transcribe_results (MODEL *pmod,
  213                    const gretl_matrix *y,
  214                    double tau,
  215                    const double *b,
  216                    const double *u,
  217                    int stage)
  218 {
  219     double SAR = 0.0;
  220     int i, s, t;
  221 
  222     if (stage == RQ_STAGE_1) {
  223     gretl_model_set_double(pmod, "tau", tau);
  224     }
  225 
  226     for (i=0; i<pmod->ncoeff; i++) {
  227     pmod->coeff[i] = b[i];
  228     if (stage == RQ_STAGE_1 || stage == RQ_LAD) {
  229         pmod->sderr[i] = NADBL; /* these will be replaced later */
  230     }
  231     }
  232 
  233     pmod->ess = 0.0;
  234 
  235     s = 0;
  236     for (t=pmod->t1; t<=pmod->t2; t++) {
  237     if (!na(pmod->yhat[t])) {
  238         pmod->uhat[t] = u[s];
  239         pmod->yhat[t] = y->val[s] - u[s];
  240         SAR += fabs(u[s]);
  241         pmod->ess += u[s] * u[s];
  242         s++;
  243     }
  244     }
  245 
  246     /* sum of absolute residuals */
  247     gretl_model_set_double(pmod, "ladsum", SAR);
  248 
  249     /* invalidate unused stats */
  250     pmod->rsq = pmod->adjrsq = NADBL;
  251     pmod->fstt = pmod->chisq = NADBL;
  252 
  253     /* LaPlace errors: equivalent of standard error is sum of
  254        absolute residuals over nobs */
  255     pmod->sigma = SAR / pmod->nobs;
  256     pmod->lnL = rq_loglik(pmod, tau);
  257     mle_criteria(pmod, 0);
  258 }
  259 
  260 /* Extract the interpolated lower and upper bounds from the matrix ci
  261    into a new matrix and attach this to the model for printing.
  262 */
  263 
  264 static int rq_attach_intervals (MODEL *pmod, struct br_info *rq,
  265                 double alpha, gretlopt opt)
  266 {
  267     gretl_matrix *ci = rq->ci;
  268     gretl_matrix *rqci;
  269     double blo, bhi;
  270     int i, p = ci->cols;
  271     int err = 0;
  272 
  273     rqci = gretl_matrix_alloc(p, 2);
  274     if (rqci == NULL) {
  275     return E_ALLOC;
  276     }
  277 
  278     for (i=0; i<p; i++) {
  279     blo = gretl_matrix_get(ci, 1, i);
  280     bhi = gretl_matrix_get(ci, 2, i);
  281     if (na(blo) || na(bhi) || (blo == 0.0 && bhi == 0.0)) {
  282         err = E_NOCONV;
  283         break;
  284     }
  285     gretl_matrix_set(rqci, i, 0, blo);
  286     gretl_matrix_set(rqci, i, 1, bhi);
  287     }
  288 
  289     if (err) {
  290     gretl_matrix_free(rqci);
  291     } else {
  292     gretl_model_set_matrix_as_data(pmod, "coeff_intervals", rqci);
  293     gretl_model_set_double(pmod, "rq_alpha", alpha);
  294     }
  295 
  296     return err;
  297 }
  298 
  299 static void correct_multi_tau_model (MODEL *pmod)
  300 {
  301     pmod->lnL = NADBL;
  302     pmod->ess = NADBL;
  303     pmod->sigma = NADBL;
  304     pmod->rsq = NADBL;
  305     pmod->adjrsq = NADBL;
  306     pmod->fstt = NADBL;
  307 
  308     pmod->criterion[C_AIC] = NADBL;
  309     pmod->criterion[C_BIC] = NADBL;
  310     pmod->criterion[C_HQC] = NADBL;
  311 
  312     free(pmod->uhat);
  313     free(pmod->yhat);
  314     free(pmod->xpx);
  315     free(pmod->vcv);
  316 
  317     pmod->uhat = NULL;
  318     pmod->yhat = NULL;
  319     pmod->xpx = NULL;
  320     pmod->vcv = NULL;
  321 }
  322 
  323 /* Attach the special results matrix generated when we estimate the
  324    model for several values of tau. In this case we retain the
  325    OLS coefficients and standard errors for comparison, but
  326    otherwise we invalidate most statistics.
  327 */
  328 
  329 static int rq_attach_multi_results (MODEL *pmod,
  330                     const gretl_vector *tauvec,
  331                     gretl_matrix *tbeta,
  332                     double alpha, gretlopt opt)
  333 {
  334     gretl_vector *tcpy;
  335 
  336 #if QDEBUG
  337     gretl_matrix_print(tauvec, "tauvec");
  338     gretl_matrix_print(tbeta, "tbeta");
  339 #endif
  340 
  341     tcpy = gretl_matrix_copy(tauvec);
  342     gretl_model_set_matrix_as_data(pmod, "rq_tauvec", tcpy);
  343     gretl_model_set_matrix_as_data(pmod, "rq_sequence", tbeta);
  344 
  345     if (alpha > 0) {
  346     gretl_model_set_double(pmod, "rq_alpha", alpha);
  347     }
  348 
  349     correct_multi_tau_model(pmod);
  350 
  351     return 0;
  352 }
  353 
  354 static int rq_interpolate_intervals (struct br_info *rq)
  355 {
  356     double c1j, c2j, c3j, c4j;
  357     double tn1, tn2, tn3, tn4;
  358     int p = gretl_matrix_cols(rq->ci);
  359     int j;
  360 
  361     for (j=0; j<p; j++) {
  362     c1j = gretl_matrix_get(rq->ci, 0, j);
  363     c2j = gretl_matrix_get(rq->ci, 1, j);
  364     c3j = gretl_matrix_get(rq->ci, 2, j);
  365     c4j = gretl_matrix_get(rq->ci, 3, j);
  366     tn1 = gretl_matrix_get(rq->tnmat, 0, j);
  367     tn2 = gretl_matrix_get(rq->tnmat, 1, j);
  368     tn3 = gretl_matrix_get(rq->tnmat, 2, j);
  369     tn4 = gretl_matrix_get(rq->tnmat, 3, j);
  370 
  371     c3j += fabs(c4j - c3j) * (rq->cut - fabs(tn3)) / fabs(tn4 - tn3);
  372     c2j -= fabs(c1j - c2j) * (rq->cut - fabs(tn2)) / fabs(tn1 - tn2);
  373 
  374     /* check for garbage */
  375     if (na(c3j) || na(c2j)) {
  376         fprintf(stderr, "rq_interpolate_intervals: coeff %d: low = %g, "
  377             "high = %g\n", j, c2j, c3j);
  378         gretl_matrix_print(rq->ci, "rq->ci");
  379         gretl_matrix_print(rq->tnmat, "rq->tnmat");
  380         gretl_errmsg_set(_("Couldn't calculate confidence intervals "
  381                    "for this model"));
  382         return E_NAN;
  383     }
  384 
  385     /* Write the (1 - alpha) intervals into rows 1 and 2
  386        of the matrix ci.
  387     */
  388     gretl_matrix_set(rq->ci, 2, j, c3j);
  389     gretl_matrix_set(rq->ci, 1, j, c2j);
  390     }
  391 
  392     return 0;
  393 }
  394 
  395 static void bad_f_count (const gretl_matrix *f)
  396 {
  397     int n = gretl_vector_get_length(f);
  398     int i, badf = 0;
  399 
  400     for (i=0; i<n; i++) {
  401     if (f->val[i] <= 0) {
  402         badf++;
  403     }
  404     }
  405 
  406     if (badf > 0) {
  407     fprintf(stderr, "Warning: %g percent of fi's <= 0\n",
  408         100 * (double) badf / n);
  409     }
  410 }
  411 
  412 /* prep for robust version of rank-inversion confidence interval
  413    calculation
  414 */
  415 
  416 static int make_nid_qn (gretl_matrix *y, gretl_matrix *X,
  417             struct br_info *rq)
  418 {
  419     int n = rq->n;
  420     int p = rq->p;
  421     double eps = rq->tol;
  422     gretl_matrix *b = NULL;
  423     gretl_matrix *f = NULL;
  424     gretl_matrix *xj = NULL;
  425     gretl_matrix *uj = NULL;
  426     gretl_matrix *Xcmp = NULL;
  427     double h, fi, ui, xik;
  428     int i, j, k, jj;
  429     int err = 0;
  430 
  431     h = hs_bandwidth(rq->tau, n, &err);
  432     if (err) {
  433     return err;
  434     }
  435 
  436     b    = gretl_matrix_alloc(p, 1);
  437     f    = gretl_matrix_alloc(n, 1);
  438     xj   = gretl_matrix_alloc(n, 1);
  439     uj   = gretl_matrix_alloc(n, 1);
  440     Xcmp = gretl_matrix_alloc(n, p - 1);
  441 
  442     if (b == NULL || f == NULL || xj == NULL ||
  443     uj == NULL || Xcmp == NULL) {
  444     err = E_ALLOC;
  445     goto bailout;
  446     }
  447 
  448     real_br_calc(y, X, rq->tau + h, rq, 0);
  449     for (j=0; j<p; j++) {
  450     b->val[j] = rq->coeff[j];
  451     }
  452 
  453     real_br_calc(y, X, rq->tau - h, rq, 0);
  454     for (j=0; j<p; j++) {
  455     b->val[j] -= rq->coeff[j];
  456     }
  457 
  458 #if QDEBUG
  459     gretl_matrix_print(b, "coeff diffs");
  460 #endif
  461 
  462     gretl_matrix_multiply(X, b, f);
  463 
  464     bad_f_count(f);
  465 
  466     for (i=0; i<n; i++) {
  467     fi = (2 * h) / (f->val[i] - eps);
  468     fi = (fi < eps)? eps : fi;
  469     f->val[i] = sqrt(fi);
  470     }
  471 
  472     /* Now set each qn[j] to the SSR from an f-weighted regression of
  473        X_j on the other regressors.
  474     */
  475 
  476     gretl_matrix_reuse(b, p - 1, 1);
  477 
  478     for (j=0; j<p; j++) {
  479     for (i=0; i<n; i++) {
  480         /* weighted dependent var */
  481         xj->val[i] = f->val[i] * gretl_matrix_get(X, i, j);
  482     }
  483 
  484     jj = 0;
  485     for (k=0; k<p; k++) {
  486         if (k != j) {
  487         /* weighted independent var */
  488         for (i=0; i<n; i++) {
  489             xik = gretl_matrix_get(X, i, k);
  490             gretl_matrix_set(Xcmp, i, jj, f->val[i] * xik);
  491         }
  492         jj++;
  493         }
  494     }
  495 
  496     err = gretl_matrix_ols(xj, Xcmp, b, NULL, uj, NULL);
  497 
  498     if (!err) {
  499         rq->qn[j] = 0.0;
  500         for (i=0; i<n; i++) {
  501         ui = uj->val[i] / f->val[i];
  502         rq->qn[j] += ui * ui;
  503         }
  504     }
  505     }
  506 
  507  bailout:
  508 
  509     gretl_matrix_free(b);
  510     gretl_matrix_free(f);
  511     gretl_matrix_free(xj);
  512     gretl_matrix_free(uj);
  513     gretl_matrix_free(Xcmp);
  514 
  515     return err;
  516 }
  517 
  518 /* Get {X'X}^{-1}, handling both the B-R case, where the X matrix
  519    has the shape you'd expect, and the F-N case, where "X" is
  520    actually X-transpose.
  521 */
  522 
  523 static gretl_matrix *get_XTX_inverse (const gretl_matrix *X, int *err)
  524 {
  525     int k = min(X->rows, X->cols);
  526     GretlMatrixMod flag1, flag2;
  527     gretl_matrix *XTX;
  528 
  529     XTX = gretl_matrix_alloc(k, k);
  530     if (XTX == NULL) {
  531     *err = E_ALLOC;
  532     return NULL;
  533     }
  534 
  535     flag1 = (X->cols == k)? GRETL_MOD_TRANSPOSE : GRETL_MOD_NONE;
  536     flag2 = (X->cols == k)? GRETL_MOD_NONE : GRETL_MOD_TRANSPOSE;
  537 
  538     *err = gretl_matrix_multiply_mod(X, flag1,
  539                      X, flag2,
  540                      XTX, GRETL_MOD_NONE);
  541 
  542     if (!*err) {
  543     *err = gretl_invert_symmetric_matrix(XTX);
  544     }
  545 
  546     if (*err) {
  547     gretl_matrix_free(XTX);
  548     XTX = NULL;
  549     }
  550 
  551     return XTX;
  552 }
  553 
  554 /* prep for IID version of rank-inversion confidence interval
  555    calculation
  556 */
  557 
  558 static int make_iid_qn (const gretl_matrix *X, double *qn)
  559 {
  560     int i, p = gretl_matrix_cols(X);
  561     gretl_matrix *XTX;
  562     int err = 0;
  563 
  564     XTX = get_XTX_inverse(X, &err);
  565 
  566     if (!err) {
  567     for (i=0; i<p; i++) {
  568         qn[i] = 1 / gretl_matrix_get(XTX, i, i);
  569     }
  570     gretl_matrix_free(XTX);
  571     }
  572 
  573 #if QDEBUG
  574     fprintf(stderr, "make_iid_qn: returning %d\n", err);
  575 #endif
  576 
  577     return err;
  578 }
  579 
  580 /* Allocate workspace to be fed to the function rqbr, and
  581    initialize various things.
  582 */
  583 
  584 static int br_info_alloc (struct br_info *rq, int n, int p,
  585               double tau, double alpha,
  586               gretlopt opt)
  587 {
  588     size_t rsize, isize;
  589 
  590     rq->rspace = NULL;
  591     rq->ispace = NULL;
  592     rq->ci = NULL;
  593     rq->tnmat = NULL;
  594 
  595     rq->n5 = n + 5;
  596     rq->p3 = p + 3;
  597     rq->p4 = p + 4;
  598     rq->nsol = rq->ndsol = 2;
  599 
  600     rsize = p + n + rq->n5 * rq->p4 + n + p + rq->nsol * rq->p3 + rq->ndsol * n;
  601     isize = n + rq->nsol * p;
  602 
  603     rq->rspace = malloc(rsize * sizeof *rq->rspace);
  604     if (rq->rspace == NULL) {
  605     return E_ALLOC;
  606     }
  607 
  608     rq->ispace = malloc(isize * sizeof *rq->ispace);
  609     if (rq->ispace == NULL) {
  610     return E_ALLOC;
  611     }
  612 
  613     if (!(opt & OPT_L)) {
  614     /* doing confidence intervals */
  615     rq->ci = gretl_matrix_alloc(4, p);
  616     rq->tnmat = gretl_matrix_alloc(4, p);
  617 
  618     if (rq->ci == NULL || rq->tnmat == NULL) {
  619         return E_ALLOC;
  620     }
  621     }
  622 
  623     /* real arrays */
  624     rq->coeff = rq->rspace;
  625     rq->resid = rq->coeff + p;
  626     rq->wa    = rq->resid + n;
  627     rq->wb    = rq->wa + rq->n5 * rq->p4;
  628     rq->qn    = rq->wb + n;
  629     rq->sol   = rq->qn + p;
  630     rq->dsol  = rq->sol + rq->nsol * rq->p3;
  631 
  632     /* integer arrays */
  633     rq->s = rq->ispace;
  634     rq->h = rq->s + n;
  635 
  636     rq->warning = 0;
  637     rq->n = n;
  638     rq->p = p;
  639     rq->tau = tau;
  640     rq->tol = calc_eps23;
  641     rq->big = DBL_MAX / 100;
  642     rq->rmax = libset_get_int(RQ_MAXITER);
  643 
  644     if (opt & OPT_L) {
  645     /* doing simple LAD */
  646     rq->cut = 0;
  647     } else if (opt & OPT_N) {
  648     /* asymptotic: no df correction */
  649     rq->cut = normal_cdf_inverse(1 - alpha/2);
  650     } else {
  651     rq->cut = student_cdf_inverse(n - p, 1 - alpha/2);
  652     }
  653 
  654     if (show_activity_func_installed()) {
  655     rq->callback = show_activity_callback;
  656     } else {
  657     rq->callback = NULL;
  658     }
  659 
  660     return 0;
  661 }
  662 
  663 /* Call Barrodale-Roberts code */
  664 
  665 static int real_br_calc (gretl_matrix *y, gretl_matrix *X,
  666              double tau, struct br_info *rq,
  667              int calc_ci)
  668 {
  669     double *ci_val, *tnmat_val;
  670     int ret, err = 0;
  671 
  672 #if QDEBUG
  673     fprintf(stderr, "real_br_calc: calling rqbr, calc_ci = %d\n", calc_ci);
  674 #endif
  675 
  676     /* these inputs are not needed if we're not computing
  677        confidence intervals */
  678     ci_val = (rq->ci == NULL)? NULL : rq->ci->val;
  679     tnmat_val = (rq->tnmat == NULL)? NULL : rq->tnmat->val;
  680 
  681     ret = rqbr_(rq->n, rq->p, X->val, y->val, tau, rq->tol,
  682         rq->coeff, rq->resid, rq->s, rq->wa, rq->wb,
  683         rq->sol, rq->dsol, rq->h, rq->qn, rq->cut,
  684         ci_val, tnmat_val,
  685         rq->big, rq->rmax, calc_ci,
  686         rq->callback);
  687 
  688 #if QDEBUG
  689     fprintf(stderr, "rqbr: ift = %d\n", ret);
  690 #endif
  691 
  692     if (ret == 1) {
  693     rq->warning = 1;
  694     fprintf(stderr, "Warning: solution may be non-unique\n");
  695     } else if (ret == 2){
  696     fprintf(stderr, "Premature end: conditioning problem in X?\n");
  697     err = E_NOCONV;
  698     } else if (ret == 3) {
  699     gretl_errmsg_sprintf("Maximum number of iterations (%d) exceeded",
  700                  (int) rq->rmax);
  701     err = E_NOCONV;
  702     }
  703 
  704     return err;
  705 }
  706 
  707 /* allocate workspace for F-N algorithm */
  708 
  709 static int fn_info_alloc (struct fn_info *rq, int n, int p,
  710               double tau, gretlopt opt)
  711 {
  712     int n10 = n * 10;
  713     int pp4, rp = p;
  714     size_t rsize;
  715 
  716     if (p == 1 && !(opt & OPT_R)) {
  717     /* just one regressor, doing "iid" standard errors:
  718        we'll need workspace for p = 2 */
  719     rp = 2;
  720     }
  721 
  722     pp4 = rp * (rp + 4);
  723     rsize = rp + n + n + n10 + pp4;
  724 
  725     rq->rspace = malloc(rsize * sizeof *rq->rspace);
  726 
  727     if (rq->rspace == NULL) {
  728     return E_ALLOC;
  729     }
  730 
  731     rq->rhs = rq->rspace;
  732     rq->d   = rq->rhs + rp;
  733     rq->u   = rq->d + n;
  734     rq->resid = rq->u + n;
  735     rq->coeff = rq->resid + n10;
  736 
  737     rq->n = n;
  738     rq->p = p;
  739     rq->tau = tau;
  740     rq->beta = .99995;
  741     rq->eps = 1.0e-7;
  742 
  743     if (show_activity_func_installed()) {
  744     rq->callback = show_activity_callback;
  745     } else {
  746     rq->callback = NULL;
  747     }
  748 
  749     return 0;
  750 }
  751 
  752 /* Initialize the tau-dependent arrays rhs and resid (rhs is
  753    also X-dependent) and set other entries to 0 or 1.
  754 */
  755 
  756 static void rq_workspace_init (struct fn_info *rq,
  757                    const gretl_matrix *XT,
  758                    double tau)
  759 {
  760     double xsum;
  761     int p = gretl_matrix_rows(XT);
  762     int n = gretl_matrix_cols(XT);
  763     int n10 = n * 10;
  764     int i, t;
  765 
  766     for (i=0; i<p; i++) {
  767     xsum = 0.0;
  768     for (t=0; t<n; t++) {
  769         xsum += gretl_matrix_get(XT, i, t);
  770     }
  771     rq->rhs[i] = tau * xsum;
  772     }
  773 
  774     for (i=0; i<n; i++) {
  775     rq->d[i] = rq->u[i] = 1.0;
  776     rq->resid[i] = tau;
  777     }
  778 
  779     for (i=n; i<n10; i++) {
  780     rq->resid[i] = 0.0;
  781     }
  782 }
  783 
  784 static int rq_call_FN (integer *n, integer *p, gretl_matrix *XT,
  785                gretl_matrix *y, struct fn_info *rq,
  786                double tau)
  787 {
  788     rq_workspace_init(rq, XT, tau);
  789 
  790     return rqfnb_(n, p, XT->val, y->val, rq->rhs,
  791           rq->d, rq->u, &rq->beta, &rq->eps,
  792           rq->resid, rq->coeff, rq->nit, &rq->info,
  793           rq->callback);
  794 }
  795 
  796 static int rq_write_variance (const gretl_matrix *V,
  797                   MODEL *pmod, double *se)
  798 {
  799     double x;
  800     int i, err = 0;
  801 
  802     if (se != NULL) {
  803     for (i=0; i<V->cols; i++) {
  804         x = gretl_matrix_get(V, i, i);
  805         if (na(x) || x < 0) {
  806         se[i] = NADBL;
  807         } else {
  808         se[i] = sqrt(x);
  809         }
  810     }
  811     } else {
  812     err = gretl_model_write_vcv(pmod, V);
  813     }
  814 
  815     return err;
  816 }
  817 
  818 /* IID version of asymptotic F-N covariance matrix */
  819 
  820 static int rq_fn_iid_VCV (MODEL *pmod, gretl_matrix *y,
  821               gretl_matrix *XT, double tau,
  822               struct fn_info *rq,
  823               double *se)
  824 {
  825     gretl_matrix *V = NULL;
  826     gretl_matrix *vx = NULL;
  827     gretl_matrix *vy = NULL;
  828     gretl_matrix *S0 = NULL;
  829     gretl_matrix *S1 = NULL;
  830     double h, sparsity, eps23;
  831     int i, df, pz = 0;
  832     integer vn, vp = 2;
  833     int err = 0;
  834 
  835     V = get_XTX_inverse(XT, &err);
  836     if (err) {
  837     return err;
  838     }
  839 
  840     h = ceil(rq->n * hs_bandwidth(tau, rq->n, NULL));
  841     h = (rq->p + 1 > h)? rq->p + 1 : h;
  842     vn = h + 1;
  843     eps23 = calc_eps23;
  844 
  845     for (i=0; i<rq->n; i++) {
  846     if (fabs(rq->resid[i]) < eps23) {
  847         pz++;
  848     }
  849     }
  850 
  851 #if QDEBUG
  852     fprintf(stderr, "rq_fn_iid_VCV: eps = %g, h = %g, pz = %d\n", eps23, h, pz);
  853 #endif
  854 
  855     vy = gretl_column_vector_alloc(vn);
  856     vx = gretl_matrix_alloc(2, vn);
  857     S0 = gretl_matrix_alloc(rq->n, 2);
  858 
  859     if (vy == NULL || vx == NULL || S0 == NULL) {
  860     err = E_ALLOC;
  861     goto bailout;
  862     }
  863 
  864     /* construct sequence of length h + 1 */
  865 
  866     for (i=0; i<vn; i++) {
  867     vy->val[i] = i + pz + 1;
  868     }
  869 
  870     /* construct artificial X (transpose) matrix, constant in
  871        first column
  872     */
  873 
  874     df = rq->n - rq->p;
  875     for (i=0; i<vn; i++) {
  876     gretl_matrix_set(vx, 0, i, 1.0);
  877     gretl_matrix_set(vx, 1, i, vy->val[i] / df);
  878     }
  879 
  880     /* sort the residuals by absolute magnitude */
  881 
  882     for (i=0; i<rq->n; i++) {
  883     gretl_matrix_set(S0, i, 0, rq->resid[i]);
  884     gretl_matrix_set(S0, i, 1, fabs(rq->resid[i]));
  885     }
  886 
  887     S1 = gretl_matrix_sort_by_column(S0, 1, &err);
  888     if (err) {
  889     goto bailout;
  890     }
  891 
  892     /* extract the desired range of residuals */
  893 
  894     for (i=0; i<vn; i++) {
  895     vy->val[i] = gretl_matrix_get(S1, pz + i, 0);
  896     }
  897 
  898     /* and re-sort, respecting sign */
  899 
  900     qsort(vy->val, vn, sizeof *vy->val, gretl_compare_doubles);
  901 
  902     /* run artificial L1 regression to get sparsity measure */
  903 
  904     err = rq_call_FN(&vn, &vp, vx, vy, rq, 0.5);
  905 
  906     if (err) {
  907     fprintf(stderr, "rq_fn_iid_VCV: rqfn: info = %d\n", rq->info);
  908     } else {
  909     /* scale X'X-inverse appropriately */
  910     sparsity = rq->coeff[1];
  911     h = sparsity * sparsity * tau * (1 - tau);
  912     gretl_matrix_multiply_by_scalar(V, h);
  913     err = rq_write_variance(V, pmod, se);
  914     }
  915 
  916  bailout:
  917 
  918     gretl_matrix_free(V);
  919     gretl_matrix_free(vy);
  920     gretl_matrix_free(vx);
  921     gretl_matrix_free(S0);
  922     gretl_matrix_free(S1);
  923 
  924     return err;
  925 }
  926 
  927 /* Robust sandwich version of F-N covariance matrix */
  928 
  929 static int rq_fn_nid_VCV (MODEL *pmod, gretl_matrix *y,
  930               gretl_matrix *XT, double tau,
  931               struct fn_info *rq,
  932               double *se)
  933 {
  934     gretl_matrix *p1 = NULL;
  935     gretl_matrix *f = NULL;
  936     gretl_matrix *fX = NULL;
  937     gretl_matrix *fXX = NULL;
  938     gretl_matrix *XTX = NULL;
  939     gretl_matrix *V = NULL;
  940     double h, x, eps23;
  941     int n = rq->n;
  942     int p = rq->p;
  943     int i, t, err = 0;
  944 
  945     h = hs_bandwidth(tau, rq->n, &err);
  946     if (err) {
  947     return err;
  948     }
  949 
  950     eps23 = calc_eps23;
  951 
  952     p1  = gretl_matrix_alloc(p, 1);
  953     f   = gretl_matrix_alloc(n, 1);
  954     fX  = gretl_matrix_alloc(p, n);
  955     fXX = gretl_matrix_alloc(p, p);
  956     XTX = gretl_matrix_alloc(p, p);
  957     V   = gretl_matrix_alloc(p, p);
  958 
  959     if (p1 == NULL || f == NULL || fX == NULL ||
  960     fXX == NULL || XTX == NULL || V == NULL) {
  961     err = E_ALLOC;
  962     goto bailout;
  963     }
  964 
  965     err = rq_call_FN(&n, &p, XT, y, rq, tau + h);
  966     if (err) {
  967     fprintf(stderr, "tau + h: info = %d\n", rq->info);
  968     goto bailout;
  969     }
  970 
  971     /* grab coeffs */
  972     for (i=0; i<p; i++) {
  973     p1->val[i] = rq->coeff[i];
  974     }
  975 
  976     err = rq_call_FN(&n, &p, XT, y, rq, tau - h);
  977     if (err) {
  978     fprintf(stderr, "tau - h: info = %d\n", rq->info);
  979     goto bailout;
  980     }
  981 
  982     /* diff coeffs */
  983     for (i=0; i<p; i++) {
  984     p1->val[i] -= rq->coeff[i];
  985     }
  986 
  987 #if QDEBUG
  988     fprintf(stderr, "rq_fn_nid_VCV: bandwidth = %g\n", h);
  989     gretl_matrix_print(p1, "coeff diffs");
  990 #endif
  991 
  992     gretl_matrix_multiply_mod(XT, GRETL_MOD_TRANSPOSE,
  993                   p1, GRETL_MOD_NONE,
  994                   f, GRETL_MOD_NONE);
  995 
  996     bad_f_count(f);
  997 
  998     for (t=0; t<n; t++) {
  999     x = (2 * h) / (f->val[t] - eps23);
 1000     f->val[t] = (x > 0)? sqrt(x) : 0;
 1001     f->val[t] = (x > 0)? x : 0;
 1002     for (i=0; i<p; i++) {
 1003         x = gretl_matrix_get(XT, i, t);
 1004         x *= f->val[t];
 1005         gretl_matrix_set(fX, i, t, x);
 1006     }
 1007     }
 1008 
 1009     gretl_matrix_multiply_mod(fX, GRETL_MOD_NONE,
 1010                   XT, GRETL_MOD_TRANSPOSE,
 1011                   fXX, GRETL_MOD_NONE);
 1012 
 1013     gretl_invert_symmetric_matrix(fXX);
 1014 
 1015     gretl_matrix_multiply_mod(XT, GRETL_MOD_NONE,
 1016                   XT, GRETL_MOD_TRANSPOSE,
 1017                   XTX, GRETL_MOD_NONE);
 1018 
 1019     gretl_matrix_qform(fXX, GRETL_MOD_NONE,
 1020                XTX, V, GRETL_MOD_NONE);
 1021 
 1022     gretl_matrix_multiply_by_scalar(V, tau * (1 - tau));
 1023 
 1024     err = rq_write_variance(V, pmod, se);
 1025 
 1026  bailout:
 1027 
 1028     gretl_matrix_free(p1);
 1029     gretl_matrix_free(f);
 1030     gretl_matrix_free(fX);
 1031     gretl_matrix_free(fXX);
 1032     gretl_matrix_free(XTX);
 1033     gretl_matrix_free(V);
 1034 
 1035     return err;
 1036 }
 1037 
 1038 static int get_ci_alpha (double *a)
 1039 {
 1040     int err = 0;
 1041     double c;
 1042 
 1043     c = get_optval_double(QUANTREG, OPT_I, &err);
 1044     if (err) {
 1045     return err;
 1046     }
 1047 
 1048     if (na(c)) {
 1049     *a = 0.1;
 1050     } else if (c > 1.0 && c < 100.0) {
 1051     *a = 1 - c / 100;
 1052     } else if (c < 0.1 || c > .999) {
 1053     gretl_errmsg_sprintf("Confidence level out of bounds");
 1054     return E_DATA;
 1055     } else {
 1056     *a = 1 - c;
 1057     }
 1058 
 1059     return 0;
 1060 }
 1061 
 1062 /* write coefficients and lower/upper c.i. values for
 1063    all parameters at a given value of tau */
 1064 
 1065 static int
 1066 write_tbeta_block_br (gretl_matrix *tbeta, int nt, double *coeff,
 1067               gretl_matrix *ci, int k)
 1068 
 1069 {
 1070     double blo, bhi;
 1071     int p = ci->cols;
 1072     int i;
 1073 
 1074     for (i=0; i<p; i++) {
 1075     blo = gretl_matrix_get(ci, 1, i);
 1076     bhi = gretl_matrix_get(ci, 2, i);
 1077     gretl_matrix_set(tbeta, k, 0, coeff[i]);
 1078     gretl_matrix_set(tbeta, k, 1, blo);
 1079     gretl_matrix_set(tbeta, k, 2, bhi);
 1080     k += nt;
 1081     }
 1082 
 1083     return 0;
 1084 }
 1085 
 1086 /* write coefficients or standard errors for
 1087    all parameters at a given value of tau */
 1088 
 1089 static int
 1090 write_tbeta_block_fn (gretl_matrix *tbeta, int nt, double *x,
 1091               integer p, int k, int j)
 1092 
 1093 {
 1094     int i;
 1095 
 1096     for (i=0; i<p; i++) {
 1097     if (na(x[i])) {
 1098         fprintf(stderr, "write_tbeta_block_fn: x[%d] = %g\n", i, x[i]);
 1099         return E_NAN;
 1100     }
 1101     gretl_matrix_set(tbeta, k, j, x[i]);
 1102     k += nt;
 1103     }
 1104 
 1105     return 0;
 1106 }
 1107 
 1108 /* Sub-driver for Barrodale-Roberts estimation, with confidence
 1109    intervals.
 1110 */
 1111 
 1112 static int rq_fit_br (gretl_matrix *y, gretl_matrix *X,
 1113               const gretl_vector *tauvec, gretlopt opt,
 1114               MODEL *pmod)
 1115 {
 1116     struct br_info rq;
 1117     gretl_matrix *tbeta = NULL;
 1118     integer n = y->rows;
 1119     integer p = X->cols;
 1120     double tau, alpha = 0;
 1121     int i, ntau;
 1122     int err = 0;
 1123 
 1124     err = get_ci_alpha(&alpha);
 1125     if (err) {
 1126     return err;
 1127     }
 1128 
 1129     ntau = gretl_vector_get_length(tauvec);
 1130     tau = gretl_vector_get(tauvec, 0);
 1131 
 1132     err = br_info_alloc(&rq, n, p, tau, alpha, opt);
 1133 
 1134     if (!err && ntau > 1) {
 1135     tbeta = gretl_zero_matrix_new(p * ntau, 3);
 1136     if (tbeta == NULL) {
 1137         err = E_ALLOC;
 1138     }
 1139 #if QDEBUG
 1140     fprintf(stderr, "p = %d, ntau = %d, alpha = %g\n", p, ntau, alpha);
 1141     fprintf(stderr, "tbeta = %d x %d\n", tbeta->rows, tbeta->cols);
 1142 #endif
 1143     }
 1144 
 1145     for (i=0; i<ntau && !err; i++) {
 1146     tau = rq.tau = gretl_vector_get(tauvec, i);
 1147 
 1148 #if QDEBUG
 1149     fprintf(stderr, "rq_fit_br: i = %d, tau = %g\n", i, tau);
 1150 #endif
 1151 
 1152     /* preliminary calculations relating to confidence intervals */
 1153     if (opt & OPT_R) {
 1154         /* robust variant */
 1155         err = make_nid_qn(y, X, &rq);
 1156     } else {
 1157         /* assuming iid errors */
 1158         err = make_iid_qn(X, rq.qn);
 1159     }
 1160 
 1161     if (!err) {
 1162         /* get the actual estimates */
 1163         err = real_br_calc(y, X, tau, &rq, 1);
 1164     }
 1165 
 1166     if (!err) {
 1167         /* post-process confidence intervals */
 1168         err = rq_interpolate_intervals(&rq);
 1169     }
 1170 
 1171     if (!err) {
 1172         if (ntau == 1) {
 1173         /* done: put intervals onto the model */
 1174         err = rq_attach_intervals(pmod, &rq, alpha, opt);
 1175         if (!err) {
 1176             rq_transcribe_results(pmod, y, tau, rq.coeff,
 1177                       rq.resid, RQ_STAGE_2);
 1178         }
 1179         } else {
 1180         /* using multiple tau values */
 1181         err = write_tbeta_block_br(tbeta, ntau, rq.coeff, rq.ci, i);
 1182         }
 1183     }
 1184     }
 1185 
 1186     if (!err && rq.warning) {
 1187     gretl_model_set_int(pmod, "nonunique", 1);
 1188     }
 1189 
 1190     if (tbeta != NULL) {
 1191     /* multiple tau values */
 1192     if (err) {
 1193         gretl_matrix_free(tbeta);
 1194     } else {
 1195         err = rq_attach_multi_results(pmod, tauvec, tbeta, alpha, opt);
 1196     }
 1197     }
 1198 
 1199     br_info_free(&rq);
 1200 
 1201     return err;
 1202 }
 1203 
 1204 /* sub-driver for Frisch-Newton interior point variant */
 1205 
 1206 static int rq_fit_fn (gretl_matrix *y, gretl_matrix *XT,
 1207               const gretl_vector *tauvec, gretlopt opt,
 1208               MODEL *pmod)
 1209 {
 1210     struct fn_info rq;
 1211     gretl_matrix *tbeta = NULL;
 1212     double *se = NULL;
 1213     integer n = y->rows;
 1214     integer p = XT->rows;
 1215     double tau;
 1216     int i, ntau;
 1217     int err = 0;
 1218 
 1219     ntau = gretl_vector_get_length(tauvec);
 1220     tau = gretl_vector_get(tauvec, 0);
 1221 
 1222     err = fn_info_alloc(&rq, n, p, tau, opt);
 1223     if (err) {
 1224     return err;
 1225     }
 1226 
 1227     if (ntau > 1) {
 1228     tbeta = gretl_zero_matrix_new(p * ntau, 2);
 1229     se = malloc(p * sizeof *se);
 1230     if (tbeta == NULL || se == NULL) {
 1231         err = E_ALLOC;
 1232     }
 1233     }
 1234 
 1235     for (i=0; i<ntau && !err; i++) {
 1236     tau = rq.tau = gretl_vector_get(tauvec, i);
 1237 
 1238 #if QDEBUG
 1239     fprintf(stderr, "rq_fit_fn: i = %d, tau = %g\n", i, tau);
 1240 #endif
 1241 
 1242     /* get coefficients and residuals */
 1243     err = rq_call_FN(&n, &p, XT, y, &rq, tau);
 1244     if (err) {
 1245         fprintf(stderr, "rqfn gave info = %d\n", rq.info);
 1246     }
 1247 
 1248     if (!err) {
 1249         if (ntau == 1) {
 1250         /* save coeffs, residuals, etc., before computing VCV */
 1251         rq_transcribe_results(pmod, y, tau, rq.coeff, rq.resid,
 1252                       RQ_STAGE_1);
 1253         } else {
 1254         /* write coeffs for this tau value */
 1255         write_tbeta_block_fn(tbeta, ntau, rq.coeff, p, i, 0);
 1256         }
 1257     }
 1258 
 1259     if (!err) {
 1260         /* compute covariance matrix */
 1261         if (opt & OPT_R) {
 1262         err = rq_fn_nid_VCV(pmod, y, XT, tau, &rq, se);
 1263         } else {
 1264         err = rq_fn_iid_VCV(pmod, y, XT, tau, &rq, se);
 1265         }
 1266     }
 1267 
 1268     if (!err && ntau > 1) {
 1269         /* write std errs for this tau */
 1270         write_tbeta_block_fn(tbeta, ntau, se, p, i, 1);
 1271     }
 1272     }
 1273 
 1274     if (tbeta != NULL) {
 1275     /* multiple tau values */
 1276     if (err) {
 1277         gretl_matrix_free(tbeta);
 1278     } else {
 1279         err = rq_attach_multi_results(pmod, tauvec, tbeta, 0, opt);
 1280     }
 1281     }
 1282 
 1283     fn_info_free(&rq);
 1284     free(se);
 1285 
 1286     return err;
 1287 }
 1288 
 1289 /* Write y and X from pmod into gretl matrices, respecting the sample
 1290    range.  Note that for use with the Frisch-Newton version of the
 1291    underlying rq function we want the X matrix in transpose form,
 1292    which is signaled by @tr = 1. Also note that @pmod was first
 1293    estimated via OLS and its sample range may contain interior
 1294    NAs, in which case we have to avoid them when filling the
 1295    data matrices.
 1296 */
 1297 
 1298 static int rq_make_matrices (MODEL *pmod,
 1299                  DATASET *dset,
 1300                  gretl_matrix **py,
 1301                  gretl_matrix **pX,
 1302                  int tr)
 1303 {
 1304     int n = pmod->nobs;
 1305     int p = pmod->ncoeff;
 1306     int yno = pmod->list[1];
 1307     gretl_matrix *X = NULL;
 1308     gretl_matrix *y = NULL;
 1309     int i, s, t, v;
 1310     int err = 0;
 1311 
 1312     y = gretl_matrix_alloc(n, 1);
 1313 
 1314     if (tr) {
 1315     X = gretl_matrix_alloc(p, n);
 1316     } else {
 1317     X = gretl_matrix_alloc(n, p);
 1318     }
 1319 
 1320     if (X == NULL || y == NULL) {
 1321     gretl_matrix_free(y);
 1322     gretl_matrix_free(X);
 1323     return E_ALLOC;
 1324     }
 1325 
 1326     s = 0;
 1327     for (t=pmod->t1; t<=pmod->t2; t++) {
 1328     if (!na(pmod->uhat[t])) {
 1329         gretl_vector_set(y, s++, dset->Z[yno][t]);
 1330     }
 1331     }
 1332 
 1333     for (i=0; i<p; i++) {
 1334     v = pmod->list[i+2];
 1335     s = 0;
 1336     for (t=pmod->t1; t<=pmod->t2; t++) {
 1337         if (!na(pmod->uhat[t])) {
 1338         if (tr) {
 1339             gretl_matrix_set(X, i, s++, dset->Z[v][t]);
 1340         } else {
 1341             gretl_matrix_set(X, s++, i, dset->Z[v][t]);
 1342         }
 1343         }
 1344     }
 1345     }
 1346 
 1347     if (err) {
 1348     gretl_matrix_free(y);
 1349     gretl_matrix_free(X);
 1350     } else {
 1351     *py = y;
 1352     *pX = X;
 1353     }
 1354 
 1355     return err;
 1356 }
 1357 
 1358 static int rq_transpose_X (gretl_matrix **pX)
 1359 {
 1360     gretl_matrix *XT = gretl_matrix_copy_transpose(*pX);
 1361     int err = 0;
 1362 
 1363     if (XT == NULL) {
 1364     err = E_ALLOC;
 1365     } else {
 1366     gretl_matrix_free(*pX);
 1367     *pX = XT;
 1368     }
 1369 
 1370     return err;
 1371 }
 1372 
 1373 static int check_user_tau (const gretl_matrix *tau, int *ntau)
 1374 {
 1375     int err = 0;
 1376 
 1377     *ntau = gretl_vector_get_length(tau);
 1378 
 1379     if (*ntau == 0) {
 1380     err = E_DATA;
 1381     } else {
 1382     double p;
 1383     int i;
 1384 
 1385     for (i=0; i<*ntau; i++) {
 1386         p = gretl_vector_get(tau, i);
 1387         if (p < .01 || p > .99) {
 1388         gretl_errmsg_sprintf("quantreg: tau must be >= .01 and <= .99");
 1389         err = E_DATA;
 1390         }
 1391     }
 1392     }
 1393 
 1394     return err;
 1395 }
 1396 
 1397 int rq_driver (const gretl_matrix *tau, MODEL *pmod,
 1398            DATASET *dset, gretlopt opt, PRN *prn)
 1399 {
 1400     gretl_matrix *y = NULL;
 1401     gretl_matrix *X = NULL;
 1402     int ntau = 0;
 1403     int err;
 1404 
 1405     err = check_user_tau(tau, &ntau);
 1406 
 1407     if (!err && (opt & OPT_I) && pmod->list[0] < 3) {
 1408     gretl_errmsg_set("quantreg: can't do confidence intervals with "
 1409              "only one regressor");
 1410     err = E_DATA;
 1411     }
 1412 
 1413     if (!err) {
 1414     if ((opt & OPT_I) && ntau == 1) {
 1415         /* doing intervals, only one tau: we'll do a first run
 1416            using F-N so the model will be equipped with standard
 1417            errors
 1418         */
 1419         err = rq_make_matrices(pmod, dset, &y, &X, 1);
 1420         if (!err) {
 1421         err = rq_fit_fn(y, X, tau, opt, pmod);
 1422         }
 1423         if (!err) {
 1424         /* flip the X matrix for use with B-R */
 1425         err = rq_transpose_X(&X);
 1426         }
 1427     } else {
 1428         int tr = !(opt & OPT_I);
 1429 
 1430         err = rq_make_matrices(pmod, dset, &y, &X, tr);
 1431     }
 1432     }
 1433 
 1434     if (!err) {
 1435     if (opt & OPT_I) {
 1436         /* doing confidence intervals -> use Barrodale-Roberts */
 1437         err = rq_fit_br(y, X, tau, opt, pmod);
 1438     } else {
 1439         /* otherwise use Frisch-Newton */
 1440         err = rq_fit_fn(y, X, tau, opt, pmod);
 1441     }
 1442     }
 1443 
 1444     if (!err) {
 1445     /* some common finishing touches */
 1446     gretl_model_add_y_median(pmod, dset->Z[pmod->list[1]]);
 1447     pmod->ci = LAD;
 1448     gretl_model_set_int(pmod, "rq", 1);
 1449     if (opt & OPT_R) {
 1450         gretl_model_set_int(pmod, "rq_nid", 1);
 1451         pmod->opt |= OPT_R;
 1452     }
 1453     gretl_model_set_vcv_info(pmod, VCV_RQ, (opt & OPT_R)?
 1454                  RQ_NID : RQ_ASY);
 1455     }
 1456 
 1457     gretl_matrix_free(y);
 1458     gretl_matrix_free(X);
 1459 
 1460     if (err && pmod->errcode == 0) {
 1461     pmod->errcode = err;
 1462     }
 1463 
 1464     return err;
 1465 }
 1466 
 1467 /* restock the y and X matrices with a bootstrap sample */
 1468 
 1469 static void rq_refill_matrices (MODEL *pmod,
 1470                 DATASET *dset,
 1471                 gretl_matrix *y,
 1472                 gretl_matrix *X,
 1473                 int *sample)
 1474 {
 1475     int n = pmod->nobs;
 1476     int p = pmod->ncoeff;
 1477     int yno = pmod->list[1];
 1478     int i, j, t, v;
 1479 
 1480     for (i=0; i<n; i++) {
 1481     t = sample[i];
 1482     gretl_vector_set(y, i, dset->Z[yno][t]);
 1483     }
 1484 
 1485     for (j=0; j<p; j++) {
 1486     v = pmod->list[j+2];
 1487     for (i=0; i<n; i++) {
 1488         t = sample[i];
 1489         gretl_matrix_set(X, i, j, dset->Z[v][t]);
 1490     }
 1491     }
 1492 }
 1493 
 1494 static int *good_observations_array (MODEL *pmod)
 1495 {
 1496     int *g = malloc(pmod->nobs * sizeof *g);
 1497 
 1498     if (g != NULL) {
 1499     int t, s = 0;
 1500 
 1501     for (t=pmod->t1; t<=pmod->t2; t++) {
 1502         if (!na(pmod->uhat[t])) {
 1503         g[s++] = t;
 1504         }
 1505     }
 1506     }
 1507 
 1508     return g;
 1509 }
 1510 
 1511 #define ITERS 500
 1512 
 1513 /* obtain bootstrap estimates of LAD covariance matrix */
 1514 
 1515 static int lad_bootstrap_vcv (MODEL *pmod, DATASET *dset,
 1516                   gretl_matrix *y, gretl_matrix *X,
 1517                   struct br_info *rq)
 1518 {
 1519     double **coeffs = NULL;
 1520     double *meanb = NULL;
 1521     int *sample = NULL;
 1522     int *goodobs = NULL;
 1523     double xi, xj;
 1524     int i, j, k;
 1525     int nc = pmod->ncoeff;
 1526     int nvcv, n = pmod->nobs;
 1527     int err = 0;
 1528 
 1529     /* note: new_vcv sets all entries to zero */
 1530     err = gretl_model_new_vcv(pmod, &nvcv);
 1531     if (err) {
 1532     return err;
 1533     }
 1534 
 1535     /* an array for each coefficient */
 1536     coeffs = doubles_array_new(nc, ITERS);
 1537 
 1538     /* a scalar for each coefficient mean */
 1539     meanb = malloc(nc * sizeof *meanb);
 1540 
 1541     /* resampling array of length pmod->nobs */
 1542     sample = malloc(n * sizeof *sample);
 1543 
 1544     if (coeffs == NULL || meanb == NULL || sample == NULL) {
 1545     err = E_ALLOC;
 1546     goto bailout;
 1547     }
 1548 
 1549     /* apparatus for handling interior NAs */
 1550     if (model_has_missing_obs(pmod)) {
 1551     goodobs = good_observations_array(pmod);
 1552     if (goodobs == NULL) {
 1553         err = E_ALLOC;
 1554         goto bailout;
 1555     }
 1556     }
 1557 
 1558     for (k=0; k<ITERS && !err; k++) {
 1559     /* create random sample index array */
 1560     for (i=0; i<n; i++) {
 1561         j = gretl_rand_int_max(n);
 1562         if (goodobs != NULL) {
 1563         sample[i] = goodobs[j];
 1564         } else {
 1565         sample[i] = pmod->t1 + j;
 1566         }
 1567     }
 1568 
 1569     rq_refill_matrices(pmod, dset, y, X, sample);
 1570 
 1571     /* re-estimate LAD model */
 1572     err = real_br_calc(y, X, 0.5, rq, 0);
 1573 
 1574     if (!err) {
 1575         for (i=0; i<nc; i++) {
 1576         coeffs[i][k] = rq->coeff[i];
 1577         }
 1578     }
 1579     }
 1580 
 1581     /* find means of coeff estimates */
 1582     for (i=0; i<nc && !err; i++) {
 1583     double bbar = 0.0;
 1584 
 1585     for (k=0; k<ITERS; k++) {
 1586        bbar += coeffs[i][k];
 1587     }
 1588     meanb[i] = bbar / ITERS;
 1589     }
 1590 
 1591     /* find variances and covariances */
 1592     for (i=0; i<nc && !err; i++) {
 1593     double vi = 0.0;
 1594 
 1595     for (k=0; k<ITERS; k++) {
 1596         xi = coeffs[i][k] - meanb[i];
 1597         vi += xi * xi;
 1598         for (j=0; j<=i; j++) {
 1599         xj = coeffs[j][k] - meanb[j];
 1600         pmod->vcv[ijton(i, j, nc)] += xi * xj;
 1601         }
 1602     }
 1603     pmod->sderr[i] = sqrt(vi / ITERS);
 1604     }
 1605 
 1606     if (!err) {
 1607     for (i=0; i<nvcv; i++) {
 1608         pmod->vcv[i] /= ITERS;
 1609     }
 1610     }
 1611 
 1612  bailout:
 1613 
 1614     free(sample);
 1615     free(meanb);
 1616     doubles_array_free(coeffs, nc);
 1617 
 1618     if (goodobs != NULL) {
 1619     free(goodobs);
 1620     }
 1621 
 1622     return err;
 1623 }
 1624 
 1625 static void lad_scrub_vcv (MODEL *pmod)
 1626 {
 1627     int i;
 1628 
 1629     free(pmod->vcv);
 1630     pmod->vcv = NULL;
 1631     free(pmod->xpx);
 1632     pmod->xpx = NULL;
 1633 
 1634     for (i=0; i<pmod->ncoeff; i++) {
 1635     pmod->sderr[i] = NADBL;
 1636     }
 1637 }
 1638 
 1639 static int lad_fit_br (MODEL *pmod, DATASET *dset,
 1640                gretl_matrix *y, gretl_matrix *X,
 1641                gretlopt opt)
 1642 {
 1643     struct br_info rq;
 1644     integer n = y->rows;
 1645     integer p = X->cols;
 1646     int err = 0;
 1647 
 1648     err = br_info_alloc(&rq, n, p, 0.5, 0.0, OPT_L);
 1649 
 1650     if (!err) {
 1651     /* get the actual estimates */
 1652     err = real_br_calc(y, X, 0.5, &rq, 0);
 1653     }
 1654 
 1655     if (!err) {
 1656     rq_transcribe_results(pmod, y, 0.5, rq.coeff,
 1657                   rq.resid, RQ_LAD);
 1658     if (rq.warning) {
 1659         gretl_model_set_int(pmod, "nonunique", 1);
 1660     }
 1661     }
 1662 
 1663     if (!err) {
 1664     if (opt & OPT_N) {
 1665         /* --no-vcv */
 1666         lad_scrub_vcv(pmod);
 1667     } else {
 1668         err = lad_bootstrap_vcv(pmod, dset, y, X, &rq);
 1669     }
 1670     }
 1671 
 1672     br_info_free(&rq);
 1673 
 1674     return err;
 1675 }
 1676 
 1677 /* Support the "lad" command: estimate with tau = 0.5, using
 1678    the Barrodale-Roberts method, and add a bootstrapped
 1679    covariance matrix.
 1680 */
 1681 
 1682 int lad_driver (MODEL *pmod, DATASET *dset, gretlopt opt)
 1683 {
 1684     gretl_matrix *y = NULL;
 1685     gretl_matrix *X = NULL;
 1686     int err;
 1687 
 1688     err = rq_make_matrices(pmod, dset, &y, &X, 0);
 1689 
 1690     if (!err) {
 1691     err = lad_fit_br(pmod, dset, y, X, opt);
 1692     }
 1693 
 1694     if (!err) {
 1695     gretl_model_add_y_median(pmod, dset->Z[pmod->list[1]]);
 1696     pmod->ci = LAD;
 1697     }
 1698 
 1699     gretl_matrix_free(y);
 1700     gretl_matrix_free(X);
 1701 
 1702     if (err && pmod->errcode == 0) {
 1703     pmod->errcode = err;
 1704     }
 1705 
 1706     return err;
 1707 }