Tối ưu hóa mã lệnh stata: củng cố các lệnh ước lượng với thư viện mata và lệnh bao bọc

Khi phát triển các công cụ ước lượng tùy chỉnh trong Stata, việc quản lý mã lệnh là một thách thức không nhỏ. Đặc biệt, nếu bạn xây dựng nhiều lệnh tương tự nhau, rủi ro trùng lặp mã lệnh sẽ rất cao. Mã trùng lặp không chỉ làm tăng công sức bảo trì mà còn dễ gây ra lỗi nếu bạn quên cập nhật ở tất cả các vị trí. Bài viết này sẽ giới thiệu cách củng cố mã lệnh Stata bằng cách sử dụng thư viện Mata để chia sẻ các hàm Mata và lệnh bao bọc (wrapper command) để chia sẻ mã ado-command. Phương pháp này giúp mã lệnh của bạn gọn gàng, dễ quản lý và mạnh mẽ hơn.

Lệnh Ado cho Mô Hình Kỳ Vọng Có Điều Kiện Hàm Mũ (ECM) và Probit (PCM)

Chúng ta hãy xem xét hai lệnh ado-command ước lượng tham số của mô hình kỳ vọng có điều kiện hàm mũ (ECM) và mô hình kỳ vọng có điều kiện Probit (PCM) bằng phương pháp bình phương nhỏ nhất phi tuyến tính (NLS). Ban đầu, các lệnh này được viết riêng lẻ.

Lệnh Mynlexp1.ado

Lệnh `mynlexp1` thực hiện ước lượng NLS cho các tham số của mô hình ECM.

1*! version 1.0.0  09May2016
2program define mynlexp1, eclass sortpreserve
3    version 14.1
4    syntax varlist [if] [in] [,  noCONStant ]
5    marksample touse
6    gettoken depvar indeps : varlist
7    tempname b V N rank
8    mata: mywork("`depvar'", "`indeps'", "`touse'", "`constant'", ///
9       "`b'", "`V'", "`N'", "`rank'" )
10    if "`constant'" == "" {
11        local indeps "`indeps' _cons"
12    }
13    matrix colnames `b' = `indeps'
14    matrix colnames `V' = `indeps'
15    matrix rownames `V' = `indeps'
16    ereturn post `b' `V', esample(`touse')
17    ereturn scalar N       = `N'
18    ereturn scalar rank    = `rank'
19    ereturn local  cmd     "mynlexp"
20    ereturn display

mata:

void MYNLExp(real scalar todo, real vector b, ///

real vector y, real matrix X, ///

val, grad, hess)

{

real vector r, f, xb

real matrix df

xb = X*b'

f = exp(xb)

r = y - f

val = -(r:^2)

df = f:*X

if (todo>=1) {

grad = r:*df

}

if (todo==2) {

hess = -1*quadcross(df, df)

}

void mywork( string scalar depvar, string scalar indeps,

string scalar touse, string scalar constant,

string scalar bname, string scalar Vname,

string scalar nname, string scalar rname)

{

real vector y, b

real matrix X, V

real scalar n, p, ssr

transmorphic S

y = st_data(., depvar, touse)

n = rows(y)

X = st_data(., indeps, touse)

if (constant == "") {

X = X,J(n, 1, 1)

}

p = cols(X)

S = optimize_init()

optimize_init_argument(S, 1, y)

optimize_init_argument(S, 2, X)

optimize_init_evaluator(S, &MYNLExp())

optimize_init_params(S, J(1, p, .01))

optimize_init_evaluatortype(S, "gf2")

optimize_init_conv_vtol(S, 1e-10)

b = optimize(S)

V = invsym(-1*optimize_result_Hessian(S))

ssr = (-1/(n-p))*optimize_result_value(S)

V = ssr*V

st_matrix(bname, b)

st_matrix(Vname, V)

st_numscalar(nname, n)

st_numscalar(rname, p)

}

end

Mã lệnh `mynlexp1` bao gồm phần ado-command và các hàm Mata định nghĩa trong cùng tệp. Hàm `MYNLExp` là hàm đánh giá được sử dụng bởi `optimize()` trong hàm `mywork`.

Lệnh Mynlprobit1.ado

Tương tự, lệnh `mynlprobit1` thực hiện ước lượng NLS cho các tham số của mô hình PCM.

1*! version 1.0.0  09May2016
2program define mynlprobit1, eclass sortpreserve
3    version 14.1
4    syntax varlist [if] [in] [,  noCONStant ]
5    marksample touse
6    gettoken depvar indeps : varlist
7    tempname b V N rank
8    mata: mywork("`depvar'", "`indeps'", "`touse'", "`constant'", ///
9       "`b'", "`V'", "`N'", "`rank'" )
10    if "`constant'" == "" {
11        local indeps "`indeps' _cons"
12    }
13    matrix colnames `b' = `indeps'
14    matrix colnames `V' = `indeps'
15    matrix rownames `V' = `indeps'
16    ereturn post `b' `V', esample(`touse')
17    ereturn scalar N       = `N'
18    ereturn scalar rank    = `rank'
19    ereturn local  cmd     "mynlexp"
20    ereturn display

