温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh
作者:宾健博 (中山大学)
邮箱:binjb@mail2.sysu.edu.cn
目录
近年来,贝叶斯估计 (Bayesian estimation) 在宏观经济分析和复杂的计量经济模型中应用越来越广泛。比如,以贝叶斯估计为基础的 VAR (向量自回归) 模型和 DSGE (动态随机一般均衡) 模型的前沿研究,就经常出现在主流经济学期刊上。因此,掌握贝叶斯方法逐渐成为了宏观计量研究者的必修课。
遗憾的是,一般对 VAR 模型和 DSGE 模型的贝叶斯估计都是在 MATLAB 中完成。当然,我们也可以利用 Stata 内置的贝叶斯估计命令处理一些较简单的模型。
基于上述考虑,本文将简要介绍贝叶斯估计和 MCMC 方法的基本概念、如何在 Stata 中使用 bayes
命令前缀、标准估计命令 bayesmh
及其自带的交互菜单。
在统计学和计量经济学中,参数估计都扮演着非常重要的角色。频率学派认为概率是在长期重复的随机实验中产生的性质。但在贝叶斯学派看来,概率是不确定性的一种表达。这使得我们可以对无法重复的事件的概率进行考察,但相应的问题是将主观性引入了概率之中。反映到参数估计上,频率学派认为模型参数是一个定值,应当使用观测数据去接近它;而贝叶斯学派认为模型参数是一个随机变量,有自己的分布,应当使用先验信息并结合观测值去修正它的分布。
因此,贝叶斯估计的核心任务就是结合先验信息和观测值信息,获得模型参数的概率分布 (后验分布,posterior)。我们将所有观测值记为
在确定了先验分布和似然函数后,可以使用贝叶斯定理 (Bayes' Law) 求解出后验分布:
其中:
在实际的贝叶斯分析中,由于
在本小节,我们将介绍最简单的贝叶斯估计实例——伯努利实验。假设现有一枚硬币,将其抛起,然后估计正面朝上的概率
根据前面的讨论,我们需要明确先验分布和似然函数。由于似然函数不具有主观性,故先予以分析。假设抛了
再来看先验分布。在实际操作中,我们往往选择似然函数的共轭先验 (conjugate prior)。由于共轭先验具有与似然函数相乘后分布族保持不变的特征,我们可以通过简单的数学推导求解出
在选定了先验分布的分布族之后,还需要调整分布的参数以符合先验信念。Beta 分布有两个参数
在 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
------------------------------------------------------------------------
可以看出,当设定先验分布为
上例中,核心代码如下:
. bayesmh y, likelihood(dbernoulli({theta})) prior({theta}, beta(1,1))
其中,likelihood()
指定了似然函数的形式,在本例中似然函数由伯努利分布导出,同时我们将参数称为 theta (在本例中只有一个参数,因此不需要进行区分,只是单纯的命名,在更复杂的例子中需要通过名字区分参数)。prior()
则指定了先验分布,本例中将 theta 的先验分布定为
上述内容的编程实现起来并不困难,但是在更复杂的问题中,根据似然函数和先验分布的不同,很多时候需要指定额外的参数或选项。此时,贝叶斯估计命令 bayesmh
并不适合,需要通过 Stata 自带的交互菜单生成。具体来看:
我们可以通过点击 Stata 上方的 Statistics 菜单,选择底部的 Bayesian analysis 选项,再根据具体问题选择相应选项,打开贝叶斯估计的交互界面。由于上文的例子不涉及回归分析,选择 General estimation and regression。
由于模型是一个简单的独立实验,每一个观测值都是一次独立实验的结果,因此在 Syntax 处选择 Univariate distributions,Dependent variable 处选择生成的变量 y。Distribution 用于确定似然函数的分布,在本例中选择伯努利分布,同时指定参数名为 theta。最后需要指定模型的先验分布,首先点击 Create 按钮,然后在弹出的菜单中选择相应的分布族和参数,便可以生成一个想要的先验分布。最终结果如下:
点击 OK 键,系统便会自动开始估计,可以得到前文中的结果。
为了说明交互界面的便利性,我们考虑一个稍微复杂一些的例子。假设样本不再是由伯努利分布生成,而是由二项分布
在这个例子中,原分布和似然函数中有两个参数
. 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 的输入框,也就是上文指定的
不仅如此,我们还可以通过键入命令 db postest
打开后验分析的交互菜单进行分析。常用的功能如 Bayesian analysis 下的 Graphical summaryies and convergence diagnostics 菜单可以绘制很多与后验分布相关的图,Hypothesis testing using model posterior probabilities 和 Interval hypothesis testing 菜单可以利用生成的后验分布进行假设检验。
在前面的例子中,可能会有读者感到疑惑,当先验分布选为共轭先验时,后验分布是有解析表达式的,这意味着其概率密度函数应该是确定的。但在结果显示的界面中,明显有模拟 (simulation)、迭代 (iteration) 和抽样 (sampling) 这样的常见于数值方法的字眼。这其实是因为 Stata 并非使用解析方法对后验分布进行求解。为了进一步深入了解现代贝叶斯方法的操作流程,我们有必要对 MCMC (Markov chain Monte Carlo,马尔科夫链蒙特卡洛) 方法进行介绍。
对于 MCMC 方法,维基百科定义为:马尔可夫链蒙特卡洛方法包括一类用于从概率分布中采样的算法。通过构建具有所需分布作为其平衡分布的马尔可夫链,可以通过记录链中的状态来获得所需分布的样本。
这一段话有些抽象,实际上从应用层面上来说,我们只需要了解两个最常用的 MCMC 方法:Metropolis–Hastings 算法 (简称 M-H 算法) 和吉布斯采样 (Gibbs sampling)。前者是 Stata 内置的对后验分布进行抽样的算法,其主要作用是通过后验分布的概率密度之比构造了一个具有接受-拒绝机制的马尔科夫链,并通过它抽取出所需的后验分布。这一算法也广泛应用于 DSGE 模型的贝叶斯估计中,但它不是本文的重点。
在宏观计量经济学中最常用的 (主要是因为其应用于 VAR 模型的贝叶斯估计) MCMC 方法是吉布斯采样。为了让大家对吉布斯采样的应用背景有一个简单的了解,下面将以线性回归为例进行介绍。
假设有如下的线性回归模型:
其中,
这里有两个待估参数,即系数矩阵
在第一种情况中,由于扰动项方差已知,此时的共轭先验是正态分布,因此假定系数矩阵
先验分布:
似然函数:
后验分布:
可以发现后验分布仍然是正态分布:
在第二种情况中,由于似然函数中的给定变量改变,其共轭先验也随之发生变化。查表知此时扰动项方差
在实际中,更加常见的是第三种情况,即两个待估参数都未知。此时,没有办法获得单一的共轭先验,一种解决办法是设定联合先验分布:
其中,
在这样的设定下,联合先验分布还是一个正态分布,和似然函数的分布族一致,这种情况称之为自然共轭先验 (natural conjugate prior)。
需要注意的是,最后解出的后验分布是一个联合后验分布
吉布斯采样的思想其实很简单,当有了两个能条件分布的抽取器时 (即第一种和第二种情况),可以通过迭代的方法获得逼近边缘分布的抽样。在本例中,具体操作为:
这一操作步骤看似复杂,但在 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.
Note:产生如下推文列表的命令为:
lianxh VAR DSGE 分布, m
安装最新版lianxh
命令:
ssc install lianxh, replace
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟,课程主页 Stata 33 讲 - 连玉君, 每讲 15 分钟. Stata 小白的取经之路 - 龙志能,时长:2 小时,课程主页 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh