近年来, 由于计算机技术的日益成熟, 分位数回归在理论和方法上都得到了广泛的应用. Koenker ^[1]首次提出了分位数回归, 如今分位数回归作为均值回归分析的稳健替代, 被广泛地用于探索响应变量与协变量之间的潜在关系.在实际应用中, 分位数回归可以刻画响应变量更多的分布特征. Koenker ^[2]发现分位数回归的结果可以提供比普通条件均值回归更丰富, 更有针对性.特别是, 它提供了探索异质性的来源与合作的响应变量一种方法, 并深入研究了分位回归模型及其估计.王新宇^[3]系统地介绍了分位数的基本模型及其扩展、分位数回归模型的经典统计推断. Tang等^[4]研究了加权复合分位数(WCQ)与随机截尾线性回归模型.在这个模型中, 提出了可变选择的自适应惩罚程序, 并证明了一致性和渐近正态性. Wang和Yin ^[5]研究了无界意义下的在线变化分位数回归算法.

分位数回归模型中的变量选择问题一直受到广泛的关注. Shows等^[6]针对一种多元线性模型, 提出了对随机删失数据的自适应Lasso加权LAD (AWLAD)变量选择方法. Wang等^[7]提出了BIC调整参数选择方法, 证明了这种方法能够辨别出真模型, 并在模拟中验证了理论的有效性. Wu等^[8]研究了惩罚分位数回归, 在一些较弱的条件下得到了SCAD和自适应Lasso惩罚分位数回归的Oracle性质. Zou ^[9]提出了分位数回归模型的自适应Lasso的变量选择方法, 也得到了其Oracle性质.吕亚召等^[10]研究部分线性单指标复合分位数回归模型, 提出了用自适应Lasso的变量选择方法, 该方法用BIC选择最优调整参数, 在随机模拟中验证了所提方法的优良性.

相对于横截面或是时间序列数据来说, 面板数据含有更多的信息, 因此, 面板数据回归模型的研究越来越受关注.巴尔塔基^[11]提出了面板数据模型及其参数的估计方法, 并给出了实际应用.李扬等^[12]提出了惩罚似然变量选择问题, 证明了面板数据的自适Lasso具有Oracle性质.在选择最优调整参数时, 模拟显示BIC和GCV的选择结果一般比AIC有优势.曲婷等^[13]对平衡纵向数据模型, 通过Lasso方法可将模型的系数压缩到0, 采用AIC和BIC准则选取最优参数, 从而达到变量选择的目的. Koenker^[14]首次提出了面板数据分位数回归模型, 用加权的形式控制分位数对效应的影响, 并加入$l_{1}$惩罚项, 既保持了线性规划形式, 又保持了结果设计矩阵的稀疏性.李翰芳等^[15]对随机效应面板数据, 通过引入条件Laplace先验, 构造了一种新的贝叶斯Lasso分位数回归法, 与一般贝叶斯分位回归法相比更有效的将异质变量的系数压缩到0, 从而起到变量选择的作用.

分位数回归对误差项的分布没有具体的限制, 对异质点或者是非正态分布的参数的估计具有一定的稳健性, 将分位数回归和面板数据模型两者结合起来, 在控制个体差异的同时, 可以分析各种变量在不同分位点之间的关系.基于面板数据的分位数回归模型, 本文提出了一种在改进的自适应Lasso的罚函数下对变量进行选择的方法, 对系数变量的值进行压缩, 使得异质变量的系数为0, 从而达到变量选择的效果, 并证明了相合性和渐近正态性, 在模拟中用验证了选择的有效性.

考虑一般的随机效应面板数据模型

$ y_{ij}=x_{ij}^{T}\beta+\alpha_{i}+u_{ij}, $

(2.1)

其中$ y_{ij}$是因变量, $x_{ij}$是自变量, $\alpha_{i}$是不可观测的时间不变效应, $u_{ij}$是误差项.写成矩阵的形式如下$ y=X^{T}\beta+Z\alpha+u, $其中$y$是$n\times1$维, $X$是$nm\times p$维, $Z$是$nm\times n$维的虚拟变量的关联矩阵, $\alpha$和$u$是独立的随机向量.

令$\rho_{\tau_{k}}(u)=u(\tau_{k}-I(u\leq0))$, $y_{ij}$的分位数函数为

