chap15 蒙特卡羅模擬和自助法

1、主要內(nèi)容隨機數(shù)的生成蒙特卡羅模擬重復抽樣自助法實驗實驗15-1： 隨機數(shù)的生成隨機數(shù)的生成實驗基本原理實驗基本原理服從在(0,1)區(qū)間上均勻分布的隨機變量的樣本值，被稱為“均勻分布隨機數(shù)”，簡稱“隨機數(shù)”（random number）。有了均勻分布的隨機數(shù)后，就幾乎可以產(chǎn)生任何分布的隨機數(shù)了。然而，由電腦產(chǎn)生的所謂“隨機數(shù)”只是“偽隨機數(shù)”（pseudo random number），因為它仍然由遞推公式產(chǎn)生。盡管如此，這些“偽隨機數(shù)”都能通過獨立性與均勻分布的統(tǒng)計檢驗，故可以作為“隨機數(shù)”來使用?！半S機數(shù)”初始值被稱為“種子”（seed）。為了使隨機樣本具有可重復性，保證以后抽樣或別人抽樣也

2、能得到完全一樣的樣本，在生成隨機數(shù)時通常需要設定“種子”。如果不給定“種子”，則Stata按照“電腦內(nèi)置鐘表”（computer clock）來自動選擇“種子”，這樣每次抽樣的樣本就會不同，模擬的結(jié)果也就不能完全復制。實驗內(nèi)容及數(shù)據(jù)來源實驗內(nèi)容及數(shù)據(jù)來源本實驗中，我們會介紹如何生成均勻分布和正態(tài)分布的隨機數(shù)，以及如何設定隨機數(shù)的種子。講解這些內(nèi)容不需要使用具體的數(shù)據(jù)文件，一些簡單的命令就可以實現(xiàn)。實驗操作指導實驗操作指導1 隨機數(shù)的生成隨機數(shù)的生成最基本而常用的隨機數(shù)序列是均勻分布隨機數(shù)序列。生成均勻分布隨機數(shù)序列的基本命令為：generate newvar=runiform()其中，gene

