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

基于極值理論的碳價格波動風險管理研究

2022-03-04 02:19:38王喜平閆麗娜
電力科學與工程 2022年2期
關(guān)鍵詞:模型

王喜平,閆麗娜

(華北電力大學 經(jīng)濟管理學院,河北 保定 071003)

0 引言

減少溫室氣體排放,盡早實現(xiàn)“碳達峰”“碳中和”已成為全球共識。碳交易作為一種可以有效控制溫室氣體排放的市場調(diào)節(jié)手段,已成為我國落實減排承諾、實現(xiàn)減排目標的重要政策工具。自2011 年,國內(nèi)多個省、市開展了碳交易試點。從試點市場的運行結(jié)果看,減排效應(yīng)日漸顯現(xiàn):截至2020 年底,全國8 個碳交易試點單位累計完成成交量超過4 億噸,累計成交額超過90 億元,單位GDP(gross domestic product)二氧化碳排放較2005 年降低約48.4%,提前超額完成下降40%~45%的目標[1]。由于受氣候、政策、能源價格以及供需機制等多種因素的影響,市場碳交易價格波動劇烈。目前,我國碳市場存在著發(fā)展時間短、政策框架不完善、金融化程度不足、碳市場作用發(fā)揮不充分等問題,這些問題造成碳市場的交易價格波動存在不確定性[2];因此,有必要對碳價格波動風險進行研究,以有效避免碳價格風險爆發(fā)對碳市場運行帶來不利影響。

近年,碳市場價格波動及風險測度成為學界關(guān)注的焦點。對于碳價格的波動,多數(shù)學者都通過建立GARCH 族模型來進行分析,并得到結(jié)論:碳價格具有尖峰厚尾、波動性聚集等特征[3-6]。有學者從碳價的影響因素入手,探討其對碳價波動的作用。文獻[7]利用向量誤差修正模型對核證減排量現(xiàn)貨價格進行了分析,認為宏觀因素和氣候因素對碳價影響最為顯著。文獻[8-13]分別運用GARCH 模型、向量自回歸模型、多尺度熵等方法研究了能源價格與碳排放權(quán)價格之間的關(guān)系,認為碳價格主要受天然氣、石油價格的影響。

在此基礎(chǔ)上,諸多學者利用不同的方法對碳價收益率的在險價值進行了度量。文獻[14]運用極值理論,對歐盟碳排放權(quán)交易的現(xiàn)貨市場和期貨市場風險進行了測度,認為運用極值理論估計VaR 比傳統(tǒng)方法更加有效。文獻[15]將經(jīng)驗?zāi)B(tài)分解算法和ARMA-EGARCH 結(jié)合,構(gòu)建了多時間尺度的VaR 模型;結(jié)果表明用該方法能夠有效降低極端事件的影響,獲得更為準確的碳市場風險值。文獻[16]用經(jīng)濟狀態(tài)依賴方法估計歐盟碳市場的風險,認為該方法在樣本外VaR 預(yù)測中的結(jié)果優(yōu)于傳統(tǒng)方法結(jié)果。

以上研究主要圍繞歐盟碳市場展開。對于我國碳交易試點市場的風險測度,文獻[17]對我國碳試點碳價收益率構(gòu)建 ARCH(autoregressive conditional heteroscedasticity)族模型并計算了各交易試點的風險價值VaR,結(jié)果顯示:不同區(qū)域的碳排放交易市場,其波動程度不盡相同,計算的 VaR 有效性也有較大差異。文獻[18]利用GARCH(generalized autoregressive conditional heteroscedasticity)模型估計和預(yù)測北京碳市場的交易價格,認為模型在95%的置信水平下能夠準確度量碳價風險。文獻[19]認為QAR-GARCH(quantile autoregression-generalized autoregressive conditional heteroscedasticity)模型比 CAVaR(conditional autoregressive value at risk)族模型更加適用于對我國碳市場風險的測度。

我國在碳交易方面剛剛起步,存在數(shù)據(jù)量小的問題;所以,關(guān)于這方面的研究主要集中在碳價格的影響因素、碳價格波動特征等方面。在為數(shù)不多的對我國碳價風險的研究中,尚未發(fā)現(xiàn)有文獻將極值理論運用到風險價值度量中。極值理論更注重對尾部風險的刻畫,提供了一種更加準確的VaR 估計方法。目前,極值理論估計VaR 已經(jīng)被廣泛應(yīng)用于金融、保險、原油市場的風險度量中,并且都得到了較為準確的估計結(jié)果[20-23]。碳價收益率的波動也具有明顯的厚尾特征,受極端事件的影響較大。我國從2013 年建立碳交易試點至今,交易時間已有8 年,交易價格和交易量的相關(guān)數(shù)據(jù)也有了一定的積累,所以可以以此為基礎(chǔ)來嘗試碳價收益率波動風險方面的研究。

