Stata：贝叶斯方法-bayes

发布时间：2021-02-17 阅读 7892

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

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

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
. ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
. help lianxh

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

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

作者：宾健博 (中山大学)
邮箱：binjb@mail2.sysu.edu.cn

1. 引言

近年来，贝叶斯估计 (Bayesian estimation) 在宏观经济分析和复杂的计量经济模型中应用越来越广泛。比如，以贝叶斯估计为基础的 VAR (向量自回归) 模型和 DSGE (动态随机一般均衡) 模型的前沿研究，就经常出现在主流经济学期刊上。因此，掌握贝叶斯方法逐渐成为了宏观计量研究者的必修课。

遗憾的是，一般对 VAR 模型和 DSGE 模型的贝叶斯估计都是在 MATLAB 中完成。当然，我们也可以利用 Stata 内置的贝叶斯估计命令处理一些较简单的模型。

基于上述考虑，本文将简要介绍贝叶斯估计和 MCMC 方法的基本概念、如何在 Stata 中使用 bayes 命令前缀、标准估计命令 bayesmh 及其自带的交互菜单。

2. 贝叶斯估计

在统计学和计量经济学中，参数估计都扮演着非常重要的角色。频率学派认为概率是在长期重复的随机实验中产生的性质。但在贝叶斯学派看来，概率是不确定性的一种表达。这使得我们可以对无法重复的事件的概率进行考察，但相应的问题是将主观性引入了概率之中。反映到参数估计上，频率学派认为模型参数是一个定值，应当使用观测数据去接近它；而贝叶斯学派认为模型参数是一个随机变量，有自己的分布，应当使用先验信息并结合观测值去修正它的分布。

因此，贝叶斯估计的核心任务就是结合先验信息和观测值信息，获得模型参数的概率分布 (后验分布，posterior)。我们将所有观测值记为 $y$ ，将模型参数记为 $θ$ 。在估计之前，需要确定先验分布 (prior) 和似然函数：

其中，先验分布 $p (θ)$ 代表了先验信念，即不依赖于观测值的主观看法。比如，在抛一枚硬币之前，我们可能会主观认为它取正面的概率是 50%，这就是一种先验信念；
似然函数 $h (y ∣ θ)$ 则代表了一种将观测值转化为概率的 “模式”，如果将 $h (θ ∣ y)$ 中的 $θ$ 和 $y$ 互换，就可以得到给定参数下观测值的概率 $p (y ∣ θ)$ 。同理，我们将 $h (θ ∣ y)$ 理解为给定观测值时参数的概率。对于同一个问题而言， $h (θ ∣ y)$ 和 $p (y ∣ θ)$ 的数学形式是一样的，区别仅为谁是给定，谁是未知。

