混合Beta分布GARCH模型的EM算法求解與實證分析

石 凱，劉洪江，孫 峰

（樂山師范學院a.數理學院；b.旅游學院，四川 樂山 614000）

0 引言

時間序列分析是考察經驗樣本數據隨時間演變的學科，由于時序樣本的走勢與波動都包含著不確定的隨機因素，因此統計的理論和方法在其中起著重要作用。單變量時序數據信息的提煉主要集中在均值的平穩性和波動的方差齊性兩個方面。關于時序均值平穩性信息的提煉可以依據差分整合自回歸移動平均（ARIMA）模型，但是ARIMA模型設定的前提條件要求方差齊性。如果時序數據資料中存在違背方差齊性的特征，則還需要進一步提取異方差信息。異方差所呈現的類型也有很多，在金融市場領域中廣泛具有的是一種條件異方差特性，即在波動偏大的時段持續偏大，波動偏小的時段持續偏小，也稱為集群效應。集群效應處理對應的方法是自回歸條件異方差（ARCH）模型，最初是由Engle（1982）[1]在分析英國通貨膨脹率序列時提出的對殘差平方進行的一種自回歸建模方法。在實際應用中，ARCH模型擬合異方差函數往往會產生較高的移動平均階數，為克服這一缺陷，Bollerslev（1986）[2]提出廣義自回歸條件異方差，即GARCH 模型。GARCH模型的提出，有效解決了異方差的長記憶性問題，成為現今普遍采用的處理時序數據條件異方差的方法。

為了拓展GARCH 模型的應用范圍，提升時序波動特征的擬合和預測精度，學者們從不同角度出發，在GARCH模型的基礎上衍生出一系列模型。然而，需要指出的是，經典GARCH族類模型中均對模型作了正態分布的假定，可是在處理一些特殊領域，尤其是金融市場時，這一假定存在嚴重的問題。Hsieh（1989）[3]在對金融市場匯率的相關研究中指出，可以用GARCH（1,1）模型解釋匯率波動的絕大部分信息，但是正態分布設定卻不能捕獲金融數據特有的“高峰厚尾”特征（即在均值附近的樣本點比正態分布多，取極端值的樣本點也比正態分布多）。此后，關于金融市場的一些特殊特征一直是學術界研究的熱點之一，如Clark 和Baccar（2018）[4]在研究信用利差波動率的變化情況時證實，利差數據除了高峰厚尾外，還具有非對稱性等特征，與傳統GARCH模型的正態性假定嚴重不符。同樣，國內金融市場也不例外，徐龍炳（2001）[5]、別曉芳（2018）[6]等均指出中國金融市場也不服從正態分布的特性；同時，白仲林等（2011）[7]還認為受我國當前經濟發展水平和金融制度的限制，金融市場的資產收益率等序列還具有“有界取值區間”的特征。

針對此類特性的刻畫，眾多文獻指出要采用非正態分布來彌補模型設定的缺陷[8—10]。這些文獻擴展了模型的應用范圍，重點提取了金融市場的高峰厚尾、非對稱性等信息，但是由于金融市場，尤其是國內金融市場的一些特殊制度和交易規則，使得時序數據的波動存在有界的取值區間。而正態分布、t 分布、GED 分布、混合高斯分布等分布類型的變量取值范圍都是整個實數域，嚴重制約了GARCH 族模型在現實環境中的應用。因此，為進一步拓展GARCH模型的應用領域，本文試著引入變量取值為有界區間（0,1）的Beta分布，并構建混合Beta分布的GARCH模型，以期能有效提取金融市場波動的復雜特性。同時，為考察混合Beta 分布對GARCH模型刻畫能力的改進，本文給出了模型設立原理與求解過程，并通過模擬仿真數據和現實金融數據，對參數估計效果進行了比較和檢驗。

1 混合Beta分布GARCH模型的構建

1.1 GARCH模型的一般形式

GARCH 族模型是一種自回歸條件異方差模型，實質是將歷史波動信息作為條件，并采用自回歸形式來刻畫波動的演變規律，基本GARCH模型的設定形式如下：

