Stata程序：Monte-Carlo-模拟之产生符合特定分布的随机数

发布时间：2020-10-08 阅读 4733

Stata 连享会主页 || 视频 || 推文

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者：陈勇吏 (上海交通大学安泰经济与管理学院)

本文讨论随机数的产生原理，以及随机数的代码实现方式。在介绍随机数之前，我们先简单介绍一下蒙特卡洛模拟和随机数的概念。

假设 $X$ 是一个随机变量，其分布函数为 $F (x)$ ，概率密度函数为 $f (x)$ 。我们可能会对 $E [g (X)]$ 感兴趣（ $g$ 是一个给定的函数）。比如，若随机变量 $X$ 代表产品的销量，则销售利润可以表示为 $g (X) = (p - c) X$ ，这也是一个随机变量。在已知随机变量 $X$ 概率分布的情况下，我们希望知道销售利润的期望值 $E [g (X)]$ 是多少。带入随机变量的期望公式得到： $E [g (X)] = \int g (x) f (x) d x$ ，可以很容易求出积分结果。

但现实生活中，我们感兴趣的 $E [g (X)]$ 中的 $X$ 往往不是单一的随机变量，而是由多个随机变量组成的随机向量 $X = (X_{1}, X_{2}, . . ., X_{n})$ 。若已知随机向量 $X$ 的联合概率密度函数 $f (x_{1}, x_{2}, . . . x_{n})$ ，则： $E [g (X)] = \int \int \cdot \cdot \cdot \int g (x_{1}, x_{2}, . . ., x_{n}) f (x_{1}, x_{2}, . . ., x_{n}) d x_{1} d x_{2} \cdot \cdot \cdot d x_{n}$ 。多数情况下，这个多重积分很难求解。此时，我们考虑另一种方式来近似计算 $E [g (X)]$ 的值。

近似计算方法： $S t e p 1$ ：产生一个联合概率密度函数为 $f (x_{1}, x_{2}, . . . x_{n})$ 的随机向量 $X^{(1)} = (X_{1}^{(1)}, X_{2}^{(1)}, . . ., X_{n}^{(1)})$ ，并计算 $Y^{(1)} = g (X^{(1)})$ 。 $S t e p 2$ ：重复 $r$ 次 $S t e p 1$ 的过程，得到独立同分布的 $Y^{(i)} = g (X^{(i)}), i = 1, . . ., r$ 。 $S t e p 3$ ：使用 $Y^{(i)}$ 的平均值（样本均值）来近似 $E [g (X)]$ （总体均值）。

lim_{n \to \infty} \frac{Y^{(1)} + Y^{(2)} + \cdot \cdot \cdot + Y^{(r)}}{r} = E [Y^{(i)}] = E [g (X)] （中心极限定理）

上述近似计算方法称为蒙特卡洛模拟方法。即通过从相同的概率分布中多次抽取独立同分布的样本，使用样本平均值近似总体平均值。将每一次抽取到的独立同分布的样本记为 $X^{(i)} = (X_{1}^{(i)}, X_{2}^{(i)}, . . ., X_{n}^{(i)})$ ，则样本中的每一个元素 $X_{k}^{(i)} (k = 1, 2, . . ., n)$ 都是一个随机变量。我们将这种从给定的概率分布中抽取出来的独立的随机变量，称为随机数。

1. 如何产生随机数

随机变量的分布形式多种多样，包括连续性随机变量、离散型随机变量。其中连续性随机变量包括均匀分布、指数分布等可以写出分布函数形式的随机变量，以及正态分布、卡方分布等只存在概率密度函数形式的随机变量。针对不同的概率分布形式，生成随机数的方式有所不同。但生成任意分布的随机数的基础，是生成服从 $U (0, 1)$ 均匀分布的随机数，其他分布的随机数是在 $U (0, 1)$ 随机数的基础上产生。

1.1 服从均匀分布的随机数

