王祥賽

武漢大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，湖北武漢430072

0 引言

在對市場的風(fēng)險研究中，波動率作為其中重要的一個度量，一直被認(rèn)為是風(fēng)險管理等方面的一個重點研究對象。在波動率研究方面，Engle［1］提出了ARCH（autoregressive conditional heteroscedasticity）模型，Bollerslev［2］在其研究的基礎(chǔ)上又提出了GARCH（generalized ARCH）模型，該模型一定程度上包含了ARCH模型所不能涵蓋的市場信息。但隨著研究的深入，很多學(xué)者發(fā)現(xiàn)，現(xiàn)存的模型對一些金融數(shù)據(jù)的集聚性、尖峰厚尾性等不能做出解釋。為解釋這些特性，Nelson［3］提出了EARCH（exponential ARCH）模型，Glosten等［4］提出了GJR（Glosten-Jagannanthan-Runkle）-GARCH模型。

AR（autoregressive）-GJR-GARCH模型是在GJR-GARCH模型的基礎(chǔ)上，進一步添加了自回歸項發(fā)展而來，對此類模型的參數(shù)估計通常是采用貝葉斯估計法［5］。但在貝葉斯估計的后驗分布中通常存在一個未知的數(shù)值積分項，這使得直接獲取后驗分布的解析形式比較困難。為解決這個問題，Alexander等［6］用Griddy-Gibbs方法對后驗分布進行抽樣，但是這種方法在每次迭代前均需對數(shù)值積分做一個估計，導(dǎo)致計算成本比較高。

Chen等［7］和Talay［8］在引入物理學(xué)中哈密爾頓系統(tǒng)的基礎(chǔ)上，提出了一種隨機梯度哈密爾頓蒙特卡洛（stochastic gradient Hamiltonian Monte Carlo，SGHMC）方法，該方法在未知數(shù)值積分項時，也可以獲取一系列服從目標(biāo)后驗分布的樣本，而且不需要引入M-H（Metropolis-Hasting）步驟，因而該方法的計算復(fù)雜度大大降低?；诖朔椒?，Kim等［9］在貝葉斯神經(jīng)網(wǎng)絡(luò)的應(yīng)用上取得了進展。

本文在張新星等［10］給出的含混合高斯項的AR-GJR-GARCH模型基礎(chǔ)上，采用SGHMC方法對后驗分布進行采樣，給出了最終的參數(shù)估計結(jié)果，并結(jié)合結(jié)果對原始數(shù)據(jù)所包含的信息進行了解釋。

1 含混合高斯項的AR-GJR-GARCH模型及貝葉斯后驗密度

AR（s）-GJR-GARCH（p，q）模型如下（其中，s，p，q表示相應(yīng)的滯后階數(shù)）

其中，yt為觀察樣本，μ為觀察變量的長期穩(wěn)定均值，T為觀察樣本總數(shù)，I為示性函數(shù)，r=max{s，p，q}。實際情況中，為了確保AR（s）過程的平穩(wěn)性，滯后算子B的多項式的零點在單位圓外。

εt為一混合高斯分布項，即

各參數(shù)的先驗分布按照下面形式選取

2 隨機梯度哈密爾頓蒙特卡洛法

隨機梯度哈密爾頓蒙特卡洛法是在哈密爾頓系統(tǒng)的基礎(chǔ)上發(fā)展而來。假定觀察樣本的集合為Ω，我們感興趣的參數(shù)θ，其基于樣本觀察值x∈Ω的后驗密度p(θ|Ω)形式如下

這里將U(θ)稱為勢能函數(shù)，根據(jù)貝葉斯定理，該函數(shù)有如下形式

其中，p(x1，x2，…，xn|θ)為在給定θ時樣本的條件概率密度，p(θ)為θ的先驗概率密度。

為了使用哈密爾頓系統(tǒng)，還需要一個關(guān)于動能的函數(shù)，這里記動量變量為r，r一般為人為設(shè)置的輔助變量且服從正態(tài)分布N(0，M)，M稱之為質(zhì)量矩陣且一般為單位矩陣I。由θ和r一起構(gòu)造如下形式的微分方程組

其中，?U(θ)為U(θ)的梯度，這樣的微分方程組形成的系統(tǒng)稱為哈密爾頓系統(tǒng)。去掉采集到的樣本r那一部分，剩下的θ那一部分即為我們需要的滿足后驗分布p(θ|Ω)的樣本。

而隨機梯度漢密爾頓蒙特卡洛法是對原有使用的?U(θ)進行調(diào)整，不再選用全部的樣本集計算，而是從Ω中隨機抽取一個子集出來，并定義如下形式的隨機梯度

由中心極限定理可知

