空间自回归模型^[1]凭借其对于数据空间自相关性的有效解释, 能够很好的弥补经典自回归模型在空间数据违背独立性假设情况下的数据解释问题。由于数据多元化时代的来临, 该模型已经在统计学领域受到了广泛的关注, 并积极的应用于经济政策分析的各类领域, 特别是实产和房地产经济、环境与资源经济等^[2]此类呈现强烈自相关性的空间数据。而部分线性空间自回归模型^[3-4]通过引入未知非参数项, 能够更加灵活的解释某些潜在协变量对于响应变量的影响, 其模型结构为：

$ \boldsymbol{Y}=\boldsymbol{\rho} \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{m}(\boldsymbol{Z})+\boldsymbol{\varepsilon}。$

(1)

其中, Y=(y₁, y₂, …y_n)^T表示响应变量的n×1维观测向量; W表示预先设定的n×n的空间权重矩阵; WY表示响应变量Y的空间滞后项; X=(x₁, x₂, …x_n)^T表示协变量的n×q的观测矩阵, 其中x_i=(x_i1, x_i2, …, x_iq)^T为q×1维列向量, i=1, 2, …, n; m(Z)=(m(z₁), m(z₂), …m(z_n))^T表示n×1的未知函数向量, 其中Z=(z₁, z₂, …z_n)^T表示某个潜在变量的n×1维观测向量, m(·)是一个未知光滑函数; ρ表示空间自回归参数且满足|ρ| < 1;β =(β₁, β₂, …β_q)^T表示q×1的未知参数向量。此外, ε=(ε₁, ε₂, …ε_n)^T是一个n×1的残差向量其中ε_i, i=1, 2, …, n相互独立同分布。

目前对于此模型未知函数m(Z)有两种主流估计方法, 一是由Fan和Gijbels^[5]提出的局部多项式拟合来对m(Z)进行非参数估计; 二是由Schumaker^[6]提出的B样条(B-spline)逼近技术对m(Z)进行非参数估计。

在对模型(1)的检验研究过程中, 本文则采用Schumaker提出的B样条逼近技术对m(Z)进行处理。对模型(1)的估计问题目前已有部分文献进行了研究, 比如, Su和Jin^[7]利用剖面最大拟似然估计方法研究了模型(1)的估计问题。Koch和Krisztin^[8]结合B样条逼近技术和遗传算法研究了模型(1)的估计问题。Su^[9]则利用广义矩估计方法(GMM)研究了模型(1)的估计问题。

但是关于模型(1)的检验问题, 目前还没有相关文献进行系统研究。为此, 本文将利用广义似然比检验方法研究模型(1)的检验问题。具体地, 本文首先结合B样条逼近方法以及矩阵的正交投影技术消除模型的非参数分量对模型参数分量检验精度的影响, 进而构造了一个关于模型参数分量的广义似然比检验统计量。另外由于模型结构的复杂性, 理论上很难得到检验统计量的渐近分布。为此, 本文利用Bootstrap^[10]自助抽样法给出了一个计算检验统计量p值的迭代算法。

1 基于B样条投影矩阵的模型参数估计

B样条是数学的子学科数值分析里样条曲线一种特殊的表示形式, 其分片多项式的构造、保持一定的连续性等性质使其拟合性能优越。具体地, 记t₀, t₁, …t_kn为区间[0, 1]内的节点, 并记π(t)=(B₁(t), B₂(t), …B_kn+l+1(t))^T为一组l阶B样条基函数, 则非参数函数m(t)可渐近表示为：

$ m(t) \approx \sum\limits_{j=1}^{k_{n}+l+1} B_{j}(t) \alpha_{j}=\pi^{\mathrm{T}}(t) \alpha, $

(2)

其中α=(α₁, α₂, …α_kn+l+1)^T为样条系数向量。结合式(2)可得:

$ \begin{aligned} \boldsymbol{m}(\boldsymbol{Z}) &=\left(\boldsymbol{\pi}\left(z_{1}\right), \pi\left(z_{2}\right), \cdots, \pi\left(z_{n}\right)\right)^{\rm{T}} \alpha \\ \triangleq \Pi \alpha。& \end{aligned} $

