999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于最優指標組合的企業信用風險預測

2021-10-10 02:05:24蘇小婷
系統管理學報 2021年5期
關鍵詞:模型

周 穎,蘇小婷

(大連理工大學 經濟管理學院,遼寧 大連 116024)

信用風險預測是指構建企業歷史數據與違約狀態之間的對應關系,揭示企業的經營發展狀況,進而對企業在未來是否會發生違約做出預判。信用風險預測結果不僅能為商業銀行提供重要的貸款決策依據,而且能幫助在股票市場、債券上的投資者做出正確的投資決策。

上市公司是中國國民經濟發展的關鍵。據Wind數據庫統計,截至2019年前三季度,中國A股上市公司總營收達36 萬億元,同比增長9.54%,占GDP比重突破50%。與此同時,出現財務虧損被證監會特別處理(ST)的上市公司數量也在持續上升,2019 年內共有86 家公司被實施ST,18家公司退市。這一數據也暴露出中國上市公司規模在不斷擴張的同時,也存在嚴重的信用風險問題。因此,建立有效的信用風險預測模型顯得極為重要。

建立信用風險預測模型主要涉及兩個問題:①信用風險預測模型中指標組合的遴選。與企業信用狀態相關聯的指標眾多,不同的指標組合違約鑒別能力不同,有必要尋找一個違約鑒別能力最強的指標組合,最大限度地區分違約客戶和非違約客戶。②非平衡樣本的處理。上市公司信用數據庫中違約客戶(ST)的數量遠小于非違約客戶(非ST)的數量,屬于不平衡樣本。在利用不平衡樣本建模時,會導致模型對違約客戶的判對率下降,第2類錯誤(Type-II error)上升,進而降低模型整體的預測精度。在信用違約預測中第2類錯誤是指將違約客戶誤判為非違約客戶,若銀行錯誤識別上市公司的信用狀況,將貸款發放給違約客戶,銀行將面臨客戶無法償還貸款帶來的巨大損失。

本文最優指標組合的遴選與文獻[1]中的研究相似。與文獻[1]的差別主要有兩點:①研究問題不同。文獻[1]中研究的是腦電信號的處理,本文研究的是信用風險預測。本文將近鄰成分分析引入信用風險領域進行指標組合遴選。②指標重要性的判別標準不同。文獻[1]中取指標權重閾值T=τmax(w),τ=0.02作為指標重要性的判別標準,通過剔除指標權重小于閾值T的指標構建指標組合。本文將τ分別設置為0.2,0.02,0.002,…,2×10-8,得到8個不同的指標組合,以指標組合的違約鑒別能力AUC最大反推指標重要性的判別標準τ和最優指標組合。

本文隨機欠采樣樣本配比與文獻[2]中的研究相仿。與文獻[2]的區別在于,文獻[2]中根據經驗將違約客戶與非違約客戶的比例設置為1∶2.5,1∶5,1∶10。本文遍歷違約客戶與非違約客戶的所有可能的樣本比例,以G-mean最大反推違約預測模型中兩類客戶的最佳配比。改變了隨機欠采樣中主觀確定采樣比例的不合理做法。

本文的貢獻有:

(1)首次將近鄰成分分析引入信用風險領域進行指標組合遴選。在近鄰成分分析算法中根據違約判別準確率最大得到馬氏距離中的指標權重向量,以基于馬氏距離的K-近鄰的違約判別誤差MSE最小為目標,確定最優指標權重向量;給定一個指標權重閾值通過剔除權重小于閾值的指標得到一個指標組合,給定不同的指標權重閾值得到不同的指標組合,以指標組合的違約判別精度AUC 最大反推最優的指標組合。拓展了信用風險領域指標組合遴選的新思路。

(2)利用隨機欠采樣將違約客戶與非違約客戶組成不同比例的樣本,以基于線性支持向量機的違約預測精度G-mean最大為標準反推違約客戶與非違約客戶的最佳比例,以確定最優的訓練樣本。

(3)通過t-m年的指標數據xt-m和t年的企業違約狀態yt,利用最優指標組合和最優訓練樣本建立了支持向量機模型,達到了運用t年的數據xt預測第t+m年企業違約狀態yt+m的預測效果。

(4)研究結果表明,本文的違約預測模型的精度高于非線性SVM、LR、DT、KNN 和LDA 等典型的大數據預測模型。

研究表明:每股收益EPS-扣除/稀釋、貨幣供應量M0(億元)和貨幣供應量M1(億元)3個指標對企業未來1~3年的短期違約狀態具有關鍵影響;當日總市值/負債總計、每股EBITDA 和固定資產周轉率3個指標對企業未來4~5年的長期違約狀態具有關鍵影響;經營活動產生的現金流量凈額/經營活動凈收益和審計意見類型2個指標,不論對于企業未來1~3年的短期、還是未來4~5年的長期違約狀態,均有關鍵影響。

1 文獻綜述

1.1 信用風險預測指標遴選的研究現狀

(1)單指標。現有研究大多為單指標遴選,Song等[3]通過fisher判別分析和F值對指標重要性進行排序,刪除所有不重要指標,降低了計算成本。Gunduz等[4]提出平衡互信息(BMI)方法檢測指標之間的非線性相關關系,進而剔除冗余指標。Lin等[5]在預測P2P 貸款違約時,使用RFE 方法篩選貸款違約行為變量。Gartner等[6]采用信息熵篩選具有鑒別能力的指標。Abbasi等[7]發現,添加行業類指標能顯著提高違約預測模型的精度。周穎等[8]通過構造Brown-Mood中位數檢驗統計量值的方法遴選出違約鑒別力顯著的指標。Raghu等[1]通過剔除權重小于閾值的指標構建指標體系。

(2)指標組合。Kozodoi等[9]以利潤最大化和指標數量最少為目標函數構建指標體系。葛興浪等[10]根據企業各類指標的偏相關系數和有序Probit回歸系數構建Wald統計量進行指標體系的遴選。Uthayakumar等[11]利用蟻群優化算法構建信用風險預測的最優指標組合。Ping等[12]不斷向指標組合中添加指標,通過指標群依賴度的變化量來反映所添加指標的重要性,剔除重要性為零的指標,最終構建信用評價指標體系。石寶峰等[13]通過Logistic回歸顯著性判別遴選對農戶違約狀態影響顯著的指標,建立了由年齡、非農收入/總收入等13個指標組成的農戶小額貸款信用評級指標體系。Jadhav等[14]提出了一種由信息增益和遺傳算法相結合的新方法篩選信用評價指標體系。Oreski等[15]基于神經網絡和遺傳算法篩選出一個包含12個指標的信用評價體系。