但該形式并不能保證對應(yīng)的哈密爾頓系統(tǒng)收斂到我們想要的平穩(wěn)分布π(θ，r)，需要添加一項進行修正，因此得到如下形式的微分方程組

為方便討論，這里假定B(θ)和θ無關(guān)，將其簡記為B。下面給出證明，該方程組可以收斂到我們想要的平穩(wěn)分布π(θ，r)。

定理1微分方程組

由Fokker-Planck方程易知，在t時刻(θ，r)的聯(lián)合概率密度pt(θ，r)有如下等式

在實際應(yīng)用中，通常不知道B的具體形式，因此要先利用樣本對B做一個估計，記為然后人為指定一項C≥，從而（11）式的第二個微分方程變?yōu)槿缦滦问?/p>

為了能盡可能多地利用觀察樣本的信息，在每構(gòu)建一個哈密爾頓系統(tǒng)前隨機抽取一個樣本子集進行隨機梯度的計算，這樣就克服了僅使用單一樣本子集信息量不足的缺陷，具體算法如下：

步驟1給定初值要采集的樣本數(shù)n，內(nèi)循環(huán)次數(shù)m，步長δ，矩陣C和

3 實例分析

在實例分析中采用中證醫(yī)藥2019年3月13日—2020年11月3日的指數(shù)收盤價數(shù)據(jù)的變化率，數(shù)據(jù)集總計有398個數(shù)據(jù)，將其拆分為2019年3月13日—2020年1月2日199個數(shù)據(jù)以及2020年1月3日—2020年11月3日199個數(shù)據(jù)分別利用上述模型進行分析。這兩部分的數(shù)據(jù)集的時間序列分別如圖1和圖2。

圖1 2019.3.13—2020.1.2收盤價變化率時間序列Fig.1 Time series of closing price change rate from March 13,2019 to January 2,2020

圖2 2020.1.3—2020.11.3收盤價變化率時間序列Fig.2 Time series of closing price change rate from January 3,2020 to November 3,2020

從圖1和圖2中可以發(fā)現(xiàn)這兩部分?jǐn)?shù)據(jù)的波動幅度隨時間變化，因此適合用異方差模型進行估計，進一步繪制了兩組數(shù)據(jù)的偏自相關(guān)函數(shù)（autocorrelation function，ACF）圖，如圖3和圖4。

圖3 2019.3.13—2020.1.2收盤價變化率偏自相關(guān)函數(shù)Fig.3 Close price change rate partial autocorrelation function from March 13,2019 to January 2,2020

圖4 2020.1.3—2020.11.3收盤價變化率偏自相關(guān)函數(shù)Fig.4 Close price change rate partial autocorrelation function from January 3,2020 to November 3,2020

可以發(fā)現(xiàn)兩組數(shù)據(jù)的偏自相關(guān)函數(shù)都是1階截尾的，因此選用AR（1）-GJR-GARCH（1，1）模型比較合適。利用MATLAB軟件使用SGHMC算法來對該模型的參數(shù)進行貝葉斯估計，每次內(nèi)循環(huán)開始前從數(shù)據(jù)集中隨機選擇99個數(shù)據(jù)作為子集，并設(shè)置y1=0。利用該算法獲取1 000個樣本，分別得到這兩組數(shù)據(jù)的參數(shù)估計結(jié)果如表1和表2。

表1 2019.3.13—2020.1.2參數(shù)估計結(jié)果表Table1 Parameter estimation result from March 13,2019 to January 2,2020

表2 2020.1.3—2020.11.3參數(shù)估計結(jié)果表Table2 Parameter estimation result from January 3,2020 to November 3,2020

綜上，我們的模型一定程度上反映了該指數(shù)在一段時期內(nèi)的背景信息，利用該模型可以對市場信息作一定的解釋。

4 結(jié)語

本文利用隨機梯度哈密爾頓蒙特卡洛方法在后驗密度解析形式未知的情況下得到了含混合高斯項的AR-GJR-GARCH模型各參數(shù)的貝葉斯估計。該方法相比于前人的Griddy-Gibbs方法優(yōu)勢在于無需在每次迭代前再對一個復(fù)雜積分項進行數(shù)值估計，降低了計算復(fù)雜度。實例分析給出了AR（1）-GJR-GARCH（1，1）模型在中證醫(yī)藥2019年3月13日—2020年11月3日指數(shù)收盤價數(shù)據(jù)變化率下的參數(shù)貝葉斯估計值，并根據(jù)得到的估計值推算市場背后信息，發(fā)現(xiàn)其較好地反映了該指數(shù)的市場波動情況與好壞消息對該指數(shù)的影響，驗證了模型和方法的合理性。