其中，yt是可觀測的時序數列。式（1）是均值方程，用以刻畫均值的信息；式（2）是對式（1）殘差項的異方差進行處理，使其轉化為一個白噪聲序列；波動信息的提煉主要體現在式（3）上，也是GARCH模型的核心，使用自回歸的方式提取時序數列波動中蘊含的自相關信息。經典GARCH 模型中假定白噪聲序列et是服從均值為0、方差為1的正態分布，即(0, 1)。則由GARCH模型的結構可以得出εt|(εt-1,εt-2,…;ht-1,ht-2,…)~N(0,ht)，以及可觀測序列yt的條件分布為N(w,ht)。由此可見，經典GARCH 模型正態分布的假定限制性太強，若序列存在高峰厚尾、非對稱、波動有界等特征，則難以進行有效的信息獲取，為此，需要放寬模型分布的設定假設，以便得出更符合現實需求的研究結果。

1.2 混合Beta分布及其性質

在概率統計中，Beta 分布是一類定義在有界區間(0,1)上的連續型概率分布，若隨機變量X服從Beta 分布，用f(?)表示Beta 分布的概率密度函數，則其對應的概率密度函數形式為：

式（2）中，Γ(?)為Gamma函數，α和β為參數，且要求滿足α>0，β>0。Beta分布對應的期望與方差分別為：

相比而言，其他分布往往形態單一，而Beta 分布的形態更具有靈活性，隨著參數α、β不同而不同。當α、β>0 時為單峰型形態；當α=β時為對稱型分布；當α≠β時為非對稱型分布，隨α、β取值的增大，峰度增加；當0<α、β<1 時呈現“U”型形態；當(α-1)(β-1)≤0 時呈現“J”型形態，即在0 或1 邊界的取值范圍對應的概率增大。有學者指出，利用多個Beta 分布的混合加權，能夠對金融市場特有的高峰厚尾、非對稱、有界取值區域等特征進行測度。因此，假設隨機變量X由K個Beta分布加權混合構成，記ωk為第k個Beta 分布的權重，k=1,2,…,K，h(?)為混合Beta 分布的概率密度函數，則X的概率密度函數為：

式（5）中，為保證f(x)作為概率密度函數所須滿足的非負性和規范性等基本性質，要求ωk≥0，且若記混合Beta分布的均值為m，方差為d2，組合中的第k個Beta分布的均值為μk，方差為σ2k，則：

1.3 混合Beta分布的GARCH模型構建

將混合Beta 分布引入GARCH 模型中，令式（2）中，其中X的分布由式（5）的形式決定，即將式（2）中的et設定為服從均值為0、方差為1 的非標準混合Beta分布。同時，根據GARCH模型的結構，εt服從均值為0、方差為ht的非標準混合Beta 分布；可觀測序列yt服從均值為w、方差為ht的非標準混合Beta分布，即有：

其中，d>0，ht>0。式（8）說明可觀測序列yt是隨機變量X（服從混合Beta 分布）的單調遞增線性變換函數，記yt的概率密度函數為g(yt)，根據概率論知識可得yt的概率密度函數為：

由此，在得到觀測值序列yt,t=1,2,…,T后，可以建立GARCH 模型的似然函數。對于金融市場而言，GARCH（1,1）模型能解釋波動的絕大部分信息，所以本文選擇GARCH（1,1）模型進行分析。對式（3）選擇滯后階數p與q均等于1，得到對應GARCH（1,1）模型的對數似然函數為：

其中，待估參數Λ=(w,γ0,η,λ,ωk,αk,βk;k=1,2,…,K)。

2 模型估計的EM算法求解

EM算法最初由Dempster等（1977）[11]提出，當似然函數有多余參數或含有隱變量時，可以通過兩步迭代進行求解，先進行E步，求均值以去掉多余參數，再進行M步，求似然函數的極大值，反復迭代E步和M步，直至收斂。而要采用EM算法，就要先構建包含隱變量的完全數據的似然函數。

2.1 完全數據的似然函數構建

假設時序樣本{yt,t=1,2,…,T}是從式（9）的非標準混合Beta分布中產生的，能觀測到的是樣本序列的取值yt，但是混合權重{ωk,k=1,2,…,K}是無法觀測的，所以對混合分布數據來說含有缺失的隱變量，也稱為不完全數據或者缺失數據。假設式（9）的生成是先以概率ωk抽取到第k個分量，然后再以這個分量的概率分布密度抽取到yt，其中該分量的概率密度函數為非標準的Beta分布，即有：

其中，fk(? |αk,βk)為式（4）Beta分布的概率密度函數，m和d分別由式（6）和式（7）決定。由此可以引入一個不能觀測的隨機隱變量{zt,t=1,2,…,T} ，zt取值為1,2,…,K，當zt=k時，說明第k類分量被抽中，對應概率即為P(zt=k)=ωk，顯然滿足：