1.1.1 产生过程

产生 $U (0, 1)$ 分布的一种方式是：将写有 $0, 1, . . ., 9$ 的完全相同的 10 个纸条放在一个盒子中，每次从中有放回地抽取一张纸条，连续抽取 $n$ 次。将抽取到的 $n$ 个数字按抽取顺序排列在小数点的后面，得到的数字即为服从 $U (0, 1)$ 均匀分布的随机数。例如：如果抽取的数字序列是 $3, 8, 7, 2, 1$ ，则得到的服从 $U (0, 1)$ 分布的随机数是 $0.38721$ 。

然而，数字计算机通常不使用这种方式来产生随机数，而是使用伪随机数来代替随机数。伪随机数的产生原理是取模。即给定一个初始值 $X_{0}$ ，通过递归公式

X_{n + 1} = (a X_{n} + c) modulo m, n \geq 0

可以得到一个序列

{X_{n}}

，其中

a, c, m

是设定的正整数。该等式表示将

(a X_{n} + c)

除以

m

的余数赋值给

X_{n + 1}

，因此序列中的每个

X_{n}

都只能在

(0, 1, . . ., m - 1)

中取值。

\frac{X_{n}}{m}

是一个在

[0, 1)

之间取值的数，可以看做近似服从

U (0, 1)

均匀分布。可以看出，当设定合适的

a, c, m

值后，得到的伪随机数序列完全取决于给定的序列初始值

X_{0}

，这个初始值称为种子。设定相同的种子，可以得到完全相同的随机序列。这也是伪随机数的优点之一：可以重复相同的模拟过程。

1.1.2 Stata 代码实现

Stata 有很多用来产生随机数的内置函数，其中生成均匀分布随机数的函数为runiform()。在生成随机数前，可以通过set seed #命令设定随机数的初始值（种子）。具体代码如下：

生成 3 个 $U (0, 1)$ 的随机数序列，每个序列包含 10 个随机数。

*===== 代码 =====
clear
set obs 10
forvalues i = 1/3 {
	gen unif`i' = runiform()
}
*===== 生成结果 =====
unif1	unif2	unif3
.1387824	.56928813	.32682486
.92432766	.7060023	.30766002
.3676381	.81376687	.21982671
.11971952	.06189653	.69186391
.91902139	.27757641	.04663675
.60349285	.41005086	.60378652
.89752506	.22520046	.52267927
.48974351	.34788238	.23261151
.21933775	.17111017	.92006581
.88946164	.72442917	.39863705
/* 三次模拟，产生三个不同的随机数序列 */

生成 3 个 $U (0, 1)$ 的随机数序列，每个随机数序列都设定相同的种子。

*===== 代码 =====
clear
set obs 10
forvalues i = 1/3 {
	set seed 10
	gen unif`i' = runiform()
}
*===== 生成结果 =====
unif1	unif2	unif3
.60128311	.60128311	.60128311
.91370436	.91370436	.91370436
.26673195	.26673195	.26673195
.60736577	.60736577	.60736577
.03607367	.03607367	.03607367
.75318004	.75318004	.75318004
.45917266	.45917266	.45917266
.90569935	.90569935	.90569935
.31343658	.31343658	.31343658
.3755541	.3755541	.3755541
/* 设定相同的随机数，三次模拟都产生相同的随机数序列。*/

1.2 连续型随机变量的随机数

对于任意给定的连续型随机变量，若分布函数 $F$ 已知，则该随机变量的随机数可以通过 $F^{- 1} (u)$ 获得，其中 $u$ 为 $U (0, 1)$ 的随机数。换言之，连续型随机变量的随机数，可以通过对 $u \sim U (0, 1)$ 求分布函数的逆运算获得。证明：记 $X = F^{- 1} (U), U \sim U (0, 1)$ ，则 $X$ 的分布函数 $F_{X} (x)$ 满足：