單指標遴選的弊端是單個指標違約鑒別能力強,組成的指標組合違約鑒別能力不一定強。現有研究關于指標組合遴選的弊端在于沒有以指標組合的違約判別精度最大為標準構建指標體系。本文首次將近鄰成分分析(NCA)引入信用風險領域進行指標組合遴選,以指標組合的違約預測精度AUC最大反推最優的指標組合。彌補了現有研究忽略指標組合整體違約鑒別能力的弊端,拓展了信用風險領域指標組合遴選的新思路。

1.2 隨機欠采樣的研究現狀

傳統的欠采樣方法是隨機在多數類樣本中抽取一定數量的樣本與少數類樣本組成比例為1∶1的平衡樣本[16-18],但是傳統的欠采樣方法會丟棄大量數據和重要信息,導致模型過擬合。Louzada等[19]基于Logistic回歸研究巴西銀行客戶的信用風險問題,發現利用隨機欠采樣獲取的平衡樣本能顯著提高先用評估模型的效果。Jian等[20]對重要違約樣本采用過采樣方法生成新違約樣本,對不重要非違約樣本采用欠采樣方法剔除非違約樣本。Perols等[21]在預測公司欺詐時,將多數類樣本劃分5個準則層,在每個準則層中隨機抽取樣本與少數類樣本組合,解決了樣本的非平衡問題。Dubey等[22]采用K-Medoids欠采樣技術對多數類樣本聚類,將聚類后的樣本與少數類樣本組成平衡樣本。Paleologo等[23]將非平衡的原始數據集分為若干個子數據集,并通過改變每個子數據集中違約和非違約兩類樣本的比例進行比較分析。

本文利用隨機欠采樣技術在非違約客戶中抽取樣本,將違約客戶與非違約客戶組成比例為1∶1,1∶2,1∶3,1∶4,1∶5,1∶6的樣本,以基于支持向量機的違約預測精度G-mean最大反推違約客戶與非違約客戶的比例。

1.3 支持向量機預測信用風險的研究現狀

不少學者利用支持向量機建立信用風險預測模型,并取得了較高的預測精度[24-26]。Schebesch等[27]對比了線性和非線性支持向量機的違約判別效果,發現非線性支持向量機在實際應用中并沒有表現出優勢。Kim 等[28]考慮企業財務指標、宏觀經濟條件和企業管理水平等因素構建支持向量機模型,為韓國中小企業提供信用風險評估模型。Danenas等[29]利用粒子群優化選取最優的線性支持向量機的懲罰系數C,進而構建有效的信用評估模型。Maldonado等[30]將利潤函數納入支持向量機模型,構建智利銀行的信用評分系統。

現有研究大多證實了在大型數據庫中線性支持向量機的分類效果不會低于非線性支持向量機[31-32]。線性支持向量機能夠輸出每個指標的權重,具有較好的可解釋性,同時,與非線性支持向量機相比,線性支持向量機復雜度更低、運算速度更快。因此,本文選用線性支持向量機作為最終的信用風險預測模型。

2 信用風險預測模型的構建

2.1 最優指標組合的遴選

在原始樣本中構建最優指標體系。

2.1.1 基于偏相關性分析的第1次指標篩選

步驟1同一準則層指標相關系數rhg的計算。

設rhg為第h個指標與第g個指標的相關系數,xhj為第h個指標第j個客戶的指標數據為第h個指標的平均值,xgj為第g個指標第j個客戶的指標數據為第g個指標的平均值。則指標h和指標g的相關系數為[32]

式(1)中,第h個指標與第g個指標的相關系數越大,表明第h個指標與第g個指標的相關性越強;反之,相關性越弱。

設R為指標h與指標g之間的相關系數rhg組成的q×q矩陣,q為準則層內指標的個數。則[33]

R的逆矩陣記為

則指標h與指標g之間的偏相關系數為[32]

式(4)中,指標h與指標g之間的偏相關系數prhg越大,表明指標h與指標g之間的相關性越強;反之,相關性越弱。

步驟2F值的計算。

設Fh為第h個指標的F值為非違約客戶中第h個指標的均值為違約客戶中第h個指標的均值為全部客戶中第h個指標的均值,n(0)為非違約客戶的個數,n(1)為違約客戶的個數,xhj為第h個指標第j個客戶的指標數據,n為客戶總數。則[33]

式(5)的經濟學含義:式(5)等號右邊的分子第1項是第h個指標中非違約客戶的均值與所有客戶均值的距離,分子第2項是第h個指標中違約客戶與所有客戶均值的距離。整個分子表示第h個指標中違約客戶均值、非違約客戶均值與全部客戶均值的距離,反映違約客戶與非違約客戶的差異。分子越大、差異越大,表明第h個指標越能區分企業的違約狀態。分母中第1項是第h個指標中非違約客戶與非違約客戶均值的方差,第2項是第h個指標中違約客戶與違約客戶均值的方差,整個分母是第h個指標中違約客戶內的方差與非違約客戶內的方差之和,反映了違約客戶、非違約客戶各自的離散程度,離散程度越小,表明違約客戶、非違約客戶內部的指標特征越集中。式(5)中的Fh表示第h個指標的違約鑒別能力,Fh越大,表明指標h的違約鑒別能力越強;反之,越弱。

偏相關性分析是在控制其他變量的線性影響的條件下分析兩變量間的線性相關性。避免了當第3個指標同時影響兩個指標時,相關性分析不能如實反映兩個指標間相關程度的弊端。

步驟3基于偏相關性分析篩選指標的標準。

基于偏相關性分析進行第1次指標篩選有兩個標準:①保留具有經濟學含義的指標。②在標準①的基礎上,計算任意兩個指標的偏相關系數,若兩個指標的偏相關系數大于0.8,說明這兩個指標高度相關,刪除F值較小的指標;若兩個高度相關的指標均具有較強的經濟學含義,則兩個指標均保留。

2.1.2 基于近鄰成分分析的第2次指標篩選

步驟4定義馬氏距離公式。

設dw(xj,xz)為客戶j與客戶z的距離,s為經過第1次指標篩選后剩余的指標個數,wi為第i個指標的權重,xij為第i個指標第j個客戶的指標數據,xiz為第i個指標第z個客戶的指標數據。則[1]

