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

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

邮箱：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到计量社群提取.另外, 建议各位圈友积极转发或点赞，这应该形成一个良好的圈子文化.

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

更多相关文章

随机推荐