Bootstrap - Standard errors (almost) equal to zero

[Adam Dvir]

Dear Statalisters,

I am trying to run a two-stage procedure of Regression Discontinuity in Time using a balanced household-level panel data, based on the method of Hausman & Rapson 2016.
In the first stage, I regress the purchased quantity on seasonal covariates using FE Poisson regression (because of the nonnegative nature of my dependent variable).
Then I predict the residuals and run the second stage regression, which is a local linear RD where the dependent variable is the residuals from the first stage and the running variable is time. According to theory, standard errors should be retrieved by a bootstrapping procedure that allows first stage variance to be reflected in the second stage.
For this purpose, I wrote a simple program:

However, when I use bootstrap to run the program I get almost zero standard errors:

After an intensive exploration of a solution, I found an old version of bootstrap that includes the option -noesample-.
When I run the program using the old version, including the -noesample- option, I get reasonable standard errors,
which are similar to the standard errors I get when I run the procedure in two separate stages without bootstrapping:

The weird thing is when I remove the -noesample- option from the syntax of the old bootstrap version,
I get the exact same standard errors as in the "new" bootstrap version::

Note that -noesample- is no longer a valid option in the "new" bootstrap syntax.
see this post which is somewhat related:
https://www.statalist.org/forums/for...help_bootstrap


Question: Could someone please help me figure out what's wrong with my code that yields these tiny standard errors, and how can it be fixed?

Many many thanks,
Adam

[Leonardo Guizzetti]

I can't help you with your main question, but I noticed that you are not using -bootstrap- in the correct way to resample on clusters. For details, see Clyde's post in #9 of the thread you have linked.

[Adam Dvir]

Dear Leonardo Guizzetti ,

Thank you for replying.
I followed Clyde's suggestion, but the results did not change.
Please see below,

Adam

[Leonardo Guizzetti]

Again, you need to change all references to panelist_id in your estimation command, so your -vce()- specification is still incorrect.

Like I said, I can't help with your actual question, but please do have a read of the FAQ (linked at the top of the page). Screenshots are not as helpful as you might think. You would be better served leaning how to copy the the code and results using the [CODE] tags and preferably a reproducible data example using -dataex- (available in fully updated versions of Stata 14.2, 15 and 16). All of this is explained in the FAQ.

[Adam Dvir]

Dear Leonardo,

Sorry about the use of screenshots. I now better understand why they are less informative.

Here is a sample of my data:

