随机限制回归模型的自适应最小绝对收缩和选择算子方法

Adaptive Least Absolute Shrinkage and Selection Operator Method for Stochastic Limit Regression Model

摘要

随机限制最小绝对收缩和选择算子（M-Lasso）方法可以在变量选择的同时使用随机的先验信息，但是该方法基于的最小绝对收缩和选择算子（Lasso）对每一个系数的惩罚是等权重的，这可能会导致某些重要的信息被过度压缩。为此，本文提出随机限制自适应Lasso（Ma-Lasso）方法。该方法赋予系数不同的权重，并且在变量选择的同时使用了随机先验信息，可提高估计的精度。通过数值实验结果分析发现，该方法在稀疏模型上表现出比其他方法更小的均方误差，并且在发现率、真实发现率以及真实模型选择次数比例方面也具有一定的优势。最后，通过将该方法应用于贵州茅台各季度财报数据和股票价格数据中，发现Ma-Lasso所构建的模型的贝叶斯信息准则（Bayesian Information Criterion，BIC）值相较于M-Lasso方法下降了约5%，进一步验证了它的优越性。

Abstract

The Mixed least absolute shrinkage and selection operator (M-Lasso) method can use random prior information at the same time as variable selection, but the minimum absolute shrinkage and selection operator (Lasso) based on this method is equally weighted for each coefficient, which may cause some important information to be overcompressed. To address this issue, this paper proposes the stochastic constrained adaptive Lasso (Ma-Lasso) method. The method assigns different weights to the coefficients and has oracle properties. It uses stochastic prior information along with variable selection, which can improve the precision of the estimation. The analysis of numerical experimental results reveals that the method exhibits a smaller mean square error on the sparse model than the other methods, and also has some advantages in terms of the discovery rate, the true discovery rate, and the ratio of the number of times the true model is selected. Finally, by applying Ma-Lasso to the quarterly financial data and stock price data of Kweichow Moutai, it is found that the BIC value of the model constructed by this method decreases by about 5% compared with M-Lasso method, which further verifies its superiority.

关键词

稀疏模型 / 变量选择 / 均方误差 / 随机先验信息

Key words

sparse model / variable selection / mean squared error / random prior information

0 引言

最小二乘法（Ordinary Least Square，OLS）是最常用的线性回归模型系数的估计方法。它通过最小化结果变量的预测值与观测值之间的误差来估计回归模型中的参数。该方法可以对当前数据集进行无偏估计，但容易导致过拟合现象。当误差项为异方差或者相关时，研究人员可以用广义最小二乘法（Generalized Least Squares，GLS）来估计参数；而如果模型参数存在线性限制，则可以使用受限最小二乘法（Restricted Least Squares，RLS）^［1］。

但是当数据间存在多重共线性时，最小二乘法不再是一个有效估计^［2-5］。为此，Hoerl和Kennard提出岭回归方法^［6］，这是能够进行共线性数据分析的一种有偏估计回归方法。它可通过缩小回归系数来减少预测误差，以缓解过度拟合，但它不能精确地将系数压缩为0，不进行协变量选择，除非调节参数\lambda趋向于\mathrm{\infty }。Tibshirani^［7］提出的最小绝对收缩和选择算子（Least Absolute Shrinkage and Selection Operator，Lasso）方法正好克服了这一点，当调节参数\lambda足够大时就会迫使一些系数估计为0，从而实现变量选择和系数估计^［8-9］。Lasso在生物信息学和经济学方面应用广泛^［10-12］，但是对于特征数大于样本个数的数据集，Lasso方法最多可以选择n个变量。Zou等^［13］将岭回归与Lasso的惩罚函数结合提出带有L1范数与L2范数惩罚项的弹性网 （Elastic net）方法，该方法可以用来解决变量个数大于样本个数的情况。Zou^［14］还在Lasso的基础上提出了自适应Lasso方法（Adaptive Lasso，A-Lasso），该方法在L1范数惩罚下引入数据自适应的权重，比Lasso算法具有更准确的变量选择和模型预测能力。

针对回归模型中随机先验限制问题，Theil^［15］提出了混合估计方法（Mixed Estimator，ME）。该方法将OLS和GLS相结合，同时使用样本的信息和额外的先验信息（如方差比等）进行估计。通过这种方式，ME方法可以提高估计的准确性和稳定性，特别是在样本容量较小的情况下，更能体现出优势。但当数据特征大于样本量的时候，传统的解决随机限制的ME模型并不能得到很好的效果。为了解决这一问题，Guler等^［16］将ME与Lasso进行结合，提出了混合Lasso（Mixed-Lasso，M-Lasso）估计。该方法很好地解决了随机限制在数据特征大于样本量的数据集下ME不能准确估计的问题。

但由于M-Lasso基于的Lasso方法对每一个系数给与的惩罚权重是等同的，这就不可避免的会导致对于重要变量的惩罚过重，对于不重要变量惩罚过轻，因而这可能就会导致模型的性能下降。为此，本文将M-Lasso基于的Lasso方法改进为A-Lasso方法。通过引入A-Lasso方法，可以更准确地为不同系数分配不同的惩罚权重，从而平衡对重要变量和不重要变量的惩罚程度，并提升模型的性能和解释能力。另外，本文还通过仿真实验验证了在稀疏模型上，与其他方法相比，该方法具有更小的均方误差。最后本文还通过实际应用，进一步验证了该方法在稀疏模型上的有效性。