$ Q_{y_{ij}}(\tau_{k}|x_{ij}, \alpha_{i})=x_{ij}^{T}\beta(\tau_{k})+\alpha_{i}, \qquad j=1, \cdots, m;i=1, \cdots, n. $

(2.2)

为了更好的估计参数, 对(2.1) 式提出加权分位数估计方法,

$ \min\limits_{\alpha, \beta} \sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}\rho_{\tau_{k}}(y_{ij}-x_{ij}^{T}\beta(\tau_{k})-\alpha_{i}), $

(2.3)

最小化(2.3) 是一个凸规划问题, 加权分位数回归估计方法可以凸优化来实现.在分位数函数(2.2) 中, $\alpha$与因变量的条件分位数相对应, 为了更好的估计截面的分位数方程, Koenker ^[14]引入了惩罚项$\lambda_{1}\sum\limits_{i=1}^{n}|\alpha_{i}|$, 代替高斯惩罚项,

$ (\widetilde{\alpha}, \widetilde{\beta})=\text{argmin}_{\alpha, \beta} \sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}\rho_{\tau_{k}}(y_{ij}-x_{ij}^{T}\beta(\tau_{k})-\alpha_{i})+\lambda_{1}\sum\limits_{i=1}^{n}|\alpha_{i}|. $

(2.4)

惩罚项$\lambda_{1}\sum\limits_{i=1}^{n}|\alpha_{i}|$与高斯惩罚项相比, 保持了结果设计矩阵的稀疏性的统计优势和线性规划计算优势.

在Koenker^[14]的定理1中给出了$\beta(\tau_{1})$ ($\tau$的个数为1) 正则化后, $\widetilde{\delta_{1}}$的渐近正态性, 本节用另一方法, 将$\beta(\tau_{k})$中$\tau$的个数推广到$q$个, 证明正则化后$\widetilde{u}$的渐近正态性.为了得到$(\alpha, \beta)$估计的渐近性质, 下面给出正则化条件

(A1) $y_{ij}$, $i=1, \cdots, n$, $j=1, \cdots, m$是独立的, 其条件分布为$F_{ij}$ (在$x_{ij}, \alpha_{i}$条件下), 连续密度函数为$f_{ij}$, $0 < f_{ij} < \infty$处处具有有界的微分$f_{ij}^{'}$.

(A2) 令$\xi_{ij}(\tau_{k})=Q_{y_{ij}}(\tau_{k}|x_{ij}, \alpha_{i})=x_{ij}^\mathrm{T}\beta(\tau_{k})+\alpha_{i}$, 当$mn\rightarrow\infty$, 式子

$ \begin{eqnarray*} &&\frac{1}{mn}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}E[f(\xi_{ij}(\tau_{k})|X)z_{ij}^\mathrm{T}z_{ij}], ~~ \frac{1}{mn}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}E[f(\xi_{ij}(\tau_{k})|X)z_{ij}^\mathrm{T}x_{ij}], \\ &&\frac{1}{mn}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}E_{X}[f(\xi_{ij}(\tau_{k})|X)x_{ij}^\mathrm{T}x_{ij}] \end{eqnarray*} $

均收敛到正定矩阵.

(A3) $ \max\limits_{1\leq i\leq n, 1\leq j\leq m}\parallel x_{ij} \parallel/\sqrt{mn}\rightarrow 0$.

(A4) 随机变量$\alpha_{i}$与$x_{ij}$独立.

考虑下面的目标函数

$ V_{mn}(\delta)=\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}[\rho_{\tau_{k}}(y_{ij}-\xi_{ij}(\tau_{k})-v_{ijk})-\rho_{\tau_{k}}(y_{ij}-\xi_{ij}(\tau_{k}))], $

(2.5)

其中$v_{ijk}=z_{ij}^{T}\delta_{0}/\sqrt{mn}+x_{ij}^{T}\delta_{k}/\sqrt{mn}$, 其中$\tilde{\delta}=(\tilde{\delta_{0}}, \tilde{u})^\mathrm{T}=(\sqrt{mn}(\tilde{\alpha}-\alpha), \tilde{\delta_{1}}, \cdots, \tilde{\delta}_{K})^\mathrm{T}$, $\tilde{\delta}_{k}=\sqrt{mn}(\tilde{\beta}(\tau_{k})-\beta(\tau_{k}))^\mathrm{T}$, $k=1, \cdots, K$, $\tilde{\delta}$是(2.5) 式的最小值点.则有如下定理.