* Example generated by -dataex-. To install: ssc install dataex
clear
input double(panelist_id month_id weight) float gweight double(advertising_cost price_ix) float(time month_int after d_indep_day kernel_weight band)
102 201411               .684  684 4037819.6100000003                  . -12 11 0 0     -.9375 .
102 201412               .554  554           874339.2                  . -11 12 0 0  -.6679688 .
102 201501              3.034 3034         2560573.57                100 -10  1 0 0   -.421875 .
102 201502               .882  882         3375242.41 100.69500000000001  -9  2 0 0 -.19921875 .
102 201503               1.59 1590          862094.97 100.58399999999999  -8  3 0 0          0 1
102 201504 .20400000000000001  204         1457067.48             99.648  -7  4 0 1  .17578125 1
102 201505               .842  842 256811.37000000002            102.953  -6  5 0 0    .328125 1
102 201506               .305  305 2404866.7399999998            100.914  -5  6 0 0  .45703125 1
102 201507  .8379999999999999  838           482087.9            101.969  -4  7 0 0      .5625 1
102 201508              1.232 1232          189708.32            101.352  -3  8 0 0   .6445313 1
102 201509              2.434 2434          419498.38            101.997  -2  9 0 0    .703125 1
102 201510               .374  374            5907.82            103.048  -1 10 0 0   .7382813 1
102 201511 .39799999999999996  398  530074.5599999999 102.74000000000001   0 11 1 0   .7382813 1
102 201512              1.092 1092         1081485.96            100.299   1 12 1 0    .703125 1
102 201601 .21200000000000002  212         1511828.59            100.119   2  1 1 0   .6445313 1
102 201602                  0    0         1696591.45             100.78   3  2 1 0      .5625 1
102 201603               .362  362         1333454.33 102.15799999999999   4  3 1 0  .45703125 1
102 201604 .44799999999999995  448         1872099.81             102.34   5  4 1 0    .328125 1
102 201605               .712  712 1329448.3199999998 101.26200000000001   6  5 1 1  .17578125 1
102 201606               .388  388 230812.27999999997            103.402   7  6 1 0          0 1
102 201607 .22000000000000003  220         2627139.71 103.04299999999998   8  7 1 0 -.19921875 .
102 201608                 .4  400  946631.7300000001 102.89200000000001   9  8 1 0   -.421875 .
102 201609                  0    0          233178.01 101.81700000000001  10  9 1 0  -.6679688 .
102 201610               .434  434          284624.85 102.66600000000001  11 10 1 0     -.9375 .
102 201611                .26  260         1217520.99 102.47200000000001  12 11 1 0 -1.2304688 .
102 201612                  0    0          685104.74 102.20700000000001  13 12 1 0  -1.546875 .
102 201701                 .4  400             365433 101.96799999999999  14  1 1 0 -1.8867188 .
102 201702                  0    0         1431956.75 101.48299999999999  15  2 1 0      -2.25 .
102 201703                .16  160          389932.85            101.174  16  3 1 0  -2.636719 .
102 201704               .724  724          177224.15  99.96600000000001  17  4 1 1  -3.046875 .
102 201705                  0    0          161195.05            100.679  18  5 1 0  -3.480469 .
102 201706                  0    0           37077.04 101.03200000000001  19  6 1 0    -3.9375 .
102 201707                  0    0         2280233.05 101.28399999999999  20  7 1 0 -4.4179688 .
102 201708  .5279999999999999  528          388680.26            101.265  21  8 1 0  -4.921875 .
102 201709                  0    0         2638988.91 101.45200000000001  22  9 1 0  -5.449219 .
102 201710                  0    0          125068.09            102.174  23 10 1 0         -6 .
102 201711                  0    0 2373066.4699999997            102.424  24 11 1 0  -6.574219 .
105 201411                  0    0 4037819.6100000003                  . -12 11 0 0     -.9375 .
105 201412                 .7  700           874339.2                  . -11 12 0 0  -.6679688 .
105 201501                  3 3000         2560573.57                100 -10  1 0 0   -.421875 .
105 201502                  0    0         3375242.41 100.69500000000001  -9  2 0 0 -.19921875 .
105 201503                  1 1000          862094.97 100.58399999999999  -8  3 0 0          0 1
105 201504                  2 2000         1457067.48             99.648  -7  4 0 1  .17578125 1
105 201505                 .7  700 256811.37000000002            102.953  -6  5 0 0    .328125 1
105 201506                 .7  700 2404866.7399999998            100.914  -5  6 0 0  .45703125 1
105 201507                 .6  600           482087.9            101.969  -4  7 0 0      .5625 1
105 201508                1.1 1100          189708.32            101.352  -3  8 0 0   .6445313 1
105 201509                 .7  700          419498.38            101.997  -2  9 0 0    .703125 1
105 201510                  0    0            5907.82            103.048  -1 10 0 0   .7382813 1
105 201511               .381  381  530074.5599999999 102.74000000000001   0 11 1 0   .7382813 1
105 201512                 .7  700         1081485.96            100.299   1 12 1 0    .703125 1
105 201601                  0    0         1511828.59            100.119   2  1 1 0   .6445313 1
105 201602                2.8 2800         1696591.45             100.78   3  2 1 0      .5625 1
105 201603                  0    0         1333454.33 102.15799999999999   4  3 1 0  .45703125 1
105 201604                 .6  600         1872099.81             102.34   5  4 1 0    .328125 1
105 201605                1.9 1900 1329448.3199999998 101.26200000000001   6  5 1 1  .17578125 1
105 201606                  1 1000 230812.27999999997            103.402   7  6 1 0          0 1
105 201607                 .6  600         2627139.71 103.04299999999998   8  7 1 0 -.19921875 .
105 201608                 .8  800  946631.7300000001 102.89200000000001   9  8 1 0   -.421875 .
105 201609                  1 1000          233178.01 101.81700000000001  10  9 1 0  -.6679688 .
105 201610                 .8  800          284624.85 102.66600000000001  11 10 1 0     -.9375 .
105 201611                1.5 1500         1217520.99 102.47200000000001  12 11 1 0 -1.2304688 .
105 201612                 .8  800          685104.74 102.20700000000001  13 12 1 0  -1.546875 .
105 201701                 .6  600             365433 101.96799999999999  14  1 1 0 -1.8867188 .
105 201702                2.3 2300         1431956.75 101.48299999999999  15  2 1 0      -2.25 .
105 201703                 .8  800          389932.85            101.174  16  3 1 0  -2.636719 .
105 201704                  3 3000          177224.15  99.96600000000001  17  4 1 1  -3.046875 .
105 201705                 .9  900          161195.05            100.679  18  5 1 0  -3.480469 .
105 201706                1.3 1300           37077.04 101.03200000000001  19  6 1 0    -3.9375 .
105 201707                 .8  800         2280233.05 101.28399999999999  20  7 1 0 -4.4179688 .
105 201708                1.7 1700          388680.26            101.265  21  8 1 0  -4.921875 .
105 201709                1.6 1600         2638988.91 101.45200000000001  22  9 1 0  -5.449219 .
105 201710                  0    0          125068.09            102.174  23 10 1 0         -6 .
105 201711                  0    0 2373066.4699999997            102.424  24 11 1 0  -6.574219 .
109 201411                2.2 2200 4037819.6100000003                  . -12 11 0 0     -.9375 .
109 201412                 .8  800           874339.2                  . -11 12 0 0  -.6679688 .
109 201501                2.2 2200         2560573.57                100 -10  1 0 0   -.421875 .
109 201502               3.45 3450         3375242.41 100.69500000000001  -9  2 0 0 -.19921875 .
109 201503                6.1 6100          862094.97 100.58399999999999  -8  3 0 0          0 1
109 201504                1.9 1900         1457067.48             99.648  -7  4 0 1  .17578125 1
109 201505                2.4 2400 256811.37000000002            102.953  -6  5 0 0    .328125 1
109 201506                1.8 1800 2404866.7399999998            100.914  -5  6 0 0  .45703125 1
109 201507                4.8 4800           482087.9            101.969  -4  7 0 0      .5625 1
109 201508                1.4 1400          189708.32            101.352  -3  8 0 0   .6445313 1
109 201509                1.7 1700          419498.38            101.997  -2  9 0 0    .703125 1
109 201510                1.4 1400            5907.82            103.048  -1 10 0 0   .7382813 1
109 201511                2.1 2100  530074.5599999999 102.74000000000001   0 11 1 0   .7382813 1
109 201512                2.9 2900         1081485.96            100.299   1 12 1 0    .703125 1
109 201601               2.88 2880         1511828.59            100.119   2  1 1 0   .6445313 1
109 201602                 .1  100         1696591.45             100.78   3  2 1 0      .5625 1
109 201603               1.75 1750         1333454.33 102.15799999999999   4  3 1 0  .45703125 1
109 201604               1.05 1050         1872099.81             102.34   5  4 1 0    .328125 1
109 201605               1.45 1450 1329448.3199999998 101.26200000000001   6  5 1 1  .17578125 1
109 201606                .15  150 230812.27999999997            103.402   7  6 1 0          0 1
109 201607                1.5 1500         2627139.71 103.04299999999998   8  7 1 0 -.19921875 .
109 201608               1.25 1250  946631.7300000001 102.89200000000001   9  8 1 0   -.421875 .
109 201609                2.8 2800          233178.01 101.81700000000001  10  9 1 0  -.6679688 .
109 201610                 .8  800          284624.85 102.66600000000001  11 10 1 0     -.9375 .
109 201611                 .4  400         1217520.99 102.47200000000001  12 11 1 0 -1.2304688 .
109 201612                  1 1000          685104.74 102.20700000000001  13 12 1 0  -1.546875 .
end
label values month_int month_int
label def month_int 1 "01", modify
label def month_int 2 "02", modify
label def month_int 3 "03", modify
label def month_int 4 "04", modify
label def month_int 5 "05", modify
label def month_int 6 "06", modify
label def month_int 7 "07", modify
label def month_int 8 "08", modify
label def month_int 9 "09", modify
label def month_int 10 "10", modify
label def month_int 11 "11", modify
label def month_int 12 "12", modify