1 模型及Ma-Lasso方法

考虑如下广义线性模型：

{\mathit{y}}^{\mathit{*}}={\mathit{X}}^{\mathit{*}}\mathit{\beta }+{\mathit{\epsilon }}^{\mathit{*}}，（1）

其中{\mathit{y}}^{\mathit{*}}={\mathit{V}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{y}，{\mathit{X}}^{\mathit{*}}={\mathit{V}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{X}，{\mathit{\epsilon }}^{\mathrm{*}}={\mathit{V}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{\epsilon }。\mathit{y}是n\times \mathrm{1}维的可观测的向量，\mathit{X}是n\times p维列满秩矩阵，\mathit{\beta }是p\times \mathrm{1}维未知参数向量，\mathit{\epsilon }是n\times \mathrm{1}维误差向量，其中\mathit{\epsilon }服从均值为0，协方差为{\sigma }^{\mathrm{2}}\mathit{V}的正态分布，\mathit{V}是给定的正定矩阵。将OLS应用于方程（1）中，得到广义最小二乘估计^［17］，具体公式如下所示：

{\widehat{\mathit{\beta }}}^{\mathrm{G}\mathrm{L}\mathrm{S}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\beta }{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{R}\mathrm{S}\mathrm{S}={\left({\mathit{X}}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}\mathit{X}\right)}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}\mathit{y}={\mathit{S}}^{-\mathrm{1}}{{\mathit{X}}^{\mathrm{*}}}^{\mathrm{T}}{\mathit{y}}^{\mathrm{*}}，（2）

其中残差平方和\mathrm{R}\mathrm{S}\mathrm{S}=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^n}(y_i^{\mathrm{*}}{-{\displaystyle \sum _{j=\mathrm{1}}^p}{x}_{ij}^{\mathrm{*}}{\beta }_j)}^{\mathrm{2}}，并且\mathit{S}=\left({\mathit{X}}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}\mathit{X}\right)。

若存在对\mathit{\beta }的先验线性限制如下所示：

\mathit{r}=\mathit{R}\mathit{\beta }，（3）

其中\mathit{R}是一个m\times p维先验矩阵，\mathit{r}是已知的m\times \mathrm{1}维常数向量，并且\mathit{R}的秩\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathit{R}\right)=m\le p。在公式（3）中给出的限制下，得到RLS模型估计如下：

{\widehat{\mathit{\beta }}}^{\mathrm{R}\mathrm{L}\mathrm{S}}={\widehat{\mathit{\beta }}}^{\mathrm{G}\mathrm{L}\mathrm{S}}-{\mathit{S}}^{-\mathrm{1}}{\mathit{R}}^{\mathrm{T}}{\left\{\mathit{R}{\mathit{S}}^{-\mathrm{1}}{\mathit{R}}^{\mathrm{T}}\right\}}^{-\mathrm{1}}\left(\mathit{R}{\widehat{\mathit{\beta }}}^{\mathrm{G}\mathrm{L}\mathrm{S}}-\mathit{r}\right)。（4）

然而得到的先验信息可能具有一些随机性，这些随机性来自先前的统计分析或不确定理论所获得的先验估计^［18］。先验信息中的不确定性模型如下：

\mathit{r}=\mathit{R}\mathit{\beta }+\mathit{\varphi }，（5）

其中\mathit{\varphi }是m\times \mathrm{1}维误差向量，其均值为0，协方差为{\mathrm{\sigma }}^{\mathrm{2}}\mathit{W}，\mathit{W}是已知的正定矩阵。

令{\mathit{y}}_{\mathit{A}}={\left(\begin{array}{cc}{\mathit{y}}^{\mathrm{T}}& {\mathit{r}}^{\mathrm{T}}\end{array}\right)}^{\mathrm{T}}，{\mathit{X}}_A={\left(\begin{array}{cc}{\mathit{X}}^{\mathrm{T}}& {\mathit{R}}^{\mathrm{T}}\end{array}\right)}^{\mathrm{T}}，\mathit{u}={\left(\begin{array}{cc}{\mathit{\epsilon }}^{\mathrm{T}}& {\mathit{\varphi }}^{\mathrm{T}}\end{array}\right)}^{\mathrm{T}}，那么线性模型就可写成如下公式：

{\mathit{y}}_A={\mathit{X}}_A\mathit{\beta }+\mathit{u}，（6）

令\mathit{T}=\left(\begin{array}{cc}\mathit{V}& \mathrm{0}\\ {\mathrm{0}}^{\mathrm{T}}& \mathit{W}\end{array}\right)是已知的正定矩阵，将GLS应用到上式中，得到ME估计，估计式如下^［18］：

