查看原文
其他

概率分布生成和蒙特卡洛模拟的实战示例

计量经济圈 计量经济圈 2019-06-30

可有偿投稿计量经济圈,计量相关则可

箱:econometrics666@sina.cn

所有计量经济圈方法论丛的do文件, 微观数据库和各种软件都放在社群里.完整版本do file到计量社群提取.

先看看这个: 编程语言中的函数什么鬼?Stata所有函数在此集结

直接放上运行程序, **后面是每个子程序的解释。

产生不同分布的随机数过程

**single draw of a uniform number
set seed 10101
scalar u = runiform ( )
display u


**1000 draws of uniform numbers
quietly set obs 1000
set seed 10101
generate x = runiform()
list x in 1/5, clean


**First three autocorrelations for the uniform draws
generate x = _n
tsset x
pwcorr x L.x L2.x L3.x , star (0.05)


**normal and uniform
clear
quietly set obs 1000
set seed 10101
generate uniform = runiform() //uniform (0, 1)
generate stnormal = rnormal( )
generate norm5and2 = rnormal(5, 2) //N (0, 1)
tabstat uniform stnormal norm5and2, stat(mean sd skew kurt min max) col(stat)


**t, chi-squared, and F with constant degrees of freedom
clear
quietly set obs 2000
set seed 10101
generate xt= rt(10) //result xt~ t(10)
generate xc = rchi2(10)/10 //result xc~ chisquared(10)
generate xfn = rchi2(10)/5 //result numerator of F(10, 5)
generate xfd = rchi2(10)/5 //result denominator of F(10, 5)
generate xf = xfn/xfd //result xf ~ F(10, 5)
summarize xt xc xf


*binomial
set seed 10101
generate p1 = runiform() //here p1~ uniform(0,1)
generate trials = ceil(10runiform()) //here # trials varies btwn 1&10
generate xbin = rbinomial(trials, p1) //draws from binomial(n, p1)
summarize p1 trials xbin
independent poisson and negative bin draws
set seed 10101
generate xb= 4 + 2runiform()
generate xg = rgamma(1,1) //draw from gamma; E(v) = 1
generate xbh = xbxg //apply multiplicative heterogeneity
generate xp = rpoisson(5)
generate xp1 = rpoisson(xb)
generate xp2 = rpoisson(xbh)
summarize xg xb xp xp1 xp2


**Example of histogram and kernel density plus graph combine
quietly twoway (histogram xc) (kdensity xc), title("Draws from chisquared(lO)")
quietly graph save mus04cdistr.gph, replace
quietly twoway (histogram xp) (kdensity xp), title("Draws from Poisson(mu) for 5<mu< 6")
quietly graph save mus04poissdistr.gph, replace
graph combine mus04cdistr.gph mus04poissdistr.gph, title("Random-number generation examples", margin(b=2) size(vlarge))


**Draw 1 sample of size 30 from uniform distribution
set obs 3000
set seed 10101
generate x=runiform( )


**Summarize x and produce a histogram
summarize x
quietly histogram x, width(0.1) xtitle("x from one sample")
**Program to draw 1 sample of size 30 from uniform and return sample mean
program onesample1, rclass
drop _all
quietly set obs 30
generate x = runiform( )
summarize x
return scalar meanforonesample = r(mean)
end


**Run program onesample once as a check----
set seed 10101
onesample1
return list


**Run program onesample 10000 times to get 10000 sample means
simulate xbar=r(meanforonesample), seed(10101) reps(10000) nodots: onesample1


**Summarize the 10000 sample means and draw histogram
summarize xbar
quietly histogram xbar, normal xtitle("xbar from many samples")


**Inverse probability transformation example: standard normal
set obs 2000
set seed 10101
generate xstn = invnormal(runiform( ))


**Inverse probability transformation example: unit exponential
generate xue = -ln(1-runiform( ))


**Inverse probability transformation example: Bernoulli (p=0.6)
generate xbernoulli = runiform() > 0.6 //Bernoulli(0.6)
summarize xstn xue xbernoulli
Draws from truncated normal x - N(mu, sigma-2) in [a,b]
quietly set obs 2000
set seed 10101
scalar a=0 //lower truncation point
scalar b=12 //upper truncation point
scalar mu=5 //mean
scalar sigma = 4 //standard deviation
generate u = runiform()
generate w=normal((a-mu)/sigma)+u(normal((b-mu)/sigma)-normal((a-mu)/sigma))
generate xtrunc = mu + sigmainvnormal()
summarize xtrunc


*Bivariate normal example:
means 10 , 20; variances 4, 9; and correlation 0 . 5
clear
quietly set obs 1000
set seed 10101
matrix MU=(10, 20)
scalar sig12 = 0.5sqrt(49)
matrix SIGMA = (4, sig12\sig12, 9) //SIGMA is 2 x 2
drawnorm y1 y2, means(MU) cov(SIGMA)
summarize y1 y2
correlate y1 y2


**Integral evaluation by Monte Carlo simulation Yith S=100
clear all
quietly set obs 100
set seed 10101
generate double y= invnormal(runiform( ))
generate double gy=exp(-exp(y))
quietly summarize gy, meanonly
scalar Egy=r(mean)
display "After 100 draYs the MC estimate of E[exp(-exp(x))] is "Egy


**Inconsistency o f OLS in errors-in-variables model (measurement error)--
clear
quietly set obs 10000
set seed 10101
matrix mu=(0 , 0 , 0)
matrix sigmasq=(9 , 0 ,0 \0 , 1 , 0\0 , 0 , 1)
drawnorm xstar u v , means(mu) cov(sigmasq)
generate y =xstar + u //DGP for y depends on xstar
generate x= xstar + v //x is mimeasured error
regress y x, noconstant