式(6)的經濟學含義:客戶i與客戶z的馬氏距離越小,表明兩個客戶的違約狀態越可能相同。保證了違約狀態相同的客戶之間距離較小,違約狀態不同的客戶之間距離較大。

步驟5客戶i與客戶z相似概率的確定。

考慮一個隨機分類準則:存在一個由n個客戶組成的樣本,對于其中一個客戶j,可以在除自身之外的n-1 個客戶中隨機選取一個客戶作為參考點,以參考點的違約狀態判別客戶j的違約狀態。

設pjz(w)為客戶j與客戶z的相似概率,由于客戶j選擇客戶z作為參考點,客戶j以客戶z的違約狀態作為客戶j自身違約狀態的理論標識,故亦可稱pjz(w)為客戶j選擇客戶z作為參考點的概率。dw(xj,xz)為客戶j、z之間的馬氏距離,n為客戶總數。則[1]

式(7)中的分子是客戶j、z距離的函數,分母是客戶j與剩余n-1個客戶距離函數的和。式(7)是以距離衡量的概率,當客戶j、z之間的距離越小時,違約狀態的相似概率越大,客戶z被客戶j選為參考點的概率Pjz越大。式(7)表示兩個客戶違約狀態相似性的度量。

當客戶j的實際違約狀態與參考點z的違約狀態相同時,則客戶j的違約狀態通過客戶z被判對;當客戶j的實際違約狀態與參考點z的違約狀態不同時,則客戶j的違約狀態通過客戶z被判錯。在n個客戶中,任意兩個客戶都可以計算其相似程度,故該相似程度pjz(w)共有n個里邊取2的組合,即pjz(w)共有個。對于任意兩個客戶的指標數據xij和xiz,都可以得到相似度pjz(w)。

步驟6客戶j與其他全部n-1個客戶相似性的總和。

設pj(w)為客戶j的違約狀態通過隨機分類準則被判對的概率,n為客戶總數,pjz(w)為客戶j選擇客戶z作為參考點的概率,yj為客戶j實際的違約狀態,yz為客戶z實際的違約狀態。則[1]

式(9)表示當客戶j實際的違約狀態yj與客戶z實際的違約狀態yz相同時,yjz=1;否則,yjz=0。式(8)是兩個客戶相似概率pjz(w)的代數和,即式(8)為第j個客戶與其他全部客戶相似的概率pj(w)。概率pj(w)越大,則說明第j個客戶與其他全部客戶的違約狀態越一致,此時用其他全部客戶的違約狀態來判斷客戶j的違約狀態就越合適。式(8)與式(7)不同,式(7)是客戶j通過參考點z將自己的違約狀態判對的概率,式(8)是客戶j通過其他全部客戶將自己的違約狀態判對的概率pj(w),j=1,2,…,n。

步驟7目標函數F(w)的構建。

設F(w)為目標函數,n為客戶總數,pj(w)為客戶j通過除自身之外的n-1個客戶將自己的違約狀態判對的概率。λ為可調節參數,s為第1次指標篩選后剩余的指標個數,wi為第i個指標的權重。則[1]

步驟8λ=λ(1)時指標權重向量的確定。

式(10)第1項是式(8)的表達式,而式(8)的表達式是用式(7)計算的,式(7)的距離函數又是由式(6)表達的。由于式(6)是權重向量w的函數,故式(10)第1項是權重向量w的函數,式(10)第2項是w的顯函數。因此,式(10)中F(w)是權重向量w的函數。當λ為常數時,通過式(10)最大可以得到一組權重向量。

本文從0 開始,以1/n(n為訓練集中樣本個數,n=2 397)為步長選取50 個點作為λ的候選值。當λ=λ(1)時,λ(1)為常數,給定一個權重向量w1,1,得到一個目標函數F(w1,1)。給定第l個權重向量,得到第l個目標函數F(w1,l),l=1,2,…,100。比較這100個目標函數值F(w1,l),l=1,2,…,100。F(w)最大所對應的權重向量,即λ=λ(1)時的最優權重向量

步驟9λ=λ(k)時指標權重向量的確定。

在步驟8中已經得到λ=λ(1)時的最優權重向量,重復步驟8,可以得到λ=λ(k),k=2,3,…,50時的最優指標權重向量

步驟10第k個距離表達式d(k)的確定。

將步驟8、9中得到的50個指標權重向量

代入式(6),得到50 個距離表達式d(1),d(2),…,d(50),如下式所示:

步驟11第1個客戶違約狀態的確定。

取步驟10中式(11a)計算第1個客戶與任意客戶之間的馬氏距離,選取與第1個客戶距離最小的客戶的違約狀態作為第1個客戶的違約狀態的預測值,若與第1個客戶距離最近的客戶的違約狀態是非違約,則第1 個客戶的違約狀態的預測值=0;若與第1個客戶距離最近的客戶的違約狀態是違約,則。將第1個客戶的違約狀態預測結果列入步驟12式(12a)等式右端第1列。

同理,取步驟10 中其他49 個表達式(11b)~(11n),運用剩余49個距離表達式得到第1個客戶的違約狀態預測值,列入下文等式(12a)右端后49列。由此可得到第1個客戶在50個距離表達式d(1),d(2),…,d(50)下的違約狀態預測值,如下式所示。

步驟12第j個客戶違約狀態的確定。

仿照步驟11,可以得到第j個客戶(j=1,2,…,n)在50 個距離表達式下的違約狀態預測值。全部n個客戶在50個距離表達式下違約狀態的預測值為:

步驟13第1個距離表達式d(1)下違約判別誤差MSE(1)的確定。

式(12a)~(12n)可以構成一個n×50的矩陣,將矩陣中第1列n個客戶在第1個距離表達式下的違約狀態的預測值代入下式,得到違約判別誤差為[1]

矩陣中的第1列是在第1個距離表達式d(1)下得到的n個客戶的違約狀態預測值。因此,MSE(1)就是在第1個距離表達式d(1)下得到的全部客戶的違約判別誤差。即步驟8中取第1個指標權重向量時得到的全部客戶的違約判別誤差。

步驟14在第k個距離表達式下違約判別誤差MSE(k)的確定。

重復步驟13,依次將式(12)構成的n×50矩陣中的第2列,3列,…,50列代入式(13),得到在第k個距離表達式下的全部客戶的違約判別誤差MSE(k)(k=1,2,…,50),故可得到由50個MSE 構成的向量MSE=(MSE(1),MSE(2),…,MSE(50))。

步驟15最優指標權重向量w*的確定。