定理3.1  假设在$(A_{1}-A_{4})$成立的条件下, $\widetilde{u}\rightarrow_{d}N(0, D^{-1}\Sigma D^{-1})$, 其中$D_{_{0}}$是B的协方差矩阵, $B=(W_{1, }, Z)^{T}$,

$ D^{-1}\Sigma D^{-1}= \begin{pmatrix} (\widetilde{D}^{-1}D_{0}\widetilde{D}^{-1})_{22}&0&\ldots&0\\ 0&(\widetilde{D}^{-1}D_{0}\widetilde{D}^{-1})_{33}&\ldots&0\\ \vdots &\vdots& &\vdots\\ 0&0&\ldots&(\widetilde{D}^{-1}D_{0}\widetilde{D}^{-1})_{KK}\\ \end{pmatrix}. $

证  由Knight ^[16]的结果, 对任意$x\neq0$, 有

$ \rho_{\tau}(x-y)=y[I(x < 0)-\tau]+\int_{0}^{y}[I(x\leq t)-I(x\leq0)]dt. $

将$V_{mn}(\delta)$分解成下面的形式

$ \begin{align*} V_{mn}(\delta)=&\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}[v_{ijk}(I(y_{ij}-\xi_{ij}(\tau_{k}) < 0)-\tau_{k})]\\ &+\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}\int_{0}^{v_{ijk}}[I(y_{ij}-\xi_{ij}(\tau_{k})\leq t)-I(y_{ij}-\xi_{ij}(\tau_{k})\leq0)]dt\\ =&\sum\limits_{k=1}^{K}Z_{n, m, k}\delta_{k}+W_{n, m, k}^{T}\delta_{0}+\sum\limits_{k=1}^{K}B_{n, m, k}, \end{align*} $

其中

$ \begin{eqnarray*} Z_{n, m, k} &=& \frac{1}{\sqrt{mn}}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}x_{ij}^{T}[I(y_{ij}-\xi_{ij}(\tau_{k}) < 0)-\tau_{k})], \\ W_{n, m, k} &=& \frac{1}{\sqrt{mn}}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}z_{ij}^{T}\sum\limits_{k=1}^{K}\omega_{k}[I(y_{ij}-\xi_{ij}(\tau_{k}) < 0)-\tau_{k}], \\ B_{n, m, k}&=& \sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}\int_{0}^{v_{ijk}}[I(y_{ij}-\xi_{ij}(\tau_{k})\leq t)-I(y_{ij}-\xi_{ij}(\tau_{k})\leq0)]dt. \end{eqnarray*} $

由于$E[I(y_{ij}-\xi_{ij}(\tau_{k})<0)-\tau_{k}]=0$, 结合中心极限定理和Cramér-Word定理, $Z_{n, m, k}$和$W_{n, m, k}$依分布收敛到$Z_{k}$和$W_{1}$, 其中$Z_{k}$是一个正态随机变量, 均值为0, $W_{1}$是一个$n$维正态向量, 均值为0.因此可以得到

$ \sum\limits_{k=1}^{K}Z_{n, m, k}\delta_{k}+W_{n, m, k}^{T}\delta_{0}\rightarrow_{d}\sum\limits_{k=1}^{K}Z_{k}\delta_{k}+W_{1}^{T}\delta_{0}. $

令

$ B_{n_{i}, m_{j}, k}=\omega_{k}\int_{0}^{v_{ijk}}[I(y_{ij}-\xi_{ij}(\tau_{k})\leq t)-I(y_{ij}-\xi_{ij}(\tau_{k})\leq0)]dt, $

则对任意$\epsilon>0$, 有$[B_{n_{i}, m_{j}, k}]^{2}=[B_{n_{i}, m_{j}, k}]^{2}I(v_{ijk}\geq\epsilon)+[B_{n_{i}, m_{j}, k}]^{2}I(v_{ijk} < \epsilon). $

一方面,

