機率分佈生成和蒙特卡洛模擬的實戰示例

可有償投稿計量經濟圈，計量相關則可

郵箱：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到計量社羣提取.另外, 建議各位圈友積極轉發或點贊，這應該造成一個良好的圈子文化.

能夠到計量經濟圈社羣進一步訪問交流各類學術問題，這年頭，咱們不能強調一我的的英雄主義，須要多多汲取他人的經驗教訓來讓本身少走彎路。