进而模型(1)可以转化为：

$ \boldsymbol{Y}=\boldsymbol{\rho} \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{\Pi} \boldsymbol{\alpha}+\boldsymbol{\varepsilon}+\boldsymbol{e}, $

(3)

其中e=m(Z)-Πα表示B样条拟合误差。令P=Π(Π^TΠ)^-1Π^T表示Π的正交投影矩阵, 则(I_n-P)Π=Π-PΠ=Π-Π=0。在模型(3)两边同乘以矩阵S=I_n-P, 则模型(3)可以转化为：

$ \boldsymbol{S} \boldsymbol{Y}=\boldsymbol{S}(\boldsymbol{\rho} \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{X} \boldsymbol{\beta})+\boldsymbol{S}(\boldsymbol{\varepsilon}+\boldsymbol{e}), $

(4)

令Q=(WY, X), θ=(ρ, β^T)^T, $\mathit{\boldsymbol{\bar Y}} = \mathit{\boldsymbol{SY}}, \mathit{\boldsymbol{\overline X}} = \mathit{\boldsymbol{SQ}}$, $\mathit{\boldsymbol{\bar \varepsilon }} = \mathit{\boldsymbol{S}}(\mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{e}})$, 则模型(4)可写成如下简约模型：

$ \overline{\boldsymbol{Y}}=\overline{\boldsymbol{X}} \boldsymbol{\theta}+\overline{\boldsymbol{\varepsilon}}, $

(5)

注意到式(5)为经典线性模型, 利用最小二乘方法可得如下估计：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\widehat \theta }} = {{\left( {{{\mathit{\boldsymbol{\overline X}} }^{\bf{T}}}\mathit{\boldsymbol{\overline X}} } \right)}^{ - 1}}{{\mathit{\boldsymbol{\overline X}} }^{\bf{T}}}\mathit{\boldsymbol{\overline Y}} }\\ {\mathit{\boldsymbol{\widehat \alpha }} = {{\left( {{{\bf{\Pi }}^{\bf{T}}}{\bf{\Pi }}} \right)}^{ - 1}}{{\bf{\Pi }}^{\bf{T}}}(\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{Q}}\mathit{\boldsymbol{\widehat \theta }})} \end{array}} \right.。$

(6)

接下来考虑σ²的估计参数, 设ε+e服从零均值同方差的正态分布, 可得关于σ²的正态拟对数似然函数为：

$ \begin{aligned} &\ln L\left(\sigma^{2}\right)=-\frac{n}{2} \mathrm{l}_{\mathrm{n}}\left(2 \pi \sigma^{2}\right)-\ln \left(\left|I_{n}-\boldsymbol{\rho} \boldsymbol{W}\right|\right)- \\ &\frac{1}{2 \sigma^{2}} \times(\boldsymbol{Y}-\hat{\rho} \boldsymbol{W} \boldsymbol{Y}-\boldsymbol{X} \hat{\boldsymbol{\beta}}-\boldsymbol{\Pi} \hat{\boldsymbol{\alpha}})^{\mathrm{T}} \times \\ &(\boldsymbol{Y}-\hat{\boldsymbol{\rho}} \boldsymbol{W} \boldsymbol{Y}-\boldsymbol{X} \hat{\boldsymbol{\beta}}-\boldsymbol{\Pi} \hat{\boldsymbol{\alpha}}), \end{aligned} $

(7)

最大化式(7)可得σ²的估计为：

$ {\hat \sigma ^2} = \frac{1}{n}{(Y - Q\hat \theta - \Pi \hat \alpha )^{\rm{T}}}(Y - Q\hat \theta - \Pi \hat \alpha )。$

(8)

2 基于似然函数的模型检验及Bootstrap实现

考虑如下参数显著性检验问题：

$ H_{0}: \boldsymbol{A \theta}=\boldsymbol{b} \text { vs } H_{1}: \boldsymbol{A \theta} \neq \boldsymbol{b}, $

(9)

其中, A为给定的d×(q+1)的列满秩矩阵, b为一个已知的d×1维向量。此检验的目的是检验空间自回归系数ρ和回归参数向量β是否显著。