{\widehat{\mathit{\beta }}}^{\mathrm{M}\mathrm{E}}={\left({\mathit{X}}_A^{\mathrm{T}}{\mathit{T}}^{-\mathrm{1}}{\mathit{X}}_A\right)}^{-\mathrm{1}}{\mathit{X}}_A^{\mathrm{T}}{\mathit{T}}^{-\mathrm{1}}{\mathit{y}}_A={\left({\mathit{X}}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}\mathit{X}+{\mathit{R}}^{\mathrm{T}}{\mathit{W}}^{-\mathrm{1}}\mathit{R}\right)}^{-\mathrm{1}}\left({\mathit{X}}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}\mathit{y}+{\mathit{R}}^{\mathrm{T}}{\mathit{W}}^{-\mathrm{1}}\mathit{r}\right)。（7）

然而上述方法并不能用于稀疏模型，Lasso估计可用于估计稀疏模型，具体估计式如下：

{\widehat{\mathit{\beta }}}^{\mathrm{L}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}}\left(\lambda \right)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\beta }{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\mathrm{R}\mathrm{S}\mathrm{S}+\lambda {\displaystyle \sum _{j=\mathrm{1}}^p}\left|{\beta }_j\right|\right\}，（8）

其中\lambda是一个非负的调优参数，用于控制回归系数压缩的程度，数值越大则惩罚力度越强。参数\lambda使用贝叶斯信息准则（Bayesian Information Criterion，BIC）方法进行选择，具体公式如下：

\mathrm{B}\mathrm{I}\mathrm{C}\left(\lambda \right)=\mathrm{l}\mathrm{o}\mathrm{g}\left({\widehat{\sigma }}^{\mathrm{2}}\right)+\widehat{df}\frac{\mathrm{l}\mathrm{o}\mathrm{g}\left(n\right)}{n}，（9）

其中{\widehat{\sigma }}^{\mathrm{2}}=\frac{\mathrm{R}\mathrm{S}\mathrm{S}}{n}，\widehat{df}是非零系数的个数，\lambda ={\mathrm{2}}^{k-\mathrm{15}}，其中k=\mathrm{1,2},\mathrm{3},\dots ,\mathrm{20}。通过给出\lambda计算BIC值，取值最小的BIC值所对应的调优参数的值即为模型的最优\lambda。

Zou^［14］通过在{\displaystyle \sum _{j=\mathrm{1}}^p}\left|{\beta }_j\right|惩罚项前增加权重，对Lasso进行改进，提出自适应Lasso方法（Adaptive-Lasso）^［14］。在该方法中，惩罚项变成了{\displaystyle \sum _{j=\mathrm{1}}^p}\widehat{\mathit{w}}\left|{\beta }_j\right|，其中\widehat{\mathit{w}}=\frac{\mathrm{1}}{{\left|{\beta }_j\right|}^{\gamma }}，\gamma是一个任意的非负数。公式如下：

\widehat{\mathit{\beta }}{}^{\mathrm{A}}{}^{-}{}^{\mathrm{L}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}}\left(\lambda \right)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\beta }{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\mathrm{R}\mathrm{S}\mathrm{S}+\lambda {\displaystyle \sum _{j=\mathrm{1}}^p}\widehat{\mathit{w}}\left|{\beta }_j\right|\right\}。（10）

采用A-Lasso 进行变量选择，能够使得对响应变量重要的变量更易进入模型， 而对其不重要的变量更易被剔除，并且在更好地实现变量选择的同时也能够有效减小系数估计的偏差。因此相比于 Lasso 方法， 自适应 Lasso方法也更适用于观测指标数p和样本量N的比值非常大的情况。

将Lasso模型中加入先验信息，可以提高估计的准确性和稳定性。Guler等提出M-Lasso^［16］，可用于具有随机限制的大数据集。他将ME和Lasso结合，以希望在选择正确模型的同时减少均方误差（Mean Squared Error， MSE）。M-Lasso模型如下所示：

\widehat{\mathit{\beta }}{}^{\mathrm{M}}{}^{-}{}^{\mathrm{L}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}}\left(\lambda \right)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\beta }{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\frac{\mathrm{1}}{\mathrm{2}}{{\mathit{u}}^{\mathrm{*}}}^{\mathrm{T}}{\mathit{u}}^{\mathrm{*}}+\lambda {\displaystyle \sum _{j=\mathrm{1}}^p}\left|{\beta }_j\right|\right\}，（11）

