一种基于N-W估计的异方差估计方法

张晓琴，郭雅静，李顺勇

（1.山西财经大学统计学院，山西 太原030006；2.河北建筑工程学院数理系，河北张家口075024；3.山西大学数学科学学院，山西太原030006）

一、引言

随机误差项的同方差性能够保证给出回归系数良好的估计[1]，但在实际问题的研究中（如截面数据），同方差性这一假设很难保证，即随机误差项具有异方差性[2]。异方差模型中，通过最小二乘（OLS）法得到的回归系数估计量仍具有良好的性质，但是，OLS估计的协方差阵不能保持无偏性，因此对于异方差的研究具有重大意义。

对异方差问题的处理，常常要先检验异方差性的存在[3，4]，再进行估计。常用的估计方法是广义最小二乘法（GLS）[5，6]，但其前提是随机误差项的协方差阵已知。当协方差阵为对角阵时，称为加权最小二乘（WLS）估计。在一般情况下，这一条件未知，因此需要先对其给出有效估计，目前常用方法有极大似然估计法[7]、Box-Cox 变换法[8]、贝叶斯估计法[9]、两阶段最小二乘估计[10]等。

由于异方差存在时，通过OLS方法得到的参数估计其协方差阵的估计量不能保证一致性，因此只要给出相应的一致估计量，就可以进行有效的假设检验。基于这一思想，White等提出了异方差一致协方差阵估计——HC0[11]，并得到了广泛的应用，但模拟发现，该估计在小样本中不再具有无偏性。对此，众多学者进行了改进，主要包括HC1、HC2、HC3、HC4、HC5、HC4m、HC5m 等，总称为 HCCMEs。其中，在HC0基础上提出的HC1、HC2和HC3三种估计都没有考虑到样本中有高杠杆点的情况[12]。基于HC4估计进行检验时，无论高杠杆点存在与否，效果都较好，但是在点估计中表现较差[13]。随后Cribari-Neto等提出的HC5估计虽然在杠杆程度较高的情况下效果更好，但在杠杆程度较低时效果不理想[14]。为解决这一问题，Cribarineto和Silva提出了HC4估计的改进形式——HC4m[15]。考虑到基于杠杆点估计的不足，李顺勇等在2016年提出了HC5m估计[16]，研究表明，无论杠杆程度如何，HC5m估计均表现良好。

利用正交表的优良性质，张晓琴等提出了一种非参数异方差估计——正交表（OR）估计[17]。随后，张晓琴等又提出了OR1估计[18]，用于修正OR估计中因变量选择以及容差选取的不足。

基于非参数方法对误差项方差的估计不依赖于模型中异方差存在形式的这一优点，Carroll在异方差结构未知时提出了一种的核估计，以达到提升估计准确性的目的[19]。本文将其引入异方差估计中，提出了一种新的估计异方差的方法，扩大了估计方法的普适性，并提升了估计方法的精确度。

二、异方差模型与两类已有的异方差估计方法

（一）异方差模型

假设解释变量（xi1，xi2，…，xip）和响应变量yi（i=1，2，…，n）满足如下线性回归模型：

写成矩阵形式为：

对模型（2）而言，其回归系数的GLS估计为：

当Σ未知时，无法得到βˆ，故需事先获得Σ的估计。

（二）两类已有的Σ估计方法

现阶段异方差估计的方法已有很多，本文仅介绍 HCCMEs中的 HC0、HC4、HC5 估计以及 OR、OR1估计。

1.HC0、HC4、HC5 估计。针对异方差模型（1），White提出的HC0估计为[11]：

Cribari-Neto提出的HC4估计为[13]：

其中，δi=min｛4，nhii/p｝，n 是设计矩阵 X 的行数，p+1 是 X 的列数，hii是“帽子矩阵”H=X（XTX）-1XT的对角线元素，hii的大小表示杠杆程度的高低，可用作自变量空间的异常值观测指标。若hii大于其均值两到三倍，则该值对应一个高杠杆点。

Cribari-Neto等提出的HC5估计为[14]：

