可使用蒙特卡洛法进行模拟，所谓“蒙特卡罗法”(Monte Carlo Methods，MC)，是通过计算机模拟，从总体抽取大量随机样本的计算方法。
代码语言：javascript复制
* PRF 和 SRF：蒙特卡罗模拟
clear
set obs 30
set seed 10101
gen x = rnormal(3, 4)
gen e = rnormal(0, 9)
gen y = 1*x e //Data Generation Process
reg y x
tw function PRF = 1 2*x, range(-5 15) || ///
scatter y x || lfit y x, lp(dash)
begin{cases}
frac{delta}{delta hat alpha} sum_{i=1}^n e_i^2 = -2sum_{i=1}^n (y_i - hat alpha - hat beta x_i) = 0 \
frac{delta}{delta hat beta} sum_{i=1}^n e_i^2 = -2 sum_{i=1}^n (y_i - hat alpha - hat beta x_i) x_i = 0
end{cases}
消去方程左边的
-2
，可得：
begin{cases}
sum_{i=1}^n (y_i - hat alpha - hat beta x_i) = 0 \
sum_{i=1}^n (y_i - hat alpha - hat beta x_i)x_i = 0
end{cases}
对上式各项分别求和，移项可得：
begin{cases}
n hat alpha hat beta sum_{i=1}^n x_i = sum_{i=1}^n y_i \
hat alphasum_{i=1}^n x_i hat beta sum_{i=1}^n x_i^2 = sum_{i=1}^n x_iy_i
end{cases}
上式为“正规方程组”。记
bar y = frac{1}{n}sum_{i=1}^n y_i, bar x = frac{1}{n} sum_{i=1}^n x_i
，则：
hat alpha = bar y - hat beta bar x
将上式带入可得：
(bar y - hat betabar x)sum_{i=1}^n x_i hat beta sum_{i=1}^n x_i^2 = sum_{i=1}^n x_iy_i
合并同类项，移项可得：
hat beta (sum_{i=1}^n x_i^2 - bar x sum_{i=1}^n x_i) = sum_{i=1}^n x_iy_i - bar y sum_{i=1}^n x_i
使用关系式
sum_{i=1}^n x_i = n bar x
，求解
hat beta
：
hat beta = frac{sum_{i=1}^n x_iy_i - n bar{xy}}{sum_{i=1}^n x_i^2-n bar x^2}
更直观的写为离差形式：
hat beta = frac{sum_{i=1}^n(x_i-bar x)(y_i-bar y)}{sum_{i=1}^n(x_i-bar x )^2}
最后，求得：
hat alpha = bar y - hat beta bar x
hat beta = frac{sum_{i=1}^n(x_i-bar x)(y_i-bar y)}{sum_{i=1}^n(x_i-bar x )^2}
sum_{i=1}^n (y_i - bar y)^2 = sum_{i=1}^n (y_i - hat y_i hat y_i - bar y)^2 = sum_{i=1}^n (e_i hat y_i - bar y)^2 \
= sum_{i=1}^n e_i^2 sum_{i=1}^n(hat y_i - bar y)^2 2 color{red}{sum_{i=1}^ne_i(hat y_i - bar y)}
frac{d}{d hat beta} sum_{i=1}^n e_i^2 = -2 sum_{i=1}^n (y_i - hat beta x_i) x_i = 0
消去方程左边
-2
，可得：
sum_{i=1}^n (y_i - hat beta x_i) x_i = 0
求解
hat beta
可得：
hat beta = frac{sum_{i=1}^n x_iy_i}{sum_{i=1}^n x_i^2}
如果回归模型无常数项，则平方和分解公式不成立，不宜使用
R^2
来度量拟合优度。
记
e_i = y_i - hat beta x_i
，则正规方程可写为：
sum_{i=1}^n x_ie_i = 0
记拟合值为
hat y_i equiv hat beta x_i
，则容易证明残差仍与拟合值正交：
sum_{i=1}^n hat y_i e_i = sum_{i=1}^n hat beta x_ie_i = hat beta sum_{i=1} x_ie_i= hat beta cdot0 = 0
仍可利用 OLS 的正交性将
sum_{i=1}^n y_i^2
分解为：
sum_{i=1} y_i^2 = sum_{i=1}^n (hat y_i e_i)^2 = sum_{i=1}^n hat y_i^2 2underbrace{sum_{i=1}^n hat y_ie_i}_{=0} sum_{i=1}^n e_i^2 = sum_{i=1}^n hat y_i^2 sum_{i=1}^n e_i^2
sum_{i=1}^n hat y_i^2
为可由模型解释的部分，而
sum_{i=1}^n e_i^2
为模型不可解释的部分。
定义非中心
R^2
：
R_{uc}^2 = frac{sum_{i=1}^n hat y_i^2}{sum_{i=1}^n y_i^2}
如果无常数项，Stata 汇报的
R^2
正是
R_{uc}^2
。
1.5 Stata 命令及实例
1.5.1 简单线性回归的 Stata 实例
代码语言：javascript复制
use ${d}/grilic.dta, clear
/*
数据说明：此数据集包括 758 位美国年轻男子的教育投资回报率数据
*/
reg lnw s
reg lnw s, noc // 无常数项回归
1.5.2 PRF 与 SRF: 蒙特卡洛法
代码语言：javascript复制
/*
PRF: y_i = 1 2x_i epsilon_i quad (i=1, ..., 30)
解释变量：$x_i ~ N(3, 2^2)$，扰动项 $epsilon_i ~ N(0,3^2)$，样本容量为 30 。
*/
clear
set obs 30
set seed 10101
gen x = rnormal(3, 4)
gen e = rnormal(0, 9)
gen y = 1*x e
reg y x
tw function PRF = 1 2*x, range(-5 15) || ///
scatter y x || lfit y x, lp(dash)