mata:

void MYNLProbit(real scalar todo, real vector b, ///

real vector y, real matrix X, ///

val, grad, hess)

{

real vector r, f, xb

real matrix df

xb = X*b'

f = normal(xb)

r = y - f

val = -(r:^2)

df = normalden(xb):*X

if (todo>=1) {

grad = r:*df

}

if (todo==2) {

hess = -1*quadcross(df, df)

}

void mywork( string scalar depvar, string scalar indeps,

string scalar touse, string scalar constant,

string scalar bname, string scalar Vname,

string scalar nname, string scalar rname)

{

real vector y, b

real matrix X, V

real scalar n, p, ssr

transmorphic S

y = st_data(., depvar, touse)

n = rows(y)

X = st_data(., indeps, touse)

if (constant == "") {

X = X,J(n, 1, 1)

}

p = cols(X)

S = optimize_init()

optimize_init_argument(S, 1, y)

optimize_init_argument(S, 2, X)

optimize_init_evaluator(S, &MYNLProbit())

optimize_init_params(S, J(1, p, .01))

optimize_init_evaluatortype(S, "gf2")

optimize_init_conv_vtol(S, 1e-10)

b = optimize(S)

V = invsym(-1*optimize_result_Hessian(S))

ssr = (-1/(n-p))*optimize_result_value(S)

V = ssr*V

st_matrix(bname, b)

st_matrix(Vname, V)

st_numscalar(nname, n)

st_numscalar(rname, p)

}

end

Như bạn có thể thấy, mã lệnh của `mynlprobit1` gần như y hệt `mynlexp1`. Sự khác biệt chủ yếu nằm ở hàm đánh giá `MYNLProbit` thay vì `MYNLExp`. Mã lệnh trùng lặp rất nguy hiểm: mỗi khi bạn muốn thêm một tính năng hoặc sửa lỗi, bạn phải thực hiện nó hai lần. Chúng ta nên tránh mã trùng lặp bằng cách viết lại các lệnh này để có một cơ sở mã Mata duy nhất và sau đó là một cơ sở mã ado-command duy nhất.

Thư Viện Mã Mata

Các hàm `mywork()` được sử dụng trong `mynlexp1` và `mynlprobit1` chỉ khác nhau ở hàm đánh giá mà chúng gọi. Các hàm Mata được định nghĩa ở cuối tệp ado-file chỉ có thể được sử dụng trong tệp ado-file đó. Để có thể chia sẻ hàm `mywork()` cho cả `mynlexp1` và `mynlprobit1`, chúng ta cần một tệp chứa các hàm Mata đã biên dịch có thể được gọi từ bất kỳ hàm Mata nào khác hoặc từ bất kỳ ado-file hay do-file nào. Loại tệp này được gọi là thư viện.

Chúng ta sẽ sử dụng tệp `mynllib.mata` để tạo thư viện `lmynllib.mlib` chứa các hàm Mata đã biên dịch `MYNLWork()`, `MYNLProbit()`, `MYNLExp()`.

1mata:
2mata clear
3void MYNLExp(real scalar todo, real vector b,   ///
4        real vector y, real matrix X,           ///
5        val, grad, hess)
6{
7        real vector  r, f
8        real matrix  df
9        f   = exp(X*b')
10        r   = y - f
11        val = -(r:^2)
12        df  = f:*X
13        if (todo>=1) {
14                grad = r:*df
15        }
16        if (todo==2) {
17                hess = -1*quadcross(df, df)
18        }
19}
20void MYNLProbit(real scalar todo, real vector b,  ///
21        real vector y, real matrix X,             ///
22        val, grad, hess)
23{
24        real vector  r, f, xb
25        real matrix  df
26        xb  = X*b'
27        f   = normal(xb)
28        r   = y - f
29        val = -(r:^2)
30        df  = normalden(xb):*X
31        if (todo>=1) {
32                grad = r:*df
33        }
34        if (todo==2) {
35                hess = -1*quadcross(df, df)
36        }
37}
38void MYNLWork( string scalar depvar,  string scalar indeps,
39             string scalar touse,   string scalar constant,
40             string scalar bname,   string scalar Vname,
41             string scalar nname,   string scalar rname,
42             string scalar model)
43{
44    real vector       y, b
45    real matrix       X, V
46    real scalar       n, p, ssr
47    string scalar     emsg
48    pointer(function) f
49    transmorphic      S
50    if (model=="expm") {
51        f = &MYNLExp()
52    }
53    else if (model=="probit") {
54        f = &MYNLProbit()
55    }
56    else {
57        emsg = "{red}model " + model + " invalid\n"
58        printf(emsg)
59        exit(error(498))
60    }
61    y = st_data(., depvar, touse)
62    n = rows(y)
63    X = st_data(., indeps, touse)
64    if (constant == "") {
65        X = X,J(n, 1, 1)
66    }
67    p = cols(X)
68    S = optimize_init()
69    optimize_init_argument(S, 1, y)
70    optimize_init_argument(S, 2, X)
71    optimize_init_evaluator(S, f)
72    optimize_init_params(S, J(1, p, .01))
73    optimize_init_evaluatortype(S, "gf2")
74    optimize_init_conv_vtol(S, 1e-10)
75    b   = optimize(S)
76    V   = invsym(-1*optimize_result_Hessian(S))
77    ssr = (-1/(n-p))*optimize_result_value(S)
78    V   = ssr*V
79    st_matrix(bname, b)
80    st_matrix(Vname, V)
81    st_numscalar(nname, n)
82    st_numscalar(rname, p)
83}
84mata mlib create lmynllib, replace
85mata mlib add    lmynllib MYNLWork() MYNLProbit() MYNLExp()

Hàm `MYNLWork()` chấp nhận một đối số thứ chín là `model`, một biến vô hướng kiểu string. Nếu `model` là "expm", nó sẽ lưu địa chỉ của hàm `MYNLExp()` vào `f` (một con trỏ hàm). Nếu `model` là "probit", nó sẽ lưu địa chỉ của hàm `MYNLProbit()` vào `f`. Con trỏ hàm cho phép `MYNLWork()` gọi đúng hàm đánh giá tương ứng với mô hình được chỉ định.

Lưu ý rằng các hàm trong thư viện Mata được đặt tên bắt đầu bằng chữ hoa (ví dụ: `MYNLExp`). Điều này giúp tránh trùng lặp tên hàm với các thư viện khác.

Trong trường hợp này, việc gộp các hàm đánh giá `MYNLExp()` và `MYNLProbit()` vào một hàm duy nhất với một đối số bổ sung có thể làm chậm quá trình tính toán, vì hàm đánh giá được gọi rất nhiều lần bởi `optimize()`. Do đó, chúng ta chấp nhận một chút trùng lặp mã ở đây để ưu tiên tốc độ.

Dòng cuối cùng của mã lệnh tạo thư viện `lmynllib.mlib` và thêm các hàm `MYNLWork()`, `MYNLProbit()`, `MYNLExp()` đã biên dịch vào đó.

Ví Dụ 1: Tạo Thư Viện Mata

1. program drop _all
2. mata: mata clear
3. quietly do mynllib.mata
4. mata: mata mlib index
5.mlib libraries to be searched are now
6    lmatabase;lmatapss;lmataado;lmatapostest;lmatafc;lmatasem;lmatapath;
7> lmatamcmc;lmatagsem;lmataopt;lmynllib;lfreduse;lpoparms;lspmat

Sau khi tạo thư viện và thêm nó vào danh sách các thư viện được biết đến của Mata bằng lệnh `mata mlib index`, chúng ta có thể sử dụng các hàm được định nghĩa trong đó.

Lệnh Mynlexp2.ado và Mynlprobit2.ado

Giờ đây, các lệnh ado-command có thể gọi trực tiếp hàm `MYNLWork()` từ thư viện.

1*! version 2.0.0  10May2016
2program define mynlexp2, eclass sortpreserve
3    version 14.1
4    syntax varlist [if] [in] [,  noCONStant ]
5    marksample touse
6    gettoken depvar indeps : varlist
7    tempname b V N rank
8    mata: MYNLWork("`depvar'", "`indeps'", "`touse'", "`constant'", ///
9       "`b'", "`V'", "`N'", "`rank'", "expm" )
10    if "`constant'" == "" {
11        local indeps "`indeps' _cons"
12    }
13    matrix colnames `b' = `indeps'
14    matrix colnames `V' = `indeps'
15    matrix rownames `V' = `indeps'
16    ereturn post `b' `V', esample(`touse')
17    ereturn scalar N       = `N'
18    ereturn scalar rank    = `rank'
19    ereturn local  cmd     "mynlexp"
20    ereturn display

1*! version 2.0.0  10May2016
2program define mynlprobit2, eclass sortpreserve
3    version 14.1
4    syntax varlist [if] [in] [,  noCONStant ]
5    marksample touse
6    gettoken depvar indeps : varlist
7    tempname b V N rank
8    mata: MYNLWork("`depvar'", "`indeps'", "`touse'", "`constant'", ///
9       "`b'", "`V'", "`N'", "`rank'", "probit" )
10    if "`constant'" == "" {
11        local indeps "`indeps' _cons"
12    }
13    matrix colnames `b' = `indeps'
14    matrix colnames `V' = `indeps'
15    matrix rownames `V' = `indeps'
16    ereturn post `b' `V', esample(`touse')
17    ereturn scalar N       = `N'
18    ereturn scalar rank    = `rank'
19    ereturn local  cmd     "mynlprobit"
20    ereturn display

