温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
作者:陈勇吏(上海交通大学)
目录
本文介绍 Stata 中做蒙特卡洛模拟的两种常用方法。第一种方法是使用 postfile
命令,第二种方法是 simulate
命令,并举了两个具体的例子,说明如何在 Stata 中做蒙特卡洛模拟。
蒙特卡洛模拟方法(MC),即从总体中抽取大量随机样本的计算方法。当根据总体的分布函数
比如,我们想知道
这是一个多重积分,大多时候很难求解。此时,我们可以使用蒙特卡洛模拟的方法,从总体中抽取大量的样本,通过样本来近似
MC 方法的原理是从总体中生成大量的样本,Stata 有两种常用的方法。第一种是使用 postfile
命令,另一种是使用 simulate
命令。
postfile
命令需要与循环语句结合使用。使用 foreach、while
等循环语句逐次生成独立的样本,并基于样本计算感兴趣的统计值。使用 postfile
命令生成 dta 文件,并将每次循环得到的数据追加进来。蒙特卡洛模拟次数由循环次数决定,最后生成的 dta 文件中,每一行代表一次蒙特卡洛模拟,每一列代表一个基于样本计算出的统计值。
postfile
命令语法查看帮助文件 help postfile
,可以看到三个配套使用的命令 postfile、post、postclose
。
postfile 命令
postfile
命令的原理是:在内存中划出一个区域,在该区域中存放 MC 中生成的数据。在这个过程中,需要给区域起一个名字,来告诉程序需要把数据存入哪一块区域中。同样地,需要给存入的数据取一个名字,方便索引到这个数据。具体语法规则如下:
postfile postname newvarlist using filename [, every(#) replace]
其中,
postname
表示内存中划出的区域的名字,filename
表示保存的数据的名字,newvarlists
表示数据中包含的变量名列表。简言之,postfile
的工作原理是:在 postname
内存区域中,生成一个 filename
数据文件,该数据包含的变量是 newvarlists
。
post 命令
post
命令的原理是,在 postfile
生成的数据文件中追加数据。具体语法规则如下:
post postname (exp) (exp) ... (exp)
其中,postname
告诉程序在哪一块内存区域追加数据(即追加到哪个内存区域对应的数据集中)。(exp) (exp) ... (exp)
为追加的数据内容,与 postfile
中的 newvarlist
变量名一一对应。
postclose 命令
postclose
表示结束追加数据,将所有蒙特卡洛模拟过程中 post
的数据写入硬盘中,可以使用 use
命令打开。
通常情况下,使用 postfile
做蒙特卡洛模拟的 Stata 代码格式如下:
tempname memhold //定义临时的内存区域的名字
tempfile results //定义临时的保存文件的名字
postfile `memhold' ... using `results' //在 `memhold' 内存区域中生成 `results' dta文件,包含 ... 中列出的变量。
forvalues i = 1/# { //循环语句,做 # 次蒙特卡洛模拟
... //根据总体分布,生成随机样本,进行相关运算
post `memhold' ... //追加感兴趣的计算值
}
postclose `memhold' //结束追加数据,完成 # 次蒙特卡洛模拟
use "`results'", clear //打开蒙特卡洛模拟生成的数据
simulate
命令的语法格式如下:
simulate [exp_list] , reps(#) [options] : command
其中,command
为一次蒙特卡洛模拟执行的程序,可以是 Stata 自带的或者外部下载安装的命令,也可以是用户自己封装的程序。[exp_list]
表示将每一次 command
命令的执行结果按 [exp_list]
格式提取出来。
exp_list 格式 |
表达式 | 举例 |
---|---|---|
(name: elist) | (scale: sd=r(sd) iqr=(r(p75)-r(p25)) | |
elist | newvar = (exp) (exp) |
mean=r(mean) r(sd) |
eexp | _b、_b[] _se、_se[] |
_b |
选项 reps(#)
表示做 # 次蒙特卡洛模拟,即执行 # 次 command
中的命令。其他选项 [options]
和功能如下:
options 作用
-----------------------------------------------------------
nodots MC 过程中不在屏幕上打印点
dots(#) 每隔 # 次模拟,在屏幕上打印一个点
noisily 将每次 MC 的结果都显示在屏幕上
trace 追踪命令运行过程
saving(filename, ...) 将 MC 结果保存在数据中
nolegend 不显示 MC 信息
verbose 显示 MC 信息
seed(#) 设定随机数种子为 #
从对数正态分布中,抽取样本量为 100 的样本,计算样本均值与方差。将上述过程重复 1000 次,查看 1000 次蒙特卡洛模拟的结果。
使用 postfile 命令:
clear
set obs 100
set seed 1234
gen b = .
tempname sim
tempfile results
postfile `sim' mean Var using "`results'", replace
quietly {
forvalues i = 1/1000 {
replace b = exp(rnormal())
summarize
scalar mean = r(mean)
scalar Var = r(Var)
post `sim' (mean) (Var)
}
}
postclose `sim'
use "`results'", clear
sum
使用 simulate 命令:
//封装代码
cap program drop lnsim
program lnsim, rclass
version 13
drop _all
set obs 100
gen z = exp(rnormal())
summarize z
return scalar mean = r(mean)
return scalar Var = r(Var)
end
set seed 1234
simulate mean=r(mean) var=r(Var), reps(1000) nodots: lnsim
describe *
sum
两个命令的结果是一样的:
. describe *
storage display value
variable name type format label variable label
-----------------------------------------------------------------------
mean float %9.0g r(mean)
var float %9.0g r(Var)
. sum
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
mean | 1,000 1.630648 .2173062 1.106372 2.612052
var | 1,000 4.60798 4.502166 .966087 70.55969
实验拟合模型
显然,当
在实验过程中,我们将维持保存参数估计值和标准误差,并尝试不同的
在这之前,生成一个真实数据 truth.dta :
drop _all
set obs 100
set seed 54321
gen x = rnormal()
gen true_y = 1 + 2*x
save truth.dta, replace
根据 truth.dta 数据,做 10000 次蒙特卡洛模拟。在每一次模拟过程中:
reg y x
命令拟合模型,提取参数估计值和标准误。显然,如果设定
使用 postfile 命令:
use truth.dta, clear
set seed 123
local c = 3
tempname sim
tempfile results
postfile `sim' b se using "`results'", replace
quietly {
forvalues i = 1/10000 {
capture drop y
gen y = true_y + (rnormal() + `c'*x)
reg y x
post `sim' (_b[x]) (_se[x])
}
}
postclose `sim'
use "`results'", clear
sum
使用 simulate 命令:
set seed 123
cap program drop hetero
program hetero
version 13
args c
capture drop y
gen y = true_y + (rnormal() + `c'*x)
regress y x
end
simulate _b _se, reps(10000): hetero 3
两个命令的结果是一样的:
===== postfile 命令的结果 ====
. sum
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
b | 1,000 5.000521 .0998493 4.748041 5.342362
se | 1,000 .0996826 .0069156 .0757347 .1218035
===== simulate 命令的结果 ====
. sum
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
_b_x | 1,000 5.000521 .0998493 4.748041 5.342362
_b_cons | 1,000 .9970645 .0972862 .6995795 1.323864
_se_x | 1,000 .0996826 .0069156 .0757347 .1218035
_se_cons | 1,000 .0997151 .0069178 .0757594 .1218432
在本例中,由于我们设定
有兴趣的读者,可重新设计一个 DGP (数据生成过程),把工具变量估计考虑进去,进而模拟一下 IV 或 2SLS 估计法是否存在偏误。如下是我建议的一个 DGP,我们将在下一篇推文中公布完整的分析过程。
*-Ref: Cameron and Trivedi (2009), pp.143 Section 4.6.5
*--------------DGP----------------
*
* y = a + b*x + u; u~N(0,1)
*
* x = z + 0.5*u; z~N(0,1)
*
* a=10; b=2; N=150
*
*---------------------------------
*-Q1: OLS 估计的性质如何?
*
*-Q2: IV 估计的性质如何?
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看,所有课程可以随时购买观看。
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
⭕ Stata数据清洗 | 游万海 | 直播, 2 小时,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD