StataSX2018
/
nwtregress_金钊_网络分析

# Stata：时变参数的网络分析-nwxtregress

<!-- TOC -->

> **作者：**
> 金钊（中山大学岭南学院）\
> **邮箱：**
> jinzh27@mail.sysu.edu.cn

> **备注：**
> 本文主要翻译自下文，特此致谢！\
> **Source：**
> Morad Zekhnini, 2021. "Network regressions in Stata," 2021 Stata Conference 17, Stata Users Group. [[pdf\]](http://fmwww.bc.edu/repec/scon2021/US21_Zekhnini.pdf)

- [Stata：时变参数的网络分析-nwxtregress](#stata时变参数的网络分析-nwxtregress)
  - [1. 问题背景](#1-问题背景)
  - [2. 理论背景](#2-理论背景)
  - [3. 命令介绍](#3-命令介绍)
    - [3.1 命令安装](#31-命令安装)
    - [3.2 语法格式](#32-语法格式)
    - [3.3 空间权重矩阵的形式](#33-空间权重矩阵的形式)
  - [5. 参考文献](#5-参考文献)

<!-- /TOC -->

## 1. 问题背景

&emsp;&emsp;随着网络理论的发展，关于微观主体行为的研究开始越来越多的强调社会网络的重要性。具体地，过去相关研究在讨论微观个体行为时，均假设个体间是独立不相关的。然而，现代社会微观个体之间的关联日益密切，相互影响。为此，经济学家们引入“网络理论”（Network Theory）试图对该类社会现象进行解释。此时，在实证研究中就将面临一个难题：如何在回归中引入社交网络？显然，传统的基于个体独立的实证方法将不再适用。空间计量方法为经济学家们提供了思路，他们运用空间计量模型来识别社交网络中的同群效应（Peer ffects）（Bramoull´e et al., 2009; Grieser et al., 2022）。本推文将详细介绍Zekhnini等（2021）提出的nwxtregress命令，该命令主要适用于非平衡面板，且网络结构或空间权重矩阵具有时变特征的数据。

## 2. 理论背景

&emsp;&emsp;我们常见的面板的空间自回归模型（Spatial Autocorrelation Model）为:

$$
Y_{it}=\rho\sum_{j\ne i}w_{ij}y_{jt}+\beta X_{it}+\varepsilon_{it}    \qquad(1)
$$

&emsp;&emsp;$w_{ij}$为空间权重矩阵$W$的$(i,j)$元素，\rho为空间自回归系数。传统的空间自回归模型的空间权重矩阵是不随时间变化的。然而，在现实中，微观个体间的相互关联是随着时间而变化的，即$w_{ij，t}。为了求解时变空间权重矩阵模型需要对空间序列进行变换，为此，我们把上式写成矩阵的形式如下：

$$
\mathsf{y}=\rho W \mathsf{y}+\beta X+\varepsilon   \qquad(2)
$$

&emsp;&emsp;为了去掉扰动项的空间自相关，可以对（2）式做准差分。此时，数据生成过程的简约模式如下：

$$
\mathsf{y}=（I-\rho W）^{-1} (\beta X+\varepsilon )  \qquad(3)
$$

需要注意的是，$\mathsf{y}$现在都在左边，但是模型中的参数是非线性的。给定$\rho $和$W$的数学约束为：

$$
（I-\rho W）^{-1} =I+\rho W+\rho^{2} W^{2} +...    \qquad(4)
$$

&emsp;&emsp;上式暗含了该模型具有几何衰减传播（Geometrically-Decaying Propagation）的特征。我们可以把该模型解读为以下各项的几何加总：

> $I$ ：表示自身的影响
> $W$ : 表示瞬时同伴效应（Immediate Peers' Effect）
> $W^{2}$ : 表示同伴的同伴效应（Peers of Peers Effect）
> ...：其他影响

此时，模型的偏导数为：

$$
\frac{\partial \mathsf{y}_i}{\partial\mathsf{x}_j}=（I-\rho W）^{-1}\beta,\forall i,j   \qquad(5)
$$

LeSage and Pace (2009)提出偏导系数的解释可以划分为直接效应和间接效应：

> 直接效应：$\frac{1}{N}\sum_{i}\frac{\partial \mathsf{y}_i}{\partial \mathsf{x}_i}$
> 间接效应：$\frac{1}{N}\sum_{i}\sum_{j \ne i}\frac{\partial \mathsf{y}_i}{\partial \mathsf{x}_i}$

&emsp;&emsp;我们可以在SAR模型中引入空间杜宾模型（SDM）：

$$
\mathsf{y}=\rho W \mathsf{y}+\beta X+WX\theta+\varepsilon   \qquad(6)
$$

$WX$ 表示来自同伴（Peers）的影响，一般是假设外生的，不影响估计。为了更直观的理解上述模型，重点放在截面形式下的似然函数为：

$$
\mathsf{f}(Y,X;\rho,\beta,\sigma^2)=|I_N-\rho W|(2\pi \sigma^2)^{-N/2}exp(-\frac{e'e}{2\sigma^2}) \\
   e=(I-\rho W)Y-X\beta\qquad(7)
$$

若$\rho$已知，那么$\beta$可以通过最大似然估计（MLE）得到。此时，最优化问题仅是关于$\rho$，运用MCMC估计上述的似然函数。

## 3. 命令介绍

### 3.1 命令安装

&emsp;&emsp;目前，空间面板命令 `spxtregress`和 `xsmle`只能适用于平衡面板和非时变空间权重矩阵。Morad Zekhnini等人开发的 `nwxtregress`命令可以运用于非平衡面板和空间权重矩阵时变的情形，更为详细的可见[nwxtregress: Network Regressions in Stata](https://janditzen.github.io/nwxtregress/)。 该非官方命令具体的安装方法如下：

```
net install nwxtregress , from(https://janditzen.github.io/nwxtregress/)
```
```
net from https://janditzen.github.io/nwxtregress/
```

### 3.2 语法格式

基本语法格式

```
Spatial Autocorrelation Model (SAR)
   nwxtregress depvar  indepvars [if],
        ivarlag(W1[, sparse timesparse mata id(string)])
        [mcmcoptions nosparse]

Spatial Durbin Model (SDM)
   nwxtregress depvar  indepvars [if],
        ivarlag(W1[, sparse timesparse mata id(string)])
        dvarlag(Ws:varlist[, sparse timesparse mata id(string)]
        [mcmcoptions nosparse]
```

其中，

+ `nwxtregress`：仅适用于SAR和SDM模型
+ `depvar`：被解释变量
+ `indepvars`：解释变量
+ `dvarlag`：空间滞后的被解释变量，可以支持多个和重复的空间权重矩阵
+ `ivarlag`：空间滞后的解释变量
+ `W1`和 `Ws`：空间权重矩阵，默认为方形矩阵格式

`ivarlag() `和 `dvarlag()`的 `Options`

+ `frame(name)` ：指定空间权重矩阵的形式，默认为方形矩阵格式，可以设定为稀疏形式或时变形式
+ `mata`：指定权重矩阵为mata矩阵
+ `sparse`：指定空间权重矩阵为稀疏矩阵
+ `timesparse`：指定空间权重矩阵为稀疏矩阵且时变
+ `id(string)`：指定时间变量、具有关联关系的个体
+ `normalize(string)`：指定空间权重矩阵的标准化形式，详细设定形式可见[spmat creat](https://www.stata.com/manuals/spspmatrixcreate.pdf)

一般 `Options`

+ `nosparse`：不将空间权重矩阵内部转换为稀疏矩阵
+ `asarray(name)`：更改估计结果和信息的名称

`MCMC`的 `Options`

+ `draws()`：griddy gibs的取样次数，默认2000
+ `gridlength()`：网格长度，默认1000
+ ` nomit()`：被删除的取样数，默认500
+ `barrypace(numlist)`：BarryPace Trick的设定，默认50和100
+ `usebp`：使用Barryface Trick而不是LUD处理矩阵的逆
+ `seed(#)`：设定抽样“种子”

计算总效应、直接效应和间接效应：
 ```
 estat impact [varlist] [, options]
  ```
 + `varlist`：指定需要计算的解释变量，默认计算所有的解释变量
  + `options`：包含 `seeds()`和`array(name)`

  估计后的预测：
   ```
   predict [type] varname [, options]
   ```
   + `xb`：计算线性预测
   + `res`：计算残差

  ### 3.3 空间权重矩阵的形式
`nwxtregress`的空间权重矩阵的形式主要有三种：
+ 非时变的$n\times n$方阵的形式，例如下列$5 \times 5$的矩阵：
  ```
        0     0.1   0.2   0
        0     0     0.1   0.2
        0.3   0.1   0     0
        0.2   0     0.2   0
    ```
+ 非时变的稀疏矩阵的形式，例如下列$v \times 3$的矩阵，$v$表示矩阵中的非0元素，第1列表示目的地，第2列表示发源地，第3列表示流入量
```
        Destination   Origin      Flow
           1            2         0.1
           1            3         0.2
           2            3         0.1
           2            4         0.2
           3            1         0.3
           3            2         0.1
           4            1         0.2
           4            3         0.2
```

+ 时变的稀疏矩阵的形式，如下所示：第1列表示时间，其他的与上述稀疏矩阵一致。我们可以看到，时变的稀疏矩阵具有两段时间。
```

         Time   Destination    Origin       Flow
           1       1              2         0.1
           1       1              3         0.2
           1       2              3         0.1
           1       2              4         0.2
           1       3              1         0.3
           1       3              2         0.1
           1       4              1         0.2
           1       4              3         0.2
                 (next time period)
           2       1              2         0.1
           2       1              3         0.4
           2       2              3         0.1
           2       2              4         0.4
           2       3              1         0.9
           2       3              2         0.1
           2       4              1         0.4
           2       4              3         0.4
```
## 4. Stata 实例
&emsp;&emsp;我们调用一份Zekhnini (2021)提供的示例数据集`IO.dta`和`VA.dta`.`IO.dta`数据集包含了1997-2019年美国每个各行业使用（USE）和制造（MAKE）的产品的数据。为了建立行业之间的联系，他们将这些产品看作是其在行业之间的流动（Flow），并以此构建空间权重矩阵。`VA.dta`为美国企业层面的数据。运用两个数据集，我们分析补偿金和净盈余对资本消耗的影响。
 ```
 * 导入数据 use https://janditzen.github.io/nwxtregress/examples/IO.dta

 /* 变量说明
 buyer_industry ：购买方
 seller_industry ：出售方
 ID1 ：购买方的id
 ID2 ：出售方的id
 sam ：交易量
  */
  ```

&emsp;&emsp; 首先，我们提取1998年的截面数据，构建一个方阵形式的空间权重矩阵`WSpmat`

```
   keep if Year == 1998 // 保留1998年的数据
   replace sam = 0 if sam < 0 // 交易量为负的变量设置为0
   replace sam = 0 if ID1==ID2 // 保证对角线为0
   keep ID1 ID2 sam // 保存数据
   reshape wide sam, i(ID1) j(ID2) // 转换数据排列形式，纵变横 spset ID1 // 指定为空间数据
   spmatrix fromdata WSpmat = sam* , replace // 根据数据集中的变量生成空间权重矩阵
```
&emsp;&emsp;接着，我们提取`VA.dta`数据，估计非时变空间权重矩阵下的SAR模型。在进行估计之前，需要先安装`moremata`命令，以保证该命令的顺畅运行，安装方式如下：
```dotnetcli
ssc install moremata
```
安装完毕，我们就可以进行估计，先看一下样本情况，可以发现该数据是一份非平衡面板。
```dotnetcli
. xtset ID Year        //设定面板数据
Panel variable: ID (unbalanced)
 Time variable: Year, 1998 to 2019, but with a gap
         Delta: 1 unit

. spset ID     //设定空间数据
  (data are not strongly balanced; run spbalance)
      Sp dataset: va.dta
Linked shapefile: <none>
            Data: Panel
 Spatial-unit ID: _ID (equal to ID)
         Time ID: Year (see xtset)
     Coordinates: <none>
```
SAR模型估计如下：

```dotnetcli
 * 导入数据 use https://janditzen.github.io/nwxtregress/examples/VA.dta

. nwxtregress cap_cons compensation net_surplus , dvarlag(WSpmat) seed(1234)
This is an alpha version! Results may change.

Initialise Grid (200)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100
Griddy Gibbs (2000)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100

                                                 Number of obs     =       1363
Panel Variable (i): ID                           Number of groups  =         63
Time Variable (t): Year                          Obs. of group:              22
                                                               min =          5
                                                               avg =         22
                                                               max =         22
                                                 R-squared         =       0.77
                                                 Adj. R-squared    =       0.77
-------------------------------------------------------------------------------
       cap_cons|     Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+---------------------------------------------------------------
   compensation| -.7198294   .0122852  -58.59    0.000     -.7673457  -.6796333
    net_surplus| -.7596348   .0142442  -53.33    0.000     -.8048452  -.7084874
          _cons|  .7214424   .0117858   61.21    0.000        .68282    .766456
---------------+---------------------------------------------------------------
WSpmat         |
       cap_cons|  .1213355   .0251213    4.83    0.000          .038       .202
---------------+---------------------------------------------------------------
       /sigma_u|   .002958    .000112                       .0026109   .0033606
-------------------------------------------------------------------------------
```
上面的结果显示空间自回归系数在$\rho$在1%水平上显著为正，我们可以计算直接效应，间接效应和总效应，需要用到下面的命令。
```
. estat impact
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100
Average Impacts                                  Number of obs =        1363
-------------------------------------------------------------------------------
       cap_cons|     dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+---------------------------------------------------------------
direct         |
   compensation| -1.913991   .4786752   -4.00    0.000     -4.064068  -.8469636
    net_surplus| -2.019859   .5054806   -4.00    0.000     -4.245732  -.8834585
---------------+---------------------------------------------------------------
indirect       |
   compensation|  1.094091   .4545737    2.41    0.016      .0853722   3.122538
    net_surplus|  1.154618   .4798387    2.41    0.016      .0890508   3.262116
---------------+---------------------------------------------------------------
total          |
   compensation| -.8198999   .0273082  -30.02    0.000     -.9415295  -.7382084
    net_surplus| -.8652407   .0296407  -29.19    0.000      -.983616  -.7706425
-------------------------------------------------------------------------------
```
&emsp;&emsp; 我们继续检验时变权重矩阵的情况。此时，我们需要构建新的空间权重矩阵，并将其命名为`IO`，具体构建命令如下：

```dotnetcli
. frame create IO \\创建时变的空间权重矩阵

. frame IO: use https://janditzen.github.io/nwxtregress/examples/IO.dta  \\运用IO数据集创建空间权重矩阵
```
运用`IO`矩阵估计时变SAR模型
```dotnetcli
. nwxtregress cap_cons compensation net_surplus ,dvarlag(sam, frame(IO) id(Year ID1 ID2)  timesparse) seed(1234)

This is an alpha version! Results may change.

Initialise Grid (200)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100
Griddy Gibbs (2000)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100

                                                 Number of obs     =       1363
Panel Variable (i): ID                           Number of groups  =         63
Time Variable (t): Year                          Obs. of group:              22
                                                               min =          5
                                                               avg =         22
                                                               max =         22
                                                 R-squared         =       0.77
                                                 Adj. R-squared    =       0.77
-------------------------------------------------------------------------------
       cap_cons|     Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+---------------------------------------------------------------
   compensation| -.7270558   .0121453  -59.86    0.000     -.7736677  -.6874203
    net_surplus| -.7661228   .0140726  -54.44    0.000     -.8100804  -.7154614
          _cons|  .7281612   .0114999   63.32    0.000      .6903187   .7720147
---------------+---------------------------------------------------------------
sam            |
       cap_cons|  .1129735   .0247994    4.56    0.000           .03       .193
---------------+---------------------------------------------------------------
       /sigma_u|  .0029633   .0001124                        .002616   .0033666
-------------------------------------------------------------------------------
```
在`dvarlag()`中，第一项是空间权重矩阵中的关联变量，第二项是空间权重矩阵，第三项是要指定空间权重矩阵中关联对象`id()`，由于是时变的矩阵，因此包含三项：时间，来源地、目的地（即生产方和消费方）。最后，矩阵是时变的需要指定` timesparse`。

&emsp;&emsp; 我们还可以把该矩阵导入mata中，构建一个mata矩阵**Wt**具体构建命令如下：
```dotnetcli
. frame IO: putmata Wt = (Year ID1 ID2 sam), replace
(1 matrix posted)

. nwxtregress cap_cons compensation net_surplus , dvarlag(Wt, mata timesparse) seed(1234)   ///估计的结果与上表一致，此处省略
```
&emsp;&emsp; 最后，我们可以在以上模型中加入`ivarlag()`，就可以估计SDM模型，具体命令如下：
```dotnetcli
. nwxtregress cap_cons compensation net_surplus , dvarlag(Wt,mata timesparse) ivarlag(Wt: compensation,mata timesparse )  seed(1234)

This is an alpha version! Results may change.

Initialise Grid (200)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100
Griddy Gibbs (2000)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100

                                                 Number of obs     =        1363
Panel Variable (i): ID                           Number of groups  =          63
Time Variable (t): Year                          Obs. of group:              22
                                                               min =           5
                                                               avg =          22
                                                               max =          22
                                                 R-squared         =        0.77
                                                 Adj. R-squared    =        0.77
------------------------------------------------------------------------ -------
       cap_cons|     Coef.   Std. Err.      z    P>|z|     [95% Conf. In terval]
---------------+---------------------------------------------------------------
   compensation| -.6952461   .0125736  -55.29    0.000     -.7427079   - .649156
    net_surplus| -.7428038   .0141965  -52.32    0.000     -.7879067  -. 6917178
          _cons|  .7684855   .0142525   53.92    0.000      .7242406   . 8139682
---------------+---------------------------------------------------------------
Wt             |
       cap_cons|   .039009   .0289965    1.35    0.179         -.048        .129
   compensation| -.0898353   .0174214   -5.16    0.000     -.1550298  -. 0253126
---------------+---------------------------------------------------------------
       /sigma_u|  .0029072   .0001131                       .0025408     .003339
------------------------------------------------------------------------------
```
&emsp;&emsp; 以上的回归中只包含了一个空间权重矩阵，我们可以构建两个时变的空间权重矩阵加入回归方程中，具体如下：
```dotnetcli
. mata: Wt2 = Wt[selectindex(Wt[.,4]:>2601.996),.] //生成Wt2矩阵

. nwxtregress cap_cons compensation net_surplus , dvarlag(Wt, mata timesparse) ivarlag(Wt: net_surplus, mata timesparse) ivarlag(Wt2: compensation, mata timesparse) seed(1234)

This is an alpha version! Results may change.

Initialise Grid (200)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100
Griddy Gibbs (2000)
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100

                                                 Number of obs     =    1363
Panel Variable (i): ID                           Number of groups  =      63
Time Variable (t): Year                          Obs. of group:           22
                                                               min =       5
                                                               avg        22
                                                               max =      22
                                                 R-squared         =    0.77
                                                 Adj. R-squared    =    0.77
-------------------------------------------------------------------------------
       cap_cons|     Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+---------------------------------------------------------------
   compensation| -.7009335   .0137442  -51.00    0.000     -.7456487  -.6560701
    net_surplus| -.7475387   .0148767  -50.25    0.000     -.7958839  -.6960453
          _cons|  .6908152   .0157555   43.85    0.000      .6272978    .742689
---------------+---------------------------------------------------------------
Wt             |
       cap_cons|   .132362   .0266373    4.97    0.000          .046       .233
    net_surplus|  .0730267   .0193825    3.77    0.000      .0029279   .1501232
---------------+---------------------------------------------------------------
Wt2            |
   compensation|  .0006659      .0071    0.09    0.925     -.0260179   . 0229685
---------------+---------------------------------------------------------------
       /sigma_u|  .0029302   .0001123                       .0025453   . 0033302
-------------------------------------------------------------------------------
```
上式中，两个解释变量的权重不一致。我们可以进一步分析直接效应、间接效应和总效应
```
. estat impact
----+--- 10 --+--- 20 --+--- 30 --+--- 40 --+--- 50   %
..................................................   50
..................................................  100
Average Impacts                                  Number of obs =       1363
-------------------------------------------------------------------------------
       cap_cons|     dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+---------------------------------------------------------------
direct         |
   compensation| -.7009335   .0137442  -51.00    0.000     -.7456487  -.6560701
    net_surplus| -.7475387   .0148767  -50.25    0.000     -.7958839  -.6960453
---------------+---------------------------------------------------------------
indirect       |
   compensation|  .6908152   .0157555   43.85    0.000      .6272978    .742689
    net_surplus|  .0006659      .0071    0.09    0.925     -.0260179   .0229685
---------------+---------------------------------------------------------------
total          |
   compensation| -.0116701   .0089471   -1.30    0.192     -.0394202   .0193187
    net_surplus| -.8616318   .0325265  -26.49    0.000     -.9834741  -.7587489
-------------------------------------------------------------------------------
```
该模型还可以用于预测，具体如下：
```dotnetcli
. predict xb  //计算线性预测
. predict residuals, residual  // 计算残差
```
预测的值会生成新的变量储存在`dta`中。

## 5. 参考文献
Bramoull´e, Y., H. Djebbari, and B. Fortin. 2009. Identification of peer effects through social networks. Journal of econometrics 150(1): 41–55.[[Link\]](https://www.sciencedirect.com/science/article/pii/S0304407609000335)\
Grieser, W., C. Hadlock, J. LeSage, and M. Zekhnini. 2022. Network Effects in Corporate Financial Policies. Journal of Financial Economics 144(1): 247-272.[[Link\]](https://www.sciencedirect.com/science/article/pii/S0304405X21002543)\
LeSage, J. P. & Pace, R. K. Introduction to Spatial Econometrics, Boca Raton, Taylor & Francis, 2009.