DSGE笔记系列4：Dynare基本应用(4)

$\text{DSGE笔记系列 4}$

此文内容为《动态随机一般均衡(DSGE)模型》的笔记，李向阳老师著，清华大学出版社出版。

我只将个人会用到的知识作了笔记，并对教材较难理解的部分做了进一步阐述。为了更易于理解，我还对教材上的一些部分（包括代码和正文）做了修改。

仅供学习参考，请勿转载，侵删！

$\S \text{ 第3章 } \S$

$\text{Dynare基本应用(4)}$

上一篇有关Dyanre的文章为大家介绍了有关确定性模拟的技术。在本次推文，我将介绍随机模拟技术。另外，我还会介绍有关参数估计的问题。

3.9 随机模拟分析：`stock_simul`

随机模拟是Dyanre最常用的计算命令之一。随机模拟首先使用扰动算法完成模型的求解任务，然后再求解的基础上计算内生变量的脉冲响应、各阶矩，最后进行随机模拟，——即根据随机抽取的外生冲击的样本，模拟出内生变量的增长路径。

3.9.1 随机模拟命令简介

在Dynare中，随机模拟的选项多达44个，常用的有：

drop = xx，计算内生变量的各阶矩时，需要首先指定丢弃的期数，默认为100
hp_filter = xx，计算内生变量的各阶矩时，首先使用HP滤波，xx代表 $\lambda$ 值
irf=xx，计算xx期的脉冲响应，默认为40
irf_shocks=(shock1, shock2,...)，针对给定的外生冲击计算脉冲响应，默认计算全部
irf_plot_threshold=xx，指定画脉冲响应图形的门槛大小，默认 $10^{-10}$
nocorr，指定不要在Matlab窗口显示相关系数矩阵
nofunctions，指定不要在Matlab窗口显示政策函数和转换方程
nomoments，指定不要在Matlab命令窗口显示各阶矩
nograph，指定不要显示和保存图形
nodisplay，指定不要显示图形，但保存到硬盘
noprint，指定不要显示任何计算信息，对循环很有帮助
order=xx，指定Taylor近似对阶数，可为1``2``3，默认为2
periods=xx，指定随机模拟的期数，同时也是决定计算理论矩还是模拟矩的关键所在
pruning，指定剪枝算法，即在迭代过程中丢弃高阶项

3.9.2 随机模拟示例

以第三章第一节的RBC模型为例，这里示例如何使用模拟结果中的水平变量进行做图分析：

/*
 * 《动态随机一般均衡(DSGE)模型》教材
 * 随机模拟示例
 * 2020年3月30日
 * @爱吃汉堡薯条
 * @侵删
 */
%----- 源代码9 -----%
var c k lab z I;    % 内生变量声明
varexo e;    % 外生变量声明
parameters beta theta delta alpha tau rho s;    % 参数声明
 
% 参数赋值
beta = 0.987;
theta = 0.357;
delta = 0.012;
alpha = 0.4;
tau = 2;
rho = 0.95;
s = 0.01;   

%----- 源代码10 -----%
model;
% 定义模型的“局部变量”，以提高编程的便利性
# mar_c = (c^theta * (1-lab)^(1-theta) ) ^ (1-tau);    % 局部变量的声明使用“#”命令
# mar_c1 = (c(+1)^theta * (1-lab(+1))^(1-theta) ) ^ (1-tau);

% (1) 欧拉方程
[name = 'Euler equation']    % 这是系统标签(tag)，可以对每个均衡条件进行标识
mar_c/c = beta * (mar_c1/c(+1)) * (1 + alpha*exp(z)*k^(alpha-1)*lab(+1)^(1-alpha) - delta);

% (2) 劳动供给方程
[name = 'Labor supply']
((1-theta)/theta) * (c/(1-lab)) = (1-alpha)*exp(z)*k(-1)^alpha*lab^(-alpha);

% (3) 资源约束方程
[name = 'Resource constraint']
c + k - (1-delta)*k(-1) = exp(z)*k(-1)^alpha*lab^(1-alpha);

% (4) 资本积累方程
[name = 'Capital accumulation']
I = k - (1-delta)*k(-1);

% (5) 外生冲击AR(1)过程
[name = 'Technology shock']
z = rho*z(-1) + s*e;
end;