鑒于此,本文在借鑒已有研究的基礎(chǔ)上,以碳交易試點市場——深圳相關(guān)數(shù)據(jù)為研究背景,將條件方差和極值理論結(jié)合,度量我國碳價收益率的波動風險,為碳價波動風險管理提供參考借鑒。

1 模型構(gòu)建及原理

1.1 VaR

在金融市場中,在險價值VaR 被廣泛應(yīng)用于對市場風險的衡量。它是指在一定的置信度和給定的時間期間內(nèi),由于價格變化引起的超過目標水平的最大損失。該指標最大的優(yōu)點在于,可用單一的數(shù)值對一定時間間隔內(nèi)的風險進行總結(jié)。對任一收益率序列的分布函數(shù)F(X),當置信水平為p時,如果用V表示該序列VaR,則有:

ES(expected shortfall)是指當投資組合的損失超過VaR 時所遭受的平均損失程度。VaR 難以有效反映序列尾部的小概率事件發(fā)生時的情況,而ES 很好地彌補了這一不足。

若將ES 用E表示,那么該收益率序列的ES可表示為:

1.2 ARMA-GARCH

以往的研究顯示,我國碳價收益率存在明顯的波動聚集性和異方差性;因此,可以用ARMA(p,q)-GARCH(m,n)來捕捉其波動特征。

ARMA 模型通常用來描述價格變化的歷史規(guī)律,AR(p)為自回歸部分,MA(q)為移動平均過程。

GARCH 模型又稱為異方差模型,它由ARCH模型演化而來,用于描述時間序列的異方差性。GARCH 模型由均值方程和條件方差方程進行刻畫。

將ARMA 過程加入到GARCH 模型的均值方程中,以更好地挖掘時間序列的殘差信息,使得殘差序列更加服從獨立同分布的假設(shè),從而更好地刻畫時間序列的波動性。模型的具體表達式如下:

式中:均值方程中,rt表示碳價格的對數(shù)收益率;εt表示殘差。條件方差方程中,為殘差εt的條件方差;It–1表示t時刻以前所有信息集;ηt表示均值為0、方差為1 的獨立同分布時間序列。

1.3 極值理論

金融時間序列往往具有肥尾分布的特征,通常不滿足正態(tài)分布的假定。根據(jù)極值理論,x值超過臨界值μ時的累積分布函數(shù)為:

這種分布被稱為 GPD(generalized pareto distribution),其中:y=(x-μ)/β,β>0,β為比例參數(shù);ξ為形狀參數(shù);Nμ表示序列中超過臨界值μ的數(shù)目;N為總體數(shù)目。

在極值理論中,衡量尾部風險有2 種常見的模型:一種是超越閾值的極值模型POT(peak over threshold);另一種為區(qū)間選取極值模型BMM(block maxima group of models),這種模型更適合對具有一定周期性或季節(jié)性的數(shù)據(jù)進行建模。

考慮到本文所研究的碳市場相關(guān)數(shù)據(jù)背景中,碳價收益率的波動不存在隨季節(jié)變化的特征;因此,本文更適合采用POT 對序列尾部進行建模。

經(jīng)過推導(dǎo)分別得到VaR 和ES 的求解公式:

用POT 模型計算VaR 時,首先要通過GPD分布對極值數(shù)據(jù)進行擬合,其次可通過圖像觀察安全閾值。一般通過平均超額函數(shù)圖法和Hill 圖法來確定閾值μ。

1.3.1 平均超額函數(shù)圖法

利用平均超額函數(shù)圖作為確定閾值μ的依據(jù)。定義隨機變量X超過閾值的平均超額值e(μ):

由于GPD 的平均超額函數(shù)是閾值的線性函數(shù),所以可以通過尋找超過某一觀測值后圖像趨于線性的這一點來確定閾值μ。

1.3.2 Hill 圖法

另一個可以確定閾值的工具是Hill 圖法。定義Hill 統(tǒng)計量Hk:

式中:k為樣本容量,k=1,2,···,n–1,n;次序統(tǒng)計量X(i)滿足X(i)>X(i–1)。