$ \begin{align*} mnE{[B_{n_{i}, m_{j}, k}]^{2}I(v_{ijk}\geq\epsilon)} &\leq mnE{\omega_{k}}^{2}(\int_{0}^{v_{ijk}}2dt)^{2} I(v_{ijk}\geq\epsilon)\\ &\leq 4mn\omega_{k}^{2}E[|v_{ijk}|^{2}I(v_{ijk}\geq\epsilon)]\\ &\rightarrow0, ~~ \text{当}~~n, m\rightarrow0. \end{align*} $

另一方面,

$ \begin{align*} &mn E{[B_{n_{i}, m_{j}, k}]^{2}I(v_{ijk} < \epsilon)}\\ \leq& mn E (\omega_{k}^{2} \int_{0}^{v_{ijk}}2dt \int_{0}^{v_{ijk}}[I(y_{ij}-\xi_{ij}(\tau_{k})\leq t)-I(y_{ij}-\xi_{ij}(\tau_{k})\leq0)]dt I(v_{ijk} < \epsilon))\\ \leq& 2mn \omega_{k}^{2} \epsilon E_{X}(\int_{0}^{v_{ijk}}[F(\xi_{ij}(\tau_{k})+t|X)-F(\xi_{ij}(\tau_{k})|X)]dt)\\ \leq&2mn \omega_{k}^{2} \epsilon E_{X}(\int_{0}^{v_{ijk}}f(\xi_{ij}(\tau_{k})|X)tdt)\\ \leq&2mn \omega_{k}^{2} \epsilon f(\xi_{ij}(\tau_{k})|X)E_{X}|z_{ij}\delta_{0}+x_{ij}^{T}\delta_{k}^{2}|\\ \rightarrow&0, ~~ \text{当}~~\epsilon\rightarrow0. \end{align*} $

因此当$ mn\rightarrow\infty $时,

$ \begin{align*} {\rm Var}(\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}B_{n_{i}, m_{j}, k}) =\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}{\rm Var}(B_{n_{i}, m_{j}, k})\leq mn E(B_{n_{i}, m_{j}, k})^{2} \rightarrow0. \end{align*} $

则有

$ \sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}[B_{n_{i}, m_{j}, k}-E(B_{n_{i}, m_{j}, k})]=o_{p}(1). $

另外, 由于

$ \begin{align*} E(B_{n, m, k})&=\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}E(B_{n_{i}, m_{j}, k})\\ &=m n E(\omega_{k}\int_{0}^{v_{ijk}}[I(y_{ij}-\xi_{ij}(\tau_{k})\leq t)-I(y_{ij}-\xi_{ij}(\tau_{k})\leq0)]dt)\\ &=mn \omega_{k} E_{X}(\int_{0}^{v_{ijk}}[F(\xi_{ij}(\tau_{k})+t|X)-F(\xi_{ij}(\tau_{k})|X)]dt)\\ &=mn \omega_{k} E_{X}(\int_{0}^{v_{ijk}}f(\xi_{ij}(\tau_{k})|X)tdt)+o_{p}(1)\\ &=\frac{1}{2}\delta_{0}^{T}G_{k}\delta_{0}+\frac{1}{2}\delta_{k}^{T}D_{k}\delta_{k}+\delta_{0}H_{k}\delta_{k}+o_{p}(1), \end{align*} $

其中

$ \begin{eqnarray*} % \nonumber to remove numbering (before each equation) G_{k} &=& \omega_{k}E_{X}[f(\xi_{ij}(\tau_{k})|X)z_{ij}^\mathrm{T}z_{ij}], ~~ H_{k} = \omega_{k}E_{X}[f(\xi_{ij}(\tau_{k})|X)z_{ij}^\mathrm{T}x_{ij}], \\ D_{k}&=& \omega_{k}E_{X}[f(\xi_{ij}(\tau_{k})|X)x_{ij}^\mathrm{T}x_{ij}]. \end{eqnarray*} $

则