Các lệnh `mynlexp2` và `mynlprobit2` đều gọi hàm `MYNLWork()` từ thư viện, chỉ khác nhau ở đối số `model` truyền vào ("expm" hoặc "probit"). Tuy nhiên, phần ado-code của chúng vẫn còn nhiều trùng lặp.

Viết Lệnh Ado Làm Việc (Work Ado-command)

Để củng cố phần ado-code, chúng ta sẽ tạo một lệnh ado-command duy nhất để thực hiện công việc chính, và các lệnh ban đầu sẽ trở thành các lệnh bao bọc.

Lệnh Mynlexp3.ado và Mynlprobit3.ado (Lệnh Bao Bọc)

1*! version 3.0.0  11May2016
2program define mynlexp3
3    version 14.1
4    mynlwork expm `0'

1*! version 3.0.0  11May2016
2program define mynlprobit3, eclass sortpreserve
3    version 14.1
4    mynlwork probit `0'

Các lệnh `mynlexp3` và `mynlprobit3` chỉ đơn giản gọi lệnh `mynlwork` và truyền đối số `model` ("expm" hoặc "probit") cùng với các đối số mà người dùng đã nhập (`0` là macro cục bộ chứa tất cả các đối số của người dùng).

Lệnh Mynlwork.ado (Lệnh Làm Việc)