可見，缺失數據yt對應的完全數據樣本序列應為{(y1,z1),(y2,z2),…,(yT,zT)}，對應概率密度函數為g(yt,zt)，雖然統計推斷只能基于缺失數據g(yt)進行，但是可以通過全概率公式建立兩者之間的關系：

所以，基于觀測數據{yt,t=1,2,…,T}的對數似然函數可改寫為如下形式：

其中，θ=(w,γ0,η,λ,αk,βk;k=1,2,…,K)，而包含隱變量的完全數據對應的對數似然為：

針對式（12）基于完全數據對數似然函數的最大化就可以用EM 算法求解，結合混合Beta 分布GARCH 模型的對數似然函數式（10），可以給出EM 算法的具體求解過程。

2.2 E步（求期望）

E 步，英文Expectation 的簡寫，即求期望。EM 算法也是一種迭代算法，以上標“*”標記上一步迭代值，記上一步參數θ的估計值為θ*，則在給定yt和θ*時，E步是計算完全數據對數似然函數關于未觀測數據zt的條件期望，該條件期望也稱為Q函數[12]。

其中，f(?)為Beta 分布的概率密度函數，將上式進一步分解可得由兩個部分構成的Q函數：

2.3 M步（最大化條件期望）

M 步，英文Maximum 的簡寫，表示最大化條件期望Q(θ|θ*)函數。根據式（13），前一部分僅與ωk有關，后一部分與參數向量θ=(w,γ0,η,λ,αk,βk;k=1,2,…,K)，即不包含ωk的部分有關，同時注意到，所以關于ωk有約束的一階條件為：

其中，τ為拉格朗日乘子，求解可得：

其中，ω*k為參數ωk的上一步迭代值，gk(yt|θk*)由式（11）決定。

對于其余參數(w,γ0,η,λ,αk,βk;k=1,2,…,K) 的估計只需最大化式（10）的后一部分。令：

其中：

可見，與參數(w,γ0,η,λ,αk,βk;k=1,2,…,K)有關的部分在于ln(?)部分，該H函數求解時無法得到顯式表達式，可以通過數值解方式求其極大值。最后，不斷重復E步和M步，直至收斂。EM算法在計算機上實現較為容易，其優點是每次迭代都增加似然函數的數值，一直到逼近最大值，因此尤其適合混合概率分布模型的參數估計求解。

2.4 EM算法實施流程

基于EM算法的混合Beta分布GARCH模型參數估計的流程可歸納如下：

（1）選取參數(w,γ0,η,λ,ωk,αk,βk;k=1,2,…,K)的初始值。

（2）依據當前模型參數值，生成混合Beta分布GARCH模型的ht序列。

（3）E步：代入參數當前值和序列ht，計算Q函數。

（4）M 步：最大化Q 函數，計算下一輪的參數迭代值(w,γ0,η,λ,ωk,αk,βk;k=1,2,…,K)。

（5）重復步驟（2）至步驟（4），直至收斂。

3 基于仿真模擬數據和現實數據的實證分析

3.1 仿真模擬數據的生成過程

假定一條GARCH（1,1）時間序列鏈生成過程為：

其中，et設置為由兩個Beta分布混合而成，具體形式為：0.85 ?Beta(38,38.5)+0.15 ?Beta(3.5,3.6)，即混合權重為0.85和0.15。采用模特卡洛仿真的方式，在計算機上模擬生成一條T=1000的序列yt，并將yt序列的走勢圖及其分布直方圖分別繪制為圖1 和圖2。從圖1 可以看出，yt序列的走勢符合GARCH模型特征，即自回歸條件異方差特征，當前序列的波動受到往期波動的影響；進一步對yt和進行Box檢驗，選擇滯后期為12期，yt序列的Q統計量對應的P 值為0.86，認為yt序列不存在自相關特征，但是序列的Q統計量對應的P值為2.2×10-16，顯著拒絕不存在序列自相關的原假設，所以Box檢驗結果也驗證了yt序列符合GARCH模型的生成機理。從圖2的yt序列分布直方圖來看，序列分布具有明顯的高峰厚尾、有限取值區域等特征，圖形中附加的分布曲線是以yt的均值和方差繪制的正態分布的概率密度曲線，對比可知正態分布的設定難以有效刻畫數據的真實特征，若不進行分布類型假設的調整，則勢必會產生模型設定上的偏誤問題。