结合正态拟对数似然函数式(7)以及Neyman-Pearson引理, 并利用与Fan等^[11]类似的广义似然比检验统计量推导方法, 对检验问题(9)本文给出一个基于B样条投影矩阵的广义似然比统计量。具体地, 当原假设H₀: Aθ=b成立时, 即在Aθ=b限制条件下最大化式(7), 可以得到θ的一个估计记为$\hat{\boldsymbol{\theta}}_{0}$=${\left({{{\mathit{\boldsymbol{\widehat \beta }}}_0}, \mathit{\boldsymbol{\widehat \beta }}_0^{\bf{T}}} \right)^{\bf{T}}}$, 再结式(6)可得参数向量α的估计记为${\mathit{\boldsymbol{\widehat \alpha }}_0}$, 并由式(8)可得σ²的估计：

$ \begin{array}{l} \widehat \sigma _0^2 = \frac{1}{n}{\left( {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{Q}}{{\mathit{\boldsymbol{\widehat \theta }}}_0} - {\bf{\Pi }}{{\mathit{\boldsymbol{\widehat \alpha }}}_0}} \right)^{\bf{T}}}\left( {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{Q}}{{\mathit{\boldsymbol{\widehat \theta }}}_0} - } \right.\\ \left. {{\bf{\Pi }}{{\mathit{\boldsymbol{\widehat \alpha }}}_0}} \right)。\end{array} $

进而可得在原假设H₀下的最大拟对数似然函数为：

$ \begin{gathered} l_{n}\left(H_{0}\right)=-\frac{n}{2}[\ln (2 \pi)+1]-\frac{n}{2} \times \\ \ln \left\{\frac { 1 } { n } ( Y - Q \hat { \theta } _ { 0 } - \Pi \hat { \alpha } _ { 0 } ) ^ { \mathrm { T } } \left(\boldsymbol{Y}-\boldsymbol{Q} \hat{\boldsymbol{\theta}}_{0}-\right.\right. \\ \left.\left.\mathit{\boldsymbol{ \boldsymbol{\varPi} }}\hat{\boldsymbol{\alpha}}_{0}\right)\right\}+\ln \left(\left|I_{n}-\hat{\rho}_{0} W\right|\right) 。\end{gathered} $

(10)

另外在备择假设H₁下, 即在没有约束条件Aθ=b的情况下, 通过第1节的估计过程可得到参数θ和α由式(6)所定义的估计${\mathit{\boldsymbol{\widehat \theta }}}$ 和${\mathit{\boldsymbol{\widehat \alpha }}}$。进而可得在备择假设H₁下的最大拟对数似然函数^[12]为：

$ \begin{aligned} &l_{n}\left(H_{1}\right)=-\frac{n}{2}[\ln (2 \pi)+1]-\frac{n}{2} \times \\ &\ln \left\{\frac{1}{n}(\boldsymbol{Y}-\boldsymbol{Q} \hat{\boldsymbol{\theta}}-\boldsymbol{\Pi} \hat{\boldsymbol{\alpha}})^{\bf{T}}(\boldsymbol{Y}-\boldsymbol{Q} \hat{\boldsymbol{\theta}}-\mathit{\boldsymbol{ \boldsymbol{\varPi}}} \hat{\boldsymbol{\alpha}})\right\}+ \\ &\ln \left(\left|I_{n}-\hat{\rho} W\right|\right), \end{aligned} $

(11)

因此结合式(10)和式(11), 并利用与Fan和Jiang^[11]类似的方法, 可定义广义似然比检验统计量为：