取向量MSE 中最小的一個違約判別誤差MSE*=min(MSE),最小的違約判別誤差MSE*對應的指標權重向量w*=就是最優的權重向量。

步驟16第1個指標組合D1及指標組合D1的違約預測精度AUC1的確定。

取τ1=0.2,則指標權重的臨界點[1]T1=τ1max(w*)=0.2max(w*),T1為常數,對比T1和步驟15中得到的每個指標的權重權重大于臨界點T1的重要指標即構成了一個指標組合D1。將訓練集中n個客戶對應指標組合D1中的指標數據xij和實際違約狀態yj代入下文式(14)~(16),構建線性支持向量機模型。對訓練集上客戶的違約狀態進行判別,得到客戶的違約狀態預測值,根據客戶違約狀態的實際值yj和預測值計算線性支持向量機的違約判別精度AUC1。

步驟17其他指標組合及違約預測精度的確定。

依次 取τ2=0.02,τ3=0.002,…,τ8=2×10-8(見表7),則指標權重的臨界點

重復步驟16,得到另外7 個指標組合D2,D3,…,D8及7 個指標組合的違約預測精度AUC2,AUC3,…,AUC8。

步驟18最優指標組合D*的確定。

比較步驟16、17中得到的8個AUC,AUC 最大所對應的指標組合即為最優的指標組合D*。以下文t-1年實證為例,AUC 最大對應的最優指標組合由表8前22行指標構成。

本文利用近鄰成分分析篩選指標與文獻[1]主要存在兩點差別:①研究問題不同。文獻[1]中研究的是腦電信號的處理,本文研究的是信用風險預測。本文首次將近鄰成分分析引入信用風險領域進行指標組合遴選。②指標重要性的判別標準不同。文獻[1]中取指標權重閾值T=τmax(w),τ=0.02作為指標重要性的判別標準,通過剔除指標權重小于閾值T的指標構建指標組合。本文將τ分別設置為0.2,0.02,0.002,…,2×10-8,得到8 個不同的指標組合,以指標組合的違約鑒別能力AUC 最大反推指標重要性的判別標準τ和最優指標組合。保證了指標組合整體的違約鑒別能力。

步驟4~8最優指標組合遴選的特色在于:首次將近鄰成分分析引入信用風險領域進行指標組合遴選,在近鄰成分分析算法中根據違約判別準確率最大得到馬氏距離中的指標權重向量,以基于馬氏距離的K-近鄰的違約判別誤差MSE 最小為目標,確定最優指標權重向量;給定一個指標權重閾值、并通過剔除權重小于閾值的指標得到一個指標組合,給定不同的指標權重閾值得到不同的指標組合,以指標組合的違約判別精度AUC 最大反推最優的指標組合。拓展了信用風險領域指標組合遴選的新思路。

2.2 違約預測模型中兩類客戶配比的選擇

本文在指標組合遴選中未對樣本比例進行處理,樣本比例的處理僅在指標組合遴選后建立違約預測模型時用到。

中國上市公司中違約樣本數量與非違約樣本數量的原始比例大概為1∶6.2,屬于非平衡樣本。本文以隨機欠采樣為基礎處理非平衡樣本。通常隨機欠采樣是從非違約客戶中隨機抽取與違約客戶數量相同的非違約客戶,使之構成1∶1的平衡樣本。對非違約客戶進行多次采樣,將全部違約客戶與采樣后的非違約客戶分別組成比例為1∶1、1∶2、1∶3、1∶4、1∶5和1∶6共6組采樣后的樣本。

以下文實證數據為例,將訓練集中全部333個違約客戶與訓練集中隨機抽取的333個非違約客戶組成樣本比例為1∶1的訓練子樣本。同理,將訓練樣本中全部333個違約客戶分別與在訓練樣本中隨機抽取的666,999,1 332,1 665,1 998個非違約客戶組成樣本比例為1∶2,1∶3,1∶4,1∶5,1∶6的訓練子樣本Q1,Q2,Q3,Q4,Q5,Q6。6組訓練樣本如表1所示。

表1 違約客戶與非違約客戶構成的6組訓練樣本

根據Q1,Q2,…,Q6共6組訓練樣本,采用3.3節的線性支持向量機模型建立6個違約預測模型,得到6個違約預測精度G-mean,以G-mean最大反推違約預測模型中兩類客戶的最佳配比。

本文處理非平衡樣本與現有研究的差別在于:現有研究[34,35]將非違約公司與違約公司的比率設定為1∶1研究企業信用風險。He等[2]根據經驗將違約客戶與非違約客戶的比例設置為1∶2.5,1∶5,1∶10。本文遍歷違約客戶與非違約客戶所有可能的樣本比例,以G-mean最大反推違約預測模型中兩類客戶的最佳配比。改變了隨機欠采樣中主觀確定采樣比例的不合理做法。

本文對非平衡樣本隨機欠采樣的處理特色是將違約客戶與非違約客戶組成不同比例的樣本,以基于線性支持向量機的違約預測精度G-mean最大為標準反推違約客戶與非違約客戶的最佳比例,以確定最優的訓練樣本。改變了現有研究中主觀設置欠采樣比例的做法。

2.3 基于線性支持向量機構建信用風險預測模型

有研究表明,線性支持向量機針對大型數據集(樣本數量大于等于2 000)具有較好的分類效果,而非線性支持向量機針對大型數據集的分類效果并不佳[35]。本文基于3 425家上市公司數據構建信用風險預測模型,選取線性支持向量機作為最終的信用判別模型。

設w為指標的權重向量,yj為客戶j的實際違約狀態(-1表示非違約,1表示違約),xj為客戶j的指標向量,對應步驟18中遴選出的最優指標組合,即表8中前22行指標構成的指標向量,b為截距,n為客戶總數。則支持向量機的目標函數為[35]:

通過引入拉格朗日乘子α,求解出指標權重向量w*和截距b*。

式(14)~(16)的經濟學含義:兩類客戶到超平面的距離之和為是兩類客戶之間距離的倒數,w越小,兩類客戶之間距離越大,違約判別精度越高。

本文構建的支持向量機與文獻[29,34]相比有兩個差別:①使用的指標組合是以違約判別誤差最小遴選出的最優指標組合;②構建線性支持向量機所使用的最佳樣本配比是根據違約判別精度Gmean最大反推得到的。