%----- 源代码11 -----%
initval;    % 初值模块，设定内生变量和外生变量的初始值
k = 0.5;
c = 0.5;
lab = 0.5;    
z = 0;
I = 0.5;
e = 0;
end;

shocks;    % 外生冲击设定，直接设定外生冲击的方差
var e = s^2;
end;

steady;    % 计算稳态
resid;    % 给定内生变量的取值（通常为稳态），计算均衡条件对应的静态方程的残差

%----- 源代码12 -----%
% 如果不指定模拟的次数“periods”，Dynare就不会进行模拟
stoch_simul(periods=1000, irf=400, order=1);

% 将模拟结果保存为“simudata.mat”文件
dynasave('simudata.mat');

则有下图：

在这里，随机模拟命令指定了3个选项：

模拟1000期
计算40期脉冲响应图
进行二阶泰勒近似求解

前面已经提到，内生变量的模拟路径结果储存于oo_.endo_simul结构数组中。其行的顺序对应着以声明顺序排序的内生变量，其列则为模拟期数。

其中，技术冲击的均值为：

>> mean(oo_.endo_simul(4,:))
ans = -7.2640e-05

3.10 参数估计简介

模型中结构参数的赋值方法有两种，分别是：

校准（Calibration），直接赋值
估计（Estimation），使用实际统计或观测数据进行估计

简单来说，Dynare对参数估计就是使用某些方法，比如贝叶斯估计、极大似然估计，从数据中发现真实取值

3.10.1 极大似然估计与贝叶斯估计的基本逻辑

本文仅简单阐述极大似然估计和贝叶斯估计的基本原理和逻辑，对背后的复杂数理知识，请参考教材和博客。

3.10.1.1 极大似然估计

极大似然估计（Maximum Likelihood Estimation, MLE）就是选择参数 $\theta$ 来最大化如下的似然函数（Likelihood），其原理比较简单：
$\begin{split} p(Y^{data}|\theta) &\equiv p(Y^{data}_1,Y^{data}_2,\dots,Y^{data}_T|\theta)\\ &\mathop{=}\limits_{条件概率} p(Y_1^{data}|\theta) \times p(Y_2^{data}|\theta)\\ &\times \dots \times p(Y_t^{data}|\theta,Y_1^{data},\dots,Y_{t-1}^{data})\\ &\times \dots \times p(Y_T^{data}|\theta,Y_1^{data},\dots,Y_{T-1}^{data}) \end{split}$
其中，条件概率 $p(Y_t^{data}|\theta,Y_1^{data},\dots,Y_{t-1}^{data})$ 是使用 $\text{Kalman}$ 滤波进行计算的。

总而言之，它是选择 $\theta$ 使已经发生时间的概率最大。

3.10.1.2 贝叶斯估计

贝叶斯估计（Bayesian Estimation）的基本原理非常直观，就是在参数先验分布（Prior）的基础上，结合数据的信息（ $Y^{data}$ ），找到参数的后验分布（Posterior）。也就是说，贝叶斯估计是一种简单的映射，依据数据信息，将先验分布映射为后验分布。

由于先验分布是依据经验而假设的分布，可能存在一定的偏误，而数据的作用其实是将数据中有用的信息融入，以纠正偏误并形成后验分布。

从逻辑上说，贝叶斯估计的背后就是贝叶斯法则，即参数的后验分布由条件概率公式：
$p(Y^{data}|\theta) \times p(\theta) = p(Y^{data},\theta)=p(\theta|Y^{data})\times p(Y^{data}) \tag{3.10.1}$
推导出：
$p(\theta|Y^{data}) = \frac{p(Y^{data}|\theta)\times p(\theta)}{p(Y^{data})} \tag{3.10.2}$
从数据使用和先验假设的角度看，贝叶斯估计位于极大似然估计和参数校准之间：

极大似然估计：直接使用数据中的信息，没有先验
贝叶斯估计：即使用数据中的信息，也具有先验信息
校准：仅使用先验信息，没有使用数据的信息

贝叶斯估计的先验分布是一般意义上的分布，即标准差不为0，不是一个点，而是一个区域。严格地说，是整个参数空间的一个子区域，作为候选区域。先验分布可以理解为施加于似然函数的“权重”或者“偏好”，使得参数的选择范围集中于均值附近。

换句话说，先验分布使得“巨大的参数选择空间”变小，使似然函数集中于“更有意义”的子区域去寻找合适的参数估计值。Dynare的贝叶斯分布可分为两步：

