Hướng dẫn chạy mô phỏng trên Stata qua lệnh simulate
Guide to Monte Carlo simulation in Stata with simulate command
Giới thiệu mô phỏng Monte Carlo
Để thực hiện mô phỏng thì đầu tiên cần viết một chương trình mô phỏng. Một chương trình trên Stata là một đoạn lệnh được bao bởi khối lệnh program define và lệnh end.
Trong mỗi chương trình các tham số đầu vào cần được khai báo cụ thể ở phần cú pháp (câu lệnh syntax) của chương trình. Nội dung quan trọng nhất của mỗi chương trình chính là các tính toán như quá trình tạo dữ liệu (cở mẫu, tạo biến), các thao tác xử lý hay ước lượng cần thực hiện. Cuối cùng là giá trị trả về từ chương trình. Có 4 tùy chọn thiết lập kết quả mà chương trình trả về gồm rclass(), eclass(), sclass() và nclass().
- rclass khai báo kết quả chạy của chương trình sẽ được cập nhật trong in r() thông qua câu lệnh return. Nếu chương trình không khai báo rõ rclass thì các thay đổi sẽ không được cập nhật trong r().
- eclass khai báo kết quả chạy của chương trình sẽ được cập nhật trong in e() thông qua câu lệnh ereturn. Nếu chương trình không khai báo rõ rclass thì các thay đổi sẽ không được cập nhật trong e().
- sclass khai báo kết quả chạy của chương trình sẽ được cập nhật trong in s() thông qua câu lệnh sreturn. Nếu chương trình không khai báo rõ rclass thì các thay đổi sẽ không được cập nhật trong s().
- nclass (mặc định) khai báo kết quả chạy của chương trình không được cập nhật trong r(), e(), hoặc s().
Thực hiện mô phỏng Monte Carlo
Sau khi đã có chương trình để chạy mô phỏng, bước tiếp theo là thực hiện mô phỏng bằng cách sử dụng câu lệnh simulate. Cú pháp cơ bản của câu lệnh simulate:
simulate [exp_list] , reps(#) [options] : command
Câu lệnh simulate sẽ chạy command (lệnh hoặc chương trình mô phỏng đã tạo) cho # lần tái lập và thu nhận các kết quả trong các biểu thức gán [exp_list].
Ở đây,
- command xác định lệnh mà nó thực hiện một mô phỏng. Đây là đoạn chương trình mô phỏng mà bạn đã tạo trước đó.
- exp_list xác định biểu thức được tính toán từ quá trình chạy của command.
- Nếu không có biểu thức nào được thiết lập, thì exp_list sẽ sử dụng giá trị mặc định mà command thay đổi. Các kết quả trong e() hoặc r().
- Nếu command có thay đổi kết quả trong e(), mặc định là _b. Nếu command thay đổi kết quả trong r() thì tất cả các đại lượng vô hướng sẽ được tổng hợp trong r().
- Sử dụng tùy chọn nodots hoặc dots(#) để thiết lập hiển thị các dấu “.” cho quá trình tái lập. Tùy chọn dots(#) hiển thị mỗi “.” cho # lần lặp.
- seed(#) thiết lập một số ngẫu nhiên khi bắt đầu chạy
- saving(filename [, suboptions]) tạo một file Stata (.dta) bao gồm các giá trị của exp_list trong một biến ở mỗi lần lặp.
Ví dụ mô phỏng Monte Carlo
Ví dụ 1: Tạo bảng thống kê
Tạo bảng thống kê giá trị trung bình của hệ số và sai số chuẩn cho một mô hình hồi quy với kỹ thuật lấy mẫu lặp 1000 lần.
Bước 1: Xây dựng chương trình tạo dữ liệu
name: vd log: D:\1 - BAI DANG\2022\simulate\vd.smcl log type: smcl opened on: 10 Jul 2022, 16:08:25 . program drop _all . program vd1, rclass 1. drop _all 2. set obs 25 3. generate x = rnormal() 4. generate y = 3*x + 1 + rnormal() 5. regress y x 6. . return scalar betax = _b[x] 7. return scalar sex = _se[x] 8. end . vd1 Number of observations (_N) was 0, now 25. Source | SS df MS Number of obs = 25 -------------+---------------------------------- F(1, 23) = 121.73 Model | 120.253974 1 120.253974 Prob > F = 0.0000 Residual | 22.7202116 23 .987835287 R-squared = 0.8411 -------------+---------------------------------- Adj R-squared = 0.8342 Total | 142.974186 24 5.95725773 Root MSE = .9939 ------------------------------------------------------------------------------ y | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- x | 2.761113 .2502515 11.03 0.000 2.243428 3.278798 _cons | .8035015 .1987886 4.04 0.001 .3922759 1.214727 ------------------------------------------------------------------------------
Bước 2: Thực hiện mô phỏng
. *B2 . etime, start . simulate coef = r(betax) se = r(sex), /// > seed(1234) reps(1000) nodots: vd1 Command: vd1 coef: r(betax) se: r(sex) . etime Elapsed time is 6 seconds
Ở bước này, bạn có thể sử dụng lệnh etime trước và sau khối lệnh để biết thời gian cụ thể chạy mô phỏng.
Bước 3: Tạo bảng thống kê
. *B3 . summ coef se Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- coef | 1,000 2.995262 .2172831 2.236966 3.704498 se | 1,000 .2110237 .0452831 .1017304 .4259755
Ví dụ 2: Lập bảng GTTB và Phương sai
Bước 1: Viết chương trình
. program drop _all . program vd2, rclass 1. version 16.0 2. drop _all 3. . set obs 100 4. generate z = exp(rnormal()) 5. summarize z 6. . return scalar GTTB = r(mean) 7. return scalar PSAI = r(Var) 8. end . vd2 Number of observations (_N) was 0, now 100. Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- z | 100 1.645352 1.748945 .1138481 8.178824
Bước 2: thực hiện mô phỏng
. etime, start . simulate TB = r(GTTB) PS = r(PSAI), /// > seed(1234) reps(1000) nodots: vd2 Command: vd2 TB: r(GTTB) PS: r(PSAI) . etime Elapsed time is 0 seconds
Bước 3: hiển thị thông tin
. summ TB PS Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- TB | 1,000 1.630648 .2173062 1.106372 2.612052 PS | 1,000 4.60798 4.502166 .966087 70.5597
Kết thúc ghi nhật ký
. log close vd name: vd log: D:\1 - BAI DANG\2022\simulate\vd.smcl log type: smcl closed on: 10 Jul 2022, 16:08:31