其中{\mathit{u}}^{\mathrm{*}}={\mathit{T}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{u}={\left(\begin{array}{cc}{\mathit{\epsilon }}^{{\mathrm{*}}^{\mathrm{T}}}& {\mathit{\varphi }}^{{\mathrm{*}}^{\mathrm{T}}}\end{array}\right)}^{\mathrm{T}}，{\mathit{T}}^{-\frac{\mathrm{1}}{\mathrm{2}}}=\left(\begin{array}{cc}{\mathit{V}}^{-\frac{\mathrm{1}}{\mathrm{2}}}& \mathrm{0}\\ {\mathrm{0}}^{\mathrm{T}}& {\mathit{W}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\end{array}\right)，\mathit{u}={\left(\begin{array}{cc}{\mathit{\epsilon }}^{\mathrm{T}}& {\mathit{\varphi }}^{\mathrm{T}}\end{array}\right)}^{\mathrm{T}}。

A-Lasso有着比Lasso更为精确的模型选择能力，且具有oracle性质。受M-Lasso的启发，本文将ME与A-Lasso结合，提出混合自适应Lasso估计（Ma-Lasso），希望其在具备随机限制的大数据集上有着比M-Lasso更为优秀的表现。为了得到Ma-Lasso估计，本文对线性模型进行变换：

{{\mathit{y}}_{\mathit{A}}}^{\mathrm{*}}={{\mathit{X}}_{\mathit{A}}}^{\mathrm{*}}\mathit{\beta }+{\mathit{u}}^{\mathrm{*}}，（12）

其中{{\mathit{y}}_{\mathit{A}}}^{\mathrm{*}}={\mathit{T}}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{y}}_{\mathit{A}}={\left(\begin{array}{cc}{\mathit{y}}^{{\mathit{*}}^{\mathrm{T}}}& {\mathit{r}}^{{\mathit{*}}^{\mathrm{T}}}\end{array}\right)}^{\mathrm{T}}，{{\mathit{X}}_{\mathit{A}}}^{\mathrm{*}}={\mathit{T}}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{X}}_{\mathit{A}}={\left(\begin{array}{cc}{\mathit{X}}^{{\mathit{*}}^{\mathrm{T}}}& {\mathit{R}}^{{\mathit{*}}^{\mathrm{T}}}\end{array}\right)}^{\mathrm{T}}。将变换后的{{\mathit{y}}_{\mathit{A}}}^{\mathrm{*}}，{{\mathit{X}}_{\mathit{A}}}^{\mathrm{*}}，{\mathit{u}}^{\mathrm{*}}代入A-Lasso模型中，得到如下估计：

\widehat{\mathit{\beta }}{}^{\mathrm{M}\mathrm{a}}{}^{-}{}^{\mathrm{L}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}}\left(\lambda \right)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\beta }{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\mathrm{R}\mathrm{S}{\mathrm{S}}^{\mathrm{*}}+\lambda {\displaystyle \sum _{j=\mathrm{1}}^p}\widehat{\mathit{w}}\left|{\mathit{\beta }}_j\right|\right\}，（13）

其中\mathrm{R}\mathrm{S}{\mathrm{S}}^{\mathrm{*}}=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^n}(y_{\mathit{A}}{{}^{\mathrm{*}}}_i{-{\displaystyle \sum _{j=\mathrm{1}}^p}{x}_A{{}^{\mathrm{*}}}_{ij}{\beta }_j)}^{\mathrm{2}}。其求解算法如下所示。

2 仿真实验

本节通过蒙特卡罗实验来比较上述所有方法的MSE。使用Tibshirani^［7］提到的仿真数据进行了四组实验，如下所示：

实验一：在这组实验中，生成样本量n=\mathrm{20}且 \mathit{\beta }={\left(\mathrm{3,1.5,0},\mathrm{0,2},\mathrm{0,0},\mathrm{0}\right)}^{\mathrm{T}}的数据。预测因子\mathit{X}服从均值为0，协方差为\mathit{\Sigma }的多元正态分布。\mathit{\Sigma }是对角元素为1，非对角元素为\mathrm{0}.{\mathrm{5}}^{\left|i-j\right|}的矩阵。误差项是由均值为0，标准差为3的正态分布生成的。

实验二：对于每一个j，有{\beta }_j=\mathrm{0.85}，其余参数与实验一相同。

实验三： {\beta }_{\mathrm{1}}=\mathrm{5}，{\beta }_j=\mathrm{0}\left(j=\mathrm{2,3},\dots ,\mathrm{8}\right)，误差项是由均值为0，标准差为2的正态分布生成的，其余参数与实验一相同。

实验四：样本量为n=\mathrm{100}且\mathit{\beta }={\left(\underset{\mathrm{10}}{\underbrace{\mathrm{0},\cdots ,\mathrm{0}}},\underset{\mathrm{10}}{\underbrace{\mathrm{2},\cdots ,\mathrm{2}}},\underset{\mathrm{10}}{\underbrace{\mathrm{0},\cdots ,\mathrm{0}}},\underset{\mathrm{10}}{\underbrace{\mathrm{2},\cdots ,\mathrm{2}}}\right)}^{\mathrm{T}}预测因子\mathit{X}由x_{ij}=z_{ij}+z_i生成，其中z_{ij}和z_i是独立的标准正态分布变量。误差项是由均值为0，标准差为15的正态分布生成。

误差项由两个不同过程产生：（1）White Noise（WN）误差：\mathit{\epsilon }~N(\mathrm{0},{\sigma }^{\mathrm{2}}\mathit{V})，\mathit{V}=I_n。（2）一阶自回归模型（AR（1））的误差：{\epsilon }_i=\mathrm{0.5}{\epsilon }_{i-\mathrm{1}}+e_i，其中e_i是由均值为0，方差为{\sigma }_e^{\mathrm{2}}={\sigma }^{\mathrm{2}}\left(\mathrm{1}-\mathrm{0}.{\mathrm{5}}^{\mathrm{2}}\right)的正态分布生成的。即\mathit{\epsilon }~N(\mathrm{0},{\sigma }^{\mathrm{2}}\mathit{V})，其中\mathit{V}=\left(\mathrm{0}.{\mathrm{5}}^{\left|i-j\right|}\right)。

随机限制

\mathit{r}=\mathit{R}\mathit{\beta }+\mathit{\varphi }；（14）

\mathit{\varphi }\sim N\left(\mathrm{0},{\sigma }^{\mathrm{2}}\mathit{W}\right)，（15）

其中\mathit{\varphi }的方差有两种：

（i）{\sigma }^{\mathrm{2}}\mathit{W}=\mathrm{0.1}{\mathit{I}}_m（Loose 先验信息），

（ii）{\sigma }^{\mathrm{2}}\mathit{W}=\mathrm{0.05}{\mathit{I}}_m（Tight先验信息）。

限制矩阵\mathit{R}有两种：

（1）　当p=\mathrm{8}时，\mathit{R}={\left[\begin{array}{cc}\left(\begin{array}{cc}\mathrm{1}& \begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\end{array}\right)& {\mathrm{0}}^{\mathrm{T}}\\ {\mathrm{0}}^{\mathrm{T}}& \left(\begin{array}{cc}\mathrm{1}& \begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\end{array}\right)\end{array}\right]}_{\mathrm{2}\times \mathrm{8}};

当p=\mathrm{40}时，\mathit{R}={\left[\begin{array}{ccc}\left(\begin{array}{cc}\mathrm{1}& \begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\end{array}\right)& \cdots & {\mathrm{0}}^{\mathrm{T}}\\ ⋮& \ddots & ⋮\\ {\mathrm{0}}^{\mathrm{T}}& \cdots & \left(\begin{array}{cc}\mathrm{1}& \begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\end{array}\right)\end{array}\right]}_{\mathrm{10}\times \mathrm{40}}\mathrm{。}

（2）　当p=\mathrm{8}时，\mathit{R}={\left[\begin{array}{cccc}\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)& {\mathrm{0}}^{\mathrm{T}}& {\mathrm{0}}^{\mathrm{T}}& {\mathrm{0}}^{\mathrm{T}}\\ {\mathrm{0}}^{\mathrm{T}}& \left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)& {\mathrm{0}}^{\mathrm{T}}& {\mathrm{0}}^{\mathrm{T}}\\ {\mathrm{0}}^{\mathrm{T}}& {\mathrm{0}}^{\mathrm{T}}& \left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)& {\mathrm{0}}^{\mathrm{T}}\\ {\mathrm{0}}^{\mathrm{T}}& {\mathrm{0}}^{\mathrm{T}}& {\mathrm{0}}^{\mathrm{T}}& \left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)\end{array}\right]}_{\mathrm{4}\times \mathrm{8}};

当p=\mathrm{40}时，\mathit{R}={\left[\begin{array}{cccc}\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)& {\mathrm{0}}^{\mathrm{T}}& \cdots & {\mathrm{0}}^{\mathrm{T}}\\ {\mathrm{0}}^{\mathrm{T}}& \left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)& \cdots & {\mathrm{0}}^{\mathrm{T}}\\ ⋮& ⋮& \ddots & {\mathrm{0}}^{\mathrm{T}}\\ {\mathrm{0}}^{\mathrm{T}}& {\mathrm{0}}^{\mathrm{T}}& \cdots & \left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)\end{array}\right]}_{\mathrm{20}\times \mathrm{40}}\mathrm{。}

其中第一种方案中m=\frac{p}{\mathrm{4}}，第二种方案中m=\frac{p}{\mathrm{2}}。

模拟过程中，参数\mathit{\beta }，\mathit{X}，\mathit{R}在仿真前设定。在每组实验中，\mathit{\epsilon }和\mathit{\varphi }由上述方案生成，\mathit{y}和\mathit{r}由公式（1）和（14）生成。通过GLS， RLS， ME， Lasso， A-Lasso， M-Lasso， Ma-Lasso方法进行参数估计。

本文通过均方误差（Mean Squared Error，MSE），发现百分比（Discovery Percentage， DP），真实发现率（True Discovery Percentage，TDP），选择真实模型的次数比例（Proportion of Times that the True model is Selected，PTTS）作为评价指标^［1］。公式如下：

\mathrm{M}\mathrm{S}\mathrm{E}\left(\widehat{\mathit{\beta }}\right)=\frac{\mathrm{1}}{\mathrm{10}\mathrm{000}}{\displaystyle \sum _{H=\mathrm{1}}^{\mathrm{10}\mathrm{000}}}{\left({\widehat{\mathit{\beta }}}_H-\mathit{\beta }\right)}^{\mathrm{T}}\left({\widehat{\mathit{\beta }}}_H-\mathit{\beta }\right)，（16）

\mathrm{D}\mathrm{P}\left(\widehat{\mathit{\beta }}\right)=\frac{\mathrm{1}}{\mathrm{10}\mathrm{000}}{\displaystyle \sum _{H=\mathrm{1}}^{\mathrm{10}\mathrm{000}}}\frac{\#({\widehat{\beta }}_{H,j}=\mathrm{0}/{\beta }_j=\mathrm{0})}{\#({\beta }_j=\mathrm{0})}\mathrm{100}，（17）

\mathrm{T}\mathrm{D}\mathrm{P}\left(\widehat{\mathit{\beta }}\right)=\frac{\mathrm{1}}{\mathrm{10}\mathrm{000}}{\displaystyle \sum _{H=\mathrm{1}}^{\mathrm{10}\mathrm{000}}}\frac{\#({\widehat{\beta }}_{H,j}\ne \mathrm{0}/{\beta }_j\ne \mathrm{0})}{\#({\beta }_j\ne \mathrm{0})}\mathrm{100}，（18）

其中{\widehat{\mathit{\beta }}}_H是第H次重复试验得到的估计值，\#\left(\cdot \right)是满足条件的次数。DP和TDP用于测量稀疏（或非稀疏）系数被正确估计的实际百分比。PTTS是当实际参数为零时，估计的模型参数为零，当实际参数不为零时，其不等于零的次数百分比，即满足所有零估计为零，所有非零估计为非零的次数百分比。

本文使用R语言（版本4.2.1）进行编程模拟，实验重复次数为10 000，研究结果如表1—表8所示。

表1—表8显示了不同方法在不同先验限制下的各种性能的结果。

表1和表2的结果对比表明，在一阶自回归模型（AR（1））的误差下的MSE明显小于White Noise（WN）误差下的MSE。Ma-Lasso在MSE方面展现出了显著的优势。其MSE不仅比M-Lasso小0.2，而且比其他方法减少了30%甚至更多。这表明在处理某些类型的数据时，Ma-Lasso具有更高的预测精度和稳定性。此外，从DP的结果来看，虽然Ma-Lasso不是最优的，但它在正确估计零系数的能力上与M-Lasso相差不大，并且明显优于其他方法。这进一步证明了Ma-Lasso在变量选择上的有效性。另外，从PTTS的结果上看，Ma-Lasso在带白噪声误差的随机限制下，其性能略低于M-Lasso，而在带有AR（1）误差的随机限制下，其性能与M-Lasso的性能相同。

表3和表4的结果对比表明，在非稀疏模型下，一阶自回归模型的误差下的MSE大于White Noise误差下的MSE。在非稀疏数据或特定条件下，M-Lasso的MSE表现最优，而Ma-Lasso虽然不是最优的，但在除M-Lasso外的其他方法中表现最好。同时，TDP的结果也显示Ma-Lasso在某些条件下表现最优，但在PTTS的结果中，M-Lasso再次优于Ma-Lasso。

表5和表6的结果对比与表1和表2的结果一致，进一步验证了稀疏模型下，一阶自回归模型的误差下的MSE要小于White Noise误差下的MSE。表5和表6的共同结果再次验证了Ma-Lasso在某些场景下的优势。其MSE比M-Lasso小0.1，并且远小于其他方法。这表明在处理稀疏数据时，Ma-Lasso仍然具有显著的优势。然而，从PTTS的结果来看，无论是带有白噪声误差还是AR（1）误差的随机限制条件下，Ma-Lasso都不再是最优的方法。

表7和表8的结果对比表明，在高维数据下，一阶自回归模型的误差下的MSE也大于White Noise误差下的MSE。根据表7和表8的共同结果，M-Lasso在MSE方面再次表现出最优的性能，而Ma-Lasso与其他方法相比存在一定的差距。从TDP的结果来看，Ma-Lasso仍然是最好的方法之一，但在PTTS的结果中，M-Lasso再次优于了Ma-Lasso。

综上所述，在系数稀疏的情况中，Ma-Lasso的效果比M-Lasso的效果好，而其在系数不是稀疏的情况中的效果没有M-Lasso好。这表明在处理具有稀疏解并且带有随机限制的线性模型时，Ma-Lasso比M-Lasso的效果好。并且由于数据的生成具有随机性，MSE的结果虽会在一定范围内波动，但总体来说Ma-Lasso在稀疏模型上的MSE优于其他模型。从DP的结果上来看，Ma-Lasso优于除M-Lasso以外的其他方法。这表明Ma-Lasso在发现真实的0系数能力上与M-Lasso相差无几，并且优于其他的估计量。从PTTS 的结果上看，M-Lasso的结果最优，Ma-Lasso的结果优势不明显但总体还是优于其他方法或与其他方法相差不大。因此，在实际应用中，我们需要根据具体的数据特征和问题背景来选择合适的模型和方法。

3 实际应用

我们将新方法应用于选择影响股票价格特征的数据集中^［19］，该数据收集了从2009年6月30日到2021年6月30日的贵州茅台各季度财报数据和股票价格数据。在主要指标中选取了营业收入（X_{\mathrm{1}}）、营业成本（X_{\mathrm{2}}）、营业利润（X_{\mathrm{3}}）、利润总额（X_{\mathrm{4}}）、所得税费用（X_{\mathrm{5}}）、净利润（X_{\mathrm{6}}）、每股收益（X_{\mathrm{7}}）、货币资金（X_{\mathrm{8}}）、存货（X_{\mathrm{9}}）、流动资产合计（X_{\mathrm{10}}）、固定资产净额（X_{\mathrm{11}}）、资产总计（X_{\mathrm{12}}）、流动负债合计（X_{\mathrm{13}}）、负债合计（X_{\mathrm{14}}）、所有者权益合计（X_{\mathrm{15}}）、期初现金及现金等价物余额（X_{\mathrm{16}}）、期末现金及现金等价物余额（X_{\mathrm{17}}）共计17个特征作为影响股票价格的截面数据特征。另外经过对数据进行探索性研究后，发现股票价格呈现一种指数型增长的趋势，而其他特征的增长模式为一般的线性增长，因此对股票价格p取对数。其所构建的模型如下式：

\mathrm{l}\mathrm{n}\left(p\right)={\beta }_{\mathrm{1}}{X}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{X}_{\mathrm{2}}+{\beta }_{\mathrm{3}}{X}_{\mathrm{3}}+\dots +{\beta }_{\mathrm{17}}{X}_{\mathrm{17}}+\epsilon。（19）

由于数据集没有先验限制，因此本文使用GLS和Lasso的估计值来形成这种先验限制。该限制为{\beta }_{\mathrm{4}}-{\beta }_{\mathrm{1}}-{\beta }_{\mathrm{6}}+{\beta }_{\mathrm{14}}=\mathrm{1}，{\beta }_{\mathrm{1}}+{\beta }_{\mathrm{4}}-{\beta }_{\mathrm{8}}+{\beta }_{\mathrm{13}}+{\beta }_{\mathrm{17}}=\mathrm{0.5}。

衡量模型选择功效的BIC准则计算公式如下：

\mathrm{B}\mathrm{I}\mathrm{C}=k\mathrm{l}\mathrm{n}\left(n\right)-\mathrm{2}\mathrm{l}\mathrm{n}\left(L\right)，（20）

其中k为模型参数的数量，n为样本个数，似然函数L如下式：

L={\left({\mathit{Y}}^{\mathrm{*}}-{\mathit{X}}^{\mathrm{*}}\mathit{\beta }\right)}^{\mathrm{T}}\left({\mathit{Y}}^{\mathrm{*}}-{\mathit{X}}^{\mathrm{*}}\mathit{\beta }\right)。（21）

根据表9的结果，各个参数估计值的符号是相似的，这表明以上方法在对该数据集进行建模时，对于参数的影响方向是一致的。另外通过估计值的结果发现，GLS， RLS， ME的估计值不是稀疏的，表明这些方法并不能很好地提供一个稀疏解。而从BIC 的值的结果来看，Ma-Lasso的BIC 值是最小的，说明Ma-Lasso可以很好地应用于所构建的模型，用于选择合适的影响股票价格的特征。从经济学角度来看，通常我们期望利润总额与股票价格之间存在正相关关系，即利润越高，股票价格越高；而所得税费用与股票价格之间通常会存在负相关关系，即所得税费用越高，股票价格越低。然而，根据GLS等的结果，这两个关系与通常的预期相反。但是从Ma-Lasso的结果来看，营业利润、利润总额以及资产总计对于股票价格有正影响，在一般情况下，我们期望营业利润、利润总额以及资产总计对股票价格具有正相关的影响，即这些指标的增加通常会伴随着股票价格的上涨。这表明Ma-Lasso所选出来的影响股票价格的特征与经济理论预期相吻合。

4 结论

受M-Lasso的启发，本文将ME和A-Lasso结合，提出Ma-Lasso方法，并通过仿真数据比较GLS， RLS， ME， Lasso， M-Lasso的MSE。结果表明，在稀疏模型中，Ma-Lasso的结果优于上述其他方法。这表明，Ma-Lasso在处理稀疏模型问题上能给出更精确的参数估计，较其他方法更有优势。最后将Ma-Lasso与其他方法一起应用到贵州茅台各季度财报数据和股票价格数据中，从所选特征和MSE的结果看，Ma-Lasso均给出一个较为理想的结果，符合经济理论。这进一步说明Ma-Lasso在处理实际问题中有其独有的优越性，可以应用在金融领域等领域。

尽管Ma-Lasso在稀疏模型上的表现优于M-Lasso，并能更好地选择影响股票价格的特征，但是在非稀疏模型上其表现不如M-Lasso。因此，实际应用中可根据需求进行综合考虑，选择合适的方法。

参考文献

[1] GULER H, GULER E O. Sparsely Restricted Penalized Estimators[J]. Commun Stat Theory Meth, 2021, 50(7): 1656-1670. DOI: 10.1080/03610926.2019.1682164 .

[2] DONOHO D L. High-dimensional Data Analysis: The Curses and Blessings of Dimensionality[C]//American Mathematical Society Math Challenges Lecture. Providence, Rhode Island: AMS, 2000, 1(2000): 1-32.

[3] FAN J Q, HAN F, LIU H. Challenges of Big Data Analysis[J]. Natl Sci Rev, 2014, 1(2): 293-314. DOI: 10.1093/nsr/nwt032 .

[4] 李欣. 高维数据的稀疏估计及其应用[D]. 杭州: 浙江大学, 2019.

[5] LI X. Sparse Estimation for High-dimensional Data with Applications[D]. Hangzhou: Zhejiang University, 2019.

[6] 苏锦霞. 基于特征选择的高维数据统计分析[D]. 兰州: 兰州大学, 2018.

[7] SU J X. Statistical Analysis of High-dimensional Data Based on Feature Selection[D]. Lanzhou: Lanzhou University, 2018.

[8] HOERL A E, KENNARD R W. Ridge Regression: Biased Estimation for Nonorthogonal Problems[J]. Technometrics, 1970, 12(1):55-67. DOI: 10.1080/00401706.1970.10488634 .

[9] TIBSHIRANI R. Regression Shrinkage and Selection via the Lasso[J]. J R Stat Soc Ser B Methodol, 1996, 58(1): 267-288. DOI: 10.1111/j.2517-6161.1996.tb02080.x .

[10] 胡蓉. 基于随机Lasso的Meta分析[D]. 北京: 北京建筑大学, 2019.

[11] HU R. Meta Analysis Based on Random Lasso[D]. Beijing: Beijing University of Civil Engineering and Architecture, 2019.

[12] 王璐, 孙聚波. Lasso回归方法在特征变量选择中的应用[J]. 吉林工程技术师范学院学报, 2021, 37(12): 109-112. DOI: 10.3969/j.issn.1009-9042.2021.12.032 .

[13] WANG L, SUN J B. Application of Lasso Regression Method in the Selection of Feature Variables[J]. J Jilin Teach Inst Eng Technol, 2021, 37(12): 109-112. DOI: 10.3969/j.issn.1009-9042.2021.12.032 .

[14] 苏宇腾, 吕世云, 谢文涵, 等. 基于LASSO回归与随机森林算法的2型糖尿病发病风险因素分析[J]. 环境卫生学杂志, 2023, 13(7): 485-495. DOI: 10.13421/j.cnki.hjwsxzz.2023.07.002 .

[15] SU Y T, LYU S Y, XIE W H, et al. A Risk Factor Analysis for Type 2 Diabetes Mellitus Based on LASSO Regression and Random Forest Algorithm[J]. J Environ Hyg, 2023, 13(7): 485-495. DOI: 10.13421/j.cnki.hjwsxzz.2023.07.002 .

[16] 邢耀. 胖(高)大数据的LASSO Logistic模型构建与应用: 以贷款偿还与糖尿病两数据集为例[D]. 济南: 山东大学, 2022.

[17] XING Y. Construction and Application of LASSO Logistic Model for Fat (High) Big Data-taking Loan Repayment and Diabetes as Examples[D]. Jinan: Shandong University, 2022.

[18] 车前子, 王晶, 白卫国, 等. 基于LASSO回归的骨质疏松肾阳虚状态辨识模型研究[J]. 中华中医药杂志, 2022, 37(10): 5928-5933.

[19] CHE Q Z, WANG J, BAI W G, et al. Study on the State Identification Model of Kidney Yang Deficiency in Osteoporosis Patients Based on LASSO Regression[J]. China J Tradit Chin Med Pharm, 2022, 37(10): 5928-5933.

[20] ZOU H, HASTIE T. Regularization and Variable Selection via the Elastic Net[J]. J R Stat Soc Ser B Stat Methodol, 2005, 67(2): 301-320. DOI: 10.1111/j.1467-9868.2005.00503.x .

[21] ZOU H. The Adaptive Lasso and Its Oracle Properties[J]. J Am Stat Assoc, 2006, 101(476): 1418-1429. DOI: 10.1198/016214506000000735 .

[22] THEIL H. On the Use of Incomplete Prior Information in Regression Analysis[J]. J Am Stat Assoc, 1963, 58(302): 401-414. DOI: 10.1080/01621459.1963.10500854 .

[23] GULER H, GULER E O. Mixed Lasso Estimator for Stochastic Restricted Regression Models[J]. J Appl Stat, 2021, 48(13/14/15): 2795-2808. DOI: 10.1080/02664763.2021.1922614 .

[24] 刘洪伟, 徐文科. 线性模型的广义最小二乘估计递推算法[J]. 哈尔滨师范大学自然科学学报, 2011, 27(3): 29-31. DOI: 10.3969/j.issn.1000-5617.2011.03.009 .

[25] LIU H W, XU W K. The Recursive Algorithm of Generalized Least Squares Estimation for Linear Model[J]. Nat Sci J Harbin Norm Univ, 2011, 27(3): 29-31. DOI: 10.3969/j.issn.1000-5617.2011.03.009 .

[26] CONWAY R N, MITTELHAMMER R C. The Theory of Mixed Estimation in Econometric Modeling[J]. Stud Econ Finance, 1986, 10(1): 79-120. DOI: 10.1108/eb028665 .

[27] 赵昊然. 基于LASSO和BP神经网络的股价价格预测研究[D]. 桂林: 广西师范大学, 2022.

[28] ZHAO H R. Research on Stock Price Forecasting Based on LASSO and BP Neural Network[D]. Guilin: Guangxi Normal University, 2022.

基金资助

山西省高校科技创新项目(2021L324)

山西省基础研究计划青年科学研究项目(202103021223282)

太原科技大学博士科研启动基金(20212014)