In addition, here is the program code (I replaced the macro with the real variables):

// Two stg estimation program - for bootstrapped std
capture program drop twostg_est

program twostg_est, rclass

xtset panelist_id
capture drop res fitted yhat meany meanyhat exp_alpha

// first stage regression
xtpoisson weight price_ix advertising_cost i.month_int d_indep_day, fe iterate(25) vce(robust)


//gen residuals after ppml regression
    predict double fitted, xb
    gen double yhat=exp(fitted)
    egen meany=mean(weight), by(panelist_id)
    egen meanyhat=mean(yhat), by(panelist_id)
    gen double exp_alpha=meany/meanyhat if meanyhat>0
    replace yhat=yhat*exp_alpha if meanyhat>0
    gen res = (weight - yhat)*1000

// second stage regression
reghdfe res i.after##c.time [aw=kernel_weight] if band == 1 ///
, absorb(panelist_id) vce(cluster panelist_id)

end

and the Bootstrap command, including Clyde's remark:

xtset, clear
bootstrap _b[1.after], reps(50) seed(1710) cluster(panelist_id) idcluster(newid): /// 
twostg_est $vars_first_stg $vars, vce(cluster panelist_id) i(newid) t(month_id)

I am not sure I understand what you meant by:

Again, you need to change all references to panelist_id in your estimation command, so your -vce()- specification is still incorrect.