使用数值求解算法，计算似然函数取最大值的众数（Model）
从众数开始，使用MCMC方法模拟后验分布，然后计算各种感兴趣的矩（Moments）

3.10.2 马尔可夫链——蒙特卡洛（MCMC）算法

由于参数的后验分布（密度函数）在大多数DSGE模型中（甚至某些RBC模型中）无法使用解析式的形式来表达，因此不能通过贝叶斯法则计算出来，只能通过随机抽样的方法，即蒙特卡洛算法（Monte Carlo）来近似逼近后验分布。

马尔可夫链-蒙特卡洛方法（Markov Chain-Monte Carlo, MCMC）是一类方法的统称，是贝叶斯估计的逻辑方法和框架，非具体估计的算法。可以参考知乎网友的文章。实际上我们不太需要了解背后的数学原理。

关于本节，我们最重要的是基本了解Dynare贝叶斯估计的核心逻辑，如图3.12所示：

3.10.3 Dynare参数估计和一个例子

我们着重于介绍如何用Dynare内置的估计命令，对一个或多个参数进行编程估计。

3.10.3.1 Dynare参数估计的基本逻辑

Dynare参数估计是使用以下模型的一阶线性化系统进行估计：
$E_t{f(y_{t+1},y_t,y_{t-1},u_t;\theta)}=0$
其中， $\theta$ 是待估计参数向量； $u_t$ 是外生冲击； $y_t$ 是内生变量组成的向量。

Dynare求解参数估计的基本逻辑是：

计算内生变量 $y_t$ 的稳态值 $y$
线性化系统均衡条件
求解线性化模型
使用Kalman滤波算法计算对数似然函数（Log-Likelihood）
找到极大似然值分布的众数（Model）
使用Metropolis-Hastings算法模拟后验分布

3.10.3.2 Dynare参数估计的语法

Dyanre参数估计的基本语法分成三部分：

声明参数估计模块
待估参数本身的声明
使用estimation命令

Dynare中声明估计参数模块的内置命令是以estimated_params;开头，以end;结束，中间输入待估计参数声明。待估计参数有两种常见的类型：

结构参数
外生参数的标准差参数

具体实例如下：

estimated_params;    % 结构参数的估计
gamma, normal_pdf, 1,0.05;
alpha, uniform_pdf, , , 0, 1;    % uniform分布需要留逗号，见下文
end;

estimated_params;    % 外生冲击标准差参数的估计
stderr e, inv_gamma_pdf, 0.01, inf;
end;

结构参数的声明是：

以待估计参数的名称开头
然后是该参数服从的先验分布的名称，该名称由Dynare预先定义，如表3.11所示：

先验分布名称	符号形式	分布取值范围
`beta_pdf`	$B(\mu,\sigma,p_3,p_4)$	$[p_3,p_4],\quad p_3>0,p_4>0$
`gamma_pdf`	$\Gamma\left(\mu, \sigma, p_{3}\right)$	$\left(p_{3}, \infty\right), \quad p_{3}>0$
`inv_gamma_pdf`	$I \Gamma(\mu, \sigma)$	$(0,+\infty)$
`normal_pdf`	$N(\mu,\sigma)$	$(-\infty,+\infty)$
`uniform_pdf`	$U(p_3,p_4)$	$(p_3,p_4)$

其中， $\mu$ 和 $\sigma$ 分别表示均值和标准差；参数 $p_3$ 、 $p_4$ 仅用于Beta分布、 Gamma分布和Uniform分布，可以为无穷大。此外，Uniform分布因为不需要声明均值和标准差，因此需要在参数声明中留出两个逗号隔开的空格，然后声明 $p_3$ 、 $p_4$ 参数。

标准差参数的声明：

需要使用关键字stderr后跟外生冲击的名称，而不是外生冲击对应的标准差的名称
然后是指定其服从的先验分布，声明语法与前面的相同

3.10.3.2 注意事项

对于先验分布的选择，确实是一个困难的问题。选择先验分布具有主观性，这成为贝叶斯分布受质疑和批评的主要原因。一般参考：

对于标准差 $\sigma$ 的估计，通常选择Inverse Gamma 分布
对于 $\text{AR(1)}$ 过程的滞后参数 $\rho$ 通常选择Beta分布
其他参数，大多使用Normal分布

在估计之前，还要注意以下几个问题：

数据问题：估计参数的前提是模型中某些内生变量是可以观测到或者由观测值经过处理得到，且观测变量的个数不能超过模型中外生冲击的个数，否则无法计算模型的似然函数。
观测数据与内生变量的匹配：在估计时，要求模型文件的变量形式和观测数据相匹配。

观测变量的声明问题：内生变量如果声明为可观测变量，则需要使用以下命令

% 语法
varobs [VARIABLE_NAME];
% 例如
varobs pi tau a;    % 声明三个内生变量为观测变量

3.10.3.3 estimation命令介绍

Dynare会在estimation命令中调取外部数据（使用datafile选项指定外部数据文件），并加载到内存中去。观测变量声明命令varobs会告诉Dynare去寻找变量名为pi、tau、a的三个序列，如果找不到，会提示错误。

Estimation`命令及其选项问题：这里仅介绍几个最重要的关键选项

datafile=xxx：可选择Matlab的.m文件、.mat矩阵或者.csv文件等
conf_sig=xxx：定义置信区间的宽度（如xxx=.95定义了“95%的置信区间”）
frst_obs=xxx：定义由load_Gali_est_data.m加载的数据，第一个开始用于对估计的数据位置，xxx必须是一个正整数，缺省为1
nobs=yy：Dynare估计中使用的观测变量序列长度（t=xxx到t=yy+xxx-1，同样的yy必须是正整数，缺省为全部观测变量
mode_check：画出后验分布众数附近的后验分布密度图
mode_compute=xxx：选择优化算法，取值0~10，缺省使用4（该命令在不断更新中）
forecast=n：在观测变量后一期之后预测 $n$ 期值
smmother：计算内生变量和外生冲击的后验分布的相关矩
mh_replic=xxx：定义计算后验分布时MCMC重复的次数，缺省为2000
mh_jscale=xxx：定义MCMC中的乘法因子。理想情况下应该选择使“接受率”落在25%~33%之间。该参数需要多次调试后获得比较合适的值，缺省为0.2
mh_nblocks=xxx：定义Metropolis-Hastings算法平行进行的马尔可夫链的个数，缺省为2
filtered_vars：计算滤波变量，结果存储在oo_.FilteredVarialbes中
moments_varendo：计算内生变量的后验分布的理论矩，结果存储在oo_.PosteriorTheoreticalMoments中

3.10.3.4 极大似然估计的例子

考虑一个新凯恩斯模型，不含投资但含粘性价格和产出缺口。

（1）模型
$U\left(C_{i}, N_{t}\right) \equiv \frac{C_{t}^{1-\sigma}-1}{1-\sigma}-\exp \left(\tau_{t}\right) \frac{N_{t}^{1+\varphi}}{1+\varphi}\tag{3.10.10}$

其中：

$\sigma>0$ 为消费跨期替代弹性的倒数
$\tau_t$ 为劳动负效用的外生冲击，服从 $\text{AR(1)}$ 过程：$\tau_{t}=\lambda \tau_{t-1}+\epsilon_{t}^{\tau}, \epsilon_{t}^{\tau} \sim i.i.Dá
$\phi>0$ 为劳动供给的 $\text{Frisch}$ 弹性的倒数

首先是新凯恩斯IS曲线（NKIS）和新凯恩斯菲利普斯曲线（NKPC）：
$x_{t}=E_{t}\left(x_{t+1}\right)-\frac{1}{\sigma}\left(i_{t}-E_{t}\left(\pi_{t+1}\right)-i_{t}^{*}\right)\tag{3.10.11}$

$\pi_{t}=\beta E_{t} \pi_{t+1}+\kappa x_{t}\tag{3.10.12}$

其次是自然水平产出和自然利率方程：
$y_{t}^{*}=\frac{1+\varphi}{\sigma+\varphi} a_{t}-\frac{1}{\sigma+\varphi} \tau_{t}\tag{3.10.13}$

$i_{t}^{*}=\sigma \frac{1+\varphi}{\sigma+\varphi}(\rho-1) a_{t}+\sigma \frac{1-\lambda}{\sigma+\varphi} \tau_{t}\tag{3.10.14}$

假设货币冲击政策为Taylor规则：
$i_{t}=\rho_{i} i_{t-1}+\left(1-\rho_{i}\right)\left(\phi_{\pi} \pi_{t}+\phi_{x} x_{t}\right)+\epsilon_{t}^{i}, \epsilon_{t}^{i} \sim \text { i.i.d }\tag{3.10.15}$
从而，模型由 $\{ x_t. i_t. i_t^*, \pi_t. a_t,\tau_t \}$ 和六个均衡条件：

NKIS曲线 $(3.10.11)$
NKPC曲线 $(3.10.12)$
自然利率方程 $(3.10.13)$
Taylor规则 $(3.10.15)$
技术冲击的 $\text{AR(1)}$ 过程
偏好冲击的 $\text{AR(1)}$ 过程

接下来，对上述模型文件进行随机模拟，模拟5000个样本，作为下一步估计的数据。然后使用上述模拟的样本作为观测数据，对模型中3个结构参数进行估计，如表3.12所示，以示例如何在Dynare中分别使用极大似然估计和贝叶斯估计估计此三个参数：

待估计参数	含义
$\rho=0.8$	技术冲击 $\text{AR(1)}$ 过程的持续参数
$\lambda=0.5$	劳动负效用外生冲击 $\text{AR(1)}$ 过程的持续参数
$\sigma_a=0.02$	技术冲击 $\text{AR(1)}$ 过程的外生冲击的标准差

随机模拟的代码如下（请首先运行）：

// 生成用于参数估计的随你模拟序列
//Written by Xiangyang Li@SCC


//declaration of the endogenous variables, log-devation form
var a // technology
pi      //CPI inflation
i       //nominal rate
istar //nature real rate
tau //labor disutility shock
x    //output gap = output - natural output
dy  //growth rate of output
;      

varexo  eps_a eps_tau eps_i;    

parameters 
sigma phi beta kappa phi_x phi_pi rho lambda rhoi theta; 

// Parameter Values  
sigma  = 1; //inverse of inter-temporal elasticity of substitution of consumption
beta    = 0.97; //discount factor of households
phi_x   = .15;  //coefficient of output gap in monetary policy rule
phi_pi = 1.5; //coefficient of inflation in monetary policy rule
rhoi   = 0.8;  //persistence of nominal interest rate;
rho     = 0.8; //persistence of technology
lambda  = 0.5; //persistence of labor disutility shock;
phi     = 1;       //inverse of Frisch elasticity of labor supply
theta   = 0.75;  //Calvo stickiness parameter
kappa   = ((1-theta)*(1-beta*theta))/theta; //simplifying parameter

// The model in log-linear
model(linear);
//(1) Philiphs Curve -  Calvo Pricing Equation
beta*pi(+1) + kappa*x = pi;         

//(2) New Keynesian IS curve
sigma*(i - pi(+1)-istar)   = x(+1) - x; 

//(3) natural real rate definition
istar = sigma*(1+phi)*(rho-1)/(sigma+phi)*a + sigma*(1-lambda)/(sigma+phi)*tau;

//(4)Taylor Rule
i= rhoi*i(-1) + (1-rhoi)*(phi_pi*pi + phi_x*x) +eps_i; 
    
//(5)Technology shock
a = rho*a(-1) + eps_a;

//(6) Labor disutility shock
tau = lambda*tau(-1) + eps_tau;

//(7) Growth rate of output - by output gap
dy  = x - x(-1) + (1+phi)/(sigma+phi)*(a-a(-1)) - (tau - tau(-1))/(sigma+phi);                    // construction of output growth

end;


shocks;
var eps_a;
stderr .02;
var eps_tau;
stderr .02;
var eps_i;
stderr 0.02;
end;

set_dynare_seed=1;
stoch_simul(periods=5000, irf=7, nograph);

gcsim = oo_.endo_simul;

save data gcsim;

（2）最大似然估计的实践

在下面的代码中，使用了estimated_param;设定了初始值，使用estimated_param_bounds;设定了上下边界，此外使用varobs定义了3个观测变量（请在运行3.10.3.4的（1）后运行）：

// Maximum Likelihood estimation of the Gali-Christiano Model
//we are going to estimate two persistence parameters, rho and lambda
// and one standard devation parameters 
//Written by Xiangyang Li@SCC@2016.12
//Modified by Xiangyang Li@2017.10.12

//declaration of the endogenous variables, log-devation form
var a // technology
pi      //CPI inflation
i       //nominal rate
istar //nature real rate
tau //labor disutility shock
x    //output gap = output - natural output
dy  //growth rate of output
;      

varexo  eps_a eps_tau eps_i;    

parameters sigma phi beta kappa phi_x phi_pi  rhoi theta; 

parameters rho lambda; //to be estimated

// Parameter Values  
sigma  = 1; //inverse of inter-temporal elasticity of substitution of consumption
beta    = 0.99; //discount factor of households
phi_x   = .15;  //coefficient of output gap in monetary policy rule
phi_pi = 1.5; //coefficient of inflation in monetary policy rule
rhoi   = 0.8;  //persistence of nominal interest rate;
//rho     = 0.8; //persistence of technology, to be estimated
//lambda  = 0.5; //persistence of labor disutility shock,to be estimated;
phi     = 1;       //inverse of Frisch elasticity of labor supply
theta   = 0.75;  //Calvo stickiness parameter
kappa   = ((1-theta)*(1-beta*theta))/theta; //simplifying parameter

// The model in log-linear
model(linear);
//(1) Philiphs Curve -  Calvo Pricing Equation
beta*pi(+1) + kappa*x = pi;         

//(2) New Keynesian IS curve
sigma*(i - pi(+1)-istar)   = x(+1) - x; 

//(3) natural real rate definition
istar = sigma*(1+phi)*(rho-1)/(sigma+phi)*a + sigma*(1-lambda)/(sigma+phi)*tau;

//(4)Taylor Rule
i= rhoi*i(-1) + (1-rhoi)*(phi_pi*pi + phi_x*x) +eps_i; 
    
//(5)Technology shock
a = rho*a(-1) + eps_a;

//(6) Labor disutility shock
tau = lambda*tau(-1) + eps_tau;

//(7) Growth rate of output - by output gap
dy  = x - x(-1) + (1+phi)/(sigma+phi)*(a-a(-1)) - (tau - tau(-1))/(sigma+phi);                    // construction of output growth

end;

// 设置冲击
shocks;
var eps_tau;    stderr 0.01;
var eps_i;      stderr 0.01;
end;

// 极大似然估计fa:
// 初始化参数值，用于优化
estimated_params;
   stderr eps_a, 0.01;    // 外生冲击参数的估计声明
   rho, .80;    // 结构参数的估计声明
   lambda, .50;  
end;

// 设置参数的边界
estimated_params_bounds;
stderr eps_a, 0.001, .2;
rho, .001,.95;
lambda, .001,.95;
end;


// 用于估计参数的观测源：观测变量
varobs pi tau a; 

// 开始估计
estimation(datafile=load_Gali_est_data,conf_sig =.95,first_obs=101,forecast =46,nobs=4000,mode_check,mode_compute=10) a pi i istar tau x dy;

运行上面的代码，意味着：

使用95%置信区间
从第101个数据开始作为观测值
在观测变量以后预测46期
使用的序列长度是4000期
画出后验分布众数附近的后验分布密度图
使用第10种优化算法（见Manual）

（3）关于观测变量的分析

在默认情况下，极大似然估计的结果令人满意，非常接近真实值；但选择不同的观测变量，会影响参数的点估计，如表3.13所示：

直观上，选择观测数据时，要尽量选择和待估计参数关系紧密的数据序列，包括使用尽可能多的数据项，以获得准确的信息。

（4）关于序列长度的分析

同样，观测变量序列的长度也会影响估计。可以看出，序列长度缩小100倍时，标准差约增大十倍。除了参数 $\lambda$ 以外，另外两个参数的估计值变化很小，说明极大似然估计具有一定的稳健型：

3.10.3.5 贝叶斯估计的例子

使用贝叶斯分布，要先指定先验分布。这里设定：

技术标准差参 $\sigma$ 数的均值和真实值相同，标准差为10
$\rho$ 和 $\lambda$ 两个参数的均值为真实值，标准差都为0.02
还设定了3个参数的上下界和初始值

如下所示（请在运行3.10.3.4的（1）后运行）：

// Maximum Likelihood estimation of the Gali-Christiano Model
//we are going to estimate two persistence parameters, rho and lambda
// and one standard devation parameters 
//Written by Xiangyang Li@SCC@2016.12
//Modified by Xiangyang Li@2017.10.12

//declaration of the endogenous variables, log-devation form
var a // technology
pi      //CPI inflation
i       //nominal rate
istar //nature real rate
tau //labor disutility shock
x    //output gap = output - natural output
dy  //growth rate of output
;      

varexo  eps_a eps_tau eps_i;    

parameters sigma phi beta kappa phi_x phi_pi  rhoi theta; 

parameters rho lambda; //to be estimated

// Parameter Values  
sigma  = 1; //inverse of inter-temporal elasticity of substitution of consumption
beta    = 0.99; //discount factor of households
phi_x   = .15;  //coefficient of output gap in monetary policy rule
phi_pi = 1.5; //coefficient of inflation in monetary policy rule
rhoi   = 0.8;  //persistence of nominal interest rate;
//rho     = 0.8; //persistence of technology, to be estimated
//lambda  = 0.5; //persistence of labor disutility shock,to be estimated;
phi     = 1;       //inverse of Frisch elasticity of labor supply
theta   = 0.75;  //Calvo stickiness parameter
kappa   = ((1-theta)*(1-beta*theta))/theta; //simplifying parameter

// The model in log-linear
model(linear);
//(1) Philiphs Curve -  Calvo Pricing Equation
beta*pi(+1) + kappa*x = pi;         

//(2) New Keynesian IS curve
sigma*(i - pi(+1)-istar)   = x(+1) - x; 

//(3) natural real rate definition
istar = sigma*(1+phi)*(rho-1)/(sigma+phi)*a + sigma*(1-lambda)/(sigma+phi)*tau;

//(4)Taylor Rule
i= rhoi*i(-1) + (1-rhoi)*(phi_pi*pi + phi_x*x) +eps_i; 
    
//(5)Technology shock
a = rho*a(-1) + eps_a;

//(6) Labor disutility shock
tau = lambda*tau(-1) + eps_tau;

//(7) Growth rate of output - by output gap
dy  = x - x(-1) + (1+phi)/(sigma+phi)*(a-a(-1)) - (tau - tau(-1))/(sigma+phi);                    // construction of output growth

end;

shocks;
var eps_tau;
stderr 0.01;
var eps_i;
stderr 0.01;
end;

// 贝叶斯估计
// 设置先验分布
estimated_params;
   // eps_a的均值为0.02，标准差为10
   stderr eps_a, inv_gamma_pdf, 0.02, 10;
   // 设置另外两个参数的先验
   rho, beta_pdf, 0.80, 0.04;
   lambda, beta_pdf, 0.80, 0.04;  
end;

estimated_params_bounds;
stderr eps_a, 0.001, .2;
rho, .001,.95;
lambda, .001,.95;
end;

// 设置初值，MLE则不需要这一模块
estimated_params_init;
stderr eps_a, 0.019;
rho, 0.82;
lambda, 0.51;
end;

// 设置观测变量
varobs pi tau a; 

estimation(datafile=load_Gali_est_data,     // 数据文件
           conf_sig =.95,                                   // 95%置信区间
           first_obs=101,                                   // 从101位开始使用数据
           forecast =46,                                        // 在观测序列后预测46期
           nobs=40,                                                 // 使用40期数据进行预测
           mode_check,                                          // 画出后验分布分布密度图
           mode_compute=4,                                  // 使用第4种优化算法
           mh_replic=1200,                                  // MCMC重复次数
           mh_jscale=1.4,                                   // MCMC的乘法因子
           mh_nblocks=2)                                        // 马尔可夫链个数
           a pi i istar tau x dy;

运行上面的代码，可以得到几个图，其中下图：

可以直观看出贝叶斯估计的先验和后验的情况。

贝叶斯估计比极大似然估计更加稳健（请看教材），所以一般使用贝叶斯估计。
ng$

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 222,104评论 6赞 515
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 94,816评论 3赞 399
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 168,697评论 0赞 360
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 59,836评论 1赞 298
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 68,851评论 6赞 397
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 52,441评论 1赞 310
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,992评论 3赞 421
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 39,899评论 0赞 276
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 46,457评论 1赞 318
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 38,529评论 3赞 341
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 40,664评论 1赞 352
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 36,346评论 5赞 350
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 42,025评论 3赞 334
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 32,511评论 0赞 24
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 33,611评论 1赞 272
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 49,081评论 3赞 377
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 45,675评论 2赞 359