在确定了先验分布和似然函数后，可以使用贝叶斯定理 (Bayes' Law) 求解出后验分布：

p (θ ∣ y) = \frac{p (y ∣ θ) p (θ)}{p (y)} = \frac{h (θ ∣ y) p (θ)}{p (y)}

其中：

p (y) = \int_{θ} p (y ∣ θ) p (θ)

在实际的贝叶斯分析中，由于 $p (y)$ 是常数，一般会忽略这一项，只考虑后验分布的分子，即：

p (θ ∣ y) \propto h (θ ∣ y) p (θ)

3. 伯努利实验

在本小节，我们将介绍最简单的贝叶斯估计实例——伯努利实验。假设现有一枚硬币，将其抛起，然后估计正面朝上的概率 $θ$ 。

根据前面的讨论，我们需要明确先验分布和似然函数。由于似然函数不具有主观性，故先予以分析。假设抛了 $n$ 次硬币，其中正面朝上的次数为 $a$ ，则似然函数为：

h (θ ∣ y) = \prod_{i = 0}^{n} P (y_{i} ∣ θ) = θ^{a} (1 - θ)^{n - a}

再来看先验分布。在实际操作中，我们往往选择似然函数的共轭先验 (conjugate prior)。由于共轭先验具有与似然函数相乘后分布族保持不变的特征，我们可以通过简单的数学推导求解出 $h (θ ∣ y) p (θ)$ 。在本例中，似然函数是由伯努利分布 (Bernoulli distribution) 推导而来，故其对应的共轭先验为 Beta 分布。

在选定了先验分布的分布族之后，还需要调整分布的参数以符合先验信念。Beta 分布有两个参数 $α$ 和 $β$ ，通过改变两者的大小，其形状会发生很大的变化。在这里，我们主要关注两种情况：

$B e t a (1, 1)$ ，即均匀分布 (uniform distribution)。均匀分布意味着参数在每一个值的概率密度相同，也称为无信息先验 (non-informative prior) 或者扁平先验 (flat prior)。如果选用这个分布，就代表完全没有参数的先验信息；
$B e t a (20, 20)$ 。由下图，可以发现当 $α = β > 1$ 时， $B e t a$ 分布会出现 “中间高，两头低” 的情况，且概率密度的最大值出现在 0.5 处。如果选用这样的分布，意味着我们认为模型参数 $θ = 0.5$ 的可能性最高，这和我们一般的常识相符 (抛硬币正面朝上的概率应该在 50% 附近)。同时还可以发现 $α$ 和 $β$ 越大，概率密度越向中间集中，因此 $α$ 和 $β$ 可以作为衡量我们先验信息确信程度的参数。这里选择的 20 并不是唯一的，代表了我们对先验信息是比较确信的。

在 Stata 中，我们可以使用上述先验分布对原问题进行估计，具体代码如下：

. clear          //清空工作区
. set seed 20    //设定随机种子
. set obs 10000  //设定观测值数量，在这里 n=10000
. gen y = runiform(0,1) > 0.5 //生成服从伯努利分布的样本，记为 y
. bayesmh y, likelihood(dbernoulli({theta})) prior({theta}, beta(1,1)) //估计命令

Burn-in ...
Simulation ...

Model summary
------------------------------------------------------------------------
Likelihood: 
  y ~ bernoulli({theta})

Prior: 
  {theta} ~ beta(1,1)
------------------------------------------------------------------------

Bayesian Bernoulli model                   MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling   Burn-in          =      2,500
                                           MCMC sample size =     10,000
                                           Number of obs    =     10,000
                                           Acceptance rate  =       .442
Log marginal likelihood = -6935.8231       Efficiency       =      .1992
 
------------------------------------------------------------------------
       |                                                Equal-tailed
       |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------+----------------------------------------------------------------
 theta |  .4991434   .0051053   .000114    .499113   .4892574   .5092233
------------------------------------------------------------------------

可以看出，当设定先验分布为 $B e t a (1, 1)$ 时，平均值的点估计结果为 0.4991434，和真实的数据生成过程中的参数 0.5 很接近。如果选择具有较强先验信息的 $B e t a (20, 20)$ ，同样的设定下点估计值达到了 0.4992197，这证明了先验信息对后验分布的影响。

4. Stata 中贝叶斯估计

上例中，核心代码如下：

. bayesmh y, likelihood(dbernoulli({theta})) prior({theta}, beta(1,1))

其中，likelihood() 指定了似然函数的形式，在本例中似然函数由伯努利分布导出，同时我们将参数称为 theta (在本例中只有一个参数，因此不需要进行区分，只是单纯的命名，在更复杂的例子中需要通过名字区分参数)。prior() 则指定了先验分布，本例中将 theta 的先验分布定为 $B e t a (1, 1)$ 。

上述内容的编程实现起来并不困难，但是在更复杂的问题中，根据似然函数和先验分布的不同，很多时候需要指定额外的参数或选项。此时，贝叶斯估计命令 bayesmh 并不适合，需要通过 Stata 自带的交互菜单生成。具体来看：

我们可以通过点击 Stata 上方的 Statistics 菜单，选择底部的 Bayesian analysis 选项，再根据具体问题选择相应选项，打开贝叶斯估计的交互界面。由于上文的例子不涉及回归分析，选择 General estimation and regression。

由于模型是一个简单的独立实验，每一个观测值都是一次独立实验的结果，因此在 Syntax 处选择 Univariate distributions，Dependent variable 处选择生成的变量 y。Distribution 用于确定似然函数的分布，在本例中选择伯努利分布，同时指定参数名为 theta。最后需要指定模型的先验分布，首先点击 Create 按钮，然后在弹出的菜单中选择相应的分布族和参数，便可以生成一个想要的先验分布。最终结果如下：

点击 OK 键，系统便会自动开始估计，可以得到前文中的结果。

为了说明交互界面的便利性，我们考虑一个稍微复杂一些的例子。假设样本不再是由伯努利分布生成，而是由二项分布 $B (10, 0.3)$ 生成，每一个观测值代表进行了 10 次成功率为 0.3 的伯努利实验，观测值序列记为 $y$ 。我们依旧可以导出其似然函数：

P (X = x) = C_{n}^{x} p^{x} (1 - p)^{n - x} h (p ∣ n, y) = \prod_{i = 1}^{n} C_{n}^{y_{i}} p^{\sum_{i = 1}^{n} y_{i}} (1 - p)^{n^{2} - \sum_{i = 1}^{n} y_{i}}

在这个例子中，原分布和似然函数中有两个参数 $n$ 和 $p$ ，其中， $n$ 代表二项分布的第一个参数， $p$ 代表二项分布的第二个参数。并且，参数 $n$ 是预先给定的，需要在使用交互界面时需要输入这个参数。沿用之前的先验分布 (二项分布的共轭先验也是 Beta 分布)，此时的代码、结果和菜单操作如下：

. clear          //清空工作区
. set seed 20    //设定随机种子
. set obs 10000  //设定观测值数量，在这里 n=10000
. gen z = rbinomial(10,0.3) //生成服从伯努利分布的样本，记为 z
. bayesmh z, likelihood(dbinomial({theta},10)) prior({theta}, normal(0,10000))

Burn-in ...
Simulation ...

Model summary
------------------------------------------------------------------------
Likelihood: 
  z ~ binomial({theta},10)

Prior: 
  {theta} ~ beta(1,1)
------------------------------------------------------------------------

Bayesian binomial model                    MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling   Burn-in          =      2,500
                                           MCMC sample size =     10,000
                                           Number of obs    =     10,000
                                           Acceptance rate  =      .4526
Log marginal likelihood = -17700.763       Efficiency       =      .2339
 
------------------------------------------------------------------------
       |                                                Equal-tailed
       |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------+----------------------------------------------------------------
 theta |  .2987419   .0014841   .000031   .2987601   .2959134   .3016467
------------------------------------------------------------------------

可以看到，当似然函数为二项分布之后，在菜单界面中多了一个 Bernoulli trials 的输入框，也就是上文指定的 $n$ 。在更复杂的问题和模型中，可能会要求输入更多的参数和其他内容，此时菜单界面可以很好地代替复杂的语句指令，免去了繁杂的语法记忆。

不仅如此，我们还可以通过键入命令 db postest 打开后验分析的交互菜单进行分析。常用的功能如 Bayesian analysis 下的 Graphical summaryies and convergence diagnostics 菜单可以绘制很多与后验分布相关的图，Hypothesis testing using model posterior probabilities 和 Interval hypothesis testing 菜单可以利用生成的后验分布进行假设检验。

5. MCMC 方法和线性回归的贝叶斯估计

在前面的例子中，可能会有读者感到疑惑，当先验分布选为共轭先验时，后验分布是有解析表达式的，这意味着其概率密度函数应该是确定的。但在结果显示的界面中，明显有模拟 (simulation)、迭代 (iteration) 和抽样 (sampling) 这样的常见于数值方法的字眼。这其实是因为 Stata 并非使用解析方法对后验分布进行求解。为了进一步深入了解现代贝叶斯方法的操作流程，我们有必要对 MCMC (Markov chain Monte Carlo，马尔科夫链蒙特卡洛) 方法进行介绍。

对于 MCMC 方法，维基百科定义为：马尔可夫链蒙特卡洛方法包括一类用于从概率分布中采样的算法。通过构建具有所需分布作为其平衡分布的马尔可夫链，可以通过记录链中的状态来获得所需分布的样本。

这一段话有些抽象，实际上从应用层面上来说，我们只需要了解两个最常用的 MCMC 方法：Metropolis–Hastings 算法 (简称 M-H 算法) 和吉布斯采样 (Gibbs sampling)。前者是 Stata 内置的对后验分布进行抽样的算法，其主要作用是通过后验分布的概率密度之比构造了一个具有接受-拒绝机制的马尔科夫链，并通过它抽取出所需的后验分布。这一算法也广泛应用于 DSGE 模型的贝叶斯估计中，但它不是本文的重点。

在宏观计量经济学中最常用的 (主要是因为其应用于 VAR 模型的贝叶斯估计) MCMC 方法是吉布斯采样。为了让大家对吉布斯采样的应用背景有一个简单的了解，下面将以线性回归为例进行介绍。

假设有如下的线性回归模型：

\begin{aligned} Y_{t} & = X_{t} B + v_{t} \\ v_{t} & \sim N (0, σ^{2}) \end{aligned}

其中， $Y_{t}$ 是一个 $T \times 1$ 的矩阵， $X_{t}$ 是一个 $T \times K$ 的矩阵， $B$ 是一个 $K \times 1$ 的矩阵。其似然函数的表达式如下：

h (B, σ^{2} ∣ Y_{t}) = {(2 π σ^{2})}^{- T / 2} \exp (- \frac{{(Y_{t} - X_{t} B)}^{'} (Y_{t} - X_{t} B)}{2 σ^{2}})

这里有两个待估参数，即系数矩阵 $B$ 和扰动项方差 $σ^{2}$ 。分别考虑以下三种情况：

第一种情况：系数矩阵 $B$ 未知，扰动项方差 $σ^{2}$ 已知；
第二种情况：系数矩阵 $B$ 已知，扰动项方差 $σ^{2}$ 未知；
第三种情况：系数矩阵 $B$ 和扰动项方差 $σ^{2}$ 均未知。

在第一种情况中，由于扰动项方差已知，此时的共轭先验是正态分布，因此假定系数矩阵 $B$ 的先验分布为 $B \sim N (B_{0}, Σ_{0})$ 。

先验分布：

\begin{aligned} f (B) & = (2 π)^{- K / 2} {| Σ_{0} |}^{- \frac{1}{2}} \exp [- 0.5 {(B - B_{0})}^{'} Σ_{0}^{- 1} (B - B_{0})] \\ \propto \exp [- 0.5 {(B - B_{0})}^{'} Σ_{0}^{- 1} (B - B_{0})] \end{aligned}

似然函数：

\begin{aligned} h (B ∣ σ^{2}, Y_{t}) & = {(2 π σ^{2})}^{- T / 2} \exp (- \frac{{(Y_{t} - X_{t} B)}^{'} (Y_{t} - X_{t} B)}{2 σ^{2}}) \\ \propto \exp (- \frac{{(Y_{t} - X_{t} B)}^{'} (Y_{t} - X_{t} B)}{2 σ^{2}}) \end{aligned}

后验分布：

p (B ∣ σ^{2}, Y_{t}) \propto \exp [- 0.5 {(B - B_{0})}^{'} Σ_{0}^{- 1} (B - B_{0})] \times \exp (- \frac{{(Y_{t} - X_{t} B)}^{'} (Y_{t} - X_{t} B)}{2 σ^{2}})

可以发现后验分布仍然是正态分布：

p (B ∣ σ^{2}, Y_{t}) \sim N (M^{*}, V^{*}) \begin{aligned} M^{*} & = {(Σ_{0}^{- 1} + \frac{1}{σ^{2}} X_{t}^{'} X_{t})}^{- 1} (Σ_{0}^{- 1} B_{0} + \frac{1}{σ^{2}} X_{t}^{'} Y_{t}) \\ V^{*} & = {(Σ_{0}^{- 1} + \frac{1}{σ^{2}} X_{t}^{'} X_{t})}^{- 1} \end{aligned}

在第二种情况中，由于似然函数中的给定变量改变，其共轭先验也随之发生变化。查表知此时扰动项方差 $σ^{2}$ 应服从倒伽马分布 (Inverse Gamma distribution)，具体运算过程和第一种情况类似，可以获得相同分布族的后验分布 $p (σ^{2} ∣ B, Y_{t})$ 。

在实际中，更加常见的是第三种情况，即两个待估参数都未知。此时，没有办法获得单一的共轭先验，一种解决办法是设定联合先验分布：

p (B, σ^{2}) = P (σ^{2}) \times P (B ∣ σ^{2})

其中， $σ^{2}$ 的分布仍是无条件的，和第二种情况中的设定一样，但 $B$ 则是有条件的，我们将其设定为：

P (B ∣ σ^{2}) \sim N (B_{0}, σ^{2} Σ_{0})

在这样的设定下，联合先验分布还是一个正态分布，和似然函数的分布族一致，这种情况称之为自然共轭先验 (natural conjugate prior)。

需要注意的是，最后解出的后验分布是一个联合后验分布 $p (σ^{2}, B ∣ Y_{t})$ ，但我们想考察的后验分布则是边缘后验分布 $p (B ∣ σ^{2}, Y_{t})$ 和 $p (σ^{2} ∣ B, Y_{t})$ 。在自然共轭先验的设定下，可以通过积分的方法求得边缘后验分布的解析形式。但这也对先验分布的选择施加了很强的限制。在较复杂的模型中，我们往往不希望限制先验分布的选择，但非自然共轭先验难以通过积分求解出边缘后验分布，因此需要应用吉布斯采样。

吉布斯采样的思想其实很简单，当有了两个能条件分布的抽取器时 (即第一种和第二种情况)，可以通过迭代的方法获得逼近边缘分布的抽样。在本例中，具体操作为：

设定系数矩阵 $B$ 的无条件先验分布和扰动项方差 $σ^{2}$ 的无条件先验分布，并从 $σ^{2}$ 的无条件先验分布中抽取一个样本；
根据第一种情况的做法，以第 1 步中抽取的 $σ^{2}$ 为条件，从获得的后验分布中抽取一个 $B$ 的样本；
根据第二种情况的做法，以第 2 步中抽取的 $B$ 为条件，从获得的后验分布中抽取一个 $σ^{2}$ 的样本；
重复第 2 步和第 3 步，获得一个经验样本，可用于估计边缘后验分布。

这一操作步骤看似复杂，但在 Stata 中的实现其实很容易，只需要在线性回归的命令前加上 bayes: 即可，但需要注意的是，此时回归命令不能缩写，并且如果要用吉布斯采样，需要使用 bayes,gibbs:。

. clear       //清空工作区
. set seed 20 //设定随机种子
. sysuse auto //使用系统自带的 auto 数据集
. bayes, gibbs: regress price mpg weight //执行吉布斯采样的贝叶斯估计

Burn-in ...
Simulation ...

Model summary
--------------------------------------------------------------------------
Likelihood: 
  price ~ normal(xb_price,{sigma2})

Priors: 
  {price:mpg weight _cons} ~ normal(0,10000)                           (1)
                  {sigma2} ~ igamma(.01,.01)
--------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_price.

Bayesian linear regression                   MCMC iterations  =     12,500
Gibbs sampling                               Burn-in          =      2,500
                                             MCMC sample size =     10,000
                                             Number of obs    =         74
                                             Acceptance rate  =          1
                                             Efficiency:  min =      .9457
                                                          avg =      .9864
Log marginal likelihood = -696.85658                      max =          1
 
--------------------------------------------------------------------------
         |                                                Equal-tailed
         |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
---------+----------------------------------------------------------------
price    |
     mpg | -5.282435   27.88942   .278894  -5.410969  -59.53671   49.12283
  weight |  2.074039    .198759   .001988   2.076242   1.682404   2.464258
   _cons |  1.928181   99.76441   .997644    2.99603  -202.2644   194.9716
---------+----------------------------------------------------------------
  sigma2 |   6448502    1099995   11311.5    6331241    4641478    8928296
--------------------------------------------------------------------------
Note: Default priors are used for model parameters.

6. 参考资料

M-H 算法和共轭先验的维基百科 -Link1- -Link2-
Blake A P, Mumtaz H. Applied Bayesian econometrics for central bankers[J]. Technical Books, 2012. -Link-
Koop G, Korobilis D. Bayesian multivariate time series methods for empirical macroeconomics[M]. Now Publishers Inc, 2010. -Link-
Bill Rising. Bayesian Analysis using Stata. 2015. -Link-

7. 相关推文

Note：产生如下推文列表的命令为：
lianxh VAR DSGE 分布, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

Stata：贝叶斯方法-bayes

1. 引言

2. 贝叶斯估计

3. 伯努利实验

4. Stata 中贝叶斯估计

5. MCMC 方法和线性回归的贝叶斯估计

6. 参考资料

7. 相关推文

相关课程

课程一览

关于我们

热门资讯