本文構建的線性支持向量機有兩個特色:①最優指標組合遴選的特色。在近鄰成分分析算法中根據違約判別準確率最大得到馬氏距離中的指標權重向量,以基于馬氏距離的K-近鄰的違約判別誤差MSE最小為目標,確定最優指標權重向量;給定一個指標權重閾值通過剔除權重小于閾值的指標得到一個指標組合,給定不同的指標權重閾值得到不同的指標組合,以指標組合的違約判別精度AUC 最大反推最優的指標組合。②最佳樣本配比的特色。利用隨機欠采樣將違約客戶與非違約客戶組成不同比例的樣本,以基于線性支持向量機的違約預測精度G-mean最大為標準,反推違約客戶與非違約客戶的最佳配比。

2.4 對比模型精度的檢驗標準

通常混淆矩陣可以用來作為違約判別模型分類效果的評價基礎。混淆矩陣如表2所示。

表2 違約判別混淆矩陣

采用精確度(Acc)、第1類錯誤(Type I error)、第2類錯誤(Type II error)、幾何平均值(G-mean)以及AUC 作為信用風險預測模型的精度檢驗標準[36]。各檢驗標準的定義及計算公式如下式所示。

準確率(Acc)是指違約客戶被判為違約的數量與非違約客戶被判為非違約的數量之和占總客戶數的比,即

第1類錯誤(Type I error)是指非違約客戶被判為違約的數量占非違約客戶總數的比,即

第2類錯誤(Type II error)是指違約客戶被判為非違約的數量占違約客戶總數的比,即

在利用不平衡數據進行客戶違約狀態判別時,通常利用G-mean作為評價判別模型好壞的標準,因為G-mean同時考慮了違約客戶和非違約客戶的判對率,即

違約客戶判對率(TPR)亦稱召回率(recall)、靈敏度(sensitivity),是指違約客戶被判為違約的數量占違約客戶總數的比,即

非違約客戶判對率(TNR)亦稱特異度(specificity),是指非違約客戶被判為非違約的數量占非違約客戶總數的比,即

以(1-特異度)為橫坐標,靈敏度為縱坐標,得到曲線ROC,ROC曲線與坐標軸圍成的面積即為AUC。通常利用AUC的值作為評價分類模型精度的標準,AUC的值越大,違約預測模型的預測效果越好[35]。

3 實證研究

3.1 樣本數據

3.1.1 原始樣本的來源

(1)指標海選。由于企業信用狀況涉及多方面因素,為使指標體系能全面反映企業現階段的信用狀況,從企業的財務狀況、非財務狀況以及宏觀經濟環境等方面海選出614個指標,其中包含342個公司財務指標[38],119個非財務指標[39-40],147個宏觀經濟指標和6個與貨幣發行量相關的指標[41-42]。初始指標體系涵蓋企業償債能力、盈利能力、營運能力、成長能力、非財務因素、企業高管基本情況、企業基本信用情況、商業信譽、社會責任以及外部宏觀因素等多個準則層。表3列舉了本文海選出的主要指標。

表3 上市公司信用風險預測初始指標體系

(2)原始樣本數據。以中國上市公司為研究對象,選取2000~2018年期間3 425家上市公司數據作為樣本,構建企業違約判別模型。所使用的數據來自Wind數據庫、國泰安數據庫、國家統計局。若上市公司連續兩年出現虧損,該公司就會被中國證監會特別處理(Specially Treated,ST)。以企業是否被標ST 為標準,將上市公司分為違約客戶和非違約客戶。其中,違約客戶476 家,非違約客戶2 949家。3 425家上市公司的原始數據uij列于表3第(1)~(3 425)列。

利用客戶第t-m(m=1,2,3,4,5)年的指標數據xij和第t年的違約狀態yj進行建模。對于違約客戶,將客戶違約當年作為第t年。對于非違約客戶,在不重復的情況下每年選取一定數量的客戶,將客戶對應的年份作為第t年。

3.1.2 指標數據標準化處理 標準化的目的是將指標數據轉化為[0,1]區間內的數值,消除數據的單位限制,便于不同單位間的數據進行加權。采用Chi等[43]的方法對原始指標uij進行標準化處理。將客戶標準化后的指標數據xij列于表3第(f)列。指標的標準化處理不是本文的主要工作,故不再贅述。

3.1.3 訓練集與測試集劃分 將樣本按照7∶3的方式劃分訓練集與測試集,如表4所示。

表4 樣本劃分

從非違約樣本與違約樣本中分別抽取70%的樣本組成訓練集Etrain,即訓練集中包含2 397個客戶,其中,違約客戶=333 個,非違約客戶2 064個。非違約樣本與違約樣本中分別剩余的30%樣本組成測試集Etest,即測試集中包含1 028個客戶,其中,違約客戶=885 個,非違約客戶=143個。

(1)指標遴選時使用的樣本。在指標遴選過程中僅使用表4列(1)的訓練樣本Etrain,不使用列(2)的測試樣本Etest。

(2)構建模型時使用的樣本。在表4列(1)訓練樣本中,采用表1中6個不同的比例,運用表3列(f)的數據xij和yj,建立式(14)~(16)的支持向量機模型。

(3)測試樣本。上文(2)使用6個樣本建模,由此建立了6個預測模型。對于6個模型中的每一個模型,都采用表4列(2)的測試樣本、表3列(f)的數據xij進行判別,由此得到理論判別狀態。通過將和測試樣本的實際違約狀態yj進行對比,得到了表2混淆矩陣的全部統計頻數,將統計頻數對應代入式(20),可以得到違約預測精度G-mean,如表11列(8)所示。

3.2 最優指標組合的確定

3.2.1 第1次指標篩選 依照步驟1~3 進行基于偏相關性分析的第1 次指標遴選。將表4 列(1)定義的訓練集樣本Etrain對應表3列(f)前614行的指標數據xij代入式(1)~(4)計算任意兩個客戶間的偏相關系數,代入式(5)計算每個指標的F統計量。挑選出具備經濟學含義且偏相關系數大于0.8的指標對,刪除指標對中經濟學含義不明顯或F統計量較小的指標。由此,614 個指標經過第1 次指標篩選剩余259 個指標,將剩余的259個指標列于表5 列(c)前259 行。表5 列(d)為訓練集Etrain中2 397個客戶的指標數據xij,列(e)為測試集Etest中1 028個客戶的指標數據xij。

3.2.2 第2次指標篩選 在第1次指標篩選剩余的259個指標的基礎上進行第2次指標篩選,目的是遴選出違約預測能力最強的指標組合。以t-1年為例(t-m年的其他年份類推,下同),依照步驟4~18,具體說明最優指標組合的遴選過程。

