缘起

• 之前写了一个stata如何做placebo test的文章，本来只是写给自己看，记录自己的学习过程，没想到浏览上千，也是有点意外加激动。
• 有些朋友说根据我写的代码，改到自己的项目里，跑不出来结果，正好今天有学弟问我这个问题，于是我把代码重写了一次，这次用stata官方给的auto数据写的，应该没有bug了，有问题欢迎留言评论。
• 上次的文章链接：

https://www.jianshu.com/p/ef4f920494d
（该文粗略讲了安慰剂的原理，并提供一版自明性相对差点的stata代码，本文在这篇文章的基础上，修改了stata代码。）

实现

假设我的数据集是auto.dta
因变量为price
自变量为rep78
控制变量为headroom weight length
回归模型为OLS，即主回归代码为reg rep78 headroom weight length
下面是1000次placebo test的代码

forvalue i=1/1000{
    sysuse auto, clear  //调入数据
    *- 思路：打乱rep78,即将rep78的全部取值拿出暂存，然后随机赋给每一个样本

    *- 打乱rep78,即将rep78的全部取值拿出暂存
    g obs_id= _n //初始样本序号
    gen random_digit= runiform() //生成随机数
    sort random_digit  //按新生成的随机数排序
    g random_id= _n  //产生随机序号
    preserve
        keep random_id rep78 //保留虚拟的rep78
        rename rep78 random_rep78
        rename random_id id //重命名为id，以备与其他变量合并（merge）
        label var id 原数据与虚拟处理变量的唯一匹配码
        save random_rep78, replace
    restore 
        drop random_digit random_id rep78 //删除原来的rep78
        rename obs_id id //重命名为id，以备与random_rep78合并（merge）
        label var id 原数据与虚拟处理变量的唯一匹配码
        save rawdata, replace 

    *- 合并，回归，提取系数
        use rawdata, clear
        merge 1:1 id using random_rep78,nogen
        reg  price random_rep78 headroom weight length
        g _b_random_rep78= _b[random_rep78]  //提取x的回归系数
        g _se_random_rep78= _se[random_rep78] //提取x的标准误
        keep _b_random_rep78 _se_random_rep78 
        duplicates drop _b_random_rep78, force
        save placebo`i', replace  //把第i次placebo检验的系数和标准误存起来
    }
    
*- 纵向合并1000次的系数和标准误 
use placebo1, clear
forvalue i=2/1000{
    append using placebo`i' //纵向合并1000次回归的系数及标准误
}  
 
gen tvalue= _b_random_rep78/ _se_random_rep78
kdensity tvalue, xtitle("t值") ytitle("分布") saving(placebo_test) 
    
*-删除临时文件
forvalue i=1/1000{
    erase  placebo`i'.dta 
} 

*-mark：后续有时间再将它封装成程序【完】

文献范例

• 周茂、陆 毅、杜 艳、 姚 星，2018：《开发区设立于地区制造业升级》，《中国工业经济》第3期。

image.png