set more off

import excel "D:\github\master-SEM\04-students\lab02-dataset-cps.xlsx", sheet("Sheet 1") firstrow  //导入数据


**** 查看数据情况 *******
describe
summarize

******任务1：处理定性变量******

**** 虚拟变量 gender****
tabulate female
generate gender_female = (female==1)
generate gender_male = (female==0)

**** 虚拟变量 married****
tabulate married
generate married_female = (married=="female")
generate married_male = (married=="male")
generate married_other = (married=="other")

**** 虚拟变量 union****
tabulate union
generate union_female = (union=="female")
generate union_male = (union=="male")
generate union_other = (union=="other")

**** 虚拟变量 formMarried****
tabulate formMarried
generate formMarried_female = (formMarried=="female")
generate formMarried_male = (formMarried=="male")
generate formMarried_other = (formMarried=="other")

**** 虚拟变量 hispanic****
tabulate hispanic
generate hispanic_y= (hispanic==1)
generate hispanic_n = (hispanic==0)

**** 虚拟变量 race****
tabulate hispanic
generate race_black = (race=="black")
generate race_ai = (race=="AI")
generate race_asian = (race=="asian")
generate race_mixed = (race=="mixed")
generate race_other = (race=="other")


******任务3-1：OLS快速回归分析******
regress lgwg edu exper exper2 gender_female union_female union_male married_female married_male formMarried_female formMarried_male hispanic_y race_black race_ai race_asian race_mixed
ereturn list

matrix b_quick = (e(b))'       // 获取斜率系数，转换为列向量
matrix list b_quick

scalar RSS_quick = e(rss)      // 获取残差平方和

matrix covvar_b_quick = e(V)       // 获取估计系数的方差协方差矩阵
matrix list covvar_b_quick

mata: st_matrix("s_b_quick", sqrt(diagonal(st_matrix("e(V)"))))    // 估计系数的标准差， 这里使用了mata语法
matrix list s_b_quick



******任务4-1:矩阵分析******
generate const = (female > -0.1) // 常数列，全部取值为1

// 以下采用mata语法体系
mata: st_view(y=., ., "lgwg")  // 构造y矩阵
mata: st_view(X=., ., ("edu", "exper", "exper2", "gender_female", "union_female", "union_male", "married_female", "married_male", "formMarried_female", "formMarried_male", "hispanic_y", "race_black", "race_ai", "race_asian", "race_mixed", "const"))  // 构造X矩阵

mata: b = (cholinv(cross(X,X)))*(cross(X,y))  // 计算斜率系数
mata: b

mata: N = rows(y)  // 样本数
mata: k = cols(X)  // 变量数
mata: printf("rows = %f k = %f", N, k)


******任务4-2:方差协方差矩阵******
mata: e = y - X*b         // 残差向量
mata: RSS = e'e           // 残差平方和
mata: sigma2 = RSS/(N-k)  // 回归误差方差 或 残差均方和MSS_RSS
mata: printf("RSS = %f sigma2 = %f", RSS, sigma2)

mata: covvar_b = cholinv(cross(X,X))*sigma2  // 估计系数的方差协方差矩阵
mata: s_b = sqrt(diagonal(covvar_b))         // 估计系数的标准差
mata: s_b


******任务5-1:异方差一致性标准误差******
// stata快速命令法, HC1使用', vce(robust)'命令； HC2使用', vce(hc2)'命令；HC3使用', vce(hc3)'命令。此外，HC0往往在实证中不使用，因为它没有进行自由度调整
regress lgwg edu exper exper2 gender_female union_female union_male married_female married_male formMarried_female formMarried_male hispanic_y race_black race_ai race_asian race_mixed, vce(robust)
mata: st_matrix("s_b_quick_robust", sqrt(diagonal(st_matrix("e(V)"))))    // 估计系数的稳健标准差， 这里使用了mata语法
matrix list s_b_quick_robust

// 手动编程计算法

mata: Xe = X:*e                // element product
mata: XXe2 = cross(Xe, Xe)     // (Xe)'*(Xe)
mata: hc0 =  cholinv(cross(X,X)) * XXe2 *  cholinv(cross(X,X))  // the HC0 formula
mata: hc1 = hc0 *N/(N-k)
mata: s_b_hc1 = sqrt(diagonal(hc1))  // 估计系数的HC1稳健标准差
mata: s_b_hc1


******以下为代码测试******


mata: XXi[1::3]
mata: (X*X')[1::5, 1::5]

mata: E = e[1::3]
mata: E
mata: A = (1,2 \ 3,4 \ 1, 2)
mata: A

mata: A:*E

mata: cross(A:*E, A:*E)

mata: A'
mata: A:*A
mata: rowsum(A:*A)
mata: B = (1 \ 2 \3)
mata: B
mata: B:*B

mata: (rowsum(A:*A)) :* (B:*B)