F_{X} (x) = P {X \leq x} = P {F^{- 1} (U) \leq x} = P {F (F^{- 1} (U)) \leq F (x)} = P {U \leq F (x)} = F (x)

证得，通过 $X = F^{- 1} (U)$ 得到的随机数服从 $F$ 分布。

Stata 代码实现

生成连续型随机变量的随机数函数包括：均匀分布 runiform(a,b)、正态分布 rnormal(m,s)、指数分布 rexponential(b)、卡方分布 rchi2(df)、t 分布 rt(df) 等。详细内容可以通过 help random 命令查看。

生成随机数的具体代码如下所示：

使用 Stata 的内置函数生成随机数

*========== 代码：生成随机数序列 ==========
* 1.生成随机数序列
  clear
  set obs 1000
 *1.1 U(1,5)均匀分布的随机数序列
  gen runif = runiform()
 *1.2 指数分布 Exp(5) 的随机数序列
  gen exp = rexponential(5)
 *1.3 标准正态分布的随机数序列
  gen rnormal = rnormal()
 *1.4 N(3,2)正态分布的随机数序列
  gen rnormal_3_2 = rnormal(3, 2)
 *1.5 自由度为 3 的 t 分布的随机数序列
  gen rt = rt(3)
 *1.6 自由度为 43 的卡方分布的随机数序列
  gen rchi2 = rchi2(43)

* 2.随机数序列的分布图
  foreach dist in runif exp rnormal rnormal_3_2 rt rchi2 {
  	histogram `dist', name(`dist', replace) nodraw //kdensity
  }
  graph combine runif exp rnormal rnormal_3_2 rt rchi2

使用分布函数求逆的方式生成随机数 指数分布 $E x p (λ)$ 的分布函数为： $F (x) = 1 - e^{- \frac{x}{λ}}$ ，反函数为 $F^{- 1} (u) = - λ l n (1 - u)$ 均匀分布 $U (a, b)$ 的分布函数为： $F (x) = \frac{x - a}{b - a}$ ，反函数为 $F^{- 1} (u) = a + (b - a) u$

由于 $u \sim U (0, 1)$ ，故 $1 - u \sim U (0, 1)$ 。所以，指数分布 $E x p (λ)$ 的随机数可以通过 $- \frac{1}{λ} l n (u)$ 得到， $U (a, b)$ 均匀分布的随机数可以通过 $a + (b - a) u$ 得到。

*========== 代码 ==========
clear
set obs 10
* 1.指数分布 Exp(5) 的随机数序列
set seed 20
gen exp_f = -5 * log(runiform()) //分布函数求逆的方式
set seed 20
gen exp = rexponential(5) //使用内置函数的方式
* 2.均匀分布 U(1,5) 的随机数序列
set seed 20
gen unif_f = 1 + (5-1)*runiform() //分布函数求逆的方式
set seed 20
gen unif = runiform(1,5) //使用内置函数的方式

*========== 生成结果 ==========
exp_f	exp	unif_f	unif
1.6727648	1.6727648	3.8626318	3.8626318
.47869168	.47869168	4.634807	4.634807
1.1330894	1.1330894	4.1889014	4.1889014
7.785189	7.785189	1.8430378	1.8430378
1.0036905	1.0036905	4.2725067	4.2725067
.02883674	.02883674	4.976997	4.976997
4.806301	4.806301	2.5296427	2.5296427
8.2741926	8.2741926	1.7644917	1.7644917
7.1020147	7.1020147	1.9664666	1.9664666
10.250074	10.250074	1.514932	1.514932
/* 可以看出，两种方式得到的随机数序列相同。 */

当感兴趣的随机变量在 Stata 中没有相应的内置函数来生成随机数时，可以通过分布函数求逆的方式计算。

1.3 离散型随机变量的随机数

对于任意给定的离散型随机变量，若概率分布为 $P {X = x_{j}} = P_{j}$ ，则该随机变量的随机数同样可以通过对 $u \sim U (0, 1)$ 求分布函数的逆运算获得。具体方式为：令 $U \sim U (0, 1)$ ，