$ \begin{align*} V_{mn}(\delta)&= \frac{1}{2}\delta_{0}^\mathrm{T}G\delta_{0}+\frac{1}{2}\sum\limits_{k=1}^{K}\delta_{k}^\mathrm{T}D\delta_{k}+\delta_{0}^\mathrm{T}\sum\limits_{k=1}^{K}H_{k}\delta_{k}+W_{1}^\mathrm{T}\delta_{0} +\sum\limits_{k=1}^{K}Z_{k}\delta_{k}+o_{p}(1)\\ &= \frac{1}{2}\delta_{0}^\mathrm{T}G\delta_{0}+\frac{1}{2}u^\mathrm{T}D^{'} u+\delta_{0}^{T}H^{T}u+W_{1}^\mathrm{T}\delta_{0}+Z^\mathrm{T}u+o_{p}(1)\\ &=\frac{1}{2}\delta^\mathrm{T}\widetilde{D}\delta+\delta^\mathrm{T}B+o_{p}(1), \end{align*} $

其中$H=(H_{1}, H_{2}\cdots H_{K})^\mathrm{T}$, $Z=(z_{1}, z_{2}\cdots z_{K})^\mathrm{T}$, $G=\sum\limits_{k=1}^{K}G_{k}$, $B=(W_{1, }, Z)^{T}$,

$ D^{'}= \begin{pmatrix} D_{1}&0&\ldots&0\\ 0&D_{2}&\ldots&0\\ \vdots &\vdots& &\vdots\\ 0&0&\ldots&D_{K}\\ \end{pmatrix}, \widetilde{D}= \left( {\begin{array}{*{20}{c}} {G}&{H^\mathrm{T}} \\ {H}&{D'} \\ \end{array}} \right). $

即$ V_{mn}(\delta)\rightarrow\frac{1}{2}\delta^\mathrm{T}\widetilde{D}\delta+\delta^\mathrm{T}B. $

由Koenker ^[14]中引理1, 可以得到$\widetilde{u}\rightarrow_{d}N(0, D^{-1}\Sigma D^{-1}). $

在对数据进行统计分析时, 人们一般会借助一些相关变量对所关心的变量进行分析, 建模, 以便得到理想的结果, 一般称这些相关的变量为协变量, 而所关心的变量为因变量.在开始建模的时候, 希望加入更多的相关变量, 来得到更真实的结果, 然而, 随着协变量的增多, 异质变量存在的可能性就越大, 于是, 希望寻找一个有效方法来选出对响应变量有显著影响的协变量.因此变量选择就是统计学中一个重要的问题.本节对上述面板数据分位数模型的变量选择进行分析, 在(4.1) 式中需要指定调节参数$\lambda_{2}$, 本文最优的调整参数$\lambda_{2}$可以通过BIC (Bayesian information criterion)准则选取.在加权分位数估计的同时, 同时希望对变量做选择, 本节选的罚函数是自适应Lasso罚函数.令

$ \begin{eqnarray} (\hat{\alpha}, \hat{\beta})&=& \text{argmin}_{\alpha, \beta}\sum\limits_{k=1}^{K} \sum\limits_{j=1}^{m} \sum\limits_{i=1}^{n} \omega_{k}[\rho_{\tau_{k}}(y_{ij}-x_{ij}^{T}\beta(\tau_{k})-\alpha_{i})]\nonumber\\ &&+\lambda_{1}\sum\limits_{i=1}^{n}|\alpha_{i}|+\lambda_{2}\sum\limits_{k=1}^{K}\omega_{k}\sum\limits_{j=1}^{p}\frac{|\beta_{j}(\tau_{k})|}{|\widetilde{\beta}_{j}(\tau_{k})|}, \end{eqnarray} $

(4.1)

令${\rm BIC}(\lambda)=\log P_{\lambda}+df_{\lambda}\cdot\log(mn)/mn, $其中

$ P_{\lambda}=\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}\rho_{\tau_{k}}(y_{ij}-x_{ij}^{T}\beta(\tau_{k})-\alpha_{i}). $

$df_{\lambda}$为$\widehat{\beta}$非零分量的个数.这样可以得到调节参数的一个适合的估计$\widehat{\lambda}=\arg\min_{\lambda}{{\rm BIC}(\lambda)}$.令$\beta^{\ast}=(\beta_{1}^{\ast T}, \beta_{2}^{\ast T})^\mathrm{T}$和$\widehat{\beta}=(\widehat{\beta}_{1}^\mathrm{T}, \widehat{\beta}_{2}^\mathrm{T})^\mathrm{T}$分别是参数向量$\beta$的真实值和加权分位数自适应Lasso惩罚估计值.不失一般性, 可以假设$\beta_{1}^{\ast}$包含$\beta^{\ast}$中的前$q$个非零部分, 同时$\beta_{2}^{\ast }$是剩下$p-q$维零向量, 则有下面的定理.