蒙特卡洛模拟回归过程

*Define program myreg to generate data and fit a linear regression-
program myreg, eclass //定义一个叫myreg的程序,rclass表示结果以r()形式储存
version 15 //设定所用程序版本为Stata 15
drop _all //删去内存中已有数据
set obs 100 //确定随机抽样的样本容量为50
generate x = rnormal() //生成服从分布的解释变量x
generate y = 3x + 1 + rnormal() //生成被解释变量
regress y x //线性回归
end


**Perform simulation----------
simulate _b _se, reps(500) seed(5762): myreg //模拟估计系数和标准差
summarize
program lnsim, rclass
version 15.0
drop _all
set obs 100
generate 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 *
summarize
set seed 1234
lnsim


*simulating a regression model----------------
drop _all
set obs 100
set seed 54321
generate x = rnormal() //产生服从正态分布的随机数x
generate true_y = 1+2x //产生true_y随机数
save truth.dta, replace
program hetero1 //定义hetero1程序
version 15.0
args c //args针对的也是这一系列macros
use truth, clear
generate y = true_y + (rnormal() + c'*x) //生成数据y regress y x //做y对x的回归 end simulate _b _se, reps(10000): hetero1 3 //模拟当c=3时x的系数和标准误差 program hetero3 version 15.0 args truey x c capture drop y generate y =truey' + (rnormal() + c'*x')
regress y x
end


simulate _b _se, reps(10000): hetero3 true_y x 3 //模拟当truey=true_y, x=x, c=3


*simulating a ratio of statistics-----
program myratio, rclass
version 15.0
args n1 mu1 sigma1 n2 mu2 sigma2
// generate the data
drop _all
local N = n1'+n2'
set obs N' tempvar y //assigns names to the specified local macro names generatey' = rnormal()
replace y' = cond(_n<=n1',mu1'+y'sigma1',mu2'+y'*sigma2')
// calculate the medians
tempname m1
summarize y' if _n<=n1', detail
scalar m1' = r(p50) summarizey' if _n>n1', detail // store the results return scalar ratio =m1' / r(p50)
end


set seed 19192
simulate ratio=r(ratio), reps(1000) nodots: myratio 5 3 1 10 3 2


**第二个程序--
program define bsse, rclass
version 15.0
drop _all
set obs 100
generate x = rnormal() //生成符合正太分布的数据x
tempfile bsfile //assigns names to the specified local macro names
bootstrap midp=r(p50), rep(100) saving(bsfile'): summarize x, detail //自助法获得中位数 usebsfile', clear
summarize midp
return scalar mean = r(mean)
return scalar sd = r(sd)
end


set seed 48901
simulate med=r(mean) bs_se=r(sd), reps(1000): bsse //自助法得到中位数的均值和标准差

完整版本do file到计量社群提取.另外, 建议各位圈友积极转发或点赞,这应该形成一个良好的圈子文化.


可以到计量经济圈社群进一步访问交流各种学术问题,这年头,我们不能强调一个人的英雄主义,需要多多汲取他人的经验教训来让自己少走弯路。


计量经济圈推荐

1.PSM-DID, DID, RDD, Stata程序百科全书式的宝典
2.RDD断点回归, Stata程序百科全书式的宝典
3.Generalized分位数回归, 新的前沿因果推断方法
4.Heckman模型out了,内生转换模型掌控大局
5.PSM倾向匹配Stata操作详细步骤和代码,干货
6.条件Logit绝对不输多项Logit,而混合模型最给力
7.广义PSM,连续政策变量因果识别的不二利器
8.自回归VAR模型操作指南针,为微观面板VAR铺基石
9.有限混合模型FMM,异质性分组分析的新筹码
10.政策评估中"中介效应"因果分析, 有趣的前沿方法
11.多期三重差分法和双重差分法的操作指南
12.多期双重差分法,政策实施时间不同的处理方法
13.随机前沿分析和包络数据分析 SFA,DEA 及操作
14.你的内生性解决方式out, ERM已一统天下而独领风骚
15.多期DID的经典文献big bad banks数据和do文件
16.面板数据里处理多重高维固定效应的神器
17.双栏模型Hurdle远超Tobit, 对于归并数据舍我其谁
18.面板数据计量方法全局脉络和程序使用指南篇

19.Clad还是Tobit, 归并最小绝对偏差, 做Tobit做不好的

20.合成控制法什么鬼? 因果推断的前沿方法指南

21.空间计量百科全书式的使用指南, 只此一份

22.实证研究中交叉项的使用和解读策略指南案例

23.各领域经济学手册全在这里, 不只做重复研究

24.事件研究法什么鬼? 从这里着手看"疫苗之王"

25.因果推断异质性什么鬼? 边际处理效应让你与众不同

计量经济圈当前有几个阵地,他们分别是如下4个matrix:

①计量经济圈社群——计量经管数据软件等资料中心,

②计量经济圈微信群——服务于计量经济圈社群群友,

③计量经济圈研究小组系列——因果推断研究小组、空间计量研究小组、面板数据库研究小组、微观计量研究小组、计量软件研究小组,

④计量经济圈QQ群——2000人大群服务于计量经济圈社群群友。


计量经济圈是中国计量第一大社区,我们致力于推动中国计量理论和实证技能的提升,圈子以海内外高校研究生和教师为主。计量经济圈六多精神:计量资料多,社会科学数据多,科研牛人多,名校人物多,热情互助多,前沿趋势多。如果你热爱计量并希望长见识,那欢迎你加入到咱们这个大家庭(戳这里),要不然你只能去其他那些Open access圈子了。注意:进去之后一定要看小鹅社群“群公告”,不然接收不了群息,也不知道怎么进入咱们独一无二的微信群和QQ群。


进去之后就能够看见这个群公告了



    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存