Stata: VAR - 模拟、估计和推断

发布时间:2020-02-05 阅读 344

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

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

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

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

编译: 许梦洁 (Frankfurt School of Finance and Management)
E-mail: m.xu@fs.de

Source: Ashish Rajbhandari[1]Vector autoregression—simulation, estimation, and inference in Stata


Stata 连享会时间序列专题:


目录


VAR 是分析多维时间序列动态变化的利器,该模型设定为由一组时间序列组成的序列是其自己滞后项的函数。

1. 模拟

首先使用如下设定模拟双变量 VAR(2) :

其中,y1,t 是在时间 t 的观测变量,μ 是一个 2×1 的截距项向量, A1 和 A2 均为 2×2 的参数矩阵,ϵt 是不与时间相关的扰动项。我假设 ϵt 服从分布 N(0,Σ) ,其中Σ 是一个 2×2 的协方差矩阵。

设定样本量为 1000 ,并在 Stata 中生成需要的变量:

 clear all

. set seed 2016

. local T = 1100

. set obs `T'
number of observations (_N) was 0, now 1,100

. gen time = _n

. tsset time
        time variable:  time, 1 to 1100
                delta:  1 unit

. generate y1 = .
(1,100 missing values generated)

. generate y2 = .
(1,100 missing values generated)

. generate eps1 = .
(1,100 missing values generated)

. generate eps2 = .
(1,100 missing values generated)

在第 1-6 行,我设定了随机种子,将样本量设为 1000,并且生成了一个时间变量 time 。接下来,我生成了变量 y1y2eps1eps2 来存放观测序列和扰动项。

2. 设定系数值

我为本文的 VAR(2) 模型选择了如下的参数值:

. mata:
------------------------------------------------- mata (type end to exit) -----
: mu = (0.1\0.4)

: A1 = (0.6,-0.3\0.4,0.2)

: A2 = (0.2,0.3\-0.1,0.1)

: Sigma = (1,0.5\0.5,1)

: end
-------------------------------------------------------------------------------

在 Mata 中,我分别创建了矩阵 μA1A2Σ 来放置参数值。在模拟得到样本之前,我检查了由这些参数值是否会得到一个平稳的 VAR(2) 模型。令:

其中 I2 是 2×2 的识别矩阵,0 是一个每个元素均为 0 的 2×2 矩阵。如果矩阵 F 的所有特征根均小于 1 ,则此 VAR(2) 过程即为一个平稳过程。下面放上计算特征根的代码:

. mata:
------------------------------------------------- mata (type end to exit) -----
: K = p = 2               // K = number of variables; p = number of lags

: F = J(K*p,K*p,0)

: F[1..2,1..2] = A1

: F[1..2,3..4] = A2

: F[3..4,1..2] = I(K)

: X = L = .

: eigensystem(F,X,L)

: L'
                              1
    +----------------------------+
  1 |                .858715598  |
  2 |  -.217760515 + .32727213i  |
  3 |  -.217760515 - .32727213i  |
  4 |                .376805431  |
    +----------------------------+

: end
------------------------------------------------------------------------------

我根据之前的设定创建了矩阵 F 并使用函数 eigensystem( ) 计算其特征根。 矩阵 X 中存储了特征向量,L 中保存特征根。所有 L 中的特征根均小于 1 ,因此该 VAR(2) 过程是平稳的。在检验过是否平稳之后,接下来生成 VAR(2) 模型的扰动项。

3. 由多维正态分布中生成扰动项

我从分布 N(0,Σ) 中生成两个随机正态变量,并将它们分别赋值给变量 eps1eps2

. mata:
--------------------------------------- mata (type end to exit) -----
: T = strtoreal(st_local("T"))

: u = rnormal(T,2,0,1)*cholesky(Sigma)

: epsmat = .

: st_view(epsmat,.,"eps1 eps2")

: epsmat[1..T,.] = u

: end
--------------------------------------------------------------------

我将样本规模(在 Stata 中定义为一个局部宏变量 T )赋值给一个 Mata 数值变量,这步将可以简化之后的工作。在 Mata 中,我使用了两个函数:st_local( )strtoreal( ) 来存储样本大小。第一个函数可以从 Stata 宏中获取文本值,第二个函数则可以将文本值转变为实数值。

第二行生成了一个 1100×2 的正态扰动项矩阵,其每个元素都服从分布 N(0,Σ) 。我使用函数 st_view( ) 将生成的正态扰动项存放到变量 eps1eps2 中。这个函数会以目前 Stata 中的数据集创建一个矩阵。首先,我创建了一个空矩阵 epsmat ,然后将变量 eps1eps2 填入到该矩阵中。最后,我使用矩阵 epsmat 的前 T 行位置来存储 u 中的数据。

4. 生成观测序列

沿用 Lütkepohl(2005) 的做法,我先生成了前两个观测值并使得它们的相关性与余下的样本一致。我假设前两项服从一个二维联合正态分布,其无条件期望为 θ=(IKA1A2)1μ,其协方差矩阵满足:

其中 vec() 是一个计算矩阵列数的算子,I16 是一个 16×16 的识别矩阵,

是一个 4×4 的矩阵。观测的前两项于是被生成为:

其中 Q 是一个 4×4 的矩阵,满足 QQ=Σyμ 是 4×1 的矩阵,其每个元素为服从标准正态分布的扰动项,Θ=(θθ) 则为 4×1 的期望矩阵。

下面的代码展示了如何生成前两个观测并给变量 y1y2 赋值的过程:

. mata:
----------------------------------- mata (type end to exit) -----
: Sigma_e = J(K*p,K*p,0)

: Sigma_e[1..K,1..K] = Sigma

: Sigma_y = luinv(I((K*p)^2)-F#F)*vec(Sigma_e)

: Sigma_y = rowshape(Sigma_y,K*p)'

: theta = luinv(I(K)-A1-A2)*mu

: Q = cholesky(Sigma_y)*rnormal(K*p,1,0,1)

: data = .

: st_view(data,.,"y1 y2")

: data[1..p,.] = ((Q[3..4],Q[1..2]):+mu)'

: end
-----------------------------------------------------------------

生成前两个观测后,我们可以使用如下的代码生成余下的序列:

. forvalues i=3/`T' {
.         qui {
. replace y1 = 0.1 + 0.6*l.y1 - 0.3*l.y2 + 0.2*l2.y1 + 0.3*l2.y2 + eps1 in `i'
. replace y2 = 0.4 + 0.4*l.y1 + 0.2*l.y2 - 0.1*l2.y1 + 0.1*l2.y2 + eps2 in `i'
.         }
. }
. drop in 1/100
(100 observations deleted)