X = {\begin{cases} x_{1}, & i f U < P_{1} \\ x_{2}, & i f P_{1} < U < P_{1} + P_{2} \\ ⋮ \\ x_{j}, & i f \sum_{1}^{j - 1} P_{j} < U < \sum_{1}^{j} P_{j} \\ ⋮ \end{cases}

可以看出 $X$ 的概率分布为 $P {X = x_{j}} = P {\sum_{1}^{j - 1} P_{j} < U < \sum_{1}^{j} P_{j}} = P_{j}$ 。

Stata 代码实现

生成离散型随机变量的随机数函数包括：二项分布 rbinomial(n,p)、超几何分布 rhypergeometric(N,K,n)、泊松分布 rpoisson(m)、负二项分布 rnbinomial(n,p) 等。详细内容可以通过 help random 命令查看。

生成随机数的具体代码如下所示：

*========== 代码：生成随机数序列 ==========
* 1.生成随机数序列
clear
set obs 1000
 *1.1 二项分布 B(10,0.5) 的随机数序列
gen rbinomial = rbinomial(10,0.5)
 *1.2 超几何分布 H(20,4,10) 的随机数序列
 /*
	其中：20 代表总体数量
		  4  代表总体中指定种类的数量
		  10 代表抽取的样本数量
 */
gen rhypergeometric = rhypergeometric(20,4,10)
 *1.3 泊松分布 Poisson(20) 的随机数序列
gen rpoisson = rpoisson(20)
 *1.4 负二项分布 NB(123,0.2) 的随机数序列
gen rnbinomial = rnbinomial(123,0.2)

* 2.随机数序列的分布图
foreach dist in rbinomial rhypergeometric rpoisson rnbinomial {
	histogram `dist', name(`dist', replace) nodraw //kdensity
}
graph combine rbinomial rhypergeometric rpoisson rnbinomial

2. 如何产生具有相关性的随机数

Stata 提供了 drawnorm 命令，可以从多元正态分布中生成随机数序列。具体语法为 drawnorm newvarlist [, options]，常用的 option 选项包括：

corr(matrix)：设定多元正态分布的协方差矩阵
mean(vector)：设定多元正态分布的均值向量，默认为 mean(0)。具体规则可以查看 help drawnorm。

2.1 相关性

Pearson 相关系数和 Spearman 相关系数是衡量变量相关性最常用的两个指标。

Pearson相关系数用来衡量变量间的线性关系。计算公式为
$ρ_{X, Y} = \frac{C o v (X, Y)}{σ_{X} σ_{Y}} = \frac{E [X Y] - E [X] E [Y]}{\sqrt{E [X^{2}] - E [X]^{2}} \sqrt{E [Y^{2}] - E [Y]^{2}}}$
Spearman相关系数等级变量的 Pearson 相关系数。将两个变量分别排序以后生成相应的等级数据，使用两个变量对应的等级数据计算 Pearson 相关系数。

Stata 举例如下：

clear
set obs 1000
* 1.生成相关系数矩阵
/*
  相关系数矩阵是对称矩阵，矩阵“维数”与待生成的“随机序列个数”一致。
  其中：
  	（第 i 个）主对角元素：代表（第 i 个）随机变量的【方差】，
  	（第 i 行 j 列 的）非主对角元素：代表（第 i,j）两个随机变量间的【协方差】
*/
matrix C = (1, 0.5 \ 0.5, 1)
matrix list C
*======= Stata 结果 ========
symmetric C[2,2]
    c1  c2
r1   1
r2  .5   1
*===========================

* 2.生成二元正态分布样本（生成两个随机数序列：均值为 0，协方差矩阵为 C）
drawnorm x y, corr(C) //means(0,0)选项设定均值向量为(0,0)，可省略不写
* 3.查看两个随机序列的相关性
 *3.1 Pearson 相关系数
corr x y
*========== Stata 结果 ===========
(obs=1,000)
             |        x        y
-------------+------------------
           x |   1.0000
           y |   0.5218   1.0000
*=================================

 *3.2 Spearman 相关系数
  *3.2.1 使用 Pearson 系数 计算 Spearman 相关系数
sort x
gen rank_x = _n
sort y
gen rank_y = _n
corr rank_x rank_y
*========== Stata 结果 ===========
(obs=1,000)
             |   rank_x   rank_y
-------------+------------------
      rank_x |   1.0000
      rank_y |   0.5011   1.0000
*=================================
  *3.2.2 使用 内置函数 计算 Spearman 相关系数
spearman x y
*========== Stata 结果 ===========
 Number of obs =    1000
Spearman's rho =       0.5011

Test of Ho: x and y are independent
    Prob > |t| =       0.0000
*=================================*
/* 两种方式得到的 Spearman 相关系数 相等，都等于 0.5011。*/

从 Pearson 和 Spearman 两个相关系数的定义及计算公式可以看出：对变量做正的线性变换，变量间的 Pearson 相关系数不变；对变量做非负的函数变换，变量间的 Spearman 相关系数保持不变。

2.2 产生具有相关性的随机序列

产生具有相关性的正态分布随机序列，可以使用 Stata 的 drawnorm 命令。产生具有相关性的其他分布的随机数序列没有现成的 Stata 命令，可以先将 drawnorm 得到的正态随机序列转换为均匀分布，然后通过分布函数求逆的方式得到想要的随机数序列。原理如下：

分布函数 $Y = F_{X} (X)$ 服从 $U (0, 1)$ 均匀分布
$F (y) = {\begin{cases} 0, y < 0 \\ P {Y \leq y} = P {F_{X} (X) \leq y} = P {X \leq F_{X}^{- 1} (y)} = F (F_{X}^{- 1} (y)) = y, 0 \leq y < 1 \\ 1, y \geq 1 \end{cases}$

若 $X$ 是正态分布的随机序列， $F_{X}$ 是正态分布的分布函数，则 $F_{X} (X)$ 是均匀分布的随机序列。

* 1.生成均匀分布随机序列
clear
set obs 10000
 *1.1 生成两个相关的正态随机序列：x,y
matrix C = (1, 0.5 \ 0.5, 1)
drawnorm x y, corr(C)
 *1.2 根据正态分布函数，生成与 x,y 对应的均匀分布序列：x_unif,y_unif
gen x_unif = normal(x)
gen y_unif = normal(y)
 *1.3 作图展示生成的随机数序列的分布情况
foreach dist in x y x_unif y_unif {
	hist `dist', name(`dist', replace) nodraw
}
graph combine x y x_unif y_unif
/* 可以看出，生成的 x_unif y_unif 两个随机序列在 (0,1) 上均匀分布 */

分布函数是单调不减函数，故 $F_{X} (X)$ 与 $X$ 具有相同的 Spearman 系数 可以验证，上面代码中生成的 x_unif 和 y_unif 之间的 Spearman 相关系数，与正态序列 x 和 y 之间的 Spearman 相关系数是一样的。

* 1.1 正态序列的 Spearman 值
spearman x y
/*============= Stata 结果 =============
 Number of obs =   10000
Spearman's rho =       0.4917
Test of Ho: x and y are independent
    Prob > |t| =       0.0000
*======================================*/
* 1.2 生成的均匀分布序列的 Spearman 值
spearman x_unif y_unif
/*============= Stata 结果 =============
 Number of obs =   10000
Spearman's rho =       0.4917
Test of Ho: x_unif and y_unif are independent
    Prob > |t| =       0.0000
*======================================*/

但 Pearson 相关系数有所不同：

* 2.1 正态序列的 Pearson 值
corr x y
/*============= Stata 结果 =============
(obs=10,000)
             |        x        y
-------------+------------------
           x |   1.0000
           y |   0.5046   1.0000
*======================================*/
* 2.2 生成的均匀分布序列的 Pearson 值
corr x_unif y_unif
/*============= Stata 结果 =============
(obs=10,000)
             |   x_unif   y_unif
-------------+------------------
      x_unif |   1.0000
      y_unif |   0.4912   1.0000
*======================================*/

/* 可以看出，两组变量的 Spearson 系数值完全一样，
            两组变量的 Pearson 系数值有些微不同   */

任意分布 $F$ 的随机序列都可以通过 $F^{- 1} (U) = F^{- 1} (F_{X} (X))$ 的方式得到 使用 Stata 内置的逆分布函数 invgamma() 和 invexponential()，带入正态分布产生的均匀分布 $U = F_{X} (X)$ ，生成 伽马 Gamma 随机序列和 指数 Exponential 随机序列：

* 1.生成 Gamma 随机序列和 Exponential 随机序列
gen x_gamma = invgammap(4, x_unif)
gen y_exponential = invexponential(1, y_unif)
* 2.做序列的分布图
hist x_gamma, name(x_gamma, replace)
hist y_exponential , name(y_exponential , replace)

分布函数的反函数 $F^{- 1} (U)$ 是单调不减函数，故 $F^{- 1} (F_{X} (X))$ 与 $X$ 具有相同的 Spearman 系数 可以验证，上面代码中生成的伽马随机序列 x_gamma 和指数随机序列 y_exponential 之间的 Spearman 相关系数，与正态序列 x 和 y 之间的 Spearman 相关系数是一样的。

spearman x y
/*============= Stata 结果 =============
 Number of obs =   10000
Spearman's rho =       0.4758
Test of Ho: x and y are independent
    Prob > |t| =       0.0000
*======================================*/
spearman x_gamma y_exponential
/*============= Stata 结果 =============
 Number of obs =   10000
Spearman's rho =       0.4758
Test of Ho: x_gamma and y_exponential are independent
    Prob > |t| =       0.0000
*======================================*/

但 Pearson 相关系数有所不同：

spearman x y
/*============= Stata 结果 =============
(obs=10,000)
             |        x        y
-------------+------------------
           x |   1.0000
           y |   0.4894   1.0000
*======================================*/
spearman x_gamma y_exponential
/*============= Stata 结果 =============
(obs=10,000)
             |  x_gamma y_expo~l
-------------+------------------
     x_gamma |   1.0000
y_exponent~l |   0.4515   1.0000
*======================================*/

如此，drawnorm 命令不仅可以生成 正态 随机序列，也可以进一步通过分布函数求逆的方式得到相关性相同的 其他分布 的随机序列。使用 drawnorm 命令生成具有相关性的 非正态分布 随机序列的 Stata 完整代码如下所示：

clear
set obs 10000
* 1.定义相关系数矩阵 
matrix C = (  1,  0.7, -0.3,  0.2 \ ///
            0.7,    1,  0.2, -0.1 \ ///
		   -0.3,  0.2,    1,  0.3 \ ///
			0.2, -0.1,  0.3,    1)
* 2.生成 4 个相关系数矩阵为 C 的正态随机向量
drawnorm x1 x2 x3 x4, corr(C)
* 3.根据 逆分布函数 方法生成四个其他分布的随机序列
/* 
    四个分布包括：
		自由度为 5 的卡方分布：chi2(5),
		均值为 5 的泊松分布：poisson(5),
		自由度为 5,5 的 F 分布：F(5,5)
		在(-5,5)之间的均匀分布：uniform(-5,5) 
*/
gen x_chi2 = invchi2(5, normal(x1))
gen x_poisson = invpoisson(5, 1-normal(x2)) // normal(x2) 和 1-normal(x2) 都是 U(0,1)均匀分布的随机序列
gen x_F = invF(5, 5, normal(x3))
gen x_uniform = -5 + normal(x4)*10
* 4.查看四个生成序列的分布图
foreach dist in x_chi2 x_poisson x_F x_uniform {
	hist `dist', name(`dist', replace) nodraw
}
graph combine x_chi2 x_poisson x_F x_uniform, imargin(large)
* 5.查看相关系数
mat list C
spearman x?
spearman x_*
corr x?
corr x_*

连享会计量方法专题……

3. 应用举例

在计量模型中，随机数可以很好的模拟扰动项的随机性，以及我们对数据特征的高质量要求。这有助于我们按照预计的分布特征去抽样，进而验证计量模型的有效性，比较计量模型的优劣。我们将引入一个简单的例子，介绍随机数如何应用于 似不相关回归（SUR） 和 OLS 的模型比较。

例子

我们考虑应用 似不相关回归（SUR） 方法来决定是否开发小行星铂金矿物资源的例子，具体内容参见 Pricing an Asteroid。

在本例中，我们尝试估计三个存在相关性的模型。三个模型具有相同的解释变量： 市场现有铂金数量 (X1)，投资者们的平均年龄 (X2)，宇宙飞船的新闻报道数量(X3)。被解释变量分别为：铂金价格 (Y1)，宇宙飞船的投资金额 (Y2)，宇宙飞船的新闻报道数量 (Y3)。

模型设定如下：

{\begin{cases} Y_{1} = β_{10} + β_{11} X_{1} + β_{12} X_{2} + β_{13} X_{3} + ε_{1} \\ Y_{2} = β_{20} + β_{21} X_{1} + β_{22} X_{2} + β_{23} X_{3} + ε_{2} \\ Y_{3} = β_{30} + β_{31} X_{1} + β_{32} X_{2} + β_{33} X_{3} + ε_{3} \end{cases}

先生成三个解释变量：

. clear
. set obs 10000
. gen x1 = 10 + rnormal()
. gen x2 = 50 + rnormal()
. gen x3 = rpoisson(10)

生成三个线性相关的扰动项

. matrix C = (12, -5,  5 \ -5, 12, -3 \ 5, -3,  6 )
. drawnorm u1 u2 u3, cov(C)

生成三个被解释变量

若模型中的系数矩阵为

B = {\begin{matrix} 20 & 3 & 0 & - 0.1 \\ 5 & 3 & 3.5 & 5.2 \\ 10 & 3 & 5.3 & - 10 \end{matrix}}

则生成被解释变量的 Stata 命令如下：

. gen y1 = 20 + 3*x1 + 0*x2 - 0.1*x3 + u1
. gen y2 = 5 + 3*x1 + 3.5*x2 + 5.2*x3 + u2
. gen y3 = 10 + 3*x1 + 5.3*x2 - 10*x3 + u3

OLS 回归

. reg y1 x1 x2 x3
. reg y2 x1 x2 x3
. reg y3 x1 x2 x3

SUR 回归

*-不加约束：
. sureg (y1 = x1 x2 x3) (y2 = x1 x2 x3) (y3 = x1 x2 x3)

*-加约束：三个模型的 x1 系数相等
. constraint 1 [y1]x1 = [y2]x1
. constraint 2 [y1]x1 = [y3]x1
. sureg (y1 = x1 x2 x3) (y2 = x1 x2 x3) (y3 = x1 x2 x3), const(1 2)

OLS 估计结果：

{\begin{cases} Y_{1} = 22.224 + 3.043 X_{1} - 0.049 X_{2} - 0.116 X_{3} + ε_{1} \\ Y_{2} = 2.264 + 2.944 X_{1} + 3.566 X_{2} + 5.203 X_{3} + ε_{2} \\ Y_{3} = 11.207 + 2.997 X_{1} + 5.278 X_{2} - 10.006 X_{3} + ε_{3} \end{cases}

SUR 估计结果（不加约束）：