实验背景
我国货币政策担负着辅助国家经济方针实施，促进经济增长和维护经济金融稳定的作用。而近年来我国经济增速放缓，经济进入结构性调整阶段，货币政策任务转向去杠杆、遏制通胀。货币政策主要分为相机抉择和货币政策规则。由于Kydland和Prescott（1977）认为相机抉择存在时间不一致性、Svensson（1997）和Woodford（1999）认为相比较前瞻货币政策规则，相机抉择存在稳定性偏差（Stabilization bias），因此国内外学者公认的货币政策规则有助于提高货币政策有效性和可信度。价格型货币政策规则又称Taylor规则，是一种盯住通胀缺口和产出缺口的利率制度，而我国目前常用的货币政策工具MLF、SLF的政策目标均是调节利率，因此研究利率制度的最优化问题具有很强的现实意义。
除了Taylor规则以外，还有以货币数量作为中介指标的McCallum规则。虽然我国目前对于最优货币政策规则的研究还没有定论，但国内许多学者认为我国货币政策应偏向于价格型货币政策规则（金成晓、马丽娟，2011；卞志村、胡恒强，2015）。
DSGE模型是目前宏观经济研究最热门的模型，其起源于Kydland&Prescott(1982)的真实商业周期理论（RBC）。Gali & Gertler（1999）引入了黏性价格等凯恩斯学派假定建立起NK-DSGE模型。DSGE的广泛应用归功于 Smets&Wouters（2002）的贡献，其建立的SW2002模型拟合效果超过了传统BVAR计量方法，验证了DSGE的存在价值。
大致说来，DSGE的建模过程如图所示：
其中，模型构建是严格按照经济学假定和理论，结合各国实际情况得到的。主要涵盖家庭、厂商、政府、金融机构、外国厂商等部分。由于本文着重研究本国货币政策，因此将模型简化为四部门：家庭、厂商、政府、央行。构建出的模型需要通过对数线性化将所有名义变量转变为与稳态值的偏差。所有稳态值作外生参数给定，冲击转化为结构冲击，DSGE转化为结构方程模型，即非线性理性预期系统。内生变量个数等于结构方程个数，因此不需要考虑求解时奇异值分解的问题。
因此，本文基于NK-DSGE-BVAR模型，讨论具有动态一致性，不存在稳定性偏误的前瞻Taylor规则的模拟分析。
模型构建
根据Negro and Schorfheide（2004）构建的新凯恩斯框架DSGE模型，简化后模型共有三个部门：代表性家庭、连续垄断竞争厂商和中央银行。
家庭部门
家庭部门采用CRRA消费效用，家庭即期效用与真实现金余额$M/P$、习惯存量消费$ C/A $和劳动时间$ h $有关。假定习惯存量等于技术水平$ A $，这项假定能够保证经济沿平衡增长路径增长。于是我们可以得到t期的家庭的终生效用为：
$$
E_{t}[\sum\limits_{s=t}^{\infty}\beta^{s-t}(\frac{(C_{s}/A_{t})^{1-\tau}-1}{1-\tau}+\chi ln\frac{M_{s}}{P_{s}}-h_{s})]
$$
其中内生变量为，$C_{t}$为当期消费，$A_{t}$为当期技术水平，$M_{t}$为当期货币量，$h_{t}$为当期劳动时间（小时数），$P_{t}$为当期物价。假定家庭是价格的接受者，定义通货膨胀率$\pi_{t}=P_{t}/P_{t-1}$。
外生参数包括：$\beta$为折现率，$\tau$为相对风险厌恶系数，$\chi$为比例系数，$E_{t}$为当期期望算子。
另外，家庭部门的总预算为当期工资收入，当期物价下的现金余额，加上债券收益，加上公司净利润分红。所有预算用于当期消费、债券投资、纳税和现金余额。于是家庭预算约束如下所示：
$$
C_{t}+\frac{B_{t}}{P_{t}}+\frac{M_{t}}{P_{t}}+\frac{T_{t}}{P_{t}}=W_{t}h_{t}+\frac{M_{t-1}}{P_{t}}+R_{t-1}\frac{B_{t-1}}{P_{t}}+D_{t}
$$
其中，上文没涉及的内生变量为，$B_{t}$是当期持有的政府债券净额，$R_{t}$表示当期利率，$T_{t}$为当期税收，$W_{t}$表示当期工资水平，$D_{t}$为当期公司净利润分红。
除此之外，家庭预算还受非蓬齐博弈（Non-Ponzi-Schemes）横截面条件的限制，以确保无穷远期的现金余额现值等于0。即
$$
\lim\limits_{t\to\infty}\sum\limits_{s=0}^{t}\beta^{s}M_{t+s}=0
$$
产品部门
产品部门是以连续统存在的许多垄断竞争厂商，不同厂商生产的产品具有完全异质性。假设产品的需求曲线来源于Dixit-Stiglitz偏好：
$$
P_{t}(j)=(\frac{X_{t}(j)}{X_{t}})^{-1/\nu}P_{t}
$$
其中内生变量$X_{t}(j)$和$P_{t}(j)$分别表示厂商j生产的产量和与之对应的价格，这个价格是由最大化利润得到。$X_{t}$和$P_{t}$分别表示加总的商品需求和加总的商品价格，这两个变量不是单个厂商可以改变的。外生参数$\nu$表示两种不同商品的固定替代弹性。
单个厂商虽然可以随意调整自己生产的产品的价格，但是如果其价格变化与社会的通货膨胀率不一致，则要承受一定的利润损失，这种损失被凯恩斯学者称为菜单成本（Menu costs）。其表达式为：
$$
MC_{t}=\frac{\phi}{2}(\frac{P_{t}(j)}{P_{t-1}(j)}-\pi^{\ast})X_{t}(j)
$$
其中，外生参数$\phi$表示经济的粘性价格程度。
对于单个厂商，我们采用简单的劳动生产函数。
$$
X_{t}(j)=A_{t}h_{t}(j)
$$
其中生产率即技术水平的变化由单位根过程给出，
$$
lnA_{t}=ln\lambda+lnA_{t-1}+\tilde{z_{t}}
$$
其中$\tilde{z_{t}}=\rho_{z}\tilde{z}_{t-1}+\epsilon_{z,t}$，$\epsilon_{z,t}$为高斯白噪声，表示技术冲击对所有企业的影响相同。
厂商通过最大化完整生命周期的利润来决定价格和劳动力投入，即最大化下式：
$$
\mathbb{E}_{t}[\sum\limits_{s=t}^{\infty}Q_{s}D_{s}(j)]
$$
其中$Q_{t}$为所有企业用于衡量未来收益大小的折现因子。这个因子随时间而变化，而由于状态依存（State-contingent）假定，在均衡条件下有：
$$
Q_{t+1}/Q_{t}=\beta(C_{t}/C_{t+1})^{\tau}
$$
这是因为家庭是企业股东，净利润决策应由家庭进行，因此其跨期替代率，即折现率如上所示。
$D_{t}(j)$表示j企业的利润，根据下式得到：
$$
D_{s}(j)=(\frac{P_{s}(j)}{P_{s}}X_{s}(j)-W_{s}h_{s}(j)-\frac{\phi}{2}(\frac{P_{s}(j)}{P_{s-1}(j)}-\pi^{\ast})^{2}X_{s}(j))
$$
中央银行
假设央行制定货币政策按照Taylor规则，则其利率应由下式决定：
$$
\frac{R_{t}}{R^{\ast}}=(\frac{R_{t-1}}{R^{\ast}})^{\rho_{R}}[(\frac{\pi_{t}}{\pi^{\ast}})^{\psi _{1}}(\frac{X_{t}}{X^{\ast}})^{\psi _{2}}]^{(1-\rho_{R})}e^{\epsilon_{R,t}}
$$
其中，$\rho_{R}$表示利率的平滑指数，$\epsilon_{R,t}$表示货币政策的预期外冲击。
前瞻Taylor规则在求解中给出。
政府部门
政府部门即财政部门，其财政支出主要是消费商品总额的一部分，其比例为$\zeta_{t}$，令$g_{t}=1/(1-\zeta_{t})$服从AR(1)过程
$$
\tilde{g}_{t}=\rho_{g}\tilde{g}_{t-1}+\epsilon_{g,t}
$$
其中$\epsilon_{g,t}$为政府消费的冲击。
政府收入由税收和债券融资获取，其中税收表示为$\frac{T_{t}}{P_{t}}$。另外加入铸币税考虑，于是有政府预算约束：
$$
\zeta_{t}X_{t}+R_{t-1}\frac{B_{t-1}}{P_{t}}+\frac{M_{t-1}}{P_{t}}=\frac{T_{t}}{P_{t}}+\frac{M_{t}}{P_{t}}+\frac{B_{t}}{P_{t}}
$$
模型求解
以上模型通过最大化条件即拉格朗日乘子法、变分法或Bellman方程方法可以化简为模型的均衡方程，也就是非线性理性预期差分方程(non-liner rational expectation difference equation)。再由对数线性化方法得到线性理性预期差分方程（LRED）。然后根据Sims(2002)，运用广义Schur分解法，最终得到线性差分方程(LDE)。
为什么要求解该模型？因为得到的求解结果是一个线性的决策方程组(Policy equations)，即当期的各个经济变量可以由过去的经济信息计算出最优的走向（或预测的走向）。
运用Uhlig的线性化方法，即例如对变量$Y_{t}$，取$\tilde{y}_{t}=lnY_{t}-lnY^{\ast}$，参考Sims(2002)和刘斌(2008)的求解过程，通过对消费、产出、工资的决策，最大化得到以下三个均衡条件：
$$
\tilde{x}_{t}=\mathbb{E}_{t}[\tilde{x}_{t+1}]-\tau^{-1}(\tilde{R}_{t}-\mathbb{E}_{t}[\tilde{\pi}_{t+1}])+(1-\rho_{g})\tilde{g}_{t}+\rho_{z}\frac{1}{\tau}\tilde{z}_{t}
$$
$$
\tilde{\pi}_{t}=\frac{\gamma}{r^{\ast}}\mathbb{E}_{t}[\tilde{x}_{t+1}]+\kappa[\tilde{x}_{t}-\tilde{g}_{t}]
$$
$$
\tilde{R}_{t}=\rho_{R}\tilde{R}_{t-1}+(1-\rho_{R})(\psi_{1}\tilde{\pi}_{t}+\psi_{2}\tilde{x}_{t})+\epsilon_{R,t}
$$
其中，(14)是消费的欧拉方程，即IS曲线，(15)为企业和政府的优化均衡，即菲利普斯曲线，又称总供给曲线AS(3)为货币政策规则。
前瞻Taylor规则下，(16)转化为
$$
\tilde{R}_{t+1}=\rho_{R}\tilde{R}_{t}+(1-\rho_{R})(\psi_{1}\mathbb{E}_{t}[\tilde{\pi}_{t+1}]+\psi_{2}\mathbb{E}_{t}[\tilde{x}_{t+1}])+\epsilon_{R,t}
$$
也就是说，根据模型的第一步求解，我们得到的线性理性预期系统（LRES）为：
$$
\tilde{x}_{t}=\mathbb{E}_{t}[\tilde{x}_{t+1}]-\tau^{-1}(\tilde{R}_{t}-\mathbb{E}_{t}[\tilde{\pi}_{t+1}])+(1-\rho_{g})\tilde{g}_{t}+\rho_{z}\frac{1}{\tau}\tilde{z}_{t}
$$
$$
\tilde{\pi}_{t}=\frac{\gamma}{r^{\ast}}\mathbb{E}_{t}[\tilde{x}_{t+1}]+\kappa[\tilde{x}_{t}-\tilde{g}_{t}]
$$
$$
\tilde{R}_{t+1}=\rho_{R}\tilde{R}_{t}+(1-\rho_{R})(\psi_{1}\mathbb{E}_{t}[\tilde{\pi}_{t+1}]+\psi_{2}\mathbb{E}_{t}[\tilde{x}_{t+1}])+\epsilon_{R,t}
$$
$$
\tilde{g}_{t}=\rho_{g}\tilde{g}_{t-1}+\epsilon_{g,t}
$$
$$
\tilde{z_{t}}=\rho_{z}\tilde{z}_{t-1}+\epsilon_{z,t}
$$
其中，内生变量为$y_{t}=[\tilde{x}_{t},\tilde{\pi}_{t},\tilde{R}_{t},\tilde{g}_{t},\tilde{z_{t}}]^{\top}$，结构冲击包括$v_{t}=[\epsilon_{R,t}\epsilon_{g,t},\epsilon_{z,t}]^{\top}$，结构参数包括
$$
\theta=[ln\gamma, ln\pi^{\ast}, ln r^{\ast}, \kappa, \tau, \psi_{1}, \psi_{2}, \rho_{R}, \rho_{g}, \rho_{z}, \sigma_{R}, \sigma_{g}, \sigma_{z}]^{\top}
$$
该线性理性预期系统由于预期算子未知，为了精确预测或制定政策，还需要运用Uhlig的Toolkit(1999)、Blachard-Kahn或Sims的gensys方法等，将LRES进一步化简为LDE。由于篇幅所限，本文仅使用最适合于BMR包的Uhlig方法，并简要介绍Sims方法。
首先，该模型可以化为以下形式：
$$
\Gamma_0 y_{t+1} =\Gamma_1 y_{t}+\Psi v_{t}+\Pi \eta_{t+1}
$$
其中，$y_{t}$即对数线性化的内生变量，$\eta_{t+1}=y_{t+1}-\mathbb{E}_{t}y_{t+1}$即为理性预期的偏离值。因此Sims方法的特点是隐含了期望算子。
具体地，
$$
\Gamma_0 =\begin{bmatrix}
1 &1/\tau & 0 &0 &0 \\ \gamma/r^ \ast & 0 &0 &0 &0 \\ (1-\rho_R)\psi_2 & (1-\rho_R)\psi_1 & -1 &0 &0 \\ 0 & 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -1 &
\end{bmatrix}
$$
$$
\Gamma_1 =\begin{bmatrix}
-1 & 0 & -1/\tau & (1-\rho_g) & \rho_z/\tau \\ \kappa & -1 & 0 & -\kappa & 0\\ 0 & 0 & \rho_R & 0 & 0\\ 0 & 0 & 0 & \rho_g & 0\\ 0 & 0 & 0 & 0 & \rho_z\end{bmatrix}
$$
$$
\Psi =\begin{bmatrix}0 &0 &0 \\ 0 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{bmatrix}
$$
$$
\Pi =\begin{bmatrix}-1 & 1/\tau\\ -\gamma/r^ \ast & 0\\ (\rho_R-1)\psi_2 & (\rho_R-1)\psi_1\\ 0 & 0\\ 0 & 0\end{bmatrix}
$$
由于MATLAB里的Dynare包支持Sims的解法，下面简单介绍一下Sims(2002)的解法。
矩阵$A$和$B$的广义Schur分解为：
$$
\Gamma_0 =QSZ^{H}, \Gamma_1 =QTZ^{H}
$$
其中$Q$和$Z$都是酉（Unitary，酉矩阵的共轭转置和它的逆矩阵相等）的，$S$和$T$都是上三角的，因此两边左乘$Q^{H}$可以得到：
$$
\begin{bmatrix}S_{11}&S_{12} \\ 0&S_{22} \end{bmatrix}\begin{bmatrix}\tilde{y}_{1t+1} \\ \tilde{y}_{2t+1}\end{bmatrix}=\begin{bmatrix}T_{11}&T_{12} \\ 0&T_{22} \end{bmatrix}\begin{bmatrix}\tilde{y}_{1t+1} \\ \tilde{y}_{2t+1}\end{bmatrix}+\begin{bmatrix}Q^{H}_{1} \\ Q^{H}_{2}\end{bmatrix}(Cv_{t+1}+D\eta_{t+1})
$$
其中，$\tilde{y}_{t}=Z^{H}y_{t}$，下半部分可以向前递归得到：
$$
\tilde{y}_{2t}=\sum\limits_{j=0}^{\infty}(T_{22}^{-1}S_{22})^{j}[T_{22}^{-1}Q_{2}^{H}(Cv_{t+1+j}+D\eta_{t+1+j})]=0
$$
得到了$\tilde{y}_{2t}$，就可以求出$\tilde{y}_{1t}$，再利用$\tilde{y}_{t}=Z^{H}y_{t}$就可以得到如下形式的解（LDE）：
$$
y_{t}=Fy_{t-1}+Gv_{t}
$$
其中，
$$
F=Z^{H}S_{11}^{-1}[T_{11}(T_{22}-\Phi T_{22})]Z
$$
$$
G=Z^{H}\begin{bmatrix}S_{11}^{-1}&-S_{11}^{-1}(S_{12}-\Phi S_{22}) \\ 0&I \end{bmatrix}\begin{bmatrix}Q_{1}^{H}-\Phi Q_{2}^{H} \\0\end{bmatrix}
$$
也就是说，经济模型最终转化为如上的线性差分方程。
以上线性差分方程用于模拟经济运行的优势在于线性差分方程可能出现的各种现象（如混沌）可以拟合经济波动和周期现象，因此用于研究中短期波动问题很有效。另外，由于得到的模型具有长期均衡的特征，因此同样可以用于长期经济发展的研究。
BMR包虽然提供了R语言运用Sims方法求解的函数gensys()，但是却没有支持该解法的后续步骤，因此仅作了粗略求解。
我们以M. D. Negro和F. Schorfheide(2004)的先验参数均值为例求解该LRES，得到LDE。代码如下：
library(BMR)
#SDSGE
# theta
mats = list()
mats$Gamma0 = matrix(0,5,5)
mats$Gamma1 = matrix(0,5,5)
mats$C = matrix(0,5,1)
mats$Psi = matrix(0,5,3)
mats$Pi = matrix(0,5,2)
lngam = 0.5
lnpi = 1
lnr =0.5
kap = 0.3
tau =0.2
psi1 = 1.5
psi2 = 0.25
rhog = 0.5
rhoz = 0.8
rhor = 0.3
sigg = 0.251
sigz = 0.63
sigr = 0.875
mats$Gamma0[1,]= -c(1, 1/tau,0,0,0)
mats$Gamma0[2,]= -c(exp(lngam-lnr),0,0,0,0)
mats$Gamma0[3,]= -c((1-rhor) * psi2,(1-rhor) * psi1,-1,0,0)
mats$Gamma0[4,]= -c(0, 0, 0,-1, 0)
mats$Gamma0[5,]= -c(0, 0, 0, 0,-1)
mats$Gamma1[1,]= c(-1,0,-1/tau,(1-rhog),rhoz/tau)
mats$Gamma1[2,]= c(kap,-1,0,-kap,0)
mats$Gamma1[3,]= c(0,0,rhor,0,0)
mats$Gamma1[4,]= c(0,0,0,rhog,0)
mats$Gamma1[5,]= c(0,0,0,0,rhoz)
mats$Pi[1,] = c(-1,1/tau)
mats$Pi[2,] = c(-exp(lngam-lnr),0)
mats$Pi[3,] = c((rhor-1) * psi2,(rhor-1) * psi1)
mats$Psi[3,1] = 1
mats$Psi[4,2] = 1
mats$Psi[5,3] = 1
dsgetest = SDSGE(mats)
print(dsgetest$G1)
print(dsgetest$impact)BMR包默认的DSGE求解方法为Uhlig方法，为了方便后续的计算，我们再运用Uhlig方法进行求解。
$$
0=A\zeta_{1,t}+B\zeta_{1,t-1},C\zeta_{2,t},Dz_t
$$
$$
0=\mathbb{E}_{t}(F\zeta_{1,t+1}+G\zeta_{1,t}+H\zeta_{1,t-1}+J\zeta_{2,t+1}+K\zeta_{2,t}+Lz_{t+1}+Mz_{t})
$$
$$
z_t=Nz_{t-1}+\varepsilon_t
$$
Uhlig方法分为两种，一种是考虑控制变量，和跳跃变量的方程，跳跃变量方程为理性预期方程，而控制变量为线性方程，但我们知道，$\mathbb{E}_{t}y_t=y_t$，因此，我们可以将模型化为：
$$
0=\mathbb{E}_{t}(F\zeta_{t+1}+G\zeta_{t}+H\zeta_{t-1}+Lz_{t+1}+Mz_{t})
$$
$$
z_t=Nz_{t-1}+\varepsilon_t
$$
其中，$\zeta_{t}=[\zeta_{1,t},\zeta_{2,t}]^{\top}$。然而，对于本模型，由于第三项包含了随机扰动，因此令$z_t=[p_{t},\tilde{g}_{t},\tilde{z}_{t}]^{\top}$，其中，$p_t=p_{t-1}+\epsilon_{R,t}$，则$
\zeta_t=[\tilde{x}_{t},\tilde{\pi}_{t},\tilde{R}_{t}]^{\top}$，$\varepsilon_t=[\epsilon_{R,t}\epsilon_{g,t},\epsilon_{z,t}]^{\top}
$。
可以得到如下待求解系数矩阵：
$$
A=B=C=D=E=H=J=K=L
$$
$$
F=\begin{bmatrix}1 &1/\tau & 0\\ \gamma/r^ \ast & 0 &0 \\ (1-\rho_R)\psi_2 & (1-\rho_R)\psi_1 & -1\end{bmatrix}
$$
$$
G=\begin{bmatrix} -1 & 0 & -1/\tau \\ \kappa & -1 & 0 \\ 0 & 0 & \rho_R \end{bmatrix}
$$
$$
M=\begin{bmatrix} 0 & (1-\rho_g) & \rho_z/\tau \\ 0 & -\kappa & 0\\ 1 & 0 & 0 \end{bmatrix}
$$
$$
N=\begin{bmatrix}0 &0 &0 \\ 0 & \rho_g & 0\\ 0 & 0 & \rho_z \end{bmatrix}
$$
进入BMR的Uhlig函数进行求解，首先运用M. D. Negro和F. Schorfheide（2004）的先验参数均值进行求解：
library(BMR)
#SDSGE
# theta
lngam = 0.5
lnpi = 1
lnr =0.5
kap = 0.3
tau =0.2
psi1 = 1.5
psi2 = 0.25
rhog = 0.5
rhoz = 0.8
rhor = 0.3
sigg = 0.251
sigz = 0.63
sigr = 0.875
mats = list()
mats$A = matrix(0,0,0)
mats$B = matrix(0,0,0)
mats$C = matrix(0,0,0)
mats$D = matrix(0,0,0)
mats$F = matrix(0,3,3)
mats$G = matrix(0,3,3)
mats$H = matrix(0,3,3)
mats$J = matrix(0,0,0)
mats$K = matrix(0,0,0)
mats$L = matrix(0,3,3)
mats$M = matrix(0,3,3)
mats$N = matrix(0,3,3)
mats$F[1,]= -c(1, 1/tau,0)
mats$F[2,]= -c(exp(lngam-lnr),0,0)
mats$F[3,]= -c((1-rhor) * psi2,(1-rhor) * psi1,-1)
mats$G[1,]= c(-1,0,-1/tau)
mats$G[2,]= c(kap,-1,0)
mats$G[3,]= c(0,0,rhor)
mats$M[1,] = c(0,(1-rhog),rhoz/tau)
mats$M[2,] = c(0,-kap,0)
mats$M[3,] = c(1,0,0)
#mats$N[1,] =
mats$N[2,] = c(0,rhog,0)
mats$N[3,] = c(0,0,rhoz)
dsgetest = uhlig(mats$A,mats$B,C=mats$C,D=mats$D,mats$F,
G=mats$G,H=mats$H,J=mats$J,K=mats$K,L=mats$L,
M=mats$M,N=mats$N)模型仿真
至此，我们可以先制定先验参数，由Uhlig方法得到DSGE的求解结果（LDE）。该方程组是一个包含外生冲击的结构式方程，因此可以运用MC方法进行仿真和分析。需要注意的是，这一步还没有引入实际变量的观测值，而所有的分析结果仅能用来说明模型本身的性质而非经济体的性质。
首先，对该结构方程可以分析各外生冲击对各内生经济变量的影响，我们以上一节求解的方程组为例，可以运用R语言BMR包中的IRF函数得到
library(grDevices)
IRF(dsgetest,
shock=1,irf.periods=10,varnames = c(
'output','inflation','interest','monetary','government','technology'),
save = TRUE)然后，根据BMR包里的DSGESim函数，我们可以得到DSGE的仿真结果，代码如下：
sim_data = DSGESim(dsgetest,shocks.cov = shocks,sim.periods = 200,burnin = 200)
par(mfrow=c(2,3))
plot(sim_data[,1],type = 'l',main='output_simu')
plot(sim_data[,2],type = 'l',main='inflat_simu')
plot(sim_data[,3],type = 'l',main='intere_simu')
plot(sim_data[,4],type = 'l',main='moneta_simu')
plot(sim_data[,5],type = 'l',main='govern_simu')
plot(sim_data[,6],type = 'l',main='techno_simu')模型估计
DSGE模型估计主要有两个主流方法。一是通过构建可观测变量和不可观测变量（内生变量）的状态空间模型，运用卡尔曼滤波得到估计参数用的对数条件似然函数，然后运用贝叶斯估计输入先验分布，得到参数的后验分布。另一种是直接将LDE简化为SVAR方法，通过简化线性差分方程为可观测变量，如运用HP滤波得到潜在变量当做稳态值，以偏离稳态的产出、价格、利率缺口建立VAR模型研究脉冲响应关系。
状态空间模型构建
为了进行实证分析，我们首先得到内生变量的度量方法：
$$
\Delta ln x_{t}=ln \gamma + \Delta \tilde{x}_{t} +\tilde{z}_{t}
$$
$$
\Delta ln P_{t}=ln \pi ^{\ast}+\tilde{\pi}_{t}
$$
$$
ln R^{a}_{t}=12[(ln r^{\ast}+ln \pi^{\ast})+\tilde{R}_{t}]
$$
其中$ln R^{a}_{t}$为年化名义利率，$\Delta ln P_{t}$为物价增长率，$\Delta ln x_{t}$为产出增长率。
令可观测变量$\mu_{t}=[\Delta ln x_{t},\Delta ln P_{t},ln R^{a}_{t}]$，该形式与求解得到的LDE模型联立，即可转化为状态空间模型。
下面，简单介绍DSGE的估计算法。
为了估计DSGE模型，首先我们可以看到最终待估形式为LDE。其中结构冲击$v_{t}$前有个系数矩阵，且所有内生变量均不可观测，不能简单由可观测变量替代。（这与传统的货币政策分析时运用潜在产出得到产出缺口，目标利率得到利率缺口的做法不同，理论上更加严谨）
其次，为了结合校准法得到的先验参数分布信息和样本的信息，一般采用Bayesian估计的技术来获得参数的估计。
因此，DSGE模型估计分为两步，第一步，建立状态空间模型，运用Kalman滤波得到Bayesian估计所需的似然函数；第二步，运用MCMC生产Markov链，Gibbs抽样得到先验参数输入，求得参数输出的后验分布。
本文的模型可以转化为如下的一般的状态空间模型：
$$
\mu_{t}=f(y_t,\varepsilon_{\mu,t}|\theta)
$$
$$
y_{t}=g(y_{t-1},\varepsilon_{y,t}|\theta)
$$
其中，$\mu_{t}$表示t时期的可观测变量样本，$y_t$表示不可观测变量即状态变量。这两个方程又称为观测方程和转移方程。根据定义，似然函数为给定参数值$\theta$下产生样本$\mu^{\top}$的概率，其中样本的观测序列为$\mu^{\top}=[\mu_1,\mu_2,…,\mu_T]$。
DSGE深参数$\theta$估计的似然函数可以表示为如下条件概率密度函数：
$$
p(\mu^{\top}|\theta)=p(\mu_{1}|\theta)\prod\limits_{t=2}^{T}p(\mu_{t}|\mu_{t-1},\theta)
$$
将$y_{t}$的条件概率密度$p(y_t|\mu_t,\theta)$用全概率公式展开可以得到：
$$
p(\mu^{\top}|\theta)=p(\mu_{1}|\theta)\prod\limits_{t=2}^{T}\int_{R}p(\mu_{t}|y_t,\theta)p(y_{t}|\mu_{t-1},\theta)dy_t
$$
其中，$p(\mu_1|\theta)=\int p(\mu_1|y_1,\theta)p(y_1)dy_1$。$p(y_1)$是内生不可观测变量的初始分布。调整初始分布，通过预测和更新两个过程，可以得到一个概率密度为$p(y_t|\mu_{t-1},\theta)$状态变量动态序列。
首先，在已知t-1期信息时，我们可以预测状态变量的分布如下：
$$
p(y_t|\mu_{t-1},\theta)=\int_{R}p(y_t|y_{t-1},\mu_{t-1},\theta)p(y_{t-1}|\mu_{t-1},\theta)dy_{t-1}
$$
这个估计结果可以得到一个新的条件概率密度函数，$p(y_t|\mu_t,\theta)$。这个密度函数是基于t期“新息”得到的，因此我们可以运用Bayes’ rule更新我们对$y_t$的估计
$$
p(y_t|\mu_t,\theta)=\frac{p(\mu_t|y_t,\theta)p(y_t|\mu_{t-1},\theta)}{\int p(\mu_t|y_t,\theta)p(y_t|\mu_{t-1},\theta)dy_t}
$$
因此，我们得到了滤波问题的预测概率密度函数。
当转移方程和观测方程都是线性的，且各自的扰动项均为高斯白噪声时，这个似然函数可以用Kalman滤波来计算。
Kalman滤波的过程较为复杂，大致就是一期一期预测，从而估计似然函数。
数据获取与处理
本文以国家统计局的GDP2006-2017季度数据、居民消费价格指数（上年同月）2006-2017月度数据和上海银行间同业拆放利率的2007年-2017年日度数据为基础，计算得到三个变量2007-2017年的月度数据。基础数据获取代码如下：
import tushare as ts
gdpq = ts.get_gdp_quarter()
gdpq.to_csv('gdpq.csv')
df_s = ts.shibor_data(2007)
for year in range(2008,2017):
df = ts.shibor_data(year).sort('date',ascending=True)
print(df.head(10))
df_s = df_s.append(df)
print(df_s.head(10))
df_s.to_csv('shibor.csv')
cpi = ts.get_cpi()
cpi.to_csv('cpi.csv')对GDP季度数据计算季度同比增长率，以季度数据代替1、4、7、10月数据，用线性插值法求得月度数据。对shibor数据按月平均，求得月度数据。
首先将GDP季度数据用Excel呈现，计算得到季度的对数差分值：
然后复制该序列到新的Excel，根据公式：
$$
Y_{2}=Y_{1}+1 \times (Y_{3}-Y_{1})/3
$$
$$
Y_{2}=Y_{1}+2 \times (Y_{3}-Y_{1})/3
$$
在Excel手动计算得到各月的差分结果。
如图所示，第一列表示1月、4月、7月和10月，第二列表示2月、5月、8月和11月，第三列3月、6月、9月和12月，通过以下程序读取名为“py.xlsx”文件，并输出名为“dlny.csv”的文件：
library(readxl)
# GDP季度转月度
dataset = read_excel('py.xlsx', col_names = FALSE)
allset = as.numeric(dataset[1,])
for(i in 2:nrow(dataset)){
allset = append(allset,as.numeric(dataset[i,]))
}
write.csv(allset,'dlny.csv', row.names=F)然后，用R语言实现shibor日数据向月度平均数据的转换。由于shibor数据包含“隔日”、“一周”、“两周”、“一月”、“一季度”、“半年”、“九个月”和“一年”共8种期限，而发布日相同，因此我们得到的日数据如下所示：
用以下思路编写R程序，输入为“shibor.xlsx”，输出为月度数据“tshibor.csv”
代码如下：
library(readxl)
# shibor数据整理
shibor = read_excel('shibor.xlsx')
# 初始数据准备
{
mon = ''
on = c()
w1 = c()
w2 = c()
m1 = c()
m3 = c()
m6 = c()
m9 = c()
y1 = c()
ton = c(shibor$'ON'[1])
tw1 = c(shibor$'1W'[1])
tw2 = c(shibor$'2W'[1])
tm1 = c(shibor$'1M'[1])
tm3 = c(shibor$'3M'[1])
tm6 = c(shibor$'6M'[1])
tm9 = c(shibor$'9M'[1])
ty1 = c(shibor$'1Y'[1])
}
for(i in 2:nrow(shibor)){
if(mon == months(shibor$date[i]) & i != nrow(shibor)){
#载入临时数据
ton = c(ton,shibor$'ON'[i])
tw1 = c(tw1,shibor$'1W'[i])
tw2 = c(tw2,shibor$'2W'[i])
tm1 = c(tm1,shibor$'1M'[i])
tm3 = c(tm3,shibor$'3M'[i])
tm6 = c(tm6,shibor$'6M'[i])
tm9 = c(tm9,shibor$'9M'[i])
ty1 = c(ty1,shibor$'1Y'[i])
}else{
# 月度数据导入
mon = months(shibor$date[i])
on = c(on,mean(ton))
w1 = c(w1,mean(tw1))
w2 = c(w2,mean(tw2))
m1 = c(m1,mean(tm1))
m3 = c(m3,mean(tm3))
m6 = c(m6,mean(tm6))
m9 = c(m9,mean(tm9))
y1 = c(y1,mean(ty1))
# 临时数据重置
ton = c(shibor$'ON'[i])
tw1 = c(shibor$'1W'[i])
tw2 = c(shibor$'2W'[i])
tm1 = c(shibor$'1M'[i])
tm3 = c(shibor$'3M'[i])
tm6 = c(shibor$'6M'[i])
tm9 = c(shibor$'9M'[i])
ty1 = c(shibor$'1Y'[i])
}
}
tshibor = data.frame(cbind(
on[-1],w1[-1],w2[-1],m1[-1],m3[-1],m6[-1],m9[-1],y1[-1]))
write.csv(tshibor,'tshi.csv', row.names=F)
shibor_ma = read_excel('~/Desktop/金融综合实验/货币政策/Rvar/shibor_ma.xlsx', col_types = c('date', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric'))
length(shibor_ma)
date = shibor_ma$date
# 选用一个利率指标，得到月度数据
da = shibor_ma$`1Y_5`
{
mon = ''
mo = c()
tmo = c(da[1])
}
for(i in 2:(length(da)-1)){
if(mon == months(date[i])){
#载入临时数据
tmo = c(tmo,da[i])
}else{
# 月度数据导入
mon = months(date[i])
mo = c(mo,mean(tmo))
# 临时数据重置
tmo = c(da[i])
}
}
tmo = c(tmo,da[i])
mo = c(mo,mean(tmo))
plot(mo,type='l')
library(tseries)
adf.test(log(mo))将数据整理为data，并保存为var_.xlsx文件，根据shibor利率的类别写出后缀。
模型估计
R语言BMR包提供了针对Uhlig求解方法得到的DSGE模型的估计函数EDSGE。
该方法即运用基于Gibb’s抽样的MCMC模拟样本，运用Kalman滤波计算似然函数，再根据贝叶斯估计得到后验参数分布。本文采用M. D. Negro和F. Schorfheide（2004）的先验参数分布和先验参数范围，运用上文得到的季度数据估计DSGE模型并对比BVAR估计结果。
首先，先验参数如下所示
先验参数范围给定如下：先验参数设置代码如下：
# 参数设置
parnames = c('lngam','lnpi','lnr','kap','tau','psi1','psi2','rhor','rhog','rhoz','sigr','sigg','sigz')
priorform = c('Normal','Normal','Gamma','Gamma','Gamma','Gamma','Gamma','Beta','Beta','Beta','IGamma','IGamma','IGamma')
priorpars = cbind(c(0.5,1,.5,.3,2,1.5,.125,.5,.8,.3,.251,.63,.875),
c(.25,.5,.25,.15,.5,.25,.1,.2,.1,.1,.139,.323,.42))
priorbounds = cbind(c(.101,.219,.132,.06,1.197,1.121,.001,.157,NA,NA,NA,NA,NA),
c(.922,1.863,.88,.513,2.788,1.91,.26,.812))由于R语言的BMR包需要给出重要函数partomat，该函数的输入参数即Deep参数，而输出为求解结果、冲击协方差矩阵、状态方程观测常数和观测误差协方差矩阵。该函数需要自己编写，代码如下：
# 该函数依赖于BMR、GRdevices等包
# partomat函数
partomat = function(lngam,lnpi,lnr,kap,tau,psi1,psi2,rhor,rhog,rhoz,sigr,sigg,sigz){
#SDSGE
# theta
mats = list()
mats$A = matrix(0,0,0)
mats$B = matrix(0,0,0)
mats$C = matrix(0,0,0)
mats$D = matrix(0,0,0)
mats$F = matrix(0,3,3)
mats$G = matrix(0,3,3)
mats$H = matrix(0,3,3)
mats$J = matrix(0,0,0)
mats$K = matrix(0,0,0)
mats$L = matrix(0,3,3)
mats$M = matrix(0,3,3)
mats$N = matrix(0,3,3)
mats$F[1,]= -c(1, 1/tau,0)
mats$F[2,]= -c(exp(lngam-lnr),0,0)
mats$F[3,]= -c((1-rhor) * psi2,(1-rhor) * psi1,-1)
mats$G[1,]= c(-1,0,-1/tau)
mats$G[2,]= c(kap,-1,0)
mats$G[3,]= c(0,0,rhor)
mats$M[1,] = c(0,(1-rhog),rhoz/tau)
mats$M[2,] = c(0,-kap,0)
mats$M[3,] = c(1,0,0)
#mats$N[1,] =
mats$N[2,] = c(0,rhog,0)
mats$N[3,] = c(0,0,rhoz)
dsgetest = uhlig(mats$A,mats$B,C=mats$C,D=mats$D,mats$F,
G=mats$G,H=mats$H,J=mats$J,K=mats$K,L=mats$L,
M=mats$M,N=mats$N)
shocks = matrix(0,3,3)
shocks[1,1] = sigr
shocks[2,2] = sigg
shocks[3,3] = sigz
ObsCons = matrix(0,3,1)
MeasErrs = matrix(0,3,3)
return=list(A=A,B=B,C=C,D=D,F=F,G=G,H=H,J=J,K=K,L=L,M=M,N=N,
shocks=shocks,ObsCons=ObsCons,MeasErrs=MeasErrs)EDSGE还有以下重要参数ObserveMat,其表示的是观测方程的系数。本文的观测矩阵为：
$$
H=\begin{bmatrix} 1& & & & & 1\\ & 1 & & & & \\ & & 12 & & & \\ & & & 0 & & \\ & & & & 0 & \\ & & & & & 0\end{bmatrix}
$$
表示的是深参数的初始值 可以直接设定为其均值：
# 参数设置
parnames = c('lngam','lnpi','lnr','kap','tau','psi1','psi2','rhor','rhog','rhoz','sigr','sigg','sigz')
priorform = c('Normal','Normal','Gamma','Gamma','Gamma','Gamma','Gamma','Beta','Beta','Beta','IGamma','IGamma','IGamma')
priorpars = cbind(c(0.5,1,.5,.3,2,1.5,.125,.5,.8,.3,.251,.63,.875),
c(.25,.5,.25,.15,.5,.25,.1,.2,.1,.1,.139,.323,.42))
parbounds = cbind(c(.101,.219,.132,.06,1.197,1.121,.001,.157,NA,NA,NA,NA,NA),
c(.922,1.863,.88,.513,2.788,1.91,.26,.812,NA,NA,NA,NA,NA))
# 该函数依赖于BMR、GRdevices等包
# partomat函数
partomat = function(parameters){
#SDSGE
lngam = parameters[1]
lnpi = parameters[2]
lnr = parameters[3]
kap = parameters[4]
tau = parameters[5]
psi1 = parameters[6]
psi2 = parameters[7]
rhor = parameters[8]
rhog = parameters[9]
rhoz = parameters[10]
sigr = parameters[11]
sigg = parameters[12]
sigz = parameters[13]
# theta
mats = list()
mats$A = matrix(0,0,0)
mats$B = matrix(0,0,0)
mats$C = matrix(0,0,0)
mats$D = matrix(0,0,0)
mats$F = matrix(0,3,3)
mats$G = matrix(0,3,3)
mats$H = matrix(0,3,3)
mats$J = matrix(0,0,0)
mats$K = matrix(0,0,0)
mats$L = matrix(0,3,3)
mats$M = matrix(0,3,3)
mats$N = matrix(0,3,3)
mats$F[1,]= -c(1, 1/tau,0)
mats$F[2,]= -c(exp(lngam-lnr),0,0)
mats$F[3,]= -c((1-rhor) * psi2,(1-rhor) * psi1,-1)
mats$G[1,]= c(-1,0,-1/tau)
mats$G[2,]= c(kap,-1,0)
mats$G[3,]= c(0,0,rhor)
mats$M[1,] = c(0,(1-rhog),rhoz/tau)
mats$M[2,] = c(0,-kap,0)
mats$M[3,] = c(1,0,0)
#mats$N[1,] =
mats$N[2,] = c(0,rhog,0)
mats$N[3,] = c(0,0,rhoz)
dsgetest = uhlig(mats$A,mats$B,C=mats$C,D=mats$D,mats$F,
G=mats$G,H=mats$H,J=mats$J,K=mats$K,L=mats$L,
M=mats$M,N=mats$N)
shocks = matrix(0,3,3)
shocks[1,1] = sigr
shocks[2,2] = sigg
shocks[3,3] = sigz
ObsCons = matrix(0,3,1)
MeasErrs = matrix(0,3,3)
ObsCons[1]=lngam
ObsCons[2]=lnpi
ObsCons[3]=12 * (lnr+lnpi)
return=list(A=mats$A,B=mats$B,C=mats$C,D=mats$D,F=mats$F,G=mats$G,H=mats$H,J=mats$J,K=mats$K,L=mats$L,M=mats$M,N=mats$N,
shocks=shocks,ObsCons=ObsCons,MeasErrs=MeasErrs)
}
# Observ
ObserveMat = matrix(0,3,6)
ObserveMat[1,1] = 1
ObserveMat[2,2] = 1
ObserveMat[3,3] = 12
ObserveMat[1,6] = 1
initialvals = c(0.5,1,.5,.3,2,1.5,.125,.5,.8,.3,.251,.63,.875)
NKest = EDSGE(bvar.data,chains = 1,cores=2,ObserveMat,initialvals,partomat,priorform,priorpars,parbounds,parnames,optimMethod = 'Nelder-Mead',optimLower=NULL,optimUpper=NULL,optimControl=list(maxit=10000),DSGEIRFs=TRUE,irf.periods=40,scalepar=0.27,keep=50000,burnin=75000)下面是运用BVAR估计的过程和结果
library(tibble)
library(readxl)
library(Rcpp)
install.packages('foreach')
install.packages('iterators')
install.packages('parallel')
install.packages('ggplot2')
library(BMR)
library(readxl)
library(tseries)
library(vars)
library(lmtest)
library(dse)
bvar.data = read_excel('var_1w.xlsx')
data(BMRVARData)
prior=c(0.9,0.95,0.95)
testbvarm = BVARM(bvar.data,NULL,p=4,constant=T,irf.periods=20,
keep=10000,burnin=5000,VType=1,
HP1=0.5,HP2=0.5,HP3=10)
plot(testbvarm,save=F)
IRF(testbvarm,save=F)
forecast(testbvarm,shocks=TRUE,backdata=10,save=FALSE)
由此可见，货币政策的利率机制还是相对平滑的，而货币政策的利率政策对通货膨胀冲击的反馈在滞后第8期左右才达到峰值。
而对产出冲击的反馈在滞后4期就达到峰值。这说明我国货币政策利率政策对产出冲击更为敏感。
参考文献
- Smets F，Wouters R. An estimated stochastic dynamic general equilibrium model for the Euro area[J]. Journal of the European Economic Association, 2003, 1: 1123-1175．
- Smets F，Wouters R. Shocks and Frictions in US Business Cycles: A Bayesian DSGE Approach[J]. American Economic Review, 2007, 97(3): 586-606.
- Del Negro M，Schorfheide F. Priors from general equilibrium models for VARs[J]. International Economic Review，2004，45(2) : 643-673．
- Christopher A. Sims Solving Linear Rational Expectations Models[J]. Computational Economics 20: 1–20, 2001.
发表回复