Help with bootstrap in obtaining a standard error.

Hi Everyone: I think my problem has nothing to do with the data set and so I'm not showing a data example. In the bootstrap, I want a standard error for the average of six estimated marginal effects but I get a missing value. I get the standard errors for the marginal effects themselves, but not the average. Clearly I don't know how to compute a function of coefficients within a bootstrap program. See my calculation of tauavg. Thanks for any hints.

Code:

. capture program drop aggregate_boot

.         
. program aggregate_boot, rclass
  1. 
. poisson y i.w#c.d4#c.f04 i.w#c.d4#c.f05 i.w#c.d4#c.f06 ///
>         i.w#c.d5#c.f05 i.w#c.d5#c.f06 ///
>         i.w#c.d6#c.f06 ///
>         i.w#c.d4#c.f04#c.x i.w#c.d4#c.f05#c.x i.w#c.d4#c.f06#c.x ///
>         i.w#c.d5#c.f05#c.x i.w#c.d5#c.f06#c.x ///
>         i.w#c.d6#c.f06#c.x ///
>         f02 f03 f04 f05 f06 ///
>         c.f02#c.x c.f03#c.x c.f04#c.x c.f05#c.x c.f06#c.x ///
>         d4 d5 d6 x c.d4#c.x c.d5#c.x c.d6#c.x, noomitted
  2. estimates store beta
  3.         
. margins, dydx(w) at(d4 = 1 d5 = 0 d6 = 0 f02 = 0 f03 = 0 f04 = 1 f05 = 0 f06 = 0) ///
>         subpop(if d4 == 1) noestimcheck post
  4. return scalar tau44 = _b[1.w]
  5. estimates restore beta
  6. margins, dydx(w) at(d4 = 1 d5 = 0 d6 = 0 f02 = 0 f03 = 0 f04 = 0 f05 = 1 f06 = 0) ///
>         subpop(if d4 == 1) noestimcheck post
  7. return scalar tau45 = _b[1.w]
  8. estimates restore beta
  9. margins, dydx(w) at(d4 = 1 d5 = 0 d6 = 0 f02 = 0 f03 = 0 f04 = 0 f05 = 0 f06 = 1) ///
>         subpop(if d4 == 1) noestimcheck post
 10. return scalar tau46 = _b[1.w]
 11. estimates restore beta
 12. margins, dydx(w) at(d4 = 0 d5 = 1 d6 = 0 f02 = 0 f03 = 0 f04 = 0 f05 = 1 f06 = 0) ///
>         subpop(if d5 == 1) noestimcheck post
 13. return scalar tau55 = _b[1.w]
 14. estimates restore beta
 15. margins, dydx(w) at(d4 = 0 d5 = 1 d6 = 0 f02 = 0 f03 = 0 f04 = 0 f05 = 0 f06 = 1) ///
>         subpop(if d5 == 1) noestimcheck post
 16. return scalar tau56 = _b[1.w]
 17. estimates restore beta
 18. margins, dydx(w) at(d4 = 0 d5 = 0 d6 = 1 f02 = 0 f03 = 0 f04 = 0 f05 = 0 f06 = 1) ///
>         subpop(if d6 == 1) noestimcheck post
 19. return scalar tau66 = _b[1.w]
 20. 
. return scalar tauavg = (tau44 + tau45 + tau46 + tau55 + tau56 + tau66)/6
 21. 
. end

. 
. bootstrap r(tau44) r(tau45) r(tau46) r(tau55) r(tau56) r(tau66) r(tauavg),  ///
>         reps(50) seed(123) cluster(id) idcluster(newid): aggregate_boot
(running aggregate_boot on estimation sample)

Bootstrap replications (50)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50

Bootstrap results                                        Number of obs = 6,000
                                                         Replications  =    50

      Command: aggregate_boot
        _bs_1: r(tau44)
        _bs_2: r(tau45)
        _bs_3: r(tau46)
        _bs_4: r(tau55)
        _bs_5: r(tau56)
        _bs_6: r(tau66)
        _bs_7: r(tauavg)

                                  (Replications based on 1,000 clusters in id)
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _bs_1 |   1.017501   1.119012     0.91   0.363    -1.175723    3.210725
       _bs_2 |    6.00713   1.978965     3.04   0.002     2.128431    9.885829
       _bs_3 |   4.569667   1.379358     3.31   0.001     1.866174    7.273159
       _bs_4 |   7.170127   2.957668     2.42   0.015     1.373203    12.96705
       _bs_5 |   7.185492   2.297489     3.13   0.002     2.682496    11.68849
       _bs_6 |   13.73294   12.83413     1.07   0.285    -11.42151    38.88738
       _bs_7 |   6.613809          .        .       .            .           .
------------------------------------------------------------------------------