步驟19λ取值范圍的確定。

如表4列(1)第3 行所示,訓練集中樣本總數n=2 397,則1/n=1/2 397=0.000 417。根據步驟8,從0開始,以0.000 417為步長,取50個點作為λ的候選值,如表6列(2)所示。

步驟20任意兩個客戶間距離表達式的確定。

將表5列(1)前259行數據xi1和列(d)剩余的2 396 列中任意一列的前259 行數據xiz代入式(6),得到第1個客戶與其他2 396個客戶關于權重w的距離表達式dw(x1,xz),其中,z=2,3,…,2 397,共計2 396個距離表達式。同理,對于表5列(d)中任意一個客戶j,均可求出其與其他2 396個客戶的2 396個距離表達式dw(xj,xz)。由于表5列(d)中有2 397個客戶,故依照上述方法可以得到2 397×2 396=5 743 212個距離表達式。由于dw(xj,xz)=dw(xz,xj),故有2 871 606個距離表達式。

步驟21任意兩個客戶間相似度pjz(w)的確定。

將步驟20中得到的第1 個客戶與其他2 396個客戶的2 396個距離表達式dw(x1,xz),z=2,3,…,2 397代入式(7),得到客戶1與其他2 396個客戶關于權重w的相似概率函數p1z(w),共計2 396個相似概率函數。同理,對于表5列(d)中任意一個客戶j,j=1,2,…,2 397,均可求出其與其他2 396個客戶關于權重w的2 396 個相似概率函數pjz(w)。

步驟22客戶j與其余2 396個客戶相似度總和pj(w)的確定。

將步驟21中得到的2 396個p1z(w)代入式(8),同時將y1z也代入式(8),得到客戶1與其余2 396個客戶相似概率的總和p1(w)。關于y1z的取值如下:對比表5列(d)最后一行中第1個客戶與其他客戶的實際違約狀態y1和yz,若客戶1與客戶z實際違約狀態同時為0或同時為1,則y1z=1;否則,y1z=0。同理,對于表5列(d)中任意一個客戶j,均可求出其與其他2 396個客戶相似概率的總和pj(w)。

步驟23權重向量的確定。

將步驟22中的Pj(w)代入式(10),同時將表6第1行第(2)列數據λ=0也代入式(10)。由此,得到一個僅由w組成的目標函數F(w),以F(w)最大為目標,計算得到=(0.003 12,…,0.001 79),列于表6第1行第(3)列。

步驟24權重向量的確定。

重復步驟23,依次取表6第(2)列第2~50行的可調節參數λk,得到對應行的權重向量其中,k=2,3,…,50,列于表6第(3)列的對應行。

步驟25兩個客戶間距離的確定。

將表5列(d)中任意兩個客戶的前259行數據xij、xiz以及表6 中第1 行第(3)列數據(0.003 12,…,0.001 79)代入式(11a),得到當w==(0.003 12,…,0.001 79)時任意兩個客戶之間的距離同理,依次取表6第(3)列第2~50行的數據代入式(11b)~(11n),將表5列(d)中任意兩個客戶的前259行數據xij和xiz也代入式(11b)~(11n),得到當w=時任意兩個客戶之間的距離(xj,xz),其中,k=2,3,…,50。

表5 基于相關性分析的指標篩選結果

表6 不同λ 下的違約判別誤差

步驟26第1個客戶違約狀態的確定。

取步驟25中第1個客戶在第1個距離表達式下與其余2 396個客戶的距離(x1,xz),其中,z=2,3,…,2 396,共計2 396 個距離。比較這2 396個距離,選取最小距離所對應的客戶。若與第1個客戶距離最小的客戶的違約狀態為0,則第1個客戶的違約狀態的預測值=0;若與第1個客戶距離最小的客戶的違約狀態為1,則第1個客戶的違約狀態的預測值=1。由此,得到第1個客戶在第1個距離表達式下的違約狀態,列于下式矩陣A的第1行第1列。同理,取步驟25中第1個客戶在第k(k=2,3,…,50)個表達式下與其余2 396個客戶的距離,得到第1個客戶在第k個距離表達式下的違約狀態列于下式矩陣A的第1行后49列。

步驟27第j個客戶違約狀態的確定。

仿照步驟26,得到第j個客戶在第k個距離表達式下的違約狀態,列于下式矩陣A第j行第k列,即

步驟28第1個距離表達式d(1)下違約判別誤差MSE(1)的確定。

將式(23)中矩陣A第1列代入式(13),同時將表5 列(d)最后一行也代入式(13),得到違約判別誤差MSE(1)=0.137 216,列于表6第(4)列第1行。

步驟29第k個距離表達式下違約判別誤差MSE(k)的確定。

重復步驟28,依次取式(23)中矩陣A第2~50列,得到違約判別誤差MSE(2),…,MSE(50),列于表6第(4)列第2~50行。

步驟30最優指標權重向量w*的確定。

比較表6第(4)列中50個違約判別誤差MSE,選取最小的違約判別誤差min MSE=0.067 167及其對應的最優指標權重向量w*=(0.000 43,…,0.002 38)。

步驟31指標權重閾值的確定。

取指標重要性的判別標準τ1=0.2,τ2=0.02,…,τ8=2×10-8,如表7第(2)列所示。

表7 不同臨界點τ 下的違約預測精度

步驟32第1個指標組合D1及指標組合D1的違約預測精度AUC1的確定。

取表7 第(2)列第1 行數據τ1=0.2 和步驟30中得到的w*=(0.000 43,…,0.002 38),則指標權重閾值T1=τ1max(w)=0.001 48。將w*中每個指標的權重與T1=0.001 48比較,保留指標權重大于0.001 48的指標。這些指標即構成了一個指標組合D1。將表4第(1)列的訓練樣本對應指標組合D1中的xij和表5第(d)列最后一行的yj代入式(14)、(15),得到線性支持向量機中的指標權重向量w*和截距b*,再將w*和b*代入式(16)構建SVM1。將表4第(1)列的訓練樣本對應指標組合D1中的xij代入SVM1,得到客戶的違約狀態預測值,根據客戶違約狀態的實際值y與預測值計算得到線性支持向量機的違約判別精度AUC1=0.982 8,列于表7 第1 行第(3)列。

步驟33第q個指標組合Dq及指標組合Dq的違約預測精度AUCq的確定。

重復步驟32,依次取表7第(2)列第2~8行數據,并得到對應行的違約預測精度,列于表7第(3)列第2~8行。