我在 replace 命令后加入了 quietly 选项来筛选掉不必要显示的输出结果。最后,我删掉了前 100 个观测来避免模拟效果受到初值的影响。

5. 估计

我使用 var 命令来拟合 VAR(2) 模型;

. var y1 y2

Vector autoregression

Sample:  103 - 1100                             Number of obs     =        998
Log likelihood =  -2693.949                     AIC               =   5.418735
FPE            =   .7733536                     HQIC              =    5.43742
Det(Sigma_ml)  =   .7580097                     SBIC              =   5.467891

Equation           Parms      RMSE     R-sq      chi2     P>chi2
----------------------------------------------------------------
y1                    5     1.14546   0.5261   1108.039   0.0000
y2                    5     .865602   0.4794   919.1433   0.0000
----------------------------------------------------------------

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y1           |
          y1 |
         L1. |   .5510793   .0324494    16.98   0.000     .4874797     .614679
         L2. |   .2749983   .0367192     7.49   0.000       .20303    .3469667
             |
          y2 |
         L1. |  -.3080881    .046611    -6.61   0.000    -.3994439   -.2167323
         L2. |   .2551285   .0425803     5.99   0.000     .1716727    .3385844
             |
       _cons |   .1285357   .0496933     2.59   0.010     .0311387    .2259327