3、rate為生成新變量的基本命令，newvar為新變量的名稱，runiform()是生成均勻分布于區(qū)間0, 1)隨機數(shù)的函數(shù)。需要注意的是，runiform()中沒有參數(shù)，但括號卻必不可少。如果要生成位于其他區(qū)間的均勻分布，我們可以進行簡單的變形。例如，要生成均勻分布于區(qū)間a, b)的隨機數(shù)，相應的函數(shù)為：a+(b-a)* runiform()要生成均勻分布于區(qū)間a, b的隨機數(shù)，相應的函數(shù)為：a+int(b-a+1)* runiform()其中，函數(shù)int()表示取整。生成標準正態(tài)分布的隨機數(shù)的函數(shù)為：invnorm(uniform()生成均值為m、標準差為s的正態(tài)分布的隨機數(shù)的函數(shù)則是：m+

4、s*invnorm(runiform()例如，我們要生成一個均值為3、方差為5且服從正態(tài)分布的序列，可輸入命令：set obs 80gen norm=3+5*invnorm(uniform()這里，第一步為設定觀測值的個數(shù)為80；第二步我們生成一個均值為3、方差為5且服從正態(tài)分布的序列，并將新生成的變量命名為norm。需要說明的一點是，如果不設定觀測值個數(shù)，則新變量的觀測值個數(shù)會與原序列的觀測值個數(shù)相同；而未打開任何數(shù)據(jù)文件時，原觀測值個數(shù)顯然為0。下面，我們看一下變量norm的描述統(tǒng)計量。輸入命令：sum norm如果我們要作圖看一下norm的分布，可輸入命令：hist norm, norm

5、al這里，hist表示做直方圖，選項normal表示畫出相應的正態(tài)分布。2 種子的設定種子的設定設定種子的基本命令為：set seed #|code其中，set seed是設定種子的基本命令。#是正整數(shù)，代表runiform()迭代的初始值， Stata的初始設置是123456789。code是字符串，代表隨機數(shù)生成器runiform()的狀態(tài)，我們可以用“c(seed)”來查看其當前值。例如，我們輸入命令：disp c(seed)這里，disp表示“顯示”，c(seed)表示種子的當前值。例如，我們先設定種子，然后生成隨機數(shù)序列，可輸入如下命令：set seed 1001set obs 60

6、gen r=runiform()這里，第一步設定種子為1001，第二步設定觀測值個數(shù)為60，第三步生成位于0, 1)隨機數(shù)序列，并將其命名為r。實驗15-2： 蒙特卡羅模擬實驗基本原理實驗基本原理通過計算機模擬從已知分布的總體中抽取大量隨機樣本的計算方法被統(tǒng)稱為“蒙特卡羅方法”（Monte Carlo Methods）。在計量經(jīng)濟學中，常使用蒙特卡羅法來確定統(tǒng)計量的小樣本性質(zhì)。我們知道，許多統(tǒng)計量的精確分布沒有解析解。一種解決方法是使用大樣本理論，用漸近分布來近似真實分布。然而，現(xiàn)實中的樣本容量常常較小。實驗內(nèi)容及數(shù)據(jù)來源實驗內(nèi)容及數(shù)據(jù)來源本實驗中，我們會介紹蒙特卡洛模擬的基本操作，以及如何編

7、寫程序并通過模擬獲得真實的顯著性水平。這里，我們不需要使用具體的數(shù)據(jù)文件，通過命令就可以實現(xiàn)相應的模擬。實驗操作指導實驗操作指導1 蒙特卡羅模擬的基本操作蒙特卡羅模擬的基本操作2 用蒙特卡羅模擬獲得真實顯著性水平用蒙特卡羅模擬獲得真實顯著性水平為此，先使用命令program定義一個叫“reschi2”程序進行一次抽樣與檢驗，然后用命令simulate來重復此程序1,000 次。程序的命令為：program reschi2, rclass （定義一個叫reschi2的程序。rclass表示結(jié)果以 r()形式儲存）version 10 （設定所用程序版本為Stata 10）drop _all （刪

8、去內(nèi)存中已有數(shù)據(jù)）set obs 50 （確定隨機抽樣的樣本容量為50）gen double x = rchi2(1) （生成服從分布的解釋變量x）gen y = 3+ 2*x + rchi2(1)-3 （生成被解釋變量）reg y x （線性回歸）return scalar t2 = (_bx-2)/_sex （計算t統(tǒng)計量的值，并將其命名為t2。這里，_bx代表x的系數(shù)，_sex代表系數(shù)的標準差，scalar表明t2是一個標量（數(shù)值）return scalar r2 = abs(return(t2)invttail(48,0.025) （是否拒絕原假設）end （程序結(jié)束）其中，命令的倒數(shù)

9、第二行的invttail(48, 0.025)表示自由度為48（樣本容量50-參數(shù)個數(shù)2）、顯著性水平為5%的雙邊檢驗的t統(tǒng)計量的臨界值。return(t2)代表我們前面模擬計算并保存的樣本t統(tǒng)計量的值t2。如果t2的絕對值大于invttail(48, 0.025)，則應拒絕原假設；這時，r2的返回值為1。此外，由于程序的首行定義結(jié)果的形式為rclass，所以，如果我們要引用計算的r2的值，其形式應該是“r(r2)”。對于上面的程序語句，我們可以將其寫入一個“Do-File Editor”，將程序保存為擴展名為“ado”的文件，并存放在Stata安裝目錄的ado文件夾下。這樣，我們后面就可以直

10、接使用該程序了。下面，我們進行1000次模擬，輸入命令：simulate reject=r(r2), reps(1000) nodots seed(101): reschi2這里，每次模擬的命令來自于程序“reschi2”，選項reps(1000)表明模擬次數(shù)為1000次，選項seed(101)表明設定種子為101，選項nodots表示不顯示模擬過程的點。表達式“reject=r(r2)”表明我們將每次模擬的返回值r(r2)保存在reject中。因為我們在程序“reschi2”中，每次拒絕原假設時r2的返回值為1；這樣，返回值為1的比例即為拒絕原假設的比例。也就是說，序列reject的均值反映

11、了真實的顯著性水平。為了獲得reject的均值，我們輸入命令：mean reject 事實上，我們在模擬之前，可以先運行一次我們所編的程序，檢查一下程序是否有問題。對于我們這個程序，可輸入如下命令：set seed 2010reschi2其中，第一步設定種子為2010，第二步為運行程序。要知道這一次模擬是否拒絕原假設，我們可以輸入命令：disp r(r2)另一種返回程序結(jié)果的命令語句為：return list實驗實驗15-3： 重復抽樣重復抽樣實驗基本原理實驗基本原理實驗內(nèi)容及數(shù)據(jù)來源實驗內(nèi)容及數(shù)據(jù)來源本書附帶光盤data文件夾下的“gender.dta”工作文件，是我們?yōu)榱酥v解重復抽樣而編制

12、的文件。該文件包括5810個觀測值，只有一個變量gender，表示性別。其取值為1時表示女性（female），取值為0時表示男性（male）。該文件中，女性觀測值個數(shù)為3418，男性觀測值個數(shù)為2392。此外，本實驗中，我們還會利用本書附帶光盤data文件夾下的數(shù)據(jù)文件“resample.dta”來介紹重復抽樣部分選項的應用。該文件也是人為構(gòu)造的，主要變量包括：group=分組變量（取值為A、B、C、D、E），strid=分層變量（取值為1或2，表明被調(diào)查者屬于哪種類型），x=觀測值。利用這些數(shù)據(jù)，我們來講解重復抽樣的基本命令和相關(guān)選項的應用。實驗操作指導實驗操作指導1 重復抽樣的基本命令重復

13、抽樣的基本命令此外，默認情況下，命令bsample會將內(nèi)存中的數(shù)據(jù)替換為抽樣的觀測值，但設定選項weight()會將抽取的樣本頻數(shù)存放在變量varname中，也就是說，這時只有varname的值改變，原數(shù)據(jù)不會改變。但選項weight()和選項idcluster()不能同時設定。另外，在bsample命令之后，選項weight(varname)中的varname可以用在Stata的其他命令中作為fweight（如果該命令允許設定fweight）。下面，我們結(jié)合“gender.dta”和“resample.dta”數(shù)據(jù)，對各選項做進一步的說明。2 簡單隨機抽樣簡單隨機抽樣對于“gender.dt

14、a”的數(shù)據(jù)，假設我們要采用簡單隨機抽樣法抽取300個樣本，可輸入命令：bsample 300下面，我們來看一下數(shù)據(jù)文件現(xiàn)在的樣本容量。輸入命令：count如果要只對男性進行簡單隨機抽樣，我們可以利用條件語句。輸入命令：bsample 200 if gender=0這里，我們對gender取值為0（表示male）的觀測值抽取了200個隨機樣本。需要注意的是，條件語句中，gender后為兩個等號。另外，需要說明的一點是，讀者若要復制這個操作，則需要在前面簡單隨機抽樣后重新打開一次數(shù)據(jù)；否則，這里的分層抽樣就會在原來抽樣結(jié)果的基礎(chǔ)上進行。后面的操作也同理。下面，我們看一下抽取的樣本的情況。輸入命令：

15、tab gender這里，命令tab表示顯示變量gender的頻數(shù)分布表。3 分層抽樣分層抽樣如果我們要令樣本中包括100個女性和100個男性，可以采取分層抽樣。輸入命令：bsample 100, strata(gender)這里，100表示每一層的樣本容量都是100，選項strata(gender)表示按變量gender的不同取值來分層。下面，我們看一下分層抽樣的樣本情況。輸入命令：tab gender此外，注意到原來的數(shù)據(jù)中約3000個為女性，2000個為男性。這樣，我們考慮按照原數(shù)據(jù)的男女比例進行抽樣；也就是說，抽取300個女性，200個男性。要做到這一點，要先生成一個新變量，令gend

16、er為female時新變量的值為300，gender為male時新變量的值為200。輸入命令：gen st = cond(gender,300,200)這里，我們將新變量命名為st。對于條件函數(shù)cond(x, a, b)，其含義為：如果x為真（或取值不是0），則返回a的值；如果x為假（或取值為0），則返回b的值。對于本例，如果變量“gender”的值為female（1），則令變量st的值為300；如果變量“gender”的值為male（0），則令變量st的值為200。下面，我們就可以利用變量st作為“exp”來進行分層抽樣。輸入命令：bsample st, strata(gender)這句命令

17、的含義為，按變量gender的值進行分層抽樣，且變量gender各個取值對應的樣本容量為st的值。也就是說，對于gender取值為female的觀測值，對應的抽樣樣本容量為300；對于gender取值為male的觀測值，對應的抽樣樣本容量為200。下面，我們看一下樣本的情況。輸入命令：tab gender4 擴展樣本容量擴展樣本容量首先要說明的一點是，如果不設定要抽取的樣本容量，則樣本容量會與原數(shù)據(jù)相同。而如果同時又設定了選項strata()，則各層中抽取的樣本個數(shù)會與原數(shù)據(jù)相同。此外，如果要抽取的樣本容量超過現(xiàn)有的觀測值個數(shù)，我們可以先對樣本容量進行擴展。輸入命令：expand 2則我們將樣

18、本容量擴展為原來的2倍（當然，我們這里也可以選擇將樣本容量擴為原來的3倍等；但對于本實驗，2倍就足夠了）。這里，需要說明的是，該命令后也可以加條件語句，選擇對分層變量的某個取值擴展樣本容量。接著，我們重新進行分層抽樣，且不設定樣本容量“exp”。輸入命令：bsample, strata(gender)這里，選項strata(gender)表示按變量gender進行分層。下面，我們來看一下抽取的樣本的頻數(shù)情況。輸入命令：tab gender5 生成頻率權(quán)重生成頻率權(quán)重在前面的幾個例子中，原樣本都被新抽取的樣本所替代。但我們也可以使用選項weight()，將抽取的樣本頻數(shù)放到一個新變量中去，并保留

19、原來的樣本。例如，我們還是要抽取200個男性樣本，并將頻數(shù)放到一個新變量中。為此，我們先生成一個變量，好用于存放頻數(shù)。輸入命令：gen fweight=.這里，我們生成一個新變量，將其命名為fweight，并令其初始值缺失（即為“.”）。下面，我們進行抽樣。輸入命令：set seed 1111bsample 200 if gender=0, weight(fweight)這里，第一步命令為設定種子，這是為了復制結(jié)果的方便。第二步，我們對gender取值為0的觀測值抽取了200個隨機樣本；選項weight(fweight)表明，我們將抽取的樣本的頻數(shù)存放在變量fweight中。下面，我們看一下樣

20、本的頻數(shù)分布情況。輸入命令：tab fweight gender這里，我們生成變量fweight和gender的二維交互表。6 分組抽樣分組抽樣對于分組數(shù)據(jù)，如果設定選項cluster()，則Stata會按組來進行抽樣。下面，我們利用“resample.dta”的數(shù)據(jù)，重點來講一下要抽取的組多于原來的組時該如何處理?！皉esample.dta”的數(shù)據(jù)只有5組，而我們想要抽取7組的樣本。顯然，我們應該先對樣本進行擴容。但需要注意的是，命令expand在這里是不恰當?shù)?。這是因為，expand命令不能為新生成的組設置其自身的標記；這樣，Stata會認為數(shù)據(jù)還是只有5組，從而，抽取7組樣本仍不可能。要

21、解決這個問題，我們可以使用分組數(shù)據(jù)的復制命令。該命令的基本格式為：expandcl = exp if in, cluster(varlist) generate(newvar)其中，expandcl為復制分組數(shù)據(jù)的基本命令（注意最后一個字母是小寫字母el，而非數(shù)字1），表達式exp用于指定復制為原來的幾份，if代表條件語句，in代表范圍語句。選項cluster()用于指定分組變量，選項generate()用于生成新的變量來保存各個組的識別標志。在對樣本進行擴容之前，我們先看一下原來的數(shù)據(jù)情況，方便與后面進行對比。輸入命令：tabstat x, stat(n mean) by(group)這里，

22、tabstat表示用表格顯示數(shù)據(jù)的描述統(tǒng)計量，x表示我們要得到變量x的統(tǒng)計量。選項stat(n mean)表示我們要獲得x的觀測值個數(shù)和均值，by(group)表示我們按組來獲得相應的統(tǒng)計量。下面，我們對“resample.dta”的數(shù)據(jù)進行擴容。輸入命令：expandcl 2, generate(newg) cluster(group)這里，我們將數(shù)據(jù)復制為原來的兩份。選項cluster(group)指定原來的分組變量為group，選項generate(newg)表示生成新變量newg用來保存復制后的各個組的識別標志。下面，我們生成一個新變量用于保存抽樣頻數(shù)，并進行分組抽樣。輸入命令：gen

23、 fweight= .set seed 1111bsample 7, cluster(newg) weight(fweight)這里，第一步生成用于保存頻數(shù)的變量fweight，并將其值設定為缺失值。第二步設定種子為1111，以方便結(jié)果的復制。第三步進行抽樣，以newg為分組變量，抽取了7組，并將抽取的頻數(shù)保存在變量fweight中。下面，我們看一下抽取的各組的頻數(shù)情況。輸入命令：tab fweight group7 分層分組抽樣分層分組抽樣對于“resample.dta”的數(shù)據(jù)，我們知道變量strid分為2層，且每層中有5組。但兩層之間的分組變量不能被唯一地識別。要在每一層中獲得7組樣本，我

24、們可以先使用expandcl命令進行擴展，并生成新變量唯一地識別每層的各個組。為了與后面的結(jié)果進行比較，我們先看一下各組各層的頻數(shù)分布情況。輸入命令：tab group strid下面，我們對數(shù)據(jù)進行擴展。輸入命令：expandcl 2, generate(newg) cluster(group strid)注意，我們這里使用了選項cluster(group strid)，表明為每一個group和strid的組合設置一個組別的標識。下面，我們來進行分層分組抽樣。輸入命令：gen fweight= .set seed 1111bsample 7, cluster(newg) strata(str

25、id) weight(fweight)其中，第三步的選項設定分組變量為newg，分層變量為strid。這樣，我們會在每一層中抽取7個組。下面，我們看一個各層分別抽到了哪些組。輸入命令：by strid: tabulate fweight group其中，“tabulate fweight group”表示繪制fweight和group的二維表，前綴“by strid”表示為各層分別繪制表格。實驗實驗15-4： 自助法自助法實驗基本原理實驗基本原理實驗內(nèi)容及數(shù)據(jù)來源實驗內(nèi)容及數(shù)據(jù)來源本書附帶光盤data文件夾下的“usaauto.dta”工作文件給出了美國的汽車業(yè)相關(guān)數(shù)據(jù)。利用該數(shù)據(jù)，我們來分析

26、各因素對每加侖油所行駛的里程數(shù)的影響，并講解自助估計的基本操作以及多種置信區(qū)間的獲得等內(nèi)容。實驗操作指導實驗操作指導1 自助估計的基本操作自助估計的基本操作進行自助法估計的基本命令為：bootstrap exp_list , options: command其中，bootstrap是自助法估計的基本命令，exp_list表示要保存的結(jié)果的表達式，command指定一次抽樣中執(zhí)行的命令，options代表其他選項。需要說明的一點是，表達式exp_list的格式與蒙特卡洛模擬處相同；如果命令command每次改變的是估計的模型的系數(shù)，則exp_list為可選項且默認值為_b。此外，表15.3列出了

27、主要的options選項。下面，我們利用“usaauto.dta”的數(shù)據(jù)來講解自助法估計的相關(guān)操作。首先，我們進行OLS回歸，獲取回歸系數(shù)和標準差，以方便和后面的對比。輸入命令：regress mpg weight gear foreign下面，我們采取自助法進行回歸。輸入命令：bootstrap, reps(100) seed(123): regress mpg weight gear foreign這里，reps(100)表示進行100次重復抽樣， seed(123)表示設置種子為123，冒號后為回歸的命令。如果假定“usaauto.dta”的數(shù)據(jù)按照變量rep78（汽車在78年的修理次數(shù)

28、）來聚類，我們可以使用cluster()選項來進行適當?shù)男拚?。另外，我們這里想獲得變量weight和gear_ratio的系數(shù)之差。這樣，輸入命令：keep if rep78 .bootstrap diff=(_bweight - _bgear), seed(123) cluster(rep78): regress mpg weight gear foreign這里，第一步是刪掉變量rep78的缺失值；這是因為命令bootstrap不允許選項cluster()中的變量存在缺失值。第二步進行自助法回歸，我們將變量weight和gear_ratio的系數(shù)之差命名為diff，并設置種子為123，聚類

29、變量為rep78。對于我們這個回歸，下面的命令也可以得到相同的結(jié)果：bootstrap diff=(_bweight - _bgear), seed(123): regress mpg weight gear foreign, vce(cluster rep78)這里，我們在回歸中設置了選項vce(cluster rep78)，表明按變量rep78聚類；這樣，bootstrap命令就會按照變量rep78的類別進行抽樣。 2 多種置信區(qū)間的獲得多種置信區(qū)間的獲得對于前面的自助回歸，我們要獲得多種置信區(qū)間。這樣，先輸入如下命令：qui bootstrap, reps(1000) saving(bs) bca seed(123): regress mpg weight gear foreign這里，選項reps()表明抽樣次數(shù)為1000；如我們前面所講，要比較可靠地獲得置信區(qū)間的估計值，抽樣次數(shù)要大一些。saving(bs)表示將自助估計的結(jié)果保存到文件“bs.dta”中。qui表明不顯示回歸結(jié)果。下面，我們來獲取多種置信區(qū)間。輸入命令：estat bootstrap, all這