其中，δi=min｛nhii/p,max｛4,nkhmax/p｝｝，hmax=max｛h11,h22,…,hnn｝，k 为常数且 0＜k＜1，当 k=0 时，该估计即为HC4估计。

2.OR、OR1估计。张晓琴等首次提出了一种基于正交表的非参数异方差估计——OR估计[17]，该方法借助正交表将模型进行扩展，使得原来的每组观测值都扩展为多个观测值，再基于扩展后的数据构建回归模型。由于正交表具有均匀整齐的性质，因此可使得原来的每组观测值对应扩展后的模型具有同方差性，在这种情况下可以直接利用OLS对回归系数进行估计，进而得到Σ的估计量。考虑到OR方法在因变量的扩展中对原始数据信息使用不充分以及容差确定的不合理性，张晓琴等对因变量的产生及容差的取法进行了改进，得到了OR1 估计[18]。

三、基于N-W估计的异方差估计方法

由于方差的非参数估计不依赖于模型中异方差的存在形式，因此估计的有效性更高。本文将在异方差估计中引入一种方差的非参数核估计方法，给出新的方差估计以及误差项协方差阵的估计量，进而用GLS估计未知参数，并对模型做出预测和拟合。

（一）N-W估计

对协方差阵Σ的估计，HCCMEs及OR估计都是从最小二乘残差平方切入进行研究，即直接用估计。在异方差形式未知时，Carroll提出的一个自适应估计可以用来提升估计的精度[19]。假设误差项的方差是yi均值的光滑函数，即=g（xiβ）=g（τ）i，其中，g（·）为未知函数，τi可以由估计得到。Carroll给出了一种的核估计——Nadaraya-Watson（N-W）估计，形式如下：

（二）新的异方差估计——KNW估计

存在异方差的情形下，需要首先对误差项的方差给出精确估计，进而才能对回归系数做出估计并对模型进行分析。N-W估计不依赖于异方差的构成形式，基于这一特点，我们认为在一般的异方差模型中，该方法更为简便并且估计的有效性更高。为了提升异方差估计方法的普适性及估计的精确度，本文在异方差模型中引入的N-W估计，然后应用GLS方法来估计未知参数，提出一种新的异方差估计，即 Kernel Nadaraya Watson（KNW）估计。

第一步，对于已得到的数据（xi1，xi2，…，xip，y）i，i=1，2，…n，先按照同方差假设得到回归系数的OLS估计，再将其带入异方差模型中得到如下回归方程：

其中，βˆOLS=（βˆ（OLS）0，βˆ（OLS）1，βˆ（OLS）2，Λ，βˆ（OLS）p）T是回归参数 β=（β0，β1，β2，Λ，βp）T的 OLS 估计。

第四步，得到协方差阵Σ的估计为：

第五步，将Σˆ带入GLS估计中，得到回归系数的GLS估计为：

第六步，将（10）式重新带入原始异方差模型（2）中，得到因变量 yi的预测值yˆi：

表1 六种核函数及函数表达式

（三）评价指标

为比较本文提出的KNW方法与已有异方差估计方法，本文采用了三个指标进行评价[22]。

第一个指标是方差的平均绝对误差，其表达式为：

第二个指标是预测值的平均绝对误差，其表达式为：

在（12）与（13）式中，m 为模拟次数，n为样本容量，代表真实的方差值，代表某种方法下方差的估计值，y代表真实的因变量值，代表因变量的预i测值。两个评价指标表示的是不同方法下估计值或预测值与真实值的接近程度，MAEΣ与MAEy的值越小说明对应的方法越有效。

第三个指标是拟合优度R2。考虑到在实际分析中，方差的真实值通常难以得到，这时拟合优度R2可以作为一个评判拟合效果的指标，表示为：

四、模拟实验及实例研究

本节从随机模拟和实例研究两个方面入手，将KNW估计与已有的异方差估计方法（包括HC0、HC4、HC5、OR、OR1估计）进行对比分析。

（一）模拟实验