圖1 yt 序列走勢圖

圖2 yt分布直方圖

3.2 模型參數的求解

根據前文給出的混合Beta分布GARCH模型的EM算法流程對yt序列進行參數估計求解，待估的參數向量包括(w,γ0,η,λ,ω1,ω2,α1,β1,α2,β2) ，對應的真值分別為（0.05, 0.0015, 0.4, 0.5, 0.85, 0.15, 38, 38.5, 3.5, 3.6），對于GARCH 模型的現實應用來說，尤其關注系數η和λ的估計，反映當前序列波動受到往期波動影響的程度。同時，為了對比不同分布類型假定下，參數求解的效果，本文將正態分布假定和混合Beta 分布假定進行了對比。傳統的GARCH 模型假定(0,1) ，根據模型結構可知(w,ht)，由此可得出樣本對數似然函數式（16），從而進行MLE求解，待估參數僅包括(w,γ0,η,λ)。

對yt序列建立GARCH(1,1)模型，基于正態分布假設的參數估計和基于混合Beta 分布的參數估計的結果如表1所示。

表1 yt 序列各參數估計結果

從表1的結果可知，混合Beta分布的參數估計與參數真實值的結果更為接近；同時，就本例而言，基于正態分布假設傾向于高估前一期條件方差的影響，低估GARCH項的影響程度，而基于混合Beta 分布正好相反。事實上，在金融市場里面，波動信息的沖擊往往具有較為持久的效應，因此采用混合Beta分布的模型設定在此類領域中更具有重要的現實意義。

3.3 現實數據分析

為了進一步體現模型的現實應用意義，本文選取從2015 年1 月1 日至2022 年6 月24 日1818 個交易日上證指數的收益率數據進行實證研究，去除掉序列自相關后的收益率序列標記為rt。序列走勢圖和其分布直方圖分別見圖3 和圖4。從圖3 的rt序列走勢圖來看，序列的波動具有GARCH模型的集群效應。從圖4的rt序列分布直方圖來看，序列分布的高峰厚尾、有限取值區域、非對稱等特征明顯，圖4中附加的分布曲線是以rt的均值和方差繪制的正態分布的概率密度曲線，對比可知正態分布的設定難以有效刻畫金融市場數據的真實特征。因此，本文運用混合Beta 分布GARCH 模型進行信息的提煉，考慮到計算量的大小，采用兩個Beta 分布的混合形式。同時，由于收益率序列rt受我國股票市場漲跌停板的限制，取值區間為[- 0.1,0.1] 的有限區域，因此通過線性變換(rt+0.1)/0.2 將rt的取值范圍轉換到Beta 分布的區間[0 ,1] 上。表2 展示了EM算法的參數求解結果。為了對比，也給出正態分布設定下的估計結果。從表2的結果來看，在正態分布假定下，對GARCH項影響的估計結果偏高，而對前一期擾動項平方影響的估計結果偏低；而在混合Beta 分布設定下，對這兩項的估計結果進行了修正，加強了前一期擾動項平方的影響，減弱了GARCH項的影響。雖然在金融市場中波動信息的沖擊具有較為持久的效應，但是相對而言，近期隨機信息沖擊對當前條件方差的影響較大，較遠時期的影響會隨著時間間隔的增加逐漸減弱。因此，混合Beta分布GARCH模型的估計結果更為合理。

表2 rt 序列各參數估計結果

圖3 rt 序列走勢圖

圖4 rt 分布直方圖

4 結束語

經典GARCH模型可以用來解決時序數據資料分析中方差齊性設定所引起的問題，通過將方差的條件自相關特性引入模型中，從而能更準確地擬合時序變量波動的變化規律。然而，在金融市場等現實領域中，還需關注的一個重要問題是模型分布類型的設定，高峰厚尾、非對稱、取值區間有限等特征難以被單一的正態分布所刻畫。有鑒于此，本文提出了基于混合Beta 分布的GARCH 模型，并通過分析給出了模型參數求解的EM 算法流程。基于模擬數據和現實數據的實證分析結果均顯示，混合Beta分布能更有效提煉波動的非正態性信息，同時，也驗證了EM 算法對模型的參數求解行之有效。需要指出的是，EM 算法也存在缺點，如對初始值的設置敏感、得到的最優解是局部最優、收斂速度較慢等，而且至今還沒有較好的解決方案。因此，在未來的進一步研究中，關于EM 算法的優化是一項值得探討的課題。