定理4.1  假设在(A1)-(A5) 成立的条件下, 如果$\sqrt{m n}\lambda_{2}\rightarrow0$, $m n\lambda_{2}\rightarrow\infty$, 则当$mn\rightarrow\infty$, $\widehat{\beta}$满足

(ⅰ)选择相合性: $P(\widehat{\beta}_{2}=0)\rightarrow0$;

(ⅱ)渐近正态性: $\sqrt{mn}(\widehat{\beta}_{1}-\beta_{1}^{\ast})\rightarrow_{d}N(0, (D_{1})^{-1}\Sigma(D_{1})^{-1})$.

证  令$\widehat{\delta}=(\widehat{\delta_{0}}, \widehat{u})^\mathrm{T}=(\sqrt{mn}(\widehat{\alpha}-\alpha), \widehat{\delta_{1}}, \cdots, \widehat{\delta}_{K})^\mathrm{T}$, $\widehat{\delta}_{k}=\sqrt{mn}(\widehat{\beta}(\tau_{k})-\beta(\tau_{k}))^\mathrm{T}$, $k=1, \cdots, K$, $\widehat{\delta}$是下面式子的最小值点.

$ \begin{eqnarray*} L_{mn}(\delta)&=& \sum\limits_{k=1}^{K} \sum\limits_{j=1}^{m} \sum\limits_{i=1}^{n} \omega_{k}[\{\rho_{\tau_{k}}(y _{ij}-\xi_{ij}(\tau_{k})-v_{ij})-\rho_{\tau_{k}}(y_{ij}-\xi_(\tau_{k}))]\\ &&+\lambda_{1}\sum\limits_{i=1}^{n}|\alpha_{i}-\delta_{0i}/\sqrt{mn}|-|\alpha_{i}|\\ &&+\lambda_{2}/|\widetilde{\beta}_{j}(\tau_{k})|\sum\limits_{k=1}^{K}\omega_{k}\sum\limits_{j=1}^{p}|\beta_{j}(\tau_{k})-\delta_{kj}/\sqrt{mn}|-|\beta_{j}(\tau_{k})|.\end{eqnarray*} $

(ⅰ)因为$L_{mn}(\delta)$是对$\delta$的分段线性函数, 在每个可微的点, 对$k=1, 2, \cdots, K$, $j=q+1, \cdots, p$取$L_{mn}(\delta)$对$\delta_{kj}$的偏导, 有

$ \begin{align*} &\frac{1}{\sqrt{mn}}\frac{\partial L_{mn}}{\partial\delta_{kj}}\\ =&-\frac{1}{\sqrt{mn}}\sum\limits_{k=1}^{K} \sum\limits_{j=1}^{m} \sum\limits_{i=1}^{n} \omega_{k}[I(y _{ij}-\xi_{ij}(\tau_{k})-v_{ijk}\leq0)-\tau_{k}]+\lambda_{2}mn\omega_{k}\frac{{\rm sign}(\delta_{kj})}{|\widetilde{\beta}(\tau_{k})|}+o_{p}(1)\\ =&\sqrt{mn}\omega_{k}\frac{{\rm sign}(\delta_{kj})}{|\widetilde{\beta}(\tau_{k})|}-\sum\limits_{k=1}^{K}(\Psi_{mn}(v_{ijk}))+o_{p}(1), \end{align*} $

其中

$ \Psi_{mn}(v_{ijk})=\sum\limits_{j=1}^{m} \sum\limits_{i=1}^{n} \omega_{k}[I(y _{ij}-\xi_{ij}(\tau_{k})-v_{ijk}\leq0)-\tau_{k}]. $

由定理3.1, 可以得到$\Psi_{mn}(0)$收敛到均值为0, 方差为$\Xi $的正态分布, $\max\limits_{1 < i < n, 1 < j < m}v_{ijk}\rightarrow0$.类似Tang ^[4]的证明可以得到$\sup_{\|\delta(k)\|}|\Psi_{mn}(v_{ijk})-\Psi_{mn}(0)+\Xi v_{ijk}|=o_{p}(1)$, 注意到$\Psi_{mn}(0)$依分布收敛到正态随机向量, $\Xi$是有限的, 当$n\rightarrow\infty$时$v_{ijk}\rightarrow0$, 因此对$k=1, 2\cdots, K$; $j=q+1, \cdots, p$, 有

$ \frac{1}{\sqrt{mn}}\frac{\partial L_{mn}(\delta)}{\partial\delta_{kj}} =\frac{\lambda_{2}\omega_{k}\sqrt{mn}{\rm sign}(\delta_{kj})}{|\widetilde{\beta}_{j}(\tau_{k})|} =\frac{\lambda_{2}\omega_{k}mn {\rm sign}(\beta_{j}(\tau_{k}))}{|\widetilde{\beta}_{j}(\tau_{k})|}. $

当$n$充分大时, $mn\lambda_{2}\rightarrow\infty$, $\omega_{k}>0$, 则$\partial L_{mn}(\delta)/\partial\delta_{kj}$的符号由$\beta_{j}(\tau_{k})$的符号决定.因此当$n$充分大时, 对于$k=1, 2\cdots, K$, $j=q+1, \cdots, p$, 如果$\beta_{j}(\tau_{k})>0$, $\partial L_{mn}(\delta)/\partial\delta_{kj}$为正; 如果$\beta_{j}(\tau_{k}) < 0$, $\partial L_{mn}(\delta)/\partial\delta_{kj}$为负.可以得到$P(\delta_{kj}=0)\rightarrow1$, $k=1, 2\cdots, K$, $j=q+1, \cdots, p$, 也就是说当$mn\rightarrow\infty$, $P(\widehat{\beta}_{2}=0)\rightarrow1$.

(ⅱ)正如定理3.1和(ⅰ), 只需建立$L_{mn}[\delta_{0}, (\beta_{1}^\mathrm{T}, \bf{0})]$的渐近形式, 因为$\sqrt{mn}\lambda_{2}\rightarrow0$, $\beta_{j}(\tau_{k})\rightarrow\beta_{j}^{\ast}(\tau(k))$, $k=1, 2\cdots, K$, $j=q+1, \cdots, p$, 则$L_{mn}[\delta_{0}, (\beta_{1}^\mathrm{T}, \bf{0})]$可以表示为

$ L_{mn}[\delta_{0}, (\beta_{1}^{\mathrm{T}}, {\bf{0}})]=\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}[\{\rho_{\tau_{k}}(y _{ij}-\xi_{ij}(\tau_{k})-v_{ijk}^{1})-\rho_{\tau_{k}}(y_{ij}-\xi_(\tau_{k}))]+O_{p}(1), $

其中$v_{ijk}^{1}=z_{ij}^{T}\delta_{0}/\sqrt{mn}+x_{ij}^{1T}\delta_{k}/\sqrt{mn}$, $x_{ij}^{1}=(x_{ij1}, x_{ij2}, \cdots, x_{ijq})$, $i=1, \cdots n$, $j=1, \cdots m$, $\widehat{u}^{1}=\sqrt{mn}(\widehat{\beta}_{1}-\beta_{1}^{\ast}).$类似定理3.1的证明, 有

$ \begin{align*} L_{mn}[\delta_{0}, (\beta_{1}^\mathrm{T}, \bf{0})] =&\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}[v_{ijk}^{1}(I(y_{ij}-\xi_{ij}(\tau_{k}) < 0)-\tau_{k})]\\ &+\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}\int_{0}^{v_{ijk}^{1}}[I(y_{ij}-\xi_{ij}(\tau_{k}) \leq t)-I(y_{ij}-\xi_{ij}(\tau_{k})\leq0)]dt+o_{p}(1)\\ =&w_{1}^\mathrm{T}\delta_{0}+\frac{1}{2}\delta_{0}^\mathrm{T}G\delta_{0}+\delta_{0}^\mathrm{T}\sum\limits_{k=1}^{K}H_{k}^{1}\delta_{k}^{1} +\frac{1}{2}\sum\limits_{k=1}^{K}\delta_{k}^\mathrm{1T}D_{k}^{1}\delta_{k}^{1}+\sum\limits_{k=1}^{K}z_{k}^\mathrm{1T}\delta_{k}^{1}+o_{p}(1)\\ =&w_{1}^\mathrm{T}\delta_{0}+\frac{1}{2}\delta_{0}^\mathrm{T}G\delta_{0}+\delta_{0}^\mathrm{T}H^{1}u^{1}+\frac{1}{2}u^\mathrm{T}D_{1}u^{1}+o_{p}(1), \end{align*} $