since I use vce(cluster panelist_id).

Thanks again,
Adam

[Leonardo Guizzetti]

When you resample on cluster, what -idcluster()- does is create a new, temporary variable which defines the newly resampled cluster. This means that any instance of the cluster variable that follows after the colon must use the -idcluster()- variable. So the two instances below in red need to be changed.

 
 bootstrap _b[1.after], reps(50) seed(1710) cluster(panelist_id) idcluster(newid): ///  twostg_est $vars_first_stg $vars, vce(cluster newid) i(newid) t(month_id)

That said, your program cannot possible be written correctly, because -twostg_est- is not aware of the options you have tried to pass to it. Your program will still run as written, but it's not doing what you think it is. Also, within that program the cluster id will be incorrectly specified as above.

[Adam Dvir]

Dear Leonardo,

I tried to change the -idcluster()- as you suggested and also tried no cluster at all. The results remained unchanged.

I've got a feeling that this issue is somewhat related to the fact that the residuals from the first stage are the dependent variable in the second stage (and not an explanatory variable as in many other applications - like 2sls and control function).
To test this hypothesis I replaced the dependent variable of the second stage regression (to be gweight) so the residuals are now just another regressor.
Although this regression is meaningless, standard errors are now far from zero.

// Two stg estimation program - for bootstrapped std
capture program drop twostg_est

program twostg_est, rclass

xtset panelist_id
capture drop res fitted yhat meany meanyhat exp_alpha

// first stage regression
xtpoisson weight $vars_first_stg, fe iterate($iter) //vce(robust)

//gen residuals after ppml regression
    predict double fitted, xb
    gen double yhat=exp(fitted)
    egen meany=mean(weight), by(panelist_id)
    egen meanyhat=mean(yhat), by(panelist_id)
    gen double exp_alpha=meany/meanyhat if meanyhat>0
    replace yhat=yhat*exp_alpha if meanyhat>0
    gen res = (weight - yhat)*1000

// second stage regression
reghdfe gweight res $vars [aw=kernel_weight] if band == 1 ///
, absorb(panelist_id) //vce(cluster new_id)

end

The bootstrap code:

bootstrap _b[1.after], reps(50) seed(1710) cluster(panelist_id) idcluster(newid): ///  
twostg_est $vars, i(newid) t(month_id)

