Stata：蒙特卡洛模擬(Monte-Carlo-Simulation)沒那麼神秘
作者：侯新爍 湘潭大學 [編譯] (知乎 | 簡書 | 碼雲)
Stata連享會現場班課程分享 資料參考來源： The Stata Blog ? Monte Carlo simulations using Stata
- Stata連享會 Stata 現場班報名中 || 精彩推文2
通過計算機模擬，從總體抽取大量隨機樣本的計算方法統稱為「蒙特卡羅法（Monte Carlo Method，簡記為MC）」，這個名字來源於摩納哥的蒙特卡洛賭場（Monte Carlo Casino），這是最早使用這個方法的一位美國物理學家的叔叔常去的賭場。蒙特卡羅方法長用來確定統計量的小樣本性質。，我們知道許多統計量的精確分布沒有解析解，通常的處理方法就是應用大樣本理論，用漸近分布來近似真實分布，但現實中的樣本容量可能不夠大，此時蒙特卡羅法的價值就更加凸顯了。
在蒙特卡羅模擬中通常需要設計 loop 語句，可能使用的命令有如 forvalues、 foreach 、while ； 也可能需要生成一些隨機數據，比如 help random_number_functions 、 help simulate 等，讀者可自行查閱進行擴展。
在本次推文中，我們將使用服從卡方分布的隨機數據進行蒙特卡洛模擬，分別展示隨機數據的均值分布，多個隨機均值的生成與保存，以及樣本標準誤差的蒙特卡洛模擬過程。
1.概覽
一個估計量的 蒙特卡羅模擬（ Monte Carlo simulation ，MCS）是基於特定的 數據生成過程 （ data-generating process， DGP ）和樣本量，使用模擬的方法去近似估計量的抽樣分布的過程。使用MCS可以考察一個估計方法對特定的數據生成過程的估計表現。在這篇推文中，我們將展示如何對 Stata 中的估計量進行 MCS 模擬以及如何解釋結果。
大樣本理論告訴我們，當真正的 DGP 是一個來自自由度為 1 的 分布的隨機樣本時（標記為
），樣本平均值是均值的一個很好的估計量。但一位朋友聲稱這個估計量對於這個 DGP 不會很好，因為
分布會產生異常值。在這篇文章中，使用 MCS 方法，看是否大樣本理論可以很好地用於本 DGP 的 500 個觀察樣本情形。
2.第一個 MCS
我首先從展示如何從 分布中抽取一個樣本量為 500 的隨機樣本，以及如何估計平均值和平均值的標準誤差開始。
示例1：模擬數據的平均值
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
設置種子 為 12345 ，設置隨機數生成器的種子可允許樣本結果是可重複的。該隨機樣本的平均樣本平均估計值為 0.91 ，估計的標準誤差為 0.055。
若有很多估計，每個估計都來自一個獨立的隨機樣本，我就可以估計估計量的抽樣分布均值和標準差。為了獲得許多估計，我需要多次重複以下過程：
- 從 DGP 中抽取
- 計算估算值
- 存儲估算值
我需要知道如何存儲這些估計結果以繼續這個過程。我還需要知道如何多次重複該過程以及如何訪問 Stata 估計值，鑒於許多讀者已經熟悉這些主題，我將這些細節信息附在了附錄 I 和 II 中，我們將專註於如何存儲如此多的抽樣結果。
我想把這些估計結果保存在某處，它們將成為隨後可分析數據集的一部分。我將使用命令 postfile ， post 和 postclose 將估計值存儲在內存中，並在完成後將所有存儲的估計值寫入一個數據集。示例 2 對有三次抽樣過程的情形進行了說明。
示例2：三次抽獎的估計平均值
set seed 12345
postfile buffer mhat using mcs, replace
forvalues i=1/3 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
post buffer (_b[y])
}
postclose buffer
use mcs, clear
list
命令 postfile buffer mhat using mcs, replace 在內存中創建了一個空間叫做 buffer ，在其中我們可以存儲將最終寫入數據集的結果。mhat 為估計變數的名稱，數據集將保存為 mcs.dta 。關鍵詞 using 將新變數名稱與新數據集的名稱分開。同時，該例中，定義選項 replace 以將 msc.dta 的任何先前版本替換為此處創建的版本。如圖可見，我們可以得到 mhat 的三個結果。
具體而言，此例子使用如下循環，重複了計算均值的過程三次。（更多信息可參閱附錄 I ）。
forvalues i=1/3 {
}
此外，命令
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
分別執行了 刪除以前的數據，從 $chi^2(1)$ 分布中抽取大小為 500 的樣本，並估算平均值的操作。命令語句前的 quietly 表示抑制輸出。
命令 post buffer (_b[y]) 將當前抽取的估計平均值存儲在緩衝區中，以便下一次觀察 mhat 。命令 postclose buffer 將存儲在緩衝區中的內容寫入文件 mcs.dta 。命令 use mcs, clear; list 從內存中刪除最後一個 樣本，讀入
msc.dta 數據集，並列出數據集。
示例 3 是示例 2 的修改版本; 增加了抽樣數量並對結果進行了描述統計。
示例3：2 000 個估計平均值的均值
set seed 12345
postfile buffer mhat using mcs, replace
forvalues i=1/2000 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
post buffer (_b[y])
}
postclose buffer
use mcs, clear
summarize
估計值的平均值是估計量的抽樣分布均值的估計值，它接近 1.0 的真值。2 000 個估計的樣本標準差是估計量的採樣分布的標準差的估計量，並且接近真實值 sqrt {sigma^2/N} = sqrt{2/500} ，約為 0.0632 ），其中 （sigma^2）是 隨機變數的方差。
3.包含標準誤
基於均值統計量報告的標準誤差是估計量的抽樣分布的標準誤差估計值。若大樣本分布在做近似估計的抽樣分布時表現良好，估計的標準誤均值 應當接近這些平均估計值的樣本標準差。
為了將估計量的 標準差 與估計的 標準誤 的平均值進行比較，修改了示例3，以存儲標準誤差信息。
示例4：2 000個標準誤差的平均值
set seed 12345
postfile buffer mhat sehat using mcs, replace
forvalues i=1/2000 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
post buffer (_b[y]) (_se[y])
}
postclose buffer
use mcs, clear
summarize
其中，命令 postfile buffer mhat sehat using mcs, replace 在緩衝區為新變數 mhat 和 sehat 騰出空間，並且使用命令 post buffer (_b[y]) (_se[y]) 存儲的每個估計的均值 mhat ，和估計的標準誤 sehat 。（如示例 3 所示，命令 postclose buffer 將存儲在內存中的內容寫入新數據文件。）
可見，這 2,000 個估計的樣本標準偏差為 0.0625 ，接近 估計的標準誤均值 0.0630。
你可能會認為我應該寫「非常接近」，但 0.0625 距離 0.0630 究竟有多近？老實說，我們不知道這兩個數字是否彼此足夠接近，因為它們之間的距離並不能自動地告訴我們結果推論有多麼可靠。
4.估計拒絕率
在頻率統計中，如果 p 值低於某指定數值，我們拒絕零假設。若大樣本分布很好的近似於有限樣本分布，對真實的零假設檢驗，其拒絕率應接近指定的數值大小。
為了比較拒絕率在 5％ 時的情形，本推文修改了例子 4 來計算和存儲是否拒絕零假設的 Wald 檢驗 相關指標。 測試與真正的零假設。（機制討論見附錄 Ⅲ ）
示例5：估計拒絕率
set seed 12345
postfile buffer mhat sehat reject using mcs, replace
forvalues i=1/2000 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
quietly test _b[y]=1
local r = (r(p)<.05)
post buffer (_b[y]) (_se[y]) (`r)
}
postclose buffer
use mcs, clear
summarize
本例中，通過 local r = (r(p)<.05) 識別了每次估計中 p 是否顯著。可見， 0.048 的拒絕率非常接近 0.05 的真實大小。
5.總結
在這篇文章中，我展示了如何在 Stata 中執行估計量的 MCS 。我討論了使用 post 命令來存儲大量估計結果的機制，並討論了如何解釋大量估計的均值和標準誤均值內涵。同時，還介紹了使用估計的拒絕率來評價給定DGP和樣本量條件下的大樣本近似抽樣分布的可用性。
示例結果證明，抽樣均值表現為與大樣笨原理預測的估計量均值一致。當然，這一結論並不意味著我朋友對異常值的擔憂完全是錯誤的。其他對異常值更穩健的估計量可能具有更好的性質。以後的推文中將嘗試為這種權衡提供一些證據。
附錄
附錄I：同一過程的多次重複
本附錄簡要介紹了局域宏以及如何使用它們進行多次重複命令操作; 相關更多詳細信息，請參見 [P] macro and [P] forvalues for more details。(可以 help macro help forvalues)
我們可以在局域宏中存儲和訪問字元串信息。如將字元串 hello 存儲在本地宏命名值中。命令為 local value "hello" ；為了訪問存儲的信息，需要修飾局域宏的名稱。具體來說，在它前面加上單引號 (實際上是英文狀態下的 Tab 按鍵上面按鍵符號)，並用右單引號跟隨它。如訪問和展示局域宏 value 中的值，可利用 display 在結果窗中展示。也可以將數字存儲為字元串。
左邊 ` 右邊 合起來 `value
local value "hello"
display "`value"
local value "2.134"
display "`value"
要多次重複某些命令，需將它們放在一個 forvalues 循環中。例如，下面的代碼重複三次 display 命令。
forvalues i=1/3 {
display "i is now `i"
}
上面的示例說明 forvalues 定義了一個局域宏，它包含了值列表中的每個值。在前面的示例中，局域宏的名稱是 i ，指定的值是 1/3 = {1,2,3 } 。
附錄II：獲取估計結果
在 Stata 估計命令之後，您可以通過鍵入 _b[y] 來訪問名為 y 的參數的點估計值，也可以通過鍵入 _se[y] 來訪問其標準誤。下面的示例說明了此過程。
示例6：訪問估計值
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
display _b[y]
display _se[y]
附錄Ⅲ：通過檢驗計算獲取p值
本附錄解釋了創建指示 Wald 檢驗是否會在一個特定的數值上拒絕零假設的指標的機理，首先是要生成一些數據，並對真實的零假設進行 Wald 檢驗。
示例7：Wald測試結果
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
test _b[y]=1
檢驗的結果被保存在 r() 中。接下來，使用 return list 來查看它們，輸入 help return list 以獲取詳細信息。
示例8：通過測試存儲的結果
return list
可以看到，p 值存儲在 r(p) 中，接下來生成一個 p 值是否小於 0.05 的 0/1 指示器指標，存儲在局域宏 r 中（局域宏的介紹可見附錄II）。本文通過展示局域宏包含的信息是 0 來完成說明。
local r = (r(p)<.05)
display "`r"
代碼匯總
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
set seed 12345
postfile buffer mhat using mcs, replace
forvalues i=1/3 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
post buffer (_b[y])
}
postclose buffer
use mcs, clear
list
set seed 12345
postfile buffer mhat using mcs, replace
forvalues i=1/2000 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
post buffer (_b[y])
}
postclose buffer
use mcs, clear
summarize
set seed 12345
postfile buffer mhat sehat using mcs, replace
forvalues i=1/2000 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
post buffer (_b[y]) (_se[y])
}
postclose buffer
use mcs, clear
summarize
set seed 12345
postfile buffer mhat sehat reject using mcs, replace
forvalues i=1/2000 {
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
quietly test _b[y]=1
local r = (r(p)<.05)
post buffer (_b[y]) (_se[y]) (`r)
}
postclose buffer
use mcs, clear
summarize
qui local value "hello"
display "`value"
qui local value "2.134"
display "`value"
forvalues i=1/3 {
display "i is now `i"
}
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
display _b[y]
display _se[y]
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
test _b[y]=1
return list
local r = (r(p)<.05)
display "`r"
關於我們
- 【Stata 連享會(公眾號：StataChina)】由中山大學連玉君老師團隊創辦，旨在定期與大家分享 Stata 應用的各種經驗和技巧。
- 公眾號推文同步發佈於 CSDN-Stata連享會 、簡書-Stata連享會 和 知乎-連玉君Stata專欄。可以在上述網站中搜索關鍵詞
Stata或Stata連享會後關注我們。 - 點擊推文底部【閱讀原文】可以查看推文中的鏈接並下載相關資料。
- Stata連享會 精彩推文1 || 精彩推文2
聯繫我們
- 歡迎賜稿： 歡迎將您的文章或筆記投稿至
Stata連享會(公眾號: StataChina)，我們會保留您的署名；錄用稿件達五篇以上，即可免費獲得 Stata 現場培訓 (初級或高級選其一) 資格。 - 意見和資料： 歡迎您的寶貴意見，您也可以來信索取推文中提及的程序和數據。
- 招募英才： 歡迎加入我們的團隊，一起學習 Stata。合作編輯或撰寫稿件五篇以上，即可免費獲得 Stata 現場培訓 (初級或高級選其一) 資格。
- 聯繫郵件： [email protected]
往期精彩推文
- Stata連享會推文列表1
- Stata連享會推文列表2
- Stata連享會 精彩推文1 || 精彩推文2
推薦閱讀：