$ \begin{aligned} &T=l_{n}\left(H_{1}\right)-l_{n}\left(H_{0}\right)=l_{n}(\hat{\theta}, \hat{\alpha})-l_{n}\left(\hat{\boldsymbol{\theta}}_{0}, \hat{\boldsymbol{\alpha}}_{0}\right) \\ &=\frac{n}{2} \ln \frac{\left(\boldsymbol{Y}-\boldsymbol{Q} \dot{\boldsymbol{\theta}}_{0}-\mathit{\boldsymbol{ \boldsymbol{\varPi}}} \hat{\boldsymbol{\alpha}}_{0}\right)^{\bf{T}}\left(\boldsymbol{Y}-\boldsymbol{Q} \dot{\boldsymbol{\theta}}_{0}-\mathit{\boldsymbol{ \boldsymbol{\varPi}}} \dot{\boldsymbol{\alpha}}_{0}\right)}{(\boldsymbol{Y}-\boldsymbol{Q} \hat{\boldsymbol{\theta}}-\mathit{\boldsymbol{ \boldsymbol{\varPi}}} \hat{\boldsymbol{\alpha}})^{\bf{T}}(\boldsymbol{Y}-\boldsymbol{Q} \dot{\boldsymbol{\theta}}-\mathit{\boldsymbol{ \boldsymbol{\varPi}}} \hat{\boldsymbol{\alpha}})}+ \\ &\ln \frac{\left|I_{n}-\hat{\boldsymbol{\rho}} \boldsymbol{W}\right|}{\left|I_{n}-\hat{\boldsymbol{\rho}}_{0} \boldsymbol{W}\right|}, \end{aligned} $

(12)

相应地, 检验p值的计算为：

$ p=P_{H_{0}}(T>t), $

(13)

其中, P_H0表示在原假设H₀下的累积概率, t表示检验统计量T的观测值。对任一给定的检验显著性水平α, 如果p < α则拒绝原假设, 否则就不能拒绝原假设。但是由于模型结构的复杂性, 在计算p值时, 检验统计量T的渐进分布很难直接的得到。因此本文采用Bootstrap自助法对检验p值的计算给出一个迭代算法, 具体步骤如下：

第一步：如本文第1节所述, 基于已知的观察样本{Y, X, Z}、预先给定的空间权重矩阵W和原假设(9)中给定的参数, 通过式(12)可以得到检验统计量T的一个观测值t^*和一个残差序列$\hat{\varepsilon}=\boldsymbol{Y}-\hat{\boldsymbol{Q} \boldsymbol{\theta}}-\hat{\mathit{\boldsymbol{ \boldsymbol{\varPi}}} \boldsymbol{\alpha}}=$ =$\left(\hat{\varepsilon}_{1}, \hat{\varepsilon}_{2}, \cdots, \hat{\varepsilon}_{n}\right)^{\mathrm{T}}$, 将$\hat{\varepsilon}$中心化可得$\hat{\boldsymbol{\varepsilon}}^{\prime}=$=$\left(\hat{\varepsilon}_{1}-\varepsilon^{\prime \prime}, \hat{\varepsilon}_{2}-\varepsilon^{\prime \prime}, \cdots, \hat{\varepsilon}_{n}-\varepsilon^{\prime \prime}\right)^{\mathrm{T}}$, 其中ε″=$n^{-1} \sum \hat{\varepsilon}_{i}$。

第二步：从$\hat{\boldsymbol{\varepsilon}}'$的经验分布函数中进行一个有放回的Bootstrap抽样得到一个Bootstrap样本$\hat{\varepsilon}^{*}=\left(\hat{\varepsilon}_{1}^{*}, \hat{\varepsilon}_{2}^{*}, \cdots, \hat{\varepsilon}_{n}^{*}\right)^{\mathrm{T}}$, 进一步可以定义一个新的因变量向量：

$ Y^{*}=\left(I_{n}-\hat{\rho}_{0} \boldsymbol{W}\right)^{-1}\left(\boldsymbol{X} \hat{\boldsymbol{\beta}}_{0}+\boldsymbol{\Pi} \hat{\boldsymbol{\alpha}}_{0}+\hat{\boldsymbol{\varepsilon}}^{*}\right)。$

(14)

第三步：由式(14)得到一个新的因变量向量Y^*, 则通过自助样本{Y^*, X, Z}可以得到一个新的检验统计量T的自助版本的检验统计量的值T^*。

第四步：将步骤2和步骤3重复m次, 可以得到m个检验统计量T的自助样本检验统计量的值, 记为T₁^*, T₂^*, …T_m^*, 进而式(13)定义的检验p值的估计可定义为：

$ \hat{p}=\frac{\#\left\{T_{i}^{*}: T_{i}^{*}>t^{*}\right\}}{m}, $

(15)

其中#{}表示集合中元素的个数。