-------------+----------------------------------------------------------------
y2           |
          y1 |
         L1. |   .3890191   .0245214    15.86   0.000      .340958    .4370801
         L2. |  -.0190324    .027748    -0.69   0.493    -.0734175    .0353527
             |
          y2 |
         L1. |   .1944531    .035223     5.52   0.000     .1254172     .263489
         L2. |   .0459445   .0321771     1.43   0.153    -.0171215    .1090106
             |
       _cons |   .4603854   .0375523    12.26   0.000     .3867843    .5339865
------------------------------------------------------------------------------

var 命令在估计模型系数时默认为两阶滞后,参数估计是显著的并且与用来生成两个序列的参数真实值接近。

6. 推断:脉冲响应函数

脉冲响应函数 (IRF) 可以用来分析 VAR 模型中的内生变量如何对扰动项的冲击进行反应。比如,在由通胀和利率组成的双变量 VAR 模型中,脉冲响应函数可以追踪来自通胀方程的外生冲击如何对利率造成影响。

考虑最开始提到的双变量模型,假如我想估计一单位 ϵt 的变动对系统中的内生变量会造成怎样的冲击。我可以将 VAR(2) 过程通过数学变换变形为 MA() 过程:

其中 Φp 是 MA() 过程系数矩阵的的第 p 阶滞后。根据Lütkepohl(2005),MA 系数矩阵与 AR 过程的系数矩阵有如下联系:

也就是说,第 i 个变量在未来 h 期对第 j 个方程在时点 t 的单位冲击作出的反应为:

该式意味着第 i 个变量在未来 h 期对第 j 个方程在时点 t 的单位冲击作出的反应为 MA() 过程系数矩阵中第 i 行 j 列的元素值。对于 VAR(2) 过程,使用估计 AR 参数得到的前若干个反应值为:

对于 t>2yt 的响应使用类似的递归即可得到,第一个方程的冲击对第一个变量的影响为向量(1,0.5510,0.4588,)。我使用 irf create 命令来得到脉冲影响的结果:

. irf create firstirf, set(myirf)
(file myirf.irf created)
(file myirf.irf now active)
(file myirf.irf updated)

这个命令估计了 firstirf 模型中的脉冲响应函数和其他的统计量,并存储在文件 myirf.irf 中。选项 set( ) 将文件名为 myirf.irf 的文件激活。我可以在表格中列示来自同一方程扰动项的外生冲击对 y1 的影响:

. irf table irf, impulse(y1) response(y1) noci

Results from firstirf

+--------------------+
|        |    (1)    |
|  step  |   irf     |
|--------+-----------|
|0       | 1         |
|1       | .551079   |
|2       | .458835   |
|3       | .42016    |
|4       | .353356   |
|5       | .305343   |
|6       | .263868   |
|7       | .227355   |
|8       | .196142   |
+--------------------+
(1) irfname = firstirf, impulse = y1, and response = y1

默认的期数为 8 期,我加入 noci 选项选择不显示置信区间。注意 Stata 算出的前几期响应结果与我之前手动算的结果非常接近。带 95% 置信带的脉冲响应函数图如下:

. irf graph irf, impulse(y1) response(y1)
image
image

由上图可知,来自同一方程的一单位冲击将使 y1 立即增加一个单位,然后在中长期随着时间逐渐衰减。

7. 正交化的脉冲响应函数

在前面的部分中,我展示了其他条件不变时,来自同一方程的一单位冲击对 y1 造成的影响。然而,如下表所示,方差-协方差矩阵 Σ^ 显示了两个等式之间较强的正相关关系。

. matrix list e(Sigma)

symmetric e(Sigma)[2,2]
           y1         y2
y1  1.3055041
y2   .4639629  .74551376

两个等式的估计协方差是正的,这意味着我不能假设一个扰动项变动而另一个扰动项保持不变。来自 y2 方程的冲击会对 y1 造成同期影响,反之亦然。

