Stata：蒙特卡洛模拟新命令-simulate2
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 五一论文班
作者：李坤 (华南师范大学)
邮箱：kunli1001@163.com
目录
1. 引言
2. 命令介绍
3. 具体案例
3.1 simulate2 命令
3.2 psimulate2 命令
3.3 simulate 与 simulate2 的区别
4. 结语
5. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 引言
蒙特卡洛模拟 (MC) 是指当根据总体的分布函数很难求出想要的数字特征时，从总体中生成大量的样本，使用样本的数字特征估计总体的数字特征。Stata 有两种常用的方法进行 MC 模拟，第一种是使用 postfile 命令，另一种是使用simulate 命令。
在本文中，我们为大家介绍一个新的高效的 MC 模拟命令 simulate2，以及另外一个可以进行并行实验的命令 psimulate2，这两个命令可以大大简化 MC 模拟的程序并提高运行的速度。
2. 命令介绍
simulate2 命令安装：
ssc install simulate2, replace
需要注意，simulate2 命令需要 Stata16 及以上才能正常使用；psimulate2 命令需要 Stata15 及以上才能正常使用，目前仅适用于 Windows 系统。
simulate2 命令语法：
simulate2 [exp_list] , reps(#) [options1 options2 ] : command
其中，[exp_list] 是指将每一次 command 的执行结果按照 [exp_list] 格式提取出来，若不是 [exp_list] 格式则会报错。command 代表每一次蒙特卡洛模拟执行的程序，可以是 Stata 自带或者外部命令，也可以是用户自己编写的程序。rep(#) 代表做 # 次蒙特卡洛模拟。
[options1] 主要选项如下：
nodots：在 MC 模拟过程中不在屏幕上打点；dots(#)：每隔 # 次模拟，在屏幕上打印一个点；noisily：将每次 MC 的结果都显示在屏幕上；trace：追踪命令运行的过程；nolegend：不显示 MC 的信息；verbose：显示 MC 信息。
[options2] 主要选项如下：
aving(filename, …)：将结果保存到一个文件夹中；seed(options)：设定随机种子数；seedsave(options)：保存用过的种子。
psimulate2 命令语法：
psimulate2 [exp_list] , reps(#) parallel(#2, options3) [options2] : command
其中，parallel(#2,()) 中的 #2 是设置并行运行 Stata 实例的数量，建议使用的实例数量不要超过可用的 CPU 内核。options3 主要选项如下：
exe(string)：设置 Stata.exe 的路径，若 Stata.exe 位于非常规文件夹中或具有非常规文件名，可能会查找失败；temppath(string)：设置保存临时文件的替代路径；processors(integer)：设置每个 Stata 实例允许的最大处理数，仅适用于 Stata MP 版本；simulate：使用simulate而不是simulate2，如果使用 Stata15 运行psimulate2，则默认使用simulate；seed(#)：设置随机种子数；seedstream(integer)：种子流，可以为第一个实例设置初始种子流编号。例如parallel(3, seedstream(4))表示实例 1 将使用种子流编号为 4、实例 2 将使用种子流编号为 5，依次类推。该选项只能在psimulate2中使用。
3. 具体案例
3.1 simulate2 命令
我们先通过 1000 次样本容量为 100 的 MC 模拟实验来估计 OLS 的参数和标准误，并记录时间。在生成的 dta 文件中，每一行的数据代表一次蒙特卡洛模拟，每一列的数据代表一个基于样本计算出的统计值，最后可以通过计算上述数据的均值来近似估计回归方程的系数和标准误。
在上述过程中，我们将种子保存到名为 seed_frame 的 frame 中。在 Stata16 之前的版本中，Stata 的内存只能存储一份数据集，这在实际中非常的不方便。frame 则可以允许在 Stata 内存中同时存储、操作多份数据，并且可以将这些数据连接起来，极大的提升了数据处理的效率。
* 定义程序 testsiuml，并以 r() 形式储存结果
program define testsimul, rclass
version 17 // 使用的Stata的版本
syntax anything
clear
set obs `anything'
gen x = rnormal(1,4)
gen y = 2 + 3*x + rnormal() // 假设真实的模型参数为 a=2 b=3
reg y x // 用 reg y x 命令来拟合模拟
matrix b = e(b) // 提取参数估计值
matrix se = e(V) // 提取标准误的估计值
ereturn clear
return scalar b = b[1,1]
return scalar V = se[1,1]
return local time "`c(current_time)'"
end
* 选取 100 个样本进行 1000 次 MC 模拟
simulate2 time = r(time) b = r(b) V = r(V), reps(1000) seedsave(seed_frame, frame): testsimul 100
现在我们可以利用保存的种子重新进行前 500 次实验：
. simulate2 time = r(time) b = r(b) V = r(V), reps(500) seed(seed_frame seed, frame): testsimul 100
对第二个 500 次重复实验，将起始种子编号设置为 501：
. simulate2 time = r(time) b = r(b) V = r(V), reps(500) ///
> seed(seed_frame seed, frame start(501)): testsimul 100
与上述的步骤相反，我们也可以先进行 500 次抽样保存种子，再进行 500 次抽样追加保存种子。最后对所有 1000 次抽样进行实验，并将比较结果保存在 frame 中。
. simulate2 time = r(time) b = r(b) V = r(V), reps(500) seed(123) ///
> seedsave(seed_frame, frame) saving(first500, frame): testsimul 100
. * 对于第二次运行，没有设置新种子，而是使用了选项 append
. simulate2 time = r(time) b = r(b) V = r(V), reps(500) ///
> seedsave(seed_frame, frame append) saving(second500, frame): testsimul 100
. simulate2 time = r(time) b = r(b) V = r(V), reps(1000) ///
> seed(seed_frame seed, frame): testsimul 100
分别比较用两种方法计算得到的前 500 次重复模拟结果与后 500 次重复模拟的结果：
. sum if _n <= 500
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
b | 500 3.001443 .024907 2.921489 3.071222
V | 500 .0006364 .0001251 .000357 .0010889
time | 0
. frame first500: sum
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
b | 500 3.001443 .024907 2.921489 3.071222
V | 500 .0006364 .0001251 .000357 .0010889
time | 0
. sum if _n > 500
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
b | 500 2.999219 .0241378 2.930022 3.078007
V | 500 .0006446 .0001321 .0003689 .0012382
time | 0
. frame second500: sum
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
b | 500 2.999219 .0241378 2.930022 3.078007
V | 500 .0006446 .0001321 .0003689 .0012382
time | 0
可以看出两个命令的结果是相同的，表明最终结果与进行实验的方式是无关的。同时观察数据可以看出，估计出的值都接 3，与我们设定的回归方程的系数非常接近。
3.2 psimulate2 命令
psimulate2 命令是 simulate2 的并行版本。它的工作原理是将复制的数据分成相等的块，每个块在单独的 Stata 实例上运行。为此 psimulate2 创建了一个 do 文件和一个批处理文件。然后使用批处理文件来启动一个运行相应 do 文件的新 Stata 实例。正在运行的实例充当父实例，另一个是子实例。
psimulate2 在 Stata 上输出与 simulate 或 simulate2 的输出不同。它显示已经运行完成的百分比、经过的时间和预期的剩余时间以及完成时间。因此，我们可以将上面的实验使用 psimulate2 进行并行化模拟。例如我们要同时运行两个实例，即实例 1 运行第一个 500 次重复，实例 2 运行第二个 500 次重复，采用下述的代码：
. psimulate2, reps(500) seed(123) parallel(2, temppath("D:\Stata17\do")): testsimul 200
. sum // 查看结果
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
V | 500 .0003177 .0000478 .0002151 .0005266
b | 500 2.999787 .0187097 2.946175 3.04934
time | 0
parallel(2, temppath("D:\Stata17\do")) 含义为设置两个并行实例，temppath() 指定保存临时文件的路径。仅当 Windows 不允许从临时文件夹运行批处理文件时，才需要此选项。运行第一个和第二个 500 次重复具有相同的随机数状态，但种子流不同，因此随机抽取也不同。可以看出，该命令最后的结果与使用 simulate2 命令计算出的结果都是接近的，但是运行速度更快，语法结构也更加简洁。
当我们需要将一个或多个 psimulate2 嵌套在循环中或按顺序调用，它们将使用相同的初始种子，可能造成最终的结果出现偏误。因此为了避免这种情况，但还能按顺序运行 psimulate2，可以使用以下代码来让 psimulate2 返回 r() 中最后一个实例的最后一个种子状态。
. psimulate2, reps(100) seed(_current) p(2) seedsave(seed, frame): testsimul 200
. set rng `r(rng_current)'
. set rngstate `r(rngstate)'
. psimulate2, reps(100) seed(_current) p(2) seedsave(seed, frame): testsimul 200
. sum
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
V | 100 .0003271 .0000516 .0002419 .0004835
b | 100 3.001315 .0181299 2.961705 3.052167
time | 0
在 psimulate2, seed() 中使用 current，可以将使用 current 父实例的种子作为所有子实例的初始种子，这样每个子实例仍将具有不同的种子流，以确保随机数抽取不同。
3.3 simulate 与 simulate2 的区别
simulate2 遵循 simulate 的语法并简化了执行蒙特卡洛模拟的编程任务。相比于 simulate，它的新颖之处在于允许其调用的程序将宏 (字符串) 返回到 r() 和 e() 来扩展命令，然后将它们保存到使用的数据集中。使用该命令时，还可以产生一个当下时间的数据集，记录该结果产生的时间。
此外，它还具有高级选项，可将随机数状态 (种子、种子流和类型) 分配给特定的模拟重复，种子可以保存到 Stata 数据集 (dta) 或 frame 中，以便无缝复制模拟结果，模拟结果也可以直接保存到 frame 中。
4. 结语
通过上面的介绍，向大家展示了 simulate2 以及 psimulate2 的具体使用方法。可以发现当需要进行 MC 模拟时，simulate2 和 psimulate2 命令拥有很大的优势，不仅语法简洁，而且使用起来也很方便。我们希望通过这篇文章，能够让大家对这两个命令有一个更加深入的认识。
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh MC 模拟, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata教程 Stata博文系列：边模拟边学计量经济学-Econometrics-By-Simulation 专题：Stata程序 Stata：蒙特卡洛模拟A-(Monte-Carlo-Simulation)没那么神秘 Stata程序：Monte-Carlo-模拟之产生符合特定分布的随机数 Stata：蒙特卡洛模拟分析 (Monte Carlo Simulation) 专题：回归分析 Stata：在线可视化模拟-OLS-的性质 专题：断点回归RDD Stata 新命令：多断点 RDD 分析 - rdmc 专题：时间序列 Stata: VAR - 模拟、估计和推断 专题：Python-R-Matlab MATLAB：Ramsey模型的求解与模拟 专题：其它 世代交叠OLG模型的求解与模拟
课程推荐：因果推断实用计量方法
主讲老师：丘嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【**百度一下：**连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。