步驟34最優指標組合D*的確定。

比較表7第(3)列8個數值,τ=2×10-7對應的AUC最大,為0.991 9,故τ=2×10-7對應的指標組合為t-1年最優的指標組合。t-1年最優的指標組合列于表8前22行。

表8 t-1年上市公司信用評價指標體系

同理,根據t-1年信用風險預測指標體系的構建流程,分別建立t-2年,t-3年,t-4年,t-5年的信用風險預測指標體系。其他年份下的最優指標組合遴選結果詳見附錄。

表9匯總了t-m(m=1,2,3,4,5)年最優指標組合中的指標,并統計了各個指標被選入最優指標組合的次數。

根據表9 列(3)可知,每股收益EPS-扣除/稀釋、貨幣供應量M0(億元)和貨幣供應量M1(億元)3個指標存在于t-1,t-2,t-3年的最優指標組合中,說明這3個指標對企業未來1~3年的短期違約狀態具有關鍵影響;當日總市值/負債總計、每股EBITDA 和固定資產周轉率3個指標存在于t-4,t-5年的最優指標組合中,說明這3個指標對企業未來4~5年的長期違約狀態具有關鍵影響;同時,經營活動產生的現金流量凈額/經營活動凈收益和審計意見類型2個指標存在于t-m(m=1,2,3,4,5)年的最優指標組合中,說明這2個指標不論對于企業未來1~3年的短期、還是未來4~5年的長期違約狀態,均有關鍵影響。

表9 t-m(m=1,2,3,4,5)年最優指標組合匯總

由此可以得出結論:每股收益EPS-扣除/稀釋、貨幣供應量M0(億元)和貨幣供應量M1(億元)3個指標對企業未來1~3年的短期違約狀態具有關鍵影響;當日總市值/負債總計、每股EBITDA 和固定資產周轉率3個指標對企業未來4~5年的長期違約狀態具有關鍵影響;經營活動產生的現金流量凈額/經營活動凈收益和審計意見類型2個指標,不論對于企業未來1~3年的短期、還是未來4~5年的長期違約狀態,均有關鍵影響。

3.3 基于不同比例的樣本構成

在訓練樣本Etrain中隨機抽取333個非違約客戶,將Etrain中全部333個違約客戶與抽取的333個非違約客戶組成比例為1∶1的訓練子樣本Q1。依照上述方法,將Etrain中全部333個違約客戶分別與在Etrain中隨機抽取的666,999,1 332,1 665,1 998個非違約客戶組成比例為1∶2,1∶3,1∶4,1∶5,1∶6的訓練子樣本Q2,Q3,Q4,Q5,Q6(見表1)。通過樣本Etrain,Q1,Q2,…,Q6構建7 個線性支持向量機,以線性支持向量機的G-mean最大反推違約客戶與非違約客戶的最佳樣本比。以訓練樣本Q2為例,具體樣本構成如表10所示。

3.4 最佳線性支持向量機的確定

3.4.1 線性支持向量機的構建 以預測期限m=1為例,利用訓練樣本Q2構建樣本比為1∶2的線性支持向量機。將表10中前22行的指標數據xij,第23行客戶的實際違約狀態yj代入式(14)、(15),估計得到線性支持向量機的指標權重向量w1,1:2和截距b1,1:2。將w1,1:2和b1,1:2代入式(16),得到預測期限m=1,違約與非違約客戶的樣本比為1∶2的線性支持向量機模型,如式(24)所示。同理,利用訓練樣本Q1,Q3,Q4,Q5,Q6,Etrain可以構建預測期限為m=1,違約與非違約客戶的樣本比為1∶1、1∶3、1∶4、1∶5、1∶6和原始比例的違約預測模型。當m=2,3,4,5時,利用訓練樣本Q1,Q3,Q4,Q5,Q6,Etrain構建違約與非違約客戶樣本比為1∶1、1∶2、1∶3、1∶4、1∶5、1∶6和原始比例的違約預測模型。共計35個線性支持向量機模型。

表10 t-1年訓練樣本Q 2

3.4.2 最佳樣本比的確定 將表4列(2)的測試樣本Etest對應表3列(f)的訓練樣本的數據xij代入t-1年構建的7個線性支持向量機,得到測試樣本中客戶違約狀態的預測值,將與實際值y對比,得到t-1年不同樣本比例下的違約預測精度,如表11第1~7行所示。同理,得到t-m(m=2,3,4,5)年下構建的其他28個模型的違約判別精度,如表11第8~35行所示。

本文以G-mean衡量違約預測模型的精度,因為G-mean同時考慮了違約客戶和非違約客戶的判對率。比較表11結果可以發現:當預測期限m=1,2,5,樣本比為1∶2時,G-mean為最大值;當m=3,4,樣本比為1∶1時,G-mean為最大值。由上述結果可以確定:當m=1,2,5時,違約客戶與非違約客戶的最佳樣本比例為1∶2;當m=3,4時,違約客戶與非違約客戶的最佳樣本比例為1∶1。通過表11還可以看出,在每一預測期限下,基于最佳樣本比構建的模型的第2 類錯誤(Type-II error)、召回率(recall)、G-mean以及AUC 均優于基于原始樣本構建的模型。

表11 基于不同樣本比例的預測結果

由此可得結論:利用違約樣本與非違約樣本的最佳樣本比能夠提高預測模型的預測能力。

根據不同年份下的最佳樣本比,可以確定不同年份下的最佳違約預測模型。t-m(m=1,2,3,4,5)年的最佳違約預測模型如式(24)~(28)所示。

式(24)~(28)中,sgn表示若方括號內的數值計算結果大于0,則sgn[]取“1”,表示第j個客戶違約;若小于0,則sgn[]取“-1”,表示第j個客戶非違約;若等于0,則sgn[]取“0”,此時不能識別第j個客戶的違約狀態。

3.5 模型精度的對比

3.5.1 指標篩選方法的對比 將本文的指標篩選方法與藍本文獻[1]中的指標篩選方法進行對比。藍本文獻將指標重要性的判別標準τ設置為0.02,本文將τ分別設置為0.2,0.02,0.002,…,2×10-8,得到8個不同的指標組合,以指標組合在訓練集中的違約鑒別能力AUC 最大反推指標重要性的判別標準τ和最優指標組合。

本文指標篩選方法記為NCA1,藍本文獻指標篩選方法記為NCA2。將對比結果列于表12。