正交脉冲响应函数 (OIRF) 通过将估计方差-协方差矩阵 Σ^ 分解为一个下三角矩阵来解决此问题。这种类型的分解可以将来自同一方程的冲击对 y1 的同期影响分离出来。然而,来自第一个方程的冲击依然会同期影响 y2 。比如,如果 y1 是通胀而 y1 是利率,这个分解意味着通胀的冲击会同时影响通胀和利率,但是来自利率的冲击则只会影响利率。

为了估计 OIRFs,令 P 表示 Σ 的 Cholesky 分解,即满足 PP=Σ 。令 ut 表示一个 2×1 的向量使得 Put=ϵt ,也即是 ut=P1ϵt 。 ut 中的误差与构造不相关,因为 E(utut)=P1E(ϵtϵt)P1=I2。这允许我们将 OIRFs 解读为来自 ut 的一个标准差的外生冲击造成的影响。

重写以 ut 向量表示的 MA() 过程:

OIRFs 则为 MA 过程的系数矩阵与下三角矩阵 P 的积:

通过以下代码得到 P^ 的估计:

. matrix Sigma_hat = e(Sigma)

. matrix P_hat = cholesky(Sigma_hat)

. matrix list P_hat

P_hat[2,2]
           y1         y2
y1  1.1425866          0
y2  .40606367  .76198823

使用这个矩阵,我可以计算前若干期的响应:

我将所有的 OIRFs 列示在一个表格中并画出了 y1 的脉冲响应图:

. irf table oirf, noci

                 Results from firstirf

+--------------------------------------------------------+
|        |    (1)    |    (2)    |    (3)    |    (4)    |
|  step  |   oirf    |   oirf    |   oirf    |   oirf    |
|--------+-----------+-----------+-----------+-----------|
|0       | 1.14259   | .406064   | 0         | .761988   |
|1       | .504552   | .523448   | -.23476   | .148171   |
|2       | .534588   | .294977   | .019384   | -.027504  |
|3       | .476019   | .279771   | -.0076    | .013468   |
|4       | .398398   | .242961   | -.010024  | -.00197   |
|5       | .346978   | .206023   | -.003571  | -.003519  |
|6       | .299284   | .178623   | -.004143  | -.001973  |
|7       | .257878   | .154023   | -.003555  | -.002089  |
|8       | .222533   | .13278    | -.002958  | -.001801  |
+--------------------------------------------------------+
(1) irfname = firstirf, impulse = y1, and response = y1
(2) irfname = firstirf, impulse = y1, and response = y2
(3) irfname = firstirf, impulse = y2, and response = y1
(4) irfname = firstirf, impulse = y2, and response = y2

irf table oirf 要求输出 OIRFs 的结果。注意前三行的估计结果与我们之前手算的结果相同。

. irf graph oirf, impulse(y1) response(y1)
image
image

上图即为 y1 对来自同一个方程的一个单位的冲击的脉冲响应图。

8. 结论

在这篇推文中,我展示了怎么模拟一个平稳的 VAR(2) 模型,使用 var 命令估计了这个模型的参数,展示了如何估计 IRFs 和 OIRFs ,其中 OIRFs 由对协方差矩阵的下三角分解得到。

9. 参考文献

  • Lütkepohl, H. 2005. New Introduction to Multiple Time Series Analysis. New York: Springer.

相关课程

连享会-直播课 上线了!
http://lianxh.duanshu.com

免费公开课:


课程一览

支持回看,所有课程可以随时购买观看。

专题 嘉宾 直播/回看视频
最新专题 DSGE, 因果推断, 空间计量等
Stata数据清洗 游万海 直播, 2 小时,已上线
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]

Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。


关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。直播间 有很多视频课程,可以随时观看。
  • 连享会-主页知乎专栏,300+ 推文,实证分析不再抓狂。
  • 公众号推文分类: 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类,主流方法介绍一目了然:DID, RDD, IV, GMM, FE, Probit 等。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会主页  lianxh.cn
连享会主页 lianxh.cn

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会学习群-常见问题解答汇总:
https://gitee.com/arlionn/WD

参考资料

[1]

Posts by Ashish Rajbhandari, Senior Econometrician: https://blog.stata.com/author/arajbhandari/