1*! version 1.0.0  11May2016
2program define mynlwork, eclass sortpreserve
3    version 14.1
4    gettoken model 0 : 0
5    if "`model'" == "expm" {
6        local cname "mynlexp"
7    }
8    else if "`model'" == "probit" {
9        local cname "mynlprobit"
10    }
11    else {
12        diplay "{red}model `model' invalid"
13        exit 498
14    }
15    syntax varlist [if] [in] [,  noCONStant ]
16    marksample touse
17    gettoken depvar indeps : varlist
18    tempname b V N rank
19    mata: MYNLWork("`depvar'", "`indeps'", "`touse'", "`constant'", ///
20       "`b'", "`V'", "`N'", "`N'", "`model'" )
21    if "`constant'" == "" {
22        local indeps "`indeps' _cons"
23    }
24    matrix colnames `b' = `indeps'
25    matrix colnames `V' = `indeps'
26    matrix rownames `V' = `indeps'
27    ereturn post `b' `V', esample(`touse')
28    ereturn scalar N       = `N'
29    ereturn scalar rank    = `rank'
30    ereturn local  cmd     "`cname'"
31    ereturn display

Lệnh `mynlwork` sử dụng `gettoken` để lấy đối số `model` đầu tiên được truyền vào. Dựa vào giá trị của `model`, nó sẽ xác định tên lệnh gọi (`cname`) và sau đó truyền `model` này cho hàm `MYNLWork()` trong Mata. `mynlwork` chứa tất cả logic xử lý chung, bao gồm phân tích cú pháp, gọi hàm Mata, và lưu kết quả.

