I am attempting to manually bootstrap standard errors using bsample, simulate and bstat commands. I have tried doing this using the standard bootstrap procedure to no avail. See previous post:

https://www.statalist.org/forums/for...and-panel-data

I am working from two example samples constructed from Stata's 'auto' dataset:

Sample 1: cross-sectional

Code:
* Example generated by -dataex-. For more info, type help dataex
clear
input int mpg float headroom int weight
14   4 4330
14 3.5 3900
21   3 4290
29 2.5 2110
16   4 3690
22 3.5 3180
22   2 3220
24   2 2750
19 3.5 3430
30   2 2120
18   4 3600
16   4 3600
17 4.5 3740
28 1.5 1800
21   2 2650
12 3.5 4840
12 2.5 4720
14 3.5 3830
22   3 2580
14 3.5 4060
15 3.5 3720
18   3 3370
14   3 4130
20 3.5 2830
21   4 4060
19   2 3310
19 4.5 3300
18   4 3690
19 4.5 3370
24   2 2730
16 3.5 4030
28   2 3260
34 2.5 1800
25   4 2200
26 1.5 2520
18   5 3330
18   4 3700
18 1.5 3470
19   2 3210
19 3.5 3200
19 3.5 3420
24   2 2690
17   3 2830
23 2.5 2070
25 2.5 2650
23 1.5 2370
35   2 2020
24 2.5 2280
21 2.5 2750
21 2.5 2130
25   3 2240
28 2.5 1760
30 3.5 1980
14 3.5 3420
26   3 1830
35 2.5 2050
18 2.5 2410
31   3 2200
18   2 2670
23 2.5 2160
41   3 2040
25   3 1930
25   2 1990
17 2.5 3170
end
Sample 2: panel

Code:
* Example generated by -dataex-. For more info, type help dataex
clear
input float(price mpg headroom t id)
     4099        22       2.5 1  1
 4186.009 26.934307  1.933958 2  1
     4749        17         3 1  2
 4843.107  18.91313 3.3818026 2  2
     3799        22         3 1  3
 3770.704  24.97907 2.9321864 2  3
     4816        20       4.5 1  4
4852.0933 22.037195 4.3978457 2  4
     7827        15         4 1  5
  7734.57 19.414215  3.253971 2  5
     5788        18         4 1  6
 5762.868 20.123154  4.837951 2  6
     4453        26         3 1  7
4521.1724   29.6178 2.0983803 2  7
     5189        20         2 1  8
 5146.077  23.62131 1.3961287 2  8
    10372        16       3.5 1  9
 10323.16 15.383166 2.9035466 2  9
     4082        19       3.5 1 10
 4041.077 16.604582  2.918597 2 10
end
The procedure I am using is as follows:

Code:
use sample1.dta, replace
qui reg weight mpg headroom

use sample2.dta, replace
capture drop weight_hat
predict weight_hat

xtset id t
xtreg price weight_hat, fe

matrix b = e(b)
scalar obs = r(N)


capture program drop example
program define example, eclass

    use sample1.dta, replace
    bsample
    qui reg weight mpg headroom
    
    use sample2.dta, replace
    capture drop weight_hat
    predict weight_hat
    
    xtset id t
    xtreg price weight_hat, fe
    
    exit
end

simulate, reps(10000) seed(12345): example
bstat, stat(b) n(`obs')
estat bootstrap, all
The results look promising in that the code runs and the coefficients are correct, but the bootstrapped standard errors are way off (too small).

Without bootstrap:

Code:
Fixed-effects (within) regression               Number of obs     =         20
Group variable: id                              Number of groups  =         10

R-squared:                                      Obs per group:
     Within  = 0.1927                                         min =          2
     Between = 0.3717                                         avg =        2.0
     Overall = 0.3019                                         max =          2

                                                F(1,9)            =       2.15
corr(u_i, Xb) = -0.5577                         Prob > F          =     0.1768

------------------------------------------------------------------------------
       price | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
  weight_hat |  -.0608991   .0415549    -1.47   0.177    -.1549027    .0331046
       _cons |   5712.816   139.0787    41.08   0.000     5398.199    6027.434
-------------+----------------------------------------------------------------
     sigma_u |  2081.0277
     sigma_e |   30.10469
         rho |  .99979077   (fraction of variance due to u_i)
------------------------------------------------------------------------------
F test that all u_i=0: F(9, 9) = 6584.95                     Prob > F = 0.0000
With bootstrap:

Code:
Bootstrap results                                           Number of obs = 20
                                                            Replications  = 50

-------------------------------------------------------------------------------
              |   Observed   Bootstrap                         Normal-based
              | coefficient  std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
_b_weight_hat |  -.0608991   .0107155    -5.68   0.000    -.0819011    -.039897
      _b_cons |   5712.816   34.90237   163.68   0.000     5644.409    5781.224
-------------------------------------------------------------------------------

. eststo: estat bootstrap, all

Bootstrap results                               Number of obs     =         20
                                                Replications      =         50

------------------------------------------------------------------------------
             |    Observed               Bootstrap
             | coefficient       Bias    std. err.  [95% conf. interval]
-------------+----------------------------------------------------------------
_b_weight_~t |  -.06089906  -.0004936   .01071552   -.0819011   -.039897   (N)
             |                                      -.0800628   -.040871   (P)
             |                                      -.0800628  -.0400041  (BC)
     _b_cons |   5712.8164   1.755529    34.90237    5644.409   5781.224   (N)
             |                                       5646.295   5768.709   (P)
             |                                       5644.807   5768.709  (BC)
------------------------------------------------------------------------------
Key:  N: Normal
      P: Percentile
     BC: Bias-corrected
Does anyone have any ideas what I am doing wrong?

Thank you!