其中$1/\sqrt{mn}\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}\omega_{k}[x_{ij}^\mathrm{1T}(I(y_{ij}-\xi_{ij}(\tau_{k}) < 0)-\tau_{k})]\rightarrow_{d}z_{k}^{1}$, $D_{k}^{1}=E_{X}(f(\xi_{ij}(\tau_{k})|X)x_{ij}^\mathrm{1T}x_{ij}^{1})$, $\delta_{k}^{1}=(\sqrt{mn}(\widehat{\beta}_{1}(\tau_{k})-\beta_{1}(\tau_{k})), \cdots, \sqrt{mn}(\widehat{\beta}_{p}(\tau_{k})-\beta_{p}(\tau_{k})))^\mathrm{T}$, $z_{k}^{1}$是均值为$0$的$q$维正态随机向量.则

$ L_{mn}[\delta_{0}, (\beta_{1}^\mathrm{T}, {\bf{0}})]\rightarrow_{d}\frac{1}{2}\delta^\mathrm{1T}\widetilde{D}_{1}\delta^{1}+\delta^\mathrm{1T}B^{1}. $

同时可以得到

$ \widehat{u}^{1}\rightarrow_{d}N(0, (D_{1})^{-1}\Sigma(D_{1})^{-1}). $

在本节给出两个例子, 比较不同的方法对参数估计值优势, 并验证自适应Lasso罚函数对变量选择的有效性.