用k表示橫坐標軸,繪制Hk的函數(shù)圖像;選取圖像中Hk趨于穩(wěn)定時所對應(yīng)的k,將X(k)確定為閾值μ。

由于上述2 種方法均是通過觀察圖像來確定閾值,故存在一定的主觀性。DuMouchel(1983)提出的一個較為簡單的量化方法:認為應(yīng)該選取序列中超過臨界值μ的數(shù)目占總體樣本的90%時所對應(yīng)的μ作為閾值。本文將結(jié)合平均超額函數(shù)圖法與Hill 圖法,并以90%時所對應(yīng)的閾值作為參考共同確定最優(yōu)的閾值μ。

1.4 動態(tài)VaR 及回測檢驗

VaR 的測度分為靜態(tài)和動態(tài)2 種方法。前文給出了靜態(tài)VaR 的計算公式。靜態(tài)VaR 值相對更加直觀。在預(yù)測收益率風險的時候,動態(tài)VaR 會與收益率序列同時發(fā)生波動,這樣可以更好地擬合序列的波動性。用Vt表示動態(tài)的VaR,則其計算公式為:

最后,需要對VaR 的有效性和準確性進行檢驗。常用的檢驗方法是Kupiec 回測檢驗。基于伯努利分布構(gòu)造LR 統(tǒng)計量,用L表示LR 統(tǒng)計量,則有:

式中:N為失敗的時長;T為實際考察的時長;N/T為失敗率;p為給定的置信水平。該假設(shè)檢驗的原假設(shè)為:p=N/T。

可通過觀察實際的失敗率與估計的失敗率是否相近,來確定模型的優(yōu)劣。如果原假設(shè)成立,則認為似然估計量LR 近似服從自由度為1 的卡方分布。因此,當置信水平給定時,就可以通過LR 的值判斷是否拒絕原假設(shè)。

2 實例分析

深圳碳交易試點是全國首個正式運行的交易試點。試點交易品種包括SZA-2013、SZA-2014等。交易主體主要為交易所內(nèi)規(guī)定的管控單位、其他機構(gòu)投資者和個人投資者,涉及電力、制造業(yè)、燃氣等多個行業(yè)。

自深圳碳試點成立以來,碳價格發(fā)生了較大幅度的變化。圖1 顯示了深圳碳配額成交均價。從圖中可以看出,深圳碳交易試點成立之初,交易價格偏高;這可能是由于碳市場起步階段受到國家政策的積極引導(dǎo),激發(fā)了投資者參與碳市場交易的熱情。2015—2017 年,價格逐漸趨于平穩(wěn)波動的狀態(tài),并且表現(xiàn)出明顯的周期性變化趨勢。2018—2019 年間,價格一度陷入低迷。

圖1 2013—2019 年深圳碳配額成交均價走勢圖Fig.1 Trend chart of average transaction price of carbon quota in Shenzhen from 2013 to 2019

2.1 數(shù)據(jù)來源

本文選取了2013 年12 月26 日至2020 年10月9 日的日度交易數(shù)據(jù)。深圳碳交易產(chǎn)品種類較多,按交易量加權(quán)平均后得到平均交易價格,然后去除交易量為0 的數(shù)據(jù)后,得到共計1 824 個數(shù)據(jù)。為了滿足數(shù)據(jù)的平穩(wěn)性要求,將碳交易價格取對數(shù)收益率,則第t個交易日的碳價格收益率rt可表示為:

數(shù)據(jù)經(jīng)處理后,得到圖2。由圖2 可見,碳價收益率均值為–0.000 477,方差為0.267 706,偏度為0.200 868,峰度為14.052 70。正的偏度值表明碳價收益率序列存在右偏概率分布性質(zhì);較大的峰度值表明碳價收益率序列存在尖峰厚尾的特征。JB 統(tǒng)計量為9 291.448,該結(jié)果顯然拒絕了收益率序列服從正態(tài)分布的假設(shè)。

圖2 深圳碳價收益率的描述性檢驗結(jié)果Fig.2 Descriptive test results of Shenzhen carbon price yield

2.2 ARMA-GARCH 模型的選擇

2.2.1 序列的平穩(wěn)性檢驗

首先,需要對碳價收益率序列進行平穩(wěn)性檢驗。根據(jù)表1 中收益率序列平穩(wěn)性檢驗結(jié)果可知,深圳碳價格收益率序列的T 統(tǒng)計量小于99%置信度時對應(yīng)的臨界值,因此拒絕序列存在單位根的假設(shè),故認為序列是平穩(wěn)的。