The results:

Any suggestion will be highly appreciated,
Adam

[Adam Dvir]

I found a solution and just wanted to post it here for a recap, as it might help someone in the future.

The source of the problem was in the program.
Thank you Leonardo Guizzetti for the hint.

I added the lines in red and now it works like magic!

capture program drop twostg_est

program twostg_est, eclass
tempname b V 

xtset panelist_id
capture drop res fitted yhat meany meanyhat exp_alpha

// first stage regression
xtpoisson weight $vars_first_stg, fe iterate($iter) vce(robust)


//gen residuals after ppml regression
    predict double fitted, xb
    gen double yhat=exp(fitted)
    egen meany=mean(weight), by(panelist_id)
    egen meanyhat=mean(yhat), by(panelist_id)
    gen double exp_alpha=meany/meanyhat if meanyhat>0
    replace yhat=yhat*exp_alpha if meanyhat>0
    gen res = (weight - yhat)*1000

// second stage regression
reghdfe res $vars [aw=kernel_weight] if band == 1 ///
, absorb(panelist_id) vce(cluster panelist_id)

matrix `b' = e(b)
ereturn post `b' 
end

Best,
Adam

[Leonardo Guizzetti]

I cannot help you further but to point out a few things. First, I don't believe your program is correctly estimating what you want. I tried to say that, but I probably wasn't clear enough.

I think your program and bootstrap call should look something more like this (emphasis on red text):

capture program drop twostg_est
program twostg_est, eclass
  syntax, clustid(varname)    // this lets your specify the name of your cluster variable (e.g., newid)

  xtset `clustid'
  capture drop res fitted yhat meany meanyhat exp_alpha

  // first stage regression
  xtpoisson weight $vars_first_stg, fe vce(robust)

  //gen residuals after ppml regression
  predict double fitted, xb
  gen double yhat=exp(fitted)
  egen meany=mean(weight), by(`clustid')
  egen meanyhat=mean(yhat), by(`clustid')
  gen double exp_alpha=meany/meanyhat if meanyhat>0
  replace yhat=yhat*exp_alpha if meanyhat>0
  gen res = (weight - yhat)*1000

  // second stage regression
  reghdfe res $vars [aw=kernel_weight] if band == 1, absorb(`clustid') vce(cluster `clustid')

  ereturn post e(b)   // this collapses two lines into one.
end

bootstrap _b[1.after], reps(50) seed(1710) cluster(panelist_id) idcluster(newid): twostg_est $vars_first_stg $vars, clustid(newid)

You can read more about this by going to the help for bootstrap. Enter -help bootstrap- and then click the link for the PDF documentation. Skip down to Example 7 and 8 which discuss clustering, and pay particular attention to the technical notes at the end of each example.

In your earlier code, you had specified 3 options to your program (vce(), i() and t() ), but your program has no knowledge that it should be expecting them. Therefore, Stata ignores them. if you want to add those options in, you must code them and then make use of them, like I have done with the clustid() option.

[Leonardo Guizzetti]

To illustrate the above, consider this toy illustration using the pig birth-weight dataset.

clear *
cls

webuse pig

cap program drop mymodel1
program define mymodel1, rclass
  version 16.1
 
  xtset id  // omit time var to demonstrate cluster resampling error
  xtgee weight
  return scalar weight = _b[_cons]
end

mymodel1
bootstrap weight=r(weight) , noisily cluster(id) idcluster(newid) reps(5) seed(17): mymodel1

The (selected) output from bootstrap is:

// use noisily to see the output from each call to my program.
. bootstrap weight=r(weight) , noisily cluster(id) idcluster(newid) reps(5) seed(17): mymodel1
bootstrap: First call to mymodel1 with data as is:

. mymodel1
       panel variable:  id (balanced)

Iteration 1: tolerance = 2.350e-15

GEE population-averaged model                   Number of obs     =        432
Group variable:                         id      Number of groups  =         48
Link:                             identity      Obs per group:
Family:                           Gaussian                    min =          9
Correlation:                  exchangeable                    avg =        9.0
                                                              max =          9
                                                Wald chi2(0)      =          .