Ví Dụ 2: Kết Quả Mynlexp3

1. mynlexp3 accidents cvalue tickets
2Iteration 0:   f(p) =  -2530.846
3Iteration 1:   f(p) = -1116.4901
4Iteration 2:   f(p) = -248.56923
5Iteration 3:   f(p) = -225.91644
6Iteration 4:   f(p) = -225.89573
7Iteration 5:   f(p) = -225.89573
8Iteration 6:   f(p) = -225.89573
9----------------------------------------------------------------------------
10           |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
11-----------+----------------------------------------------------------------
12    cvalue |   .1759434   .0323911     5.43   0.000     .1124581    .2394287
13   tickets |   1.447672   .0333599    43.40   0.000     1.382287    1.513056
14     _cons |  -7.660608   .2355725   -32.52   0.000    -8.122322   -7.198894
15----------------------------------------------------------------------------

Ví Dụ 3: Kết Quả Mynlprobit3

1. mynlprobit3 hadaccident cvalue tickets
2Iteration 0:   f(p) = -132.90997
3Iteration 1:   f(p) = -16.917203
4Iteration 2:   f(p) = -10.995001
5Iteration 3:   f(p) = -10.437501
6Iteration 4:   f(p) = -10.427738
7Iteration 5:   f(p) = -10.427156
8Iteration 6:   f(p) = -10.427123
9Iteration 7:   f(p) = -10.427121
10Iteration 8:   f(p) =  -10.42712
11Iteration 9:   f(p) =  -10.42712
12Iteration 10:  f(p) =  -10.42712
13----------------------------------------------------------------------------
14           |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
15-----------+----------------------------------------------------------------
16    cvalue |   .3616322   .0918214     3.94   0.000     .1816656    .5415988
17   tickets |   2.177509   .1974173    11.03   0.000     1.790578     2.56444
18     _cons |  -10.95166    1.05565   -10.37   0.000    -13.02069   -8.882622
19----------------------------------------------------------------------------

Các ví dụ trên chứng minh rằng `mynlexp3` và `mynlprobit3` tạo ra kết quả chính xác, tương ứng với kết quả của lệnh `nl` gốc.

✨ **Giá Trị Đắt Giá**

Việc chia sẻ mã lệnh luôn tốt hơn việc lặp lại mã. Mã dùng chung chỉ cần được thay đổi ở một nơi duy nhất để thêm tính năng hoặc sửa lỗi, trong khi mã lặp lại phải được thay đổi ở mọi nơi. Việc sử dụng thư viện Mata để chia sẻ các hàm Mata giữa các ado-command và sử dụng các lệnh bao bọc để chia sẻ ado-code là những kỹ thuật hiệu quả giúp tối ưu hóa quá trình phát triển, giảm thiểu lỗi và nâng cao khả năng bảo trì mã.

Câu Hỏi Tư Duy hoặc Bài Tập Ứng Dụng

Hãy thử tưởng tượng bạn cần thêm một tính năng mới vào các lệnh `mynlexp3` và `mynlprobit3`, ví dụ như hỗ trợ ước lượng hiệp phương sai của tham số (VCE) mạnh mẽ (robust VCE). Bạn sẽ cần thay đổi những tệp nào và ở những vị trí nào để thực hiện tính năng này theo cấu trúc mã đã được củng cố? So sánh với việc phải thực hiện thay đổi trên các tệp `mynlexp1` và `mynlprobit1` ban đầu.