王祥賽

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

0 引言

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

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

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

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

1 含混合高斯項(xiàng)的AR-GJR-GARCH模型及貝葉斯后驗(yàn)密度

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

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

εt為一混合高斯分布項(xiàng)，即

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

2 隨機(jī)梯度哈密爾頓蒙特卡洛法

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

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

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

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

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

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

由中心極限定理可知

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

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

定理1微分方程組

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

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

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

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

3 實(shí)例分析

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

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

圖2 2020.1.3—2020.11.3收盤價(jià)變化率時(shí)間序列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ù)的波動(dòng)幅度隨時(shí)間變化，因此適合用異方差模型進(jìn)行估計(jì)，進(jìn)一步繪制了兩組數(shù)據(jù)的偏自相關(guān)函數(shù)（autocorrelation function，ACF）圖，如圖3和圖4。

圖3 2019.3.13—2020.1.2收盤價(jià)變化率偏自相關(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收盤價(jià)變化率偏自相關(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算法來對(duì)該模型的參數(shù)進(jìn)行貝葉斯估計(jì)，每次內(nèi)循環(huán)開始前從數(shù)據(jù)集中隨機(jī)選擇99個(gè)數(shù)據(jù)作為子集，并設(shè)置y1=0。利用該算法獲取1 000個(gè)樣本，分別得到這兩組數(shù)據(jù)的參數(shù)估計(jì)結(jié)果如表1和表2。

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

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

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

4 結(jié)語

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