表1 收益率序列的平穩(wěn)性檢驗結(jié)果Tab.1 Stationarity test results of the rate of return series

2.2.2 ARMA-GARCH 模型階數(shù)的確定

在收益率平穩(wěn)的基礎(chǔ)上,需要對碳價收益率序列的自相關(guān)性進行擬合。大量關(guān)于風險度量的研究發(fā)現(xiàn),一階或二階的ARMA-GARCH 模型能夠較為充分地刻畫時間序列的特征。結(jié)合AIC(akaike information criterion)和SC(schwarz criterion)原則,可以基本確定最優(yōu)模型。通過比較擬合結(jié)果,在系數(shù)均顯著的條件下,認為ARMA(2,0)模型能夠更加充分地捕捉深圳碳價收益率的自相關(guān)特征。

對模型擬合后的殘差序列進行ARCH 效應(yīng)檢驗,檢驗結(jié)果如表2 所示。序列在99%的置信度下拒絕原假設(shè),認為序列確實存在ARCH 效應(yīng)。因此,需要進一步引入GARCH 模型。

表2 殘差序列的ARCH 效應(yīng)檢驗結(jié)果Tab.2 ARCH effect test results of residual series

表3 示出了不同階數(shù)下ARMA-GARCH 模型的擬合結(jié)果。通過對比AIC 和SC 的值,發(fā)現(xiàn)最好的擬合模型是ARMA(2,0)-GARCH(1,1)模型。

表3 不同階數(shù)ARMA-GARCH 模型擬合結(jié)果Tab.3 Fitting results of ARMA-GARCH models with different orders

表4 示出了ARMA(2,0)-GARCH(1,1)模型的擬合結(jié)果。從表4 的估計結(jié)果可以看出,AR(1)和AR(2)的系數(shù)均通過了顯著性檢驗,說明深圳碳價收益率序列存在自相關(guān)性。ARCH(1)和GARCH(1)系數(shù)均顯著且大于零,滿足系數(shù)和小于1 的條件。系數(shù)和為0.99,接近于1,表明收益率序列存在波動性集聚的現(xiàn)象,即:前一時期的收益率波動會影響下一期的收益率波動。當市場受到?jīng)_擊時,波動雖然會逐漸減弱,但持續(xù)時間較長。另外,由于ARCH 項的系數(shù)明顯小于GARCH項系數(shù),說明收益率波動受到自身記憶性的影響要強于市場沖擊所帶來的影響。

表4 ARMA(2,0)-GARCH(1,1)模型的擬合結(jié)果Tab.4 Fitting results of ARMA(2,0)-GARCH(1,1) model

對ARMA(2,0)-GARCH(1,1)模型擬合的殘差序列進行ARCH-LM 檢驗,結(jié)果如表5 所示。從表5 可以看出,F(xiàn) 統(tǒng)計量的p值為0.1586,明顯大于0.05 的顯著性水平。因此,可以認為ARMA(2,0)-GARCH(1,1)基本消除了碳價收益率序列的異方差性。根據(jù)估計結(jié)果可以得到ARMA(2,0)-GARCH(1,1)模型的表達式為:

表5 殘差項的ARCH-LM 檢驗結(jié)果Tab.5 ARCH-LM test results of residuals

2.3 標準殘差序列的GPD 擬合

本文通過R 語言的evir 工具包,運用極值理論中的POT 模型對ARMA(2,0)-GARCH(1,1)模型得到的標準殘差序列進行GDP 分布的擬合。由于POT 模型擬合時需要數(shù)據(jù)均為正值,本文將深圳碳價收益率的標準殘差序列分為正、負2 部分,分別進行擬合。對于殘差正序列,首先通過繪制平均超額函數(shù)圖和Hill 圖確定合適的閾值。通過觀察圖3可以看到:平均超額函數(shù)在2 之后,數(shù)據(jù)變得稀松;在1.5~2 之間存在一定的線性關(guān)系。因此,可以初步判定閾值應(yīng)該大于1.5。結(jié)合圖4 中Hill 圖可以看出,圖像在殘差序列的第76 個之后逐漸平穩(wěn),對應(yīng)的閾值在應(yīng)小于1.9。另外根據(jù)DuMouchel 提出的量化方法,認為超出閾值的數(shù)目應(yīng)該小于10%。因此,最終確定的閾值為1.811 015。