Scale parameter:                  276.2861      Prob > chi2       =          .

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   50.40509   .5646629    89.27   0.000     49.29837    51.51181
------------------------------------------------------------------------------

Bootstrap replications (5)
. mymodel1
       panel variable:  id (unbalanced)

GEE population-averaged model                   Number of obs     =        432
Group variable:                         id      Number of groups  =         30
Link:                             identity      Obs per group:
Family:                           Gaussian                    min =          9
Correlation:                  exchangeable                    avg =       14.4
                                                              max =         27
                                                Wald chi2(0)      =          .
Scale parameter:                  259.4306      Prob > chi2       =          .

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   49.35638   .6131599    80.50   0.000     48.15461    50.55816
------------------------------------------------------------------------------
. mymodel1
       panel variable:  id (unbalanced)

GEE population-averaged model                   Number of obs     =        432
Group variable:                         id      Number of groups  =         31
Link:                             identity      Obs per group:
Family:                           Gaussian                    min =          9
Correlation:                  exchangeable                    avg =       13.9
                                                              max =         27
                                                Wald chi2(0)      =          .
Scale parameter:                  278.6297      Prob > chi2       =          .

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   51.27121   .6103018    84.01   0.000     50.07504    52.46738
------------------------------------------------------------------------------
. mymodel1
       panel variable:  id (unbalanced)

GEE population-averaged model                   Number of obs     =        432
Group variable:                         id      Number of groups  =         31
Link:                             identity      Obs per group:
Family:                           Gaussian                    min =          9
Correlation:                  exchangeable                    avg =       13.9
                                                              max =         36
                                                Wald chi2(0)      =          .
Scale parameter:                  276.2892      Prob > chi2       =          .

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   50.44557    .886703    56.89   0.000     48.70767    52.18348
------------------------------------------------------------------------------
. mymodel1
       panel variable:  id (unbalanced)

GEE population-averaged model                   Number of obs     =        432
Group variable:                         id      Number of groups  =         30
Link:                             identity      Obs per group:
Family:                           Gaussian                    min =          9
Correlation:                  exchangeable                    avg =       14.4
                                                              max =         36
                                                Wald chi2(0)      =          .
Scale parameter:                  272.7414      Prob > chi2       =          .

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   49.69621   .7738632    64.22   0.000     48.17946    51.21295
------------------------------------------------------------------------------
. mymodel1
       panel variable:  id (unbalanced)

GEE population-averaged model                   Number of obs     =        432
Group variable:                         id      Number of groups  =         25
Link:                             identity      Obs per group:
Family:                           Gaussian                    min =          9
Correlation:                  exchangeable                    avg =       17.3
                                                              max =         36
                                                Wald chi2(0)      =          .
Scale parameter:                  255.2569      Prob > chi2       =          .

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |    50.0423   .5496361    91.05   0.000     48.96503    51.11956
------------------------------------------------------------------------------

Every call after the first one now has clusters that don't make sense. The data contain weight data for 48 pigs, so I expect to have resampled data from exactly 48 pigs (clusters). Above, we get a range from 25 to 31 pigs. (In fact, I had to deliberately -xtset- the data incorrectly by omitting the time variable to get this to run to demonstrate the point. If I had included it, bootstrap would have immediately failed and issues a warning about repeated time values within panels.)

Compare this to the following program, which passes the new cluster id variable through syntax so that we may properly resample on clusters. Here I have also added in the time variable when i xtset the data.

cap program drop mymodel2
program define mymodel2, rclass
  version 16.1
  syntax , clustid(varname)
 
  xtset `clustid' week
  xtgee weight
  return scalar weight = _b[_cons]
end

mymodel2, clustid(id)
bootstrap weight=r(weight) , noisily cluster(id) group(id) idcluster(newid) reps(5) seed(17): mymodel2, clustid(newid)

Every call to mymodel2 now has exactly 48 cluster with exactly 9 observations per cluster, as expected. I've omitted the output for brevity, but you can easily check this yourself.