在模拟实验中，考虑模型（1），令p=5，从已有研究对回归参数取值的经验可知，不同的参数对结果并无显著影响，本文取 β=（β0，β1，β2，β3，β4，β5）T=（0.2，1，1，1，1，1）T，σ2i=xi12/100，模型中自变量xij为取自均匀分布U（0，1）的随机数，样本量分别取50、100、1500。

下面分别用本文提出的 KNW与 HC0、HC4、HC5、OR、OR1方法对方差进行估计。其中，OR、OR1方法均借助正交表L2（556）对观测值进行展开，并且每个因子x（ijj=1，2，3，4，5）的容差 Δj均可表示为 Δj=min｛sj，1/sj，0.01｝，j=1，2，3，4，5。其中，sj为原始数据中第j列的标准差。KNW方法分别采用表1中给出的六种核函数，并固定窗宽h为。每种情形下模拟100次，对比结果见表2至表5。

从表2至表4可见，无论模拟实验中样本容量如何改变，无论选取哪种核函数，KNW估计与HC0、HC4、HC5、OR、OR1 估计相比，都有最小的MAEΣ值和更大的R2值，但MAEy值差别不大。这说明，KNW估计能更精确地估计协方差阵Σ，对模型有更高的拟合程度。此外，六种核函数下各种估计结果表明，方差的估计精度与核函数选取没有太大联系。

表5给出了应用高斯核函数时，在不同样本量下各种估计方法的运行时间。对比发现，无论样本容量如何变化，KNW估计的运行时间均显著少于OR估计方法，但与HCCME方法相比，KNW估计的时间成本要高一些。

表2 五个自变量模型中xij取自U（0，1）时方差的平均绝对误差值MAE∑

表3 五个自变量模型中xij取自U（0，1）时因变量预测值的平均绝对误差值MAEy

（续表3）

表4 六种核函数下不同估计方法的拟合优度R2

表5 高斯核函数下不同估计方法的运行时间 （单位：秒）

（二）实例研究

采用本文提出的KNW估计对1998—2017年间北京市GDP和第一产业、第二产业、第三产业的数据（文献[23]中的例6.2.2）进行分析研究。

首先，基于数据建立回归模型，yi=β0+β1xi1+β2xi2+εi，i=1，2，…，n。

其次，利用Breusch-Pagan检验法得到的p值小于2e-5，比а=0.05小，证实了该模型中异方差的存在性。用前15年的数据样本作为训练集，后5年的数据做测试集。用训练集对参数进行估计，进而计算出对应的拟合值，再用估计出的模型计算出测试集中的预测值。表7是不同估计预测均方根误差（预测RMSE）和拟合均方根误差（拟合RMSE）的具体数值。其中：

当yˆi取因变量的预测值时称为预测 RMSE，当yˆi取因变量的拟合值时称为拟合RMSE。OR和OR1方法分别借助正交表L9（34）和L4（23）对观测值进行展开，并且采用改进后的容差，KNW估计中采用高斯核函数，并固定窗宽h为。数据结果显示，在两种评价方法下，KNW估计对应的数值都要比其他方法小，即 KNW 估计比 HC0、HC4、HC5、OR、OR1估计更有效。

表6 1998—2017年北京市GDP值以及一二三产业增加值 （单位：亿元）

表7 不同估计方法的预测RMSE和拟合RMSE

五、结论

在异方差的构成形式未知时，可以利用非参数方法对方差做出更加精确的估计，进而估计未知参数，最后对模型做出预测和拟合。本文在异方差估计中引入了N-W估计方法，并通过模拟实验和实例分析将KNW方法与已有的异方差估计方法进行比较。在模拟实验中，通过样本容量及核函数的变化说明了该方法的优越性，同时发现，不同的核函数在方差的核估计及模型的预测和拟合方面，效果差别不大。实例分析结果表明，本文所提出的KNW方法较之原有估计方法对于随机误差项协方差阵的估计更加准确，并且相比OR估计运行时间更短。

由于随机误差项的N-W估计是由核函数K（·）和窗宽h共同决定的，而本文只考虑了固定窗宽后不同核函数对核估计效果的影响，因此后期还可以从选择最优窗宽方面对KNW方法进行深入的探索性研究。