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

石 凱，劉洪江，孫 峰

（樂山師范學(xué)院a.數(shù)理學(xué)院；b.旅游學(xué)院，四川 樂山 614000）

0 引言

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

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

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

1 混合Beta分布GARCH模型的構(gòu)建

1.1 GARCH模型的一般形式

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

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

1.2 混合Beta分布及其性質(zhì)

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

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

相比而言，其他分布往往形態(tài)單一，而Beta 分布的形態(tài)更具有靈活性，隨著參數(shù)α、β不同而不同。當(dāng)α、β>0 時為單峰型形態(tài)；當(dāng)α=β時為對稱型分布；當(dāng)α≠β時為非對稱型分布，隨α、β取值的增大，峰度增加；當(dāng)0<α、β<1 時呈現(xiàn)“U”型形態(tài)；當(dāng)(α-1)(β-1)≤0 時呈現(xiàn)“J”型形態(tài)，即在0 或1 邊界的取值范圍對應(yīng)的概率增大。有學(xué)者指出，利用多個Beta 分布的混合加權(quán)，能夠?qū)鹑谑袌鎏赜械母叻搴裎?、非對稱、有界取值區(qū)域等特征進行測度。因此，假設(shè)隨機變量X由K個Beta分布加權(quán)混合構(gòu)成，記ωk為第k個Beta 分布的權(quán)重，k=1,2,…,K，h(?)為混合Beta 分布的概率密度函數(shù)，則X的概率密度函數(shù)為：

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

1.3 混合Beta分布的GARCH模型構(gòu)建

將混合Beta 分布引入GARCH 模型中，令式（2）中，其中X的分布由式（5）的形式?jīng)Q定，即將式（2）中的et設(shè)定為服從均值為0、方差為1 的非標(biāo)準(zhǔn)混合Beta分布。同時，根據(jù)GARCH模型的結(jié)構(gòu)，εt服從均值為0、方差為ht的非標(biāo)準(zhǔn)混合Beta 分布；可觀測序列yt服從均值為w、方差為ht的非標(biāo)準(zhǔn)混合Beta分布，即有：

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

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

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

2 模型估計的EM算法求解

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

2.1 完全數(shù)據(jù)的似然函數(shù)構(gòu)建

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

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

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

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

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

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

2.2 E步（求期望）

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

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

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

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

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

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

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

其中：

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

2.4 EM算法實施流程

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

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

（2）依據(jù)當(dāng)前模型參數(shù)值，生成混合Beta分布GARCH模型的ht序列。

（3）E步：代入?yún)?shù)當(dāng)前值和序列ht，計算Q函數(shù)。

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

（5）重復(fù)步驟（2）至步驟（4），直至收斂。

3 基于仿真模擬數(shù)據(jù)和現(xiàn)實數(shù)據(jù)的實證分析

3.1 仿真模擬數(shù)據(jù)的生成過程

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

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

圖1 yt 序列走勢圖

圖2 yt分布直方圖

3.2 模型參數(shù)的求解

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

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

表1 yt 序列各參數(shù)估計結(jié)果

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

3.3 現(xiàn)實數(shù)據(jù)分析

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

表2 rt 序列各參數(shù)估計結(jié)果

圖3 rt 序列走勢圖

圖4 rt 分布直方圖

4 結(jié)束語

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