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