由表12可以看出,在每一個預測期限下,利用本文的指標篩選方法(NCA1)構建最優指標組合,線性支持向量機的第2 類錯誤(Type-II error)、召回率(recall)、G-mean以及AUC 均為最優值。

表12 指標篩選方法的精度對比

由此可得結論:基于最優的指標權重閾值遴選指標組合能夠提高指標組合的違約鑒別能力。在t-m(m=1,2,3,4,5)年,基于近鄰成分分析構建的指標體系,使得模型預測精度AUC 達到0.9以上,說明本文將近鄰成分分析(NCA)方法引入信用風險領域進行指標組合遴選是有效的。

3.5.2 信用風險預測模型的對比 將本文構建的違約預測模型 (Neighborhood Component Analysis-Undersampling-Support Vector Machine,NCA-US-SVM)與其他經典違約預測模型如非線性支持向量機[44]、邏輯回歸(LR)[45]、決策樹(DT)[44]、K-近 鄰(KNN)[46]以及線性判別(LDA)[47]進行對比,結果列于表13。

通過表13可以看出,不同的預測期限下,本文構建的NCA-US-SVM 模型的第2類錯誤(Type-II error)、召回率(recall)均優于其他模型。當m=2,3,4,5時,NCA-US-SVM 模型的G-mean優于其他模型;當m=1,2,3,4 時,NCA-US-SVM 模型的AUC優于其他模型。總體而言,本文構建的違約預測模型的違約預測能力優于其他經典模型。

4 結論

本文主要結論如下:

(1)在鄰近成分分析的馬氏距離中,以違約判別誤差最小確定馬氏距離中的指標權重向量,通過淘汰權重低于臨界點的指標,得到一組AUC 最大的最優指標組合。通過對比不同預測期限下的最優指標組合可知:每股收益EPS-扣除/稀釋、貨幣供應量M0(億元)和貨幣供應量M1(億元)3個指標對企業未來1~3年的短期違約狀態具有關鍵影響;當日總市值/負債總計、每股EBITDA 和固定資產周轉率3個指標對企業未來4~5年的長期違約狀態具有關鍵影響;經營活動產生的現金流量凈額/經營活動凈收益和審計意見類型2個指標,不論對于企業未來1~3年的短期、還是未來4~5年的長期違約狀態,均有關鍵影響。

(2)利用隨機欠采樣方法將兩類客戶組成不同比例的樣本,以G-mean最大為標準得到不同年份下線性支持向量機的最佳樣本配比。以本文實證為例:對于中國上市公司,t-1,t-2,t-3年違約客戶與非違約客戶的最佳樣本比為1∶2;t-4,t-5年違約客戶與非違約客戶的最佳樣本比為1∶1。

本文主要創新點:

(1)在近鄰成分分析算法中根據違約判別準確率最大得到馬氏距離中的指標權重向量,以基于馬氏距離的K-近鄰的違約判別誤差MSE 最小為目標,確定最優指標權重向量;給定一個指標權重閾值、并通過剔除權重小于閾值的指標得到一個指標組合,給定不同的指標權重閾值得到不同的指標組合,以指標組合的違約判別精度AUC 最大反推最優的指標組合。拓展了信用風險領域指標組合遴選的新思路。

(2)利用隨機欠采樣將違約客戶與非違約客戶組成不同比例的樣本,以基于線性支持向量機的違約預測精度G-mean最大為標準反推違約客戶與非違約客戶的最佳比例,以確定最優的訓練樣本。

(3)通過t-m年的指標數據xt-m和t年的企業違約狀態yt,利用最優指標組合和最優訓練樣本建立了支持向量機模型,達到了運用t年的數據xt預測第t+m年企業違約狀態yt+m的預測效果。

(4)本文的違約預測模型的精度高于非線性SVM、LR、DT、KNN 和LDA 等典型的大數據預測模型。

附表1 t-2年上市公司信用評價指標體系

附表2 t-3年上市公司信用評價指標體系

附表3 t-4年上市公司信用評價指標體系

附表4 t-5年上市公司信用評價指標體系

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产成人精品男人的天堂| 欧美一级在线看| 91色老久久精品偷偷蜜臀| 欧美日韩亚洲国产| 欧美亚洲欧美| 欧美人人干| 日本成人一区| 青草91视频免费观看| 久久这里只精品国产99热8| 亚洲天堂精品视频| 自偷自拍三级全三级视频| 国产肉感大码AV无码| 国产视频一区二区在线观看| 国内99精品激情视频精品| 成人在线视频一区| 草草影院国产第一页| 国产青青操| 永久免费无码日韩视频| 最新加勒比隔壁人妻| 亚洲成AV人手机在线观看网站| 国产成人乱无码视频| 国产草草影院18成年视频| 国产91精品调教在线播放| 欧美成人一区午夜福利在线| 亚洲人成在线免费观看| 亚洲水蜜桃久久综合网站| 久久综合色88| 国产免费a级片| 国产成人精品亚洲77美色| 六月婷婷综合| 无码'专区第一页| 国产在线第二页| 亚洲无卡视频| 午夜国产理论| 无码不卡的中文字幕视频| 爆操波多野结衣| 青草免费在线观看| 91成人精品视频| julia中文字幕久久亚洲| 午夜爽爽视频| 国产麻豆91网在线看| 国产黄视频网站| 亚洲无码视频喷水| 亚洲精品爱草草视频在线| 国产va在线观看| 亚洲日韩精品综合在线一区二区| 国产精品短篇二区| 国产素人在线| 欧洲亚洲一区| 中文字幕亚洲综久久2021| 人与鲁专区| 黄片一区二区三区| 国产一二视频| 色婷婷色丁香| 大乳丰满人妻中文字幕日本| 国产一级妓女av网站| 狠狠色狠狠综合久久| 天天色综网| 在线国产资源| 91在线一9|永久视频在线| 日本久久网站| 久久久久亚洲AV成人网站软件| 奇米影视狠狠精品7777| 国产乱子伦精品视频| 一本大道AV人久久综合| 国产偷倩视频| 欧美精品伊人久久| 国产日韩精品一区在线不卡| 4虎影视国产在线观看精品| 美女国内精品自产拍在线播放| 老司机精品一区在线视频| 日本精品影院| 国产第一页第二页| 91青草视频| 日韩在线播放中文字幕| 国产成年女人特黄特色毛片免| 日韩在线视频网| 欧美日韩专区| 国产精品免费p区| 99热国产在线精品99| 成人免费网站在线观看| 国产人成在线视频|