Calculate quantiles with bootstrapping in Stata

Question

I have data on income and want to calculate percentiles using bootstrapping as well as population weights. What I did so far (and what worked):

collapse    (p10) p10 = income ///
            (p90) p90 = income ///
            (p95) p95 = income ///
            (p99) p99 = income ///
            [iweight = wg], by(year)

How can I add bootstrapping to this code?

Thanks!

Answer 1

There is a community-contributed command qrprocess that can do fast multiple quantile regression with weights and bootstrapping:

. sysuse auto, clear
(1978 automobile data)

. drop if missing(rep78)
(5 observations deleted)

. set seed 10312023

. ssc install moremata
checking moremata consistency and verifying not already installed...
all files already exist and are up to date.

. ssc install qrprocess
checking qrprocess consistency and verifying not already installed...
all files already exist and are up to date.

. /* M1: use weighted quantile regression + margins get the quantiles */
. qrprocess price i.rep78 [iw=mpg], q(.10 .90 .95 .99) vce(bootstrap, reps(100))
(bootstrapping .x.......x......x.......x................................x.................x.........x.x....x..x..............)

Quantile regression
No. of obs.        69      
Algorithm:         qreg.
Variance:          empirical bootstrap.

------------------------------------------------------------------------------
      price  |      Coef.   Std. Err.      t     P>|t|    [95% Conf. Interval]
-------------+----------------------------------------------------------------
Quant. 0.1   |
    2.rep78  |       -528   460.1462   -1.15    0.255    -1447.248    391.2477
    3.rep78  |       -300   423.3904   -0.71    0.481     -1145.82    545.8195
    4.rep78  |       -200   458.3257   -0.44    0.664    -1115.611    715.6108
    5.rep78  |       -447   464.7086   -0.96    0.340    -1375.362    481.3622
      _cons  |       4195   333.4818   12.58    0.000     3528.794    4861.206
------------------------------------------------------------------------------
Quant. 0.9   |
    2.rep78  |       1408   4199.176    0.34    0.738    -6980.819    9796.819
    3.rep78  |       6563   2324.565    2.82    0.006     1919.148    11206.85
    4.rep78  |       3880    992.606    3.91    0.000     1897.042    5862.958
    5.rep78  |       4756   2737.756    1.74    0.087    -713.2964     10225.3
      _cons  |       4934   333.4818   14.80    0.000     4267.794    5600.206
------------------------------------------------------------------------------
Quant. 0.95  |
    2.rep78  |       9566   4110.508    2.33    0.023     1354.317    17777.68
    3.rep78  |       8660   2004.501    4.32    0.000     4655.548    12664.45
    4.rep78  |       4801   747.2811    6.42    0.000     3308.134    6293.866
    5.rep78  |       7061   2410.464    2.93    0.005     2245.545    11876.46
      _cons  |       4934   333.4818   14.80    0.000     4267.794    5600.206
------------------------------------------------------------------------------
Quant. 0.99  |
    2.rep78  |       9566   4032.689    2.37    0.021     1509.778    17622.22
    3.rep78  |      10972   1289.523    8.51    0.000     8395.882    13548.12
    4.rep78  |       4801   717.3175    6.69    0.000     3367.993    6234.007
    5.rep78  |       7061   2024.685    3.49    0.001     3016.226    11105.77
      _cons  |       4934   333.4818   14.80    0.000     4267.794    5600.206
------------------------------------------------------------------------------

. margins rep78 [iw=mpg], predict(equation(q1)) predict(equation(q2)) predict(equation(q3)) predict(equation(q4)) 
warning: cannot perform check for estimable functions.

Adjusted predictions                                        Number of obs = 69
Model VCE: BOOTSTRAP

1._predict: 0..1 QR fitted values, predict(equation(q1))
2._predict: 0..9 QR fitted values, predict(equation(q2))
3._predict: 0..95 QR fitted values, predict(equation(q3))
4._predict: 0..99 QR fitted values, predict(equation(q4))

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
_predict#rep78 |
          1 1  |       4195   333.4818    12.58   0.000     3541.388    4848.612
          1 2  |       3667   303.8006    12.07   0.000     3071.562    4262.438
          1 3  |       3895   319.6076    12.19   0.000     3268.581    4521.419
          1 4  |       3995   333.1012    11.99   0.000     3342.134    4647.866
          1 5  |       3748   272.9058    13.73   0.000     3213.115    4282.885
          2 1  |       4934   333.4818    14.80   0.000     4280.388    5587.612
          2 2  |       6342   4197.791     1.51   0.131    -1885.519    14569.52
          2 3  |      11497    2298.64     5.00   0.000     6991.747    16002.25
          2 4  |       8814   889.7568     9.91   0.000     7070.109    10557.89
          2 5  |       9690   2660.416     3.64   0.000     4475.681    14904.32
          3 1  |       4934   333.4818    14.80   0.000     4280.388    5587.612
          3 2  |      14500   4100.996     3.54   0.000     6462.196     22537.8
          3 3  |      13594   1975.281     6.88   0.000     9722.521    17465.48
          3 4  |       9735   673.6134    14.45   0.000     8414.742    11055.26
          3 5  |      11995   2362.492     5.08   0.000       7364.6     16625.4
          4 1  |       4934   333.4818    14.80   0.000     4280.388    5587.612
          4 2  |      14500   4021.214     3.61   0.000     6618.566    22381.43
          4 3  |      15906   1240.245    12.82   0.000     13475.16    18336.84
          4 4  |       9735   639.6546    15.22   0.000       8481.3     10988.7
          4 5  |      11995   1983.466     6.05   0.000     8107.478    15882.52
--------------------------------------------------------------------------------

. /* M2: confirm output of margins matches your code */
. collapse ///
>  (p10) p10 = price ///
>  (p90) p90 = price ///
>  (p95) p95 = price ///
>  (p99) p99 = price ///
>  [iweight = mpg], by(rep78)

. list, noobs

  +------------------------------------------+
  | rep78     p10      p90      p95      p99 |
  |------------------------------------------|
  |     1   4,195    4,934    4,934    4,934 |
  |     2   3,667    6,342   14,500   14,500 |
  |     3   3,895   11,497   13,594   15,906 |
  |     4   3,995    8,814    9,735    9,735 |
  |     5   3,748    9,690   11,995   11,995 |
  +------------------------------------------+