例1  考虑$n=50, m=5, p=1$, 响应变量由下面的模型生成

$ y_{ij}=\alpha_{i}+x_{ij}\beta+u_{ij}, $

其中$\beta=1$, $\alpha_{i}$和$u_{ij}$服从标准正态分布, $\omega=(0. 25, 0. 5, 0. 25)$在三个分位点$\tau=(0. 25, 0. 5, 0. 75)$, $x_{ij}$由高斯分布生成

$ x_{ij}=\gamma_{i}+v_{ij}, $

(5.1)

$\gamma_{i}$和$v_{ij}$独立同分布, 相应的组内相关系数,

$ \rho_{x}=\sigma_{\gamma}^{2}/(\sigma_{\gamma}^{2}+\sigma_{u}^{2}), $

(5.2)

就是$x_{ij}$和$x_{ik}$之间的相关系数, 当$j\neq k$时, 在的模拟中, 都令$\rho_{x}=0. 5$.而$\lambda_{1}$选择位置参数比$\sigma_{u}/\sigma_{\alpha}$, $\lambda_{2}$的选择由上一节BIC得到, $\alpha$和$u_{ij}$分两种情况.

1.都来自于标准正态;

2.都来自于自由度为3的$t$分布.

这样可以得到分别在分位数回归的估计方法(QR)、分位数效应罚函数估计(PQR)、分位数回归自适应罚函数估计(LPQR), 对$\beta$的估计, 如表 1, 可以看出在$\alpha$和$u_{ij}$的两种情况PQR和LPQR都比QR估计更优.

例2  令$m=5, n=50, p=8$, 响应变量来自下面的模型

$ y_{ij}=\alpha_{i}+x_{ij}^\mathrm{T}\beta+u_{ij}, $

(5.3)

$\beta=(3, 1.5, 0, 0, 0, 0, 2, 0)$, $x_{ij}$由(5.1), (5.2) 式生成, $\alpha_{i}$和$u_{ij}$同样分两种情况.

1.都来自于标准正态;

2.都来自于自由度为3的$t$分布.

表 2是分位数罚估计(PQR)分别对上面两种情形下$\beta$的估计, 表 3是分位数自适应Lasso罚函数(LPQR)对参数的估计, 通过模拟可以看出PQR可以对参数做近似估计, 但对异质变量不能做选择, 而LPQR在参数估计的同时对变量做了选择, 0参数都选择出来了, 不管是参数估计还是变量选择都比PQR有优势.