如何在stata中用bootstrap方法比较组间差异？请教高手！

下载了的同学，搞明白的，用系统自带的数据演示一下，发上来大家互相学习哦，不能只下载不分享心得啊，呵呵！

这个命令目前只适用于估计固定效应模型，因此只能用在面板数据中。
一般化的方法很容易做，无论是 Logit，probit 还是其他模型，都可以通过简单编程来实现。
下面是Stata高级视频教程 B9_BS_MC.do 一讲中的例子，贴出来，希望能有助于诸位了解这个方法。
等有时间了，我编写一个一般化的命令，来实现这一方法。

*-- 组间系数差异显著性检验
*-------------------------

  *- 问题的来源:
  *
  *  模型：  y_i = x_i*b1 + e_i  (group1: i=1,2,...,N1)
  *          y_i = x_i*b2 + e_i  (group2: i=1,2,...,N2)
  *  
  *  Ho: b1 = b2  如何检验？
  *
  * 传统方法：Chow检验，F 检验，
  *           都需要对样本进行混合，生成虚拟变量；
  *           需要较强的假设条件。
  
  * Bootstrap 法：
  *    在 Ho 成立的情况下，
  *    将 N1 和 N2 个样本混合后得到的估计系数 b 与 b1(=b2)
  *    应该不存在显著差异。
  *   
  * 步骤：参见 Efron(1993,Chp16,p.221)
  *
  * 1. 混合两个子样本组的观察值，得到一个样本数为 N1+N2 的“混合样本”；
  * 2. 从“混合样本”中可重复地抽取(N1+N2)个观察值，
  *    将前   N1 个观察值定义为 group1,
  *    剩下的 N2 个观察值定义为 group2；
  * 3. 分别针对两个样本组进行估计，得到
  *    Diff(bs_j) = b1 - b2, (j=1,2,...,1000)
  * 4. 计算“实证P值”：
  *      P_bs = #{Diff(bs_j) >= Diff(0)} / 1000
  *      即，BS得到的系数差异大于真实系数差异的次数占抽样次数的比例。
  *    其中，Diff(0)=b1-b2，即两组系数的真实差异。
  
  *- 例：
  *  模型： price = a0 + a1*weight + a2*length + a3*mpg + e
  *  问题：在国产车和进口车两个子样本中，汽车重量(weight)对价格影响相同吗？
  
     * 观察到的实际差异：
       sysuse auto, clear
       reg price weight length mpg if foreign==1  /*进口车*/
         local b1 = _b[weight]   
       reg price weight length mpg if foreign==0  /*国产车*/
         local b2 = _b[weight]
       scalar diff0 = `b1' - `b2'
       dis "diff0_weight = " scalar(diff0)
     
     * Bootstrap 检验 Ho: b1 = b2
     *--------------------------------------------------
       local reps = 1000
       mat D = J(`reps', 3, .)   /*存储结果的矩阵*/
       forvalues j = 1/`reps'{
          qui sysuse auto, clear
          bsample
          qui reg price wei len mpg in 1/22  /*前22个观察值，视为进口车*/
          local b1 = _b[weight]
          qui reg price wei len mpg in 23/74 /*后52个观察值，视为国产车*/
          local b2 = _b[weight]
          local diff = `b1' -`b2'
          mat D[`j',1] = (`b1', `b2', `diff')
       }   
     *--------------------------------------------------  
       svmat D, names(d)
       rename d1 b1
       rename d2 b2
       rename d3 diff
      
       * 计算“实证P值”
         count if diff > diff0 | diff == diff0
         local p = `r(N)'/`reps'
         dis "实证P值 = " `p'   
         
       * 图示
           local diff0 = scalar(diff0)
         histogram diff, xline(`diff0',lw(thick))
           local diff0 = scalar(diff0)
         kdensity diff, lw(thick) xline(`diff0',lw(thick))         
      
       * Boostrap 参数估计值 b1 和 b2 的特征
         sum b1 b2
         * 实际估计值
         sysuse auto, clear
         reg price weight length mpg if foreign==1,noheader  /*进口车*/
         reg price weight length mpg if foreign==0,noheader  /*国产车*/        
      
  
  *------------------------------------------------
  *- 扩展 I：如何计算模型中所有系数的“实证P值”？
  *------------------------------------------------
  * 解决思路：将系数差异存放于矩阵中，而非暂元中（如上例）  
     sysuse auto, clear
     reg price wei len mpg
     eret list
     mat list e(b)
     
     * 记录真实系数差异
       sysuse auto, clear
       qui reg price weight length mpg if foreign==1  /*进口车*/
         mat b1 = e(b)   
       qui reg price weight length mpg if foreign==0  /*国产车*/
         mat b2 = e(b)
       mat D0 = b1 - b2
       mat list D0, title(系数估计值的真实差异)
  
     * 采用Bootstrap得到实证差异和实证P值
     *--------------------------------------------------------------------
       local reps = 1000
       mat D = J(`reps', 4, .)   /*存储结果的矩阵*/
       forvalues j = 1/`reps'{
          qui sysuse auto, clear
          bsample
          qui reg price wei len mpg in 1/22  /*前22个观察值，视为进口车*/
          matrix b1 = e(b)
          qui reg price wei len mpg in 23/74 /*后52个观察值，视为国产车*/
          matrix b2 = e(b)
          matrix diff = b1 - b2
          mat D[`j',1] = diff
       }  
     *--------------------------------------------------------------------
      
       svmat D, names(diff)  /*diff1-diff4 对应 weight,length,mpg和常数项*/
      
       mat P = J(4,2,.)  /*记录系数真实差异和实证P值的矩阵*/
       forvalues j = 1/4{
          local diff0_`j' = D0[1,`j']
          qui count if (diff`j'>`diff0_`j'' | diff`j'==`diff0_`j'')
          local p = `r(N)'/1000
          mat P[`j',1] = (`diff0_`j'' , `p')
       }
      
       mat colnames P = 系数真实差异  实证P值
       mat rownames P = weight length mpg cons
       mat list P
      
     * 结论：各变量的组间系数估计值并不存在显著差异。
     
     * 图示：mpg变量的BS系数差异的分布
       histogram diff3, xline(-161.1735, lw(thick))
     
     
  *------------------------------------------------
  *- 扩展 II：如何得到面板数据模型的“实证P值”？
  *------------------------------------------------  
  * 面临的主要问题：如何保证面板中每个截面内部的时序特征？
  * 解决思路：按个体进行抽样(每个面板中的截面都是一个个体)
  * see [XT] p.215 Technical Note
  
  *- 例：成长机会(Tobin)对资本结构的影响
  *      在大规模公司和小规模公司中是否相同？
  *      Ho: b1 = b2  (b_large = b_small)
  
        use xtcs.dta, clear  /*中国上市公司资本结构数据*/
        tsset code year
        egen size_mean = mean(size), by(code) /*保证每家公司都分在同一组*/
        list code year size size_mean in 1/30
        
        sort size_mean
        gen group = group(3)
        tab group
        drop if group == 2   /*去掉中间组，保证最大和最小两个组的差异足够显著*/
        xtreg tl size ndts tang npr tobin if group==1,fe /*小规模公司*/
          est store small
        xtreg tl size ndts tang npr tobin if group==3,fe /*大规模公司*/      
          est store large
        esttab small large, mtitle(small large)
        * 可见，在不同规模的公司中，Tobin对TL的影响是不同的。
        * 下面采用BS进行检验。
        
     * 记录真实系数差异
        xtreg tl size ndts tang npr tobin if group==1,fe /*小规模公司*/
          mat b1 = e(b)
        xtreg tl size ndts tang npr tobin if group==3,fe /*大规模公司*/      
          mat b3 = e(b)
        mat D0 = b1 - b3
        mat list D0, title(系数估计值的真实差异)   
     
     
     * 采用Bootstrap得到实证差异和实证P值
      
       * 解析抽样过程：
         use xtcs.dta, clear
         set seed 12345
         bsample, cluster(code) idcluster(code_bs)
         tsset code_bs year
         list code_bs code year in 1/35 if year>2001
         * 可见：公司000002, 000012 都被抽中了2次
         * 新的截面变量 code_bs 可以唯一标示公司
     
     *--------------------------------------------------------------------
       local reps = 100
       local varlist "tl size ndts tang npr tobin"
       mat D = J(`reps', 6, .)   /*存储结果的矩阵*/
       forvalues j = 1/`reps'{
          use xtcs.dta, clear
          bsample, cluster(code) idcluster(code_bs)
          gen gg = group(3)      /*将样本等分为三组*/
          qui tsset code_bs year

          qui xtreg `varlist' if gg==1,fe  /*第一组，视为小规模公司*/
          matrix b1 = e(b)
          qui xtreg `varlist' if gg==3,fe  /*第三组，视为大规模公司*/
          matrix b3 = e(b)
          matrix diff = b1 - b3
          mat D[`j',1] = diff
       }  
     *--------------------------------------------------------------------
      
       svmat D, names(diff)  
       * diff1-diff6 对应 size ndts tang npr tobin和常数项
      
       mat P = J(6,2,.)  /*记录系数真实差异和实证P值的矩阵*/
       forvalues j = 1/6{
          local diff0_`j' = D0[1,`j']
          qui count if (diff`j'>`diff0_`j'' | diff`j'==`diff0_`j'')&diff`j'!= .
          local p = `r(N)'/`reps'
          mat P[`j',1] = (`diff0_`j'' , `p')
       }
      
       mat colnames P = 系数真实差异  实证P值
       mat rownames P = size ndts tang npr tobin cons
       mat list P
      
     * 结论：除SIZE外，其它变量的组间系数估计值并不存在显著差异。

     * Reps = 5000 次的结果(大约需要20分钟)
     *----------------------------------------------
     *                系数真实差异        实证P值
     *----------------------------------------------
     * size         .0392026        .0032
     * ndts         .11700291        .7070
     * tang           -.04144458        .4664
     * npr        -.08150638        .7198
     * tobin          -.01324722        .1788
     * cons        -.62137931        .9848
     *----------------------------------------------

全是好淫，希望论坛上多些这些提问的人。问题说的很清楚，方便别人理解。向楼主致敬！！！

先把你的“全是好淫”，改成“全是好人”，真是无语，不知你是故意的，还是。。。对不起连帮主百忙之中热情解答，也对不起我发起提问！

感谢连帮主，现在国内学校很缺连老师这样的人才，更重要的是他能百忙之中抽出时间和我们无偿的交流，为人师表啊！敬佩他！

gushiydw 发表于 2011-10-9 19:17
先把你的“全是好淫”，改成“全是好人”，真是无语，不知你是故意的，还是。。。对不起连帮主百忙之中热情 ...

那娃娃显然是故意的，无论是在何种输入法下都无法直接敲出这四个字。
不过，这娃比较幽默，也算是帮大家缓和一下严肃的气氛吧。
这个发音是哪里人氏？东北银？

呵呵，感谢连帮主一直关注此贴！不跟他计较，看他的级别还是个“小学生”。。。。连帮主是版主级别，我是教授，哈哈！

不必计较这些虚衔，呵呵，我们只是在论坛里“浪费”了更多的时间而已，呵呵
只要大家开心，有所得就好。

连帮主好！我们用这种方法也要像您在12楼那样编写很多程序吗？如何用你的命令bdiff来做啊，编写程序自己还是有点困难。。。我利用一个平衡面板数据自己演示，得不到结果，系统提示两种情况：有时“varlist required   或者有时“factor variables and time-series operators not allowed”。不知到怎么回事，花了很多时间还是没搞清楚，求教！还有您举得例子中，将规模分3组，比较group1和group3的差异。那我如果做高科技行业样本组和非高科行业样本组之间的差异？比如我有一个18家公司6年的平衡面板数据（共108个观测值，其中高科技公司9家，共9*6=54个观测值，非高科技企业也是9家，共9*6=54个观测值），我要分析企业研发投资的业绩敏感性在高科技行业组和非高科行业组是否存在显著差异，数据中id是公司代码（1 2 3 .....数值型），year是年度（2005-2010，数值型），roe（企业业绩，是被解释变量），rd(研发投资，解释变量)，控制变量假定只有一个lnasset(企业规模)，假定回归模型如下，且是正确的，用indus表示行业，高科技行业取1，非高科技行业取0，高科技组和非高科技组都是9家公司共54个观测值。
    tsset id year
    xtreg roe rd lnasset if indus==1,fe
    xtreg roe rd lnasset if indus==0,fe

     我的问题是如何使用命令bdiff来检验高科技行业样本组研发投资的业绩敏感性与非高科行业研发投资业绩敏感性存在显著差异？自己水平太差，编写稍微长一点的程序难以胜任啊，求教连帮主！谢谢！