部分线性变系数空间自回归固定效应面板数据模型的估计

Estimation of Partially linear Varying Coefficient Spatial Autoregressive Panel Data Models with Fixed Effects

摘要: 为了提高线性空间自回归面板数据模型的适用性和灵活度，提出了一类固定效应部分线性变系数空间自回归面板数据模型。基于局部线性技术，利用Profile似然方法构造了模型中未知空间滞后参数、回归系数以及非参数系数函数的估计。利用数值模拟考察了所提估计方法在有限样本下的表现。

Abstract: In order to increase the model adaptability and flexibility of the standard linear spatial autoregressive panel data models, a fixed effects partially linear varying coefficient spatial autoregressive panel data model is proposed. A Profile maximum likelihood approach based on the local-linear method is introduced to estimate the unknown spatial lag parameter, regression coefficients and nonparametric coefficient functions. Some simulations are conducted to examine the performance of the proposed procedures with finite sample.

1. 引言

近年来，空间面板数据模型得到了统计学家和计量经济学家的重视，在理论上得到了深入的研究，并被广泛地应用到各类实际问题分析中，详细介绍可参考综述文章Elhorst (2003) [1] 和Anselin (2008) [2]，著作Elhorst (2014) [3] 系统介绍了空间面板数据模型的设定、估计和检验问题。然而，上述文献大都假定了因变量与自变量的关系为参数结构特别是线性形式上。众所周知，数据分析中模型形式的错误设定会导致后面推断结果偏离实际情况。为了更好地探求因变量与自变量之间蕴含的复杂关系，提高模型的灵活性，多种非参数和半参数建模方法在近二十年来相继被提出，并被应用到各类复杂数据分析中。自然，利用半参数建模方法分析空间面板数据得到了关注，最近多类半参数空间面板数据模型被提出和研究。基于部分线性模型形式，Hu等(2014) [4] 提出了一类部分线性空间误差自回归固定效应面板数据模型，Zhang和Yang (2015) [5] 则构造了部分线性空间滞后回归固定效应面板数据模型。Bai等(2015) [6] 研究了一类部分线性变系数空间误差自回归面板数据模型，Zhang和Shen (2015) [7] 则提出了部分线性变系数空间滞后随机效应面板数据模型。上面几篇文献都是基于级数(样条)方法构造了模型的估计。

在众多的半参数模型中，部分线性变系数模型作为部分线性模型和变系数模型的推广得到了广泛的研究，其中Zhang和Shen (2015) [7] 基于这一建模方法提出了部分线性变系数空间滞后随机效应面板数据模型，并且基于样条方法和工具变量估计方法给出了模型中未知参数分量和非参数分量的估计。不同于Zhang和Shen (2015) [7] 所设定的随机效应模型，本文考虑如下的部分线性变系数空间自回归面板数据固定效应模型

$y_{it}=\rho {\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}{y}_{jt}}+{\alpha }_i+v_{it}^{\text{T}}\theta \left(u_{it}\right)+x_{it}^{\text{T}}\beta +{\epsilon }_{it},i=1,2,\cdots ,n;t=1,2,\cdots ,T$ (1)

其中 $y_{it}$ 是第i个空间观测单元在时刻的因变量观测值， $u_{it}$ 和是对应的自变量观测值，也为对应的自变量观测值，不失一般性，我们假设其为一维。 $W_n=\left(w_{ij}\right),1\le i,j\le n$ 是 $n\times n$ 的空间权重矩阵， $\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_p\right)}^{\text{T}}$ 和 $\rho$ 为未知参数， $\theta (\cdot )={\left({\theta }_1(\cdot ),{\theta }_2(\cdot ),\cdots ,{\theta }_q(\cdot )\right)}^{\text{T}}$ 是q维未知函数。误差 ${\epsilon }_{it}$ 相互独立，有 $E\left({\epsilon }_{it}\right)=0,Var\left({\epsilon }_{it}\right)={\sigma }^2$，且与自变量无关。固定效应 ${\alpha }_i$ 描述个体之间的差异性，为了模型可被识别，对固定效应附加如下条件

${\displaystyle \underset{i=1}{\overset{n}{\sum }}{\alpha }_i=0}$

显然，模型(1)包含几类常见的模型。如果 $\rho =0$，那么模型变化一个没有考虑空间效应的部分线性变系数固定效应面板数据模型。如果 $\theta (\cdot )=0$，那么模型则为空间自回归固定效应面板数据模型，相关研究可参考Lee和Yu (2010) [8] 和Debarsy和Ertur (2010) [9]。如果对于每个观测单元只观察一次，即 $T=1$，模型则变为Su和Jin (2010) [10] 针对截面数据提出的部分线性空间自回归模型。

对于模型(1)，我们重点研究其估计问题，不同于Zhang和Shen (2015) [7] 采用的工具变量方法，本文采用Su和Jin (2010) [10] 中的Profile极大似然估计方法。

2. 模型的Profile极大似然估计

为了方便叙述，我们定义

$Y=\left[\begin{array}{c}{y}_{11}\\ y_{12}\\ ⋮\\ y_{nT}\end{array}\right]$， $X=\left[\begin{array}{c}{x}_{11}^{\text{T}}\\ x_{12}^{\text{T}}\\ ⋮\\ x_{nT}^{\text{T}}\end{array}\right]$， $M=\left[\begin{array}{c}{v}_{11}^{\text{T}}\theta \left(u_{11}\right)\\ v_{12}^{\text{T}}\theta \left(u_{12}\right)\\ ⋮\\ v_{nT}^{\text{T}}\theta \left(u_{nT}\right)\end{array}\right]$， $\epsilon =\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{12}\\ ⋮\\ {\epsilon }_{nT}\end{array}\right]$

因此模型(1)可以写为

$Y=\rho \left(I_T\otimes W_n\right)Y+Z\gamma +X\beta +M+\epsilon$ (2)

其中， $Z={\left(z_{11},z_{12},\cdots ,z_{nT}\right)}^{\text{T}}=\left(I_n\otimes 1_T\right){\left(-1_{n-1}{I}_{n-1}\right)}^{\text{T}}$, $\gamma ={\left({\alpha }_2,\cdots ,{\alpha }_n\right)}^{\text{T}}$。 $I_s$ 为维数为s的单位阵， $1_k$ 为元素全为1的k维列向量。

假设 $\epsilon ~\left(0,{\sigma }^2{I}_{nT}\right)$，定义 $\varphi =\left({\beta }^{\text{T}},{\gamma }^{\text{T}}\rho ,{\sigma }^2\right)$， $A\left(\rho \right)=I_{nT}-\rho \left(I_T\otimes W_n\right)$， $\bar{Y}=Y-\rho \left(I_T\otimes W_n\right)Y=A\left(\rho \right)Y$，那么模型(2)的对数似然函数为：

$\ln L\left(\varphi ,\theta (\cdot )\right)=-\frac{nT}{2}\ln \left(2\pi \right)-\frac{nT}{2}\ln \left({\sigma }^2\right)+\ln \left|A\left(\rho \right)\right|-\frac{1}{2{\sigma }^2}{\left(\bar{Y}-Z\gamma -X\beta -M\right)}^{\text{T}}\left(\bar{Y}-Z\gamma -X\beta -M\right)$ (3)

显然，上述对数似然函数中含有未知的非参数函数，因此无法直接求出未知参数的极大似然估计。接下来，将使用Profile似然方法来估计其中参数 $\varphi$ 和未知非参数函数 $\theta (\cdot )$。如果参数 $\beta ,\gamma$ 和 $\rho$ 是已知的，那么模型(1)可以为

$y_{it}^{*}=v_{it}^{\text{T}}\theta \left(u_{it}\right)+{\epsilon }_{it},i=1,2,\cdots ,n;t=1,2,\cdots ,T$ (4)

其中 $y_{it}^{*}=y_{it}-{\alpha }_i-\rho {\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}{y}_{jt}}-x_{it}^{\text{T}}\beta$。显然，模型(4)是一个标准的变系数模型。下面将采用局部线性方法

对模型(4)进行估计。假设 $\{{\theta }_j(\cdot ),j=1,2,\cdots ,q\}$ 存在连续二阶导数。那么，对u小邻域内的 $u_0$，由Taylor展开可得

${\theta }_j\left(u\right)\approx {\theta }_j\left(u_0\right)+{{\theta }^{\prime }}_j\left(u_0\right)\left(u-u_0\right),j=1,2,\cdots ,q$

其中 ${{\theta }^{\prime }}_j\left(u\right)=\partial {\theta }_j\left(u\right)/\partial u$。

由局部加权最小二乘方法， $\{\left({\theta }_j\left(u_0\right),{{\theta }^{\prime }}_j\left(u_0\right)\right),j=1,2,\cdots ,q\}$ 的估计可通过使得

${\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{t=1}{\overset{T}{\sum }}{[y_{it}^{*}-{\displaystyle \underset{j=1}{\overset{q}{\sum }}\{{\theta }_j\left(u_0\right)+{{\theta }^{\prime }}_j\left(u_0\right)\left(u_{it}-u_0\right)v_{itj}\}}]}^2{K}_h\left(u_{it}-u_0\right)}}$ (5)

达到最小而得到，其中，Ｋ是核函数，h是窗宽，且 $K_h(\cdot )=K\left(\cdot /h\right)/h$。

令

$D_{z_0}=\left[\begin{array}{cc}{v}_{11}^{\text{T}}& \left(u_{11}-u_0\right)v_{11}^{\text{T}}\\ v_{12}^{\text{T}}& \left(u_{12}-u_0\right)v_{12}^{\text{T}}\\ ⋮& ⋮\\ v_{nT}^{\text{T}}& \left(u_{nT}-u_0\right)v_{nT}^{\text{T}}\end{array}\right]$， $S=\left[\begin{array}{c}\left(v_{11}^{\text{T}}\text{0}\right){\left\{D_{u_{11}}^{\text{T}}{K}_{u_{11}}{D}_{u_{11}}\right\}}^{-1}{D}_{u_{11}}^{\text{T}}{K}_{u_{11}}\\ \left(v_{1\text{2}}^{\text{T}}\text{0}\right){\left\{D_{u_{1\text{2}}}^{\text{T}}{K}_{u_{1\text{2}}}{D}_{u\text{12}}\right\}}^{-1}{D}_{u_{1\text{2}}}^{\text{T}}{K}_{u_{1\text{2}}}\\ ⋮\\ \left(v_{nT}^{\text{T}}\text{0}\right){\left\{D_{u_{nT}}^{\text{T}}{K}_{u_{nT}}{D}_{u_{nT}}\right\}}^{-1}{D}_{u_{nT}}^{\text{T}}{K}_{u_{nT}}\end{array}\right]$。

则基于(5)可得

${\left[{\bar{\theta }}_1\left(u_0\right),\cdots ,{\bar{\theta }}_q\left(u_0\right),{{\bar{\theta }}^{\prime }}_1\left(u_0\right),\cdots ,{{\bar{\theta }}^{\prime }}_q\left(u_0\right)\right]}^{\text{T}}={\left\{D_{u_0}^{\text{T}}{K}_{u_0}{D}_{u_0}\right\}}^{-1}{D}_{u_0}^{\text{T}}{K}_{u_0}\left(\bar{Y}-Z\gamma -X\beta \right)$ (6)

将 $u_{11},u_{12},\cdots ,u_{nT}$ 代替 $u_0$，可得 $\theta \left(u_{it}\right)$ 的估计，从而可定义M的估计为

$\bar{M}={\left[v_{11}^{\text{T}}\bar{\theta }\left(u_{11}\right),v_{12}^{\text{T}}\bar{\theta }\left(u_{12}\right),\cdots ,v_{nT}^{\text{T}}\bar{\theta }\left(u_{nT}\right)\right]}^{\text{T}}=S\left(\bar{Y}-Z\gamma -X\beta \right)$ (7)

将 $\bar{M}$ 代替(3)中的M，令 $H=\left(ZX\right)$， $\phi =\left({\gamma }^{\text{T}}{\beta }^{\text{T}}\right)$，可得如下Profile对数似然函数

$\begin{array}{l}\log L\left(\phi \right)=-\frac{nT}{2}\ln \left(2\pi \right)-\frac{nT}{2}\ln \left({\sigma }^2\right)+\ln \left|A\left(\rho \right)\right|-\frac{1}{2{\sigma }^2}{\left(\bar{Y}-Z\gamma -X\beta -\bar{M}\right)}^{\text{T}}\left(\bar{Y}-Z\gamma -X\beta -\bar{M}\right)\\ =-\frac{nT}{2}\ln \left(2\pi \right)-\frac{nT}{2}\ln \left({\sigma }^2\right)+\ln \left|A\left(\rho \right)\right|-\frac{1}{2{\sigma }^2}{\left[\left(I_{nT}-S\right)\left(\bar{Y}-Z\gamma -X\beta \right)\right]}^{\text{T}}\left[\left(I_{nT}-S\right)\left(\bar{Y}-Z\gamma -X\beta \right)\right]\\ =-\frac{nT}{2}\ln \left(2\pi \right)-\frac{nT}{2}\ln \left({\sigma }^2\right)+\ln \left|A\left(\rho \right)\right|-\frac{1}{2{\sigma }^2}{\left[\left(I_{nT}-S\right)\left(\bar{Y}-H\phi \right)\right]}^{\text{T}}\left[\left(I_{nT}-S\right)\left(\bar{Y}-H\phi \right)\right]\end{array}$ (8)

给定 $\rho$，使 $\log L\left(\varphi \right)$ 对 $\phi ,{\sigma }^2$ 分别求偏导，得到如下方程组

$\{\begin{array}{l}\frac{\partial \log L\left(\varphi \right)}{\partial \phi }=0\Rightarrow V^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}\left(I_{nT}-S\right)H\phi =H^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}\left(I_{nT}-S\right)\bar{Y}\\ \frac{\partial \log L\left(\varphi \right)}{\partial {\sigma }^2}=0\Rightarrow n{\sigma }^2={\left[\left(I_{nT}-S\right)\left(\bar{Y}-H\phi \right)\right]}^{\text{T}}\left[\left(I_{nT}-S\right)\left(\bar{Y}-H\phi \right)\right]\end{array}$

解得上述方程组，得 $\phi$ 的Profile极大似然估计为

$\hat{\phi }\left(\rho \right)={\left[H^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}\left(I_{nT}-S\right)V\right]}^{-1}{H}^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}\left(I_{nT}-S\right)\bar{Y}$ (9)

同时可得到的Profile极大似然估计为

$\begin{array}{c}{\hat{\sigma }}^2\left(\rho \right)=\frac{1}{nT}{\left\{\left(I_{nT}-S\right)\left[\bar{Y}-H\hat{\phi }\left(\rho \right)\right]\right\}}^{\text{T}}\left\{\left(I_{nT}-S\right)\left[\bar{Y}-H\hat{\phi }\left(\rho \right)\right]\right\}\\ =\frac{1}{nT}{\bar{Y}}^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}D\left(I_{nT}-S\right)\bar{Y}\end{array}$ (10)

其中 $D=I_{nT}-\left(I_{nT}-S\right)H{\left[H^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}\left(I_{nT}-S\right)H\right]}^{-1}{H}^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}$。进一步，关于 $\rho$ 的极大似然函数为

$\log L\left(\rho \right)=-\frac{nT}{2}\left[\log \left(2\pi \right)+1\right]-\frac{nT}{2}\log \left[{\hat{\sigma }}^2\left(\rho \right)\right]+\log \left|A\left(\rho \right)\right|$ (11)

利用优化算法通过使 $\log L\left(\rho \right)$ 最大化，得到 $\rho$ 的估计 $\hat{\rho }$。从而，得到 $\phi$ 和 ${\sigma }^2$ 的最终估计，即 $\hat{\phi }=\hat{\phi }\left(\hat{\rho }\right)$， ${\hat{\sigma }}^2={\hat{\sigma }}^2\left(\hat{\rho }\right)$。因此，得到 $\beta$ 和 $\gamma$ 的最终估计，分别为：

$\hat{\beta }={\left[X^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}P\left(I_{nT}-S\right)X\right]}^{-1}{X}^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}P\left(I_{nT}-S\right)\left[I_{nT}-\hat{\rho }\left(I_T\otimes W_n\right)\right]Y$ (12)

$\hat{\gamma }=\left({\hat{\alpha }}_2,\cdots ,{\hat{\alpha }}_n\right)={\left[Z^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}\left(I_{nT}-S\right)Z\right]}^{-1}{Z}^{\text{T}}{\left(I_{nT}-S\right)}^{\text{T}}\left(I_{nT}-S\right)\left[Y-X\hat{\beta }\right]$ (13)

其中 $P=I_{nT}-Z^{*}{\left(Z^{*\text{T}}{Z}^{*}\right)}^{-1}{Z}^{*\text{T}}$， $Z^{*}=\left(I_{nT}-S\right)Z$， ${\alpha }_1$ 的估计为 $-{\displaystyle {\sum }_{i=2}^n{\alpha }_i}$。

最后，将 $\hat{\beta },\hat{\gamma },\hat{\rho }$ 代入(6)，我们可以得到 $\theta \left(u_0\right)$ 的最终估计为

$\hat{\theta }\left(u_0\right)=\left(I_q\text{0}\right){\left\{D_{u_0}^{\text{T}}{K}_{u_0}{D}_{u_0}\right\}}^{-1}{D}_{u_0}^{\text{T}}{K}_{u_0}\left[Y-\hat{\rho }\left(I_T\otimes W_n\right)Y-Z\hat{\gamma }-X\hat{\beta }\right]$ (14)

其中 $I_q$ 为维数为q的单位阵， $0_q$ 是 $q\times q$ 的零矩阵。最后得到残差向量的估计为

$\hat{\epsilon }=Y-\hat{\rho }\left(I_T\otimes W_n\right)Y-Z\hat{\gamma }-X\hat{\beta }-\hat{M}$ (15)

其中 $\hat{M}={\left[v_{11}^{\text{T}}\hat{\theta }\left(u_{11}\right),v_{12}^T\hat{\theta }\left(u_{12}\right),\cdots ,v_{nT}^{\text{T}}\hat{\theta }\left(u_{nT}\right)\right]}^{\text{T}}=S\left(\bar{Y}-Z\gamma -X\beta \right)$。

3. 数值模拟

假定研究区域是边长为 $m-1$ 个距离单位的正方形，n个空间观测单元正好在 $m\times m$ 个格子点上，权重矩阵W设定为rook矩阵和queen矩阵两种形式。考虑如下部分线性空间自回归面板数据模型

$y_{it}=\rho {\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}{y}_{jt}}+{\alpha }_i+v_{it}\theta \left(u_{it}\right)+x_{it}^{\text{T}}\beta +{\epsilon }_{it},i=1,2,\cdots ,n;t=1,2,\cdots ,T$

假设 $x_{it}~N\left(1,1\right)$， $v_{it}~U\left(-2,2\right)$， $u_{it}~U\left(0,1\right)$， $\theta \left(u_{it}\right)=2\cos \left(2\pi u_i\right)+1$， $\beta =2$，固定效应 ${\alpha }_i~U\left(0,1\right)$, $i=2,\cdots ,n$。上一节模型估计过程中假定了误差服从正态分布，但这只是为了写出似然函数，为了考察误差分布对模型估计效果的影响，设定如下三种形式的模型误差，都满足均值为0，方差为0.25。

1) ${\epsilon }_{it}~U\left(0,{0.5}^2\right)$ ; 2) ${\epsilon }_{it}~U\left(-\sqrt{3}/2,\sqrt{3}/2\right)$ ; 3) ${\epsilon }_{it}~\left(\frac{1}{8}{\chi }_8^2-1\right)$，计算过程中选择高斯核函数 $K\left(z\right)=\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{z^2}{2}\right)$，平滑指数取 $s_u{\left(nT\right)}^{-1/5}$，其中 $s_u$ 是 $u_{11},u_{12},\cdots ,u_{nT}$ 的标准差。

样本量分别取 $n=64,100$，令时间维度 $T=3$， $\rho =\left(\pm 0.9,\pm 0.5,\pm 0.1,0\right)$，针对每种设定，模拟重复500次，每次都利用前面介绍的估计方法给出参数分量和非参数分量的估计。对于参数分量 $\rho ,\beta$，以“MEAN”和“SD”分别表示其500次估计值的均值和标准差。运用Matlab软件， $\rho ,\beta$ 的模拟结果见表1和表2。

Table 1. Estimation of ρ , β under sample size n = 64 of system resulting data of standard experiment

表1. 标准试验系统结果数据样本量n = 64下 $\rho ,\beta$ 的估计

Table 2. Estimation of ρ , β under sample size n = 100 of system resulting data of standard experiment

表2. 标准试验系统结果数据样本量n = 100下 $\rho ,\beta$ 的估计

对于非参数系数函数 $\theta (\cdot )$，利用如下的根均方误差(RASE)考察其表现，

$\text{RASE}\left(\theta \right)={\left\{\frac{1}{500\ast n\ast T}{\displaystyle \underset{l=1}{\overset{500}{\sum }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{T=1}{\overset{T}{\sum }}{\left[{\hat{\theta }}^l\left(u_{it}\right)-\theta \left(u_{it}\right)\right]}^2}}}\right\}}^{1/2}$,

其中 ${\hat{\theta }}^l\left(u_{it}\right),l=1,2,\cdots ,500$ 是 $\theta \left(u_{it}\right)$ 在第次重复时的估计值。模拟结果见表3和表4。

Table 3. RASE of coefficient function estimation when n = 64

表3. n = 64时系数函数估计的RASE

Table 4. RASE of coefficient function estimation when n = 100

表4. n = 100时系数函数估计的RASE

从模拟结果可以看出，第一，对于参数分量，估计值的均值非常接近真实值，而且随着样本量增加，Mean更靠近对应的真实值，SD变小。模型误差和空间权重矩阵W的设定对估计结果影响不大。第二，对于非参数系数函数的估计，RASE值针对不同的误差分布和W的设定变化不大，随着样本量的增加，RASE值总体上在变小。总之，本文所提方法对于参数分量和非参数分量的估计在有限样本下表现良好。

4. 总结

本文在传统的空间面板数据模型基础上提出了一类固定效应部分线性变系数空间自回归面板数据模型，并且给出了模型的估计。由于模型本身的复杂性，本文没有研究所提估计量的渐近性质，而是利用数值模拟考察了所提估计量在有限样本情况下的表现，模拟结果表明所提估计方法有效。此外，论文重点研究了模型的估计，关于空间效应的检验没有涉及。这些问题将作为我们后续研究的目标。

基金项目

国家社科基金项目(21BTJ005)。