圖3 上尾的平均超額函數(shù)圖Fig.3 Average excess function of the upper tail

圖4 上尾的Hill 圖Fig.4 Hill of the upper tail

在確定合適的閾值后,利用極大似然估計得到殘差序列尾部擬合GPD 分布的參數(shù)ξ值為0.004 850 581,β值為0.608 867 825。

由繪制出的圖5 超額分布和和圖6 擬合的GPD 分布,可以看出GPD 曲線對數(shù)據(jù)進行了非常好的擬合。

圖5 上尾的超值分布圖Fig.5 Premium distribution chart of the upper tail

圖6 上尾的內(nèi)在分布尾部圖Fig.6 Intrinsic distribution of the upper tail

在圖7 上尾的殘差散點圖和圖8 的殘差Q-Q圖中,殘差均勻散落,極少數(shù)點沒有落在直線上,基本滿足GPD 分布的假設(shè)條件。可見,標準殘差上尾對GPD 分布的擬合程度較好。

圖7 上尾的殘差散點圖Fig.7 Residual scatter plot of the upper tail

圖8 上尾的殘差Q-Q 圖Fig.8 The residual Q-Q plot of the upper tail

根據(jù)靜態(tài)VaR 的計算公式計算結(jié)果可以得到:在99%的置信水平下,深圳碳價收益率的上尾靜態(tài)VaR 值為 3.214 282;這意味著有99%的可能性,由市場價格變化帶來的碳價收益率變化不會超過3.214 282。ES 值為 3.832 957,99%的分位數(shù)為3.214 282;這意味著如果3.214 282 的收益率被超過了,則預(yù)期損失為3.832 957。

圖9 示出了99%置信水平下的VAR 與ES。第一條垂直虛線和對應(yīng)的輪廓似然曲線顯示了估計的VaR,第二條垂直虛直線和對應(yīng)的輪廓似然虛線顯示了估計的ES。圖中,輪廓似然虛線與水平虛直線的2 個交點分別是估計出的置信區(qū)間邊界。

圖9 上尾99%的VaR 與ES 估計Fig.9 VaR and ES estimation at 99%confidence level of the upper tail

用同樣的方法對負的標準殘差序列進行GPD尾部擬合,結(jié)果如圖10 所示。從圖10 可以看出,平均超額函數(shù)圖中橫坐標值大于2 之后,數(shù)據(jù)變得稀松,在1.5~2 之間存在一定的線性關(guān)系:由此可以初步判定閾值應(yīng)該大于1.5。結(jié)合圖11 中的Hill 圖可以看出,圖像曲線在殘差序列第86 個之后逐漸平穩(wěn),對應(yīng)的閾值應(yīng)小于1.8,最終確定的閾值為–1.611 475。通過極大似然估計,得到參數(shù)ξ=0.008 154 222,β=0.525 983 924。從輸出圖12 可以看出,下尾收益率序列的擬合效果也比較好。根據(jù)靜態(tài)VaR 的計算公式,可以得到:碳價收益率的下尾靜態(tài)VaR 值為–2.826 280;這表明有99%的可能性,由市場價格變化帶來的碳價收益率變化不會低于–2.826 280。ES 值為3.366 575,99%的分位數(shù)為–2.826 280;這表示如果超過了–2.826 280 的收益率,則預(yù)期損失為–3.366 575。

圖10 下尾的平均超額函數(shù)Fig.10 Average excess function of the lower tail

圖11 下尾的Hill 圖Fig.11 Hill of the lower tail

圖12 下尾的分布輸出值Fig.12 Distribution output of the lower tail

圖13 示出了99%置信水平下下尾的VaR 與ESFig.13 VaR and ES estimation at 99% confidence level of the lower tail

3 回測檢驗

由于動態(tài)VaR 可以更好地擬合收益率的波動性,因此將ARMA(2,0)-GARCH(1,1)得到的殘差序列的條件均值和條件方差代入動態(tài)VaR 的計算公式,可以得到動態(tài)VaR 序列。

圖14 為上漲的動態(tài)VaR,圖15 為下跌的動態(tài)VaR。無論是上漲還是下跌,動態(tài)VaR 都很好地包絡(luò)了收益率曲線。在VaR 的回測檢驗中發(fā)現(xiàn):在上漲的動態(tài)VaR 中,共有847 個觀測值,VaR失效的個數(shù)共有4 個,失敗率為0.004 72;下跌的動態(tài)VaR 中,共有974 個觀測值,VaR 失效的個數(shù)共有9 個,失敗率為0.009 24。

圖14 上尾的動態(tài)VaR 與收益率Fig.14 Dynamic VaR and rate of return of the upper tail

圖15 下尾的動態(tài)VaR 與收益率Fig.15 Dynamic VaR and rate of return of the lower tail

通過計算99%置信水平下的LR 統(tǒng)計量,得到上漲的LR 值為1.286 336,下跌的LR 值為0.025 306,均小于99%置信水平下的臨界值6.635:因此,認為該模型對深圳碳價收益率VaR的估計是有效的。

4 結(jié)論

本文以碳排放權(quán)交易價格為研究對象,在考慮了極端情況碳價波動風險的條件下,將條件方差和極值理論結(jié)合起來,對碳價收益率進行了建模。對比結(jié)果顯示,ARMA-GARCH-EVT 模型能夠更加準確地衡量極端條件下的碳價波動風險。回測檢驗結(jié)果表明,模型的檢驗失敗率較低,所以認為模型在99%的置信水平下是有效的。

從得到的VaR 估計結(jié)果來看,深圳碳市場的價格波動較大,雖然下尾的靜態(tài)VaR 值小于上尾的靜態(tài)VaR 值,但下尾的動態(tài)VaR 值在2019 年4 月達到–3.216 51。當投資者面臨這種風險時,會遭受巨大的損失,引起整個碳市場的波動。因此,建議對碳市場進行必要的風險管控。政府應(yīng)在穩(wěn)定碳價方面發(fā)揮重要作用,建立相應(yīng)的市場穩(wěn)定儲備機制。金融機構(gòu)也應(yīng)積極參與碳市場建設(shè),為碳市場設(shè)計更加多樣的風險控制工具,從而幫助控排企業(yè)更好地管理風險。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 91成人在线免费观看| 国产成人精品免费av| 欧美国产精品不卡在线观看| 韩日午夜在线资源一区二区| 亚洲成人一区在线| 国产精品久久久久久久伊一| 亚洲精品日产精品乱码不卡| 国产精品妖精视频| 四虎成人在线视频| 国产精品白浆在线播放| 亚洲av色吊丝无码| 欧美日本在线观看| 亚洲中文无码h在线观看| AV熟女乱| 国产午夜在线观看视频| 国内黄色精品| 亚洲第一色视频| 激情视频综合网| 国产女人在线| 色香蕉影院| 欧洲高清无码在线| 国产成人无码播放| 成年人国产网站| 在线欧美一区| 麻豆精品久久久久久久99蜜桃| 国产精品一区二区不卡的视频| 久久毛片免费基地| 国内精品久久九九国产精品| 狠狠ⅴ日韩v欧美v天堂| 欧美中文字幕在线视频| 亚洲v日韩v欧美在线观看| 一区二区三区国产| 中文字幕在线播放不卡| 精品一区二区三区水蜜桃| 丰满人妻一区二区三区视频| 在线观看视频一区二区| 国产日本一线在线观看免费| 波多野结衣第一页| 久久精品国产电影| 国产va在线| 小说 亚洲 无码 精品| 亚洲娇小与黑人巨大交| 91久久夜色精品国产网站| 国产视频只有无码精品| 亚洲网综合| 亚洲av无码成人专区| 免费三A级毛片视频| 日本一本正道综合久久dvd | 国产精品国产三级国产专业不| 色婷婷在线影院| 国产91蝌蚪窝| 性视频一区| 制服丝袜一区| 亚洲三级成人| 国产精品成人AⅤ在线一二三四| 亚洲欧美人成人让影院| 亚洲天堂在线免费| 91福利片| 亚洲AⅤ波多系列中文字幕| 亚洲资源站av无码网址| 免费毛片视频| 国产一级二级在线观看| 久久国产精品嫖妓| 精品视频在线一区| 久久夜色撩人精品国产| 国产成人综合亚洲欧美在| 国产97视频在线观看| 国产正在播放| 成人综合网址| 夜色爽爽影院18禁妓女影院| 日韩在线视频网| 美臀人妻中出中文字幕在线| 久久99精品久久久久久不卡| 在线无码九区| 色亚洲成人| 97人妻精品专区久久久久| 国产精品一区二区在线播放| 一本色道久久88亚洲综合| 日本一区高清| av在线5g无码天天| 99热这里只有精品久久免费| 国产成人成人一区二区|