兩維異質性面板分位數模型的雙懲罰回歸方法

任燕燕，李東霖，王文悅

（山東大學 經濟學院，濟南 250100）

0 引言

在現代經濟問題實證研究中應用的樣本數據逐漸呈現大數據特征，由于數據來源具有多樣性、復雜性，經濟大數據集表現出顯著的異質性特點，原有基于同質性假設的單一化等傳統建模方法無法全面揭示經濟運行的客觀規律。因此，異質性分析是計量經濟學和統計學中的一個研究熱點，如何處理經濟大數據中的異質性并分析經濟變量之間的結構關系或變化是該研究的關鍵所在。利用面板數據模型進行異質性分析的研究大多基于均值回歸假設進行，回歸結果只能刻畫均值水平上的異質性結構，未能全面考慮數據所蘊含的信息。而基于面板分位數模型的異質性分析方法能夠從被解釋變量條件分布的角度，更加全面地反映不同分位點處變量之間的回歸關系，并識別其異質性結構。其中，一部分學者基于個體維度的異質性，研究組群結構的識別方法。Zhang 等（2019）[1]考慮個體維度上的組群結構，提出了面板分位數模型的聚類方法，并借鑒機器學習中的K-means 算法開發了用于識別某一分位點處或不同分位點間組群結構的迭代算法；Zhang 等（2023）[2]提出了一種關于面板分位數模型組群結構的非參數估計方法，使用兩兩融合的懲罰項估計了組群結構的數量，通過ADMM 迭代算法提高了運算效率。另外一部分學者關注時間維度上的異質性，試圖尋找未知的結構變點。Wang 等（2009）[3]基于回歸系數是關于時間的平滑非參數函數的假設，構建了變系數分位數回歸模型，并且使用基函數近似法進行參數估計；部分學者提出了具有結構變點的分位數自回歸模型，根據非對稱Laplace 分布構建似然函數，實現了對結構變點的估計；Lee等（2018）[4]構建了具有結構變點的分位數回歸模型，利用L1懲罰項結合加權絕對偏差項構建目標函數，估計了該模型的回歸系數和結構變點，同時證明了估計量具有Oracle屬性。上述針對異質性面板分位數模型的研究，普遍基于個體或時間的單一維度進行分析，然而隨著異質性分析方法的不斷發展，一些實證研究結果表明，變量之間的關系可能同時存在個體和時間兩個維度上的異質性，因此有必要假設模型存在兩維的異質性結構。

兩維異質性面板分位數模型存在嚴重的冗余參數問題，因此減少待估參數的數量成為相關研究的關鍵。一種方法是將兩維異質性系數分解為個體維度的異質性系數、時間維度的異質性系數和共同系數三項之和，稱為維數縮減法[5]。此類方法存在模型識別問題，系數分解需要滿足個體和時間維度的異質性系數均值為零的假設，重點考察通過系數分解得到的共同系數。另一種方法假設系數在個體和時間兩個維度上存在異質性結構，根據不同問題的研究背景，設定兩維異質性系數具有不同的稀疏結構，基于某一維度的異質性結構關注另一個維度上是否存在組群結構或結構變點。Okui 和Wang（2021）[6]在基于個體維度識別出組群結構的基礎上，考察不同組群時間維度上的結構變點，提出利用Adaptive Lasso進行組群融合的方法；Lumsdaine 等（2020）[7]關注面板數據模型在時間維度上的某一結構變點，基于K-means算法識別該結構變點附近個體維度上的組群結構是否發生變化。此外，部分研究假設模型系數在個體或時間中某一個維度上存在稀疏結構，而在另一個維度上存在完全異質性。Baltagi等（2016）[8]允許不同個體的斜率系數存在差異，但斜率系數具有相同的結構變點；Su 等（2017）[9]假定不同個體之間存在組群結構，并且具有時間維度上的時變非參數形式的異質性系數，利用篩法或B-樣條將個體維度和時間維度進行分離，考慮分離出的個體維度并識別其組群結構。迄今為止，尚未有文獻研究面板分位數模型的兩維異質性問題。本文考慮系數具有兩維異質性結構的面板分位數模型，基于SCAD懲罰函數和MCP懲罰函數提出一種能夠同時進行參數估計和兩維異質性結構識別的雙懲罰回歸方法。

1 模型設定和估計方法

1.1 兩維異質性面板分位數模型

本文考慮如下的兩維異質性面板數據模型：

其中，yit表示一維被解釋變量，μit表示隨個體和時間同時變化的固定效應項，zit=(zit,1,zit,2,…,zit,P-1)Τ表示P-1 維 解 釋 變 量 向 量，δit=(δit,1,δit,2,…,δit,P-1)Τ表 示P-1 維異質性解釋變量系數，εit表示一維隨機誤差項。定義和。在分位點τ處，可構建如下的條件分位數函數：

模型存在冗余參數問題，為了減少待估參數，本文假設斜率系數βit(τ)具有如下的兩維異質性結構：

其中，B(τ)={B1(τ),B2(τ),…,BL(τ)(τ)}表示兩維異質性結構，L(τ)表示其數量。

1.2 雙懲罰最小加權絕對偏差目標函數

其中：ρτ(u)=u?(τ-I(u<0))表示分段線性分位數損失函數；I(?)表示示性函數，當括號內不等式成立時取值為1，否則取值為0。因此ρτ(u)的分段函數為：

pλ(τ)(?)和pγ(τ)(?)表示成對融合懲罰函數，分別基于個體維度和時間維度對兩維異質性結構進行融合，λ(τ)和γ(τ)表示控制懲罰力度的調節參數。常用的Lasso懲罰函數pλ(t)=λt對所有的或使用相同的閾值，因此過度收縮了較大的模型參數，導致估計量存在偏差并且可能無法正確識別模型的異質性結構。因此，本文使用能夠產生無偏估計的懲罰函數，包括Fan和Li（2001）[10]提出的SCAD懲罰函數：

Zhang（2010）[11]提出的MCP懲罰函數：

其中：(?)+=max(0,?)；a表示控制懲罰函數凹性的參數，在本文中視作一個固定常數；t表示基于個體維度或時間維度模型參數的差異，即‖βit(τ)-βjt(τ) ‖或‖βit(τ)-βit′(τ) ‖。SCAD 懲罰函數和MCP 懲罰函數均為漸進無偏的，能夠更加精確地識別模型參數的異質性結構。特別地，當a→+∞時，SCAD懲罰函數和MCP懲罰函數收斂到Lasso懲罰函數。

對于給定的λ(τ)和γ(τ)，參數估計量可通過最小化目標函數給出：

1.3 參數估計的ADMM算法

考慮到最小化目標函數S(β(τ);γ(τ),λ(τ))這一無約束優化問題的顯式解不存在，本文擬使用求解優化問題的交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）估計分位點τ處的斜率系數β(τ)。基于目標函數S(β(τ);γ(τ),λ(τ))，考慮如下約束條件：

其中，ηij,t(τ)表示在分位點τ處和相同的時間點t不同個體i和j之間斜率系數的差異，θi,tt′(τ)表示在分位點τ處同一個體i在不同時間點t和t′之間斜率系數的差異。因此，可以將原有的無約束優化問題轉化成有約束的優化問題：

其中，η(τ)=(t(τ),i

本文使用ADMM 算法，從給定的初始值開始逐步更新β(τ)、η(τ)、θ(τ)、v(τ)和w(τ)的迭代值。假設第s步的迭代值β(s)(τ)、v(s)(τ)和w(s)(τ)是已知的，第s+1 步使用ADMM 算 法 獲 得 迭 代 值β(s+1)(τ) 、η(s+1)(τ) 、θ(s+1)(τ) 、v(s+1)(τ)和w(s+1)(τ)的具體步驟為：（1）根據β(s)(τ)、v(s)(τ)和w(s)(τ)，更新迭代值η(s+1)(τ)和θ(s+1)(τ)；（2）根據v(s)(τ)和w(s)(τ)以及第一步迭代得到的η(s+1)(τ)和θ(s+1)(τ)，更新斜率系數的迭代值β(s+1)(τ)；（3）根據相關參數，更新迭代值v(s+1)(τ)和w(s+1)(τ)。基于以上步驟，更新迭代值的詳細過程為：

首先，更新迭代值η(s+1)(τ)和θ(s+1)(τ)的子優化問題為：

使用不同的懲罰函數，即使在目標函數不可導的情況下，依然能使用凸分析中的次微分得到和的迭代公式。

如果參數滿足a>max(1/φ+1,1/?+1)并且使用SCAD懲罰函數，那么迭代公式為：

其中，ST(z,t)=(1-t/‖ω‖ )+ω表示軟閾值算子（Soft Thresholding Operator），并且(?)+=max(0,?)。

如果參數滿足a>max(1/φ,1/?)并且使用MCP懲罰函數，那么迭代公式為：

其次，模型斜率系數的迭代值β(s+1)(τ)可以通過如下無約束最優化問題計算：

其中，η(s+1)(τ)、θ(s+1)(τ)、v(s)(τ)和w(s)(τ)均已知。目標函數等價于：

最后，使用對偶上升法（Dual Ascent）計算對偶變量v(τ)和w(τ)的迭代值，迭代方向取增廣拉格朗日目標函數L對對偶變量的次微分，迭代步長選擇固定的懲罰參數φ和?。因此，(τ)和(τ)的計算公式為：

2 蒙特卡洛模擬

2.1 數據生成過程

本文使用蒙特卡洛模擬驗證提出的基于雙懲罰最小加權絕對偏差目標函數系數估計量的有限樣本性質。每輪蒙特卡洛模擬使用的數據{yit,xit}（i=1,…,N；t=1,…,T）基于如式（26）所示的數據生成過程得到。不失一般性，考慮二變量面板數據模型：

其中，N=10,20，T=10,20，固定效應項μit和解釋變量系數ηit具有相同的兩維異質性結構。假設兩維異質性結構的數量為L=2，并且每一個兩維異質性結構系數的真實值為α1=(-2,3)和α2=(3,5)。具體兩維異質性結構 可 任 意 設 定，例 如N=20 和T=20 ，當6 ≤i≤8 且8 ≤t≤12，以及9 ≤i≤13 且6 ≤t≤15 時，(i,t)屬于兩維異質性結構B2；否則，(i,t)屬于兩維異質性結構B1。B1和B2元素數量的比例 |B1|:|B2|≈5:1。假設xit的生成過程滿足：

2.2 蒙特卡洛模擬結果

本文在蒙特卡洛模擬的過程中選擇重復次數n=200 。調節參數λ和γ的網格搜索篩選范圍為[0.1,1.9]，搜索步長為0.2。為了驗證估計量的有限樣本性質，本文將基于SCAD 懲罰函數和MCP 懲罰函數的估計結果與Post 估計量（見式（28））和Oracle 估計量（見式（29））進行比較。

其中，B?表示兩維異質性結構的估計值，L?表示兩維異質性結構數量的估計值，B0表示兩維異質性結構的真實值，L0表示兩維異質性結構數量的真實值。根據參數α的估計量和相應的兩維異質性結構，能夠得到參數β的Post估計量和Oracle估計量。

為了評估模型參數估計結果的精確度，下頁表1和表2分別匯報了基于同方差假設和異方差假設的中位點處、0.25 分位點處和0.75 分位點處系數估計值的均方誤差和平均偏差。結果顯示，在各分位點處，本文提出的應用SCAD 懲罰函數和MCP 懲罰函數的參數估計方法的估計結果相似，并且Post 估計量的RMSE 與Bias 和Oracle 估計量差別不大。

表1 同方差假設下參數估計的RMSE和Bias

表2 異方差假設下參數估計的RMSE和Bias

為了檢驗本文提出的兩維異質性結構識別方法的準確性，下頁表3和表4分別匯報了同方差假設下和異方差假設下不同分位點處兩維異質性結構估計的準確率和Rand 指數。結果表明，本文提出的方法在所有情況下都能準確地識別各分位點處的兩維異質性結構，而且基于SCAD 懲罰函數和MCP 懲罰函數的估計方法在兩維異質性結構識別方面的表現相近。

表3 同方差假設下兩維異質性結構估計的Acc和RI

表4 異方差假設下兩維異質性結構估計的Acc和RI

3 應用

本文將提出的雙懲罰最小加權絕對偏差估計方法應用于2009—2019 年我國31 個省份（不含港澳臺）的GDP數據，以識別不同省份和不同年份經濟發展的兩維異質性結構并進行參數估計，研究使用的數據來源于國家統計局官網。考慮如下兩維異質性面板分位數模型：

其中，l n(GDPit)表示省份i第t年人均GDP的自然對數值，t=0 對應于2009 年，βit(τ)=(βit,0(τ),βit,1(τ),βit,2(τ))Τ表示具有未知兩維異質性結構的異質性斜率系數，(i,t)?Bl，{B1,B2,…,BL}表示未知的兩維異質性結構，εit(τ)表示隨機誤差項。調節參數λ(τ)和γ(τ)的選擇范圍為[0.1，1.9]，步長為0.2，考慮τ=0.25，τ=0.5 和τ=0.75 三個分位點。應用本文方法估計上述模型，估計結果表明，各分位點處的兩維異質性結構數量L?(τ)均為2。

表5列出了模型的參數估計結果，發現各分位點處識別出不同兩維異質性結構的參數存在顯著差異。首先，不同兩維異質性結構的截距項在各分位點處均存在顯著差異；其次，GDP 一階滯后項的系數即自回歸系數在不同兩維異質性結構下和不同分位點處均顯著但差距較小；最后，各分位點處不同兩維異質性結構之間時間斜率系數的差異主要體現在其顯著性上，其中一個兩維異質性結構時間的斜率系數顯著，而另一個兩維異質性結構的系數不顯著。上述結果表明，我國各省份在研究期間的經濟發展狀況存在兩維異質性，本文提出的方法實現了對兩維異質性結構的估計。

表5 各分位點處不同兩維異質性結構的參數估計結果

4 結論

本文研究了具有兩維異質性結構的面板分位數模型及其參數估計問題。利用SCAD 懲罰函數和MCP 懲罰函數在個體維度和時間維度上對異質性參數進行融合，構建了雙懲罰最小加權絕對偏差目標函數，并通過設計ADMM迭代算法求解該目標函數，實現了參數估計和兩維異質性結構識別。根據蒙特卡洛模擬驗證了雙懲罰估計量的有限樣本性質，無論基于同方差假設還是異方差假設，本文提出的方法均能準確地識別各分位點處的兩維異質性結構，并且Post估計量的RMSE和Bias接近于假設真實兩維異質性結構已知的Oracle估計量。進一步地，將本文方法應用于中國省級GDP回歸的估計，估計結果表明，我國各省份在研究期間的經濟增長存在兩維異質性，本文提出的方法能夠識別這種兩維異質性結構。