含混合高斯项的AR-GJR-GARCH模型的贝叶斯估计

Bayesian Estimation of AR-GJR-GARCH Model with Mixed Gaussian Terms

摘要

基于不需要后验密度解析形式的随机梯度哈密尔顿蒙特卡洛（stochastic gradient Hamiltonian Monte Carlo,SGHMC）方法对AR-GJR-GARCH模型的参数进行了贝叶斯估计。以2019.3.13—2020.1.2和2020.1.3—2020.11.3两个时间段的中证医药指数的数据为例，对本文提出的方法进行了检验。结果显示，所得的参数估计值反映了与该指数的波动性相关的市场背景信息。

Abstract

Based on the stochastic gradient Hamiltonian Monte Carlo(SGHMC) method which does not require a specific posterior density form, the parameters of AR（autoregressive）-GJR（Glosten-Jagannanthan-Runkle）-GARCH(generalized autoregressive conditional heteroscedasticity) model were estimated by Bayesian method. The method proposed in this paper was tested by taking the data of Chinese Medicine Index from March 13, 2019 to January 2, 2020 and January 3, 2020 to November 3, 2020 as examples. The results show that the estimated parameters reflect the market background information related to the volatility of the index.

Graphical abstract

关键词

AR-GJR-GARCH模型 / 贝叶斯估计 / 哈密尔顿系统 / 随机梯度哈密尔顿蒙特卡洛

Key words

AR-GJR-GARCH model / Bayesian estimation / Hamiltonian system / SGHMC(stochastic gradient Hamiltonian Monte Carlo)

0　 引 言

在对市场的风险研究中，波动率作为其中重要的一个度量，一直被认为是风险管理等方面的一个重点研究对象。在波动率研究方面，Engle^［1］提出了ARCH（autoregressive conditional heteroscedasticity）模型，Bollerslev^［2］在其研究的基础上又提出了GARCH（generalized ARCH）模型，该模型一定程度上包含了ARCH模型所不能涵盖的市场信息。但随着研究的深入，很多学者发现，现存的模型对一些金融数据的集聚性、尖峰厚尾性等不能做出解释。为解释这些特性，Nelson^［3］提出了EARCH（exponential ARCH）模型，Glosten等^［4］提出了GJR（Glosten-Jagannanthan-Runkle）-GARCH模型。

AR（autoregressive）-GJR-GARCH模型是在GJR-GARCH模型的基础上，进一步添加了自回归项发展而来，对此类模型的参数估计通常是采用贝叶斯估计法^［5］。但在贝叶斯估计的后验分布中通常存在一个未知的数值积分项，这使得直接获取后验分布的解析形式比较困难。为解决这个问题，Alexander等^［6］用Griddy-Gibbs方法对后验分布进行抽样，但是这种方法在每次迭代前均需对数值积分做一个估计，导致计算成本比较高。

Chen等^［7］和Talay^［8］在引入物理学中哈密尔顿系统的基础上，提出了一种随机梯度哈密尔顿蒙特卡洛（stochastic gradient Hamiltonian Monte Carlo，SGHMC）方法，该方法在未知数值积分项时，也可以获取一系列服从目标后验分布的样本，而且不需要引入M-H（Metropolis-Hasting）步骤，因而该方法的计算复杂度大大降低。基于此方法，Kim等^［9］在贝叶斯神经网络的应用上取得了进展。

本文在张新星等^［10］给出的含混合高斯项的AR-GJR-GARCH模型基础上，采用SGHMC方法对后验分布进行采样，给出了最终的参数估计结果，并结合结果对原始数据所包含的信息进行了解释。

1　 含混合高斯项的AR-GJR-GARCH模型及贝叶斯后验密度

AR（s）-GJR-GARCH（p，q）模型如下（其中，s，p，q表示相应的滞后阶数）

（1）

其中，y_t为观察样本，\mu为观察变量的长期稳定均值，T为观察样本总数，I为示性函数，r=\mathrm{m}\mathrm{a}\mathrm{x}\{s,p,q\}。实际情况中，为了确保AR（s）过程的平稳性，滞后算子\mathrm{B}的多项式\mathrm{1}-{\displaystyle \sum _{i=\mathrm{1}}^s}{\gamma }_i{\mathrm{B}}^i的零点在单位圆外。

{\epsilon }_t为一混合高斯分布项，即

{\epsilon }_t~\left\{\begin{array}{l}N(\mathrm{0},h_t{\sigma }^{\mathrm{2}}),\mathrm{概}\mathrm{率}\mathrm{为}\mathrm{\rho }\\ N(\mathrm{0},\frac{h_t{\sigma }^{\mathrm{2}}}{l}),\mathrm{概}\mathrm{率}\mathrm{为}\mathrm{1}-\rho \end{array}\right.,t=s+\mathrm{1},s+\mathrm{2},\cdots ,T

并假定{\epsilon }_t=\mathrm{0}(t=\mathrm{1,2},\cdots ,s)，。为了辅助参数估计，引入一个服从两点分布的辅助变量z_t，具体形式如下

z_t=\left\{\begin{array}{l}\mathrm{0},\mathrm{概}\mathrm{率}\mathrm{为}\mathrm{\rho }\\ \mathrm{1},\mathrm{概}\mathrm{率}\mathrm{为}\mathrm{1}-\rho \end{array}\right.,t=s+\mathrm{1},s+\mathrm{2},\cdots ,T

并假定z_t=\mathrm{0}(t=\mathrm{1,2},\cdots ,s)，长期平均方差项\lambda >\mathrm{0}，h_t=\mathrm{0}(t=\mathrm{1,2},\cdots ,r)，{\alpha }_{\mathrm{1}j}\ge \mathrm{0},{\alpha }_{\mathrm{2}j}\ge \mathrm{0}(j=\mathrm{1,2},\cdots ,p) （{\alpha }_{\mathrm{2}j}>{\alpha }_{\mathrm{1}j}时表明负面信息对观察变量的影响高于正面信息的影响，时则反之，{\alpha }_{\mathrm{2}j}={\alpha }_{\mathrm{1}j}时则表明二者对观察变量的影响相同），{\beta }_k>\mathrm{0}(k=\mathrm{1,2},\cdots ,q)。其中，为了确保GARCH模型是稳定的，通常要求。为方便接下来的推导，记\vec{y}=(y_{\mathrm{1}},y_{\mathrm{2}},\cdots ,y_T),\vec{\gamma }=({\gamma }_{\mathrm{1}},{\gamma }_{\mathrm{2}},\cdots ,{\gamma }_s),\vec{z}=(z_{\mathrm{1}},z_{\mathrm{2}},\cdots ,z_T),{\vec{\alpha }}_{\mathrm{1}}=({\alpha }_{\mathrm{11}},{\alpha }_{\mathrm{12}},\cdots ,{\alpha }_{\mathrm{1}p}),{\vec{\alpha }}_{\mathrm{2}}=({\alpha }_{\mathrm{21}},{\alpha }_{\mathrm{22}},\cdots ,{\alpha }_{\mathrm{2}p}),\vec{\beta }=({\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},\cdots ,{\beta }_q),\vec{\theta }=(\rho ,l,\lambda ,\mu ,\vec{\gamma },{\vec{\alpha }}_{\mathrm{1}},{\vec{\alpha }}_{\mathrm{2}},\vec{\beta })。

接下来给出参数的贝叶斯后验分布表达式。给定观察样本\vec{y}以及辅助变量\vec{z}，参数\vec{\theta }的后验分布密度\pi (\vec{\theta }|\vec{y},\vec{z})由贝叶斯定理有如下形式

\pi (\vec{\theta }|\vec{y},\vec{z})\propto L(\vec{y},\vec{z}|\vec{\theta })p(\vec{\theta })（2）

其中，L(\vec{y},\vec{z}|\vec{\theta })为似然函数，p(\vec{\theta })为\vec{\theta }的先验概率密度。

对于似然函数L(\vec{y},\vec{z}|\vec{\theta })，在给定辅助变量z_t时，可得y_t的条件分布如下

y_t|h_t,y_{t-\mathrm{1}},y_{t-\mathrm{2}},\cdots ,y_{t-s},z_t~\left\{\begin{array}{l}N(\mu +{\displaystyle \sum _{i=\mathrm{1}}^s}{\gamma }_i(y_{t-i}-\mu ),h_t{\sigma }^{\mathrm{2}}),z_t=\mathrm{0}\\ N(\mu +{\displaystyle \sum _{i=\mathrm{1}}^s}{\gamma }_i(y_{t-i}-\mu ),\frac{h_t{\sigma }^{\mathrm{2}}}{l}),z_t=\mathrm{1}\end{array}\right.

由此可得似然函数L(\vec{y},\vec{z}|\vec{\theta })如下

\begin{array}{l}L(\vec{y},\vec{z}|\vec{\theta })\propto L(\vec{y}|\vec{z},\vec{\theta })p_Z(\vec{z}|\vec{\theta })\propto \prod _{z_t=\mathrm{0}}[\rho ({\sigma }^{\mathrm{2}}{h}_t{)}^{-\frac{\mathrm{1}}{\mathrm{2}}}exp\{-\frac{[y_t{-\mu -{\displaystyle \sum _{i=\mathrm{1}}^s}{\gamma }_i(y_{t-i}-\mu )]}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}{h}_t}\}]\times \\ \prod _{z_t=\mathrm{1}}[(\mathrm{1}-\rho ){(\frac{{\sigma }^{\mathrm{2}}{h}_t}{l})}^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{e}\mathrm{x}\mathrm{p}\{-\frac{l[y_t{-\mu -{\displaystyle \sum _{i=\mathrm{1}}^s}{\gamma }_i(y_{t-i}-\mu )]}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}{h}_t}\}]\end{array}（3）

各参数的先验分布按照下面形式选取

\begin{array}{l}\rho ~U(\mathrm{0.5,1}),l~U(\mathrm{0,1}),\lambda ~U(\mathrm{0},{\widehat{\sigma }}_y^{\mathrm{2}}),\mu ~U(\overline{y}-\frac{\mathrm{5}{\widehat{\sigma }}_y}{\sqrt[]{T}},\overline{y}+\frac{\mathrm{5}{\widehat{\sigma }}_y}{\sqrt[]{T}}),\gamma ~U(-\mathrm{1,1})\\ {\vec{\alpha }}_{\mathrm{1}}~U{(\mathrm{0,1})}^p,{\vec{\alpha }}_{\mathrm{2}}~U{(\mathrm{0,1})}^p,\vec{\beta }~U{(\mathrm{0,1})}^q\end{array}

其中，\overline{y},{\widehat{\sigma }}_y分别为样本观察值的均值与标准差，并且假定各参数相互独立。将以上式子组合在一起即得\vec{\theta }的后验分布概率密度。

对于z_t的后验分布密度，由其先验分布以及关于y_t和\vec{\theta }的条件分布可得

\begin{array}{l}p(z_t=\mathrm{0}|y_{\mathrm{1}:t},\vec{\theta })=\frac{p_{\mathrm{1}}}{p_{\mathrm{1}}+p_{\mathrm{2}}}\\ p(z_t=\mathrm{1}|y_{\mathrm{1}:t},\vec{\theta })=\frac{p_{\mathrm{2}}}{p_{\mathrm{1}}+p_{\mathrm{2}}}\end{array}

其中，p_{\mathrm{1}}=\rho \mathrm{e}\mathrm{x}\mathrm{p}\{-\frac{[y_t{-\mu -{\displaystyle \sum _{i=\mathrm{1}}^s}{\gamma }_i(y_{t-i}-\mu )]}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}{h}_t}\},p_{\mathrm{2}}=(\mathrm{1}-\rho )l^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{e}\mathrm{x}\mathrm{p}\{-\frac{l[y_t{-\mu -{\displaystyle \sum _{i=\mathrm{1}}^s}{\gamma }_i(y_{t-i}-\mu )]}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}{h}_t}\}。

2　 随机梯度哈密尔顿蒙特卡洛法

随机梯度哈密尔顿蒙特卡洛法是在哈密尔顿系统的基础上发展而来。假定观察样本的集合为\Omega，我们感兴趣的参数\theta，其基于样本观察值x\in \Omega的后验密度p(\theta |\Omega )形式如下

p(\theta |\Omega )\propto {\mathrm{e}}^{-U(\theta )}（4）

这里将U(\theta )称为势能函数，根据贝叶斯定理，该函数有如下形式

U(\theta )=-\mathrm{l}\mathrm{o}\mathrm{g}p(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n|\theta )-\mathrm{l}\mathrm{o}\mathrm{g}p(\theta ),x_i\in \mathrm{\Omega }（5）

其中，p(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n|\theta )为在给定\theta时样本的条件概率密度，p(\theta )为\theta的先验概率密度。

为了使用哈密尔顿系统，还需要一个关于动能的函数，这里记动量变量为r，r一般为人为设置的辅助变量且服从正态分布N(\mathrm{0},\mathit{M})，\mathit{M}称之为质量矩阵且一般为单位矩阵\mathit{I}。由\theta和r一起构造如下形式的微分方程组

\left\{\begin{array}{l}\mathrm{d}\theta ={\mathit{M}}^{-\mathrm{1}}r\mathrm{d}t\\ \mathrm{d}r=-\nabla U(\theta )\mathrm{d}t\end{array}\right.（6）

其中，\nabla U(\theta )为U(\theta )的梯度，这样的微分方程组形成的系统称为哈密尔顿系统。去掉采集到的样本r那一部分，剩下的\theta那一部分即为我们需要的满足后验分布p(\theta |\Omega )的样本。

而随机梯度汉密尔顿蒙特卡洛法是对原有使用的\nabla U(\theta )进行调整，不再选用全部的样本集计算，而是从\Omega中随机抽取一个子集\tilde{\Omega }出来，并定义如下形式的随机梯度

\nabla \tilde{U}(\theta )=-\frac{\left|\Omega \right|}{\left|\tilde{\Omega }\right|}\nabla \mathrm{l}\mathrm{o}\mathrm{g}p(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_m|\theta )-\nabla \mathrm{l}\mathrm{o}\mathrm{g}p(\theta ),x_j\in \tilde{\Omega }（7）

由中心极限定理可知

\nabla \tilde{U}(\theta )\approx \nabla U(\theta )+N(\mathrm{0},\mathit{V}(\theta ))（8）

其中，\mathit{V}(\theta )为一未知的协方差矩阵。在实际应用哈密尔顿系统时通常会设置一个迭代步长\delta，由此可知辅助变量r的变化量\mathrm{\Delta }r在\nabla \tilde{U}(\theta )作用时有\mathrm{\Delta }r=-\delta \nabla \tilde{U}(\theta )=-\delta \nabla U(\theta )+N(\mathrm{0},{\delta }^{\mathrm{2}}\mathit{V}(\theta ))，并记B(\theta )=\frac{\mathrm{1}}{\mathrm{2}}\delta \mathit{V}(\theta )，则原微分方程组变为如下形式

\left\{\begin{array}{l}\mathrm{d}\theta ={\mathit{M}}^{-\mathrm{1}}r\mathrm{d}t\\ \mathrm{d}r=-\nabla U(\theta )\mathrm{d}t+N(\mathrm{0,2}B(\theta )\mathrm{d}t)\end{array}\right.（9）

但该形式并不能保证对应的哈密尔顿系统收敛到我们想要的平稳分布\pi (\theta ,r)，需要添加一项进行修正，因此得到如下形式的微分方程组

\left\{\begin{array}{l}\mathrm{d}\theta ={\mathit{M}}^{-\mathrm{1}}r\mathrm{d}t\\ \mathrm{d}r=-\nabla U(\theta )\mathrm{d}t-B(\theta ){\mathit{M}}^{-\mathrm{1}}r\mathrm{d}t+N(\mathrm{0,2}B(\theta )\mathrm{d}t)\end{array}\right.（10）

为方便讨论，这里假定B(\theta )和\theta无关，将其简记为B。下面给出证明，该方程组可以收敛到我们想要的平稳分布\pi (\theta ,r)。

定理1 微分方程组

\left\{\begin{array}{l}\mathrm{d}\theta ={\mathit{M}}^{-\mathrm{1}}r\mathrm{d}t\\ \mathrm{d}r=-\nabla U(\theta )\mathrm{d}t-B{\mathit{M}}^{-\mathrm{1}}r\mathrm{d}t+N(\mathrm{0,2}B\mathrm{d}t)\end{array}\right.（11）

当其满足：任意t时刻的联合概率密度p_t(\theta ,r)关于\theta，r的二阶混合偏导数连续时，可以收敛到\pi (\theta ,r)\propto \mathrm{e}\mathrm{x}\mathrm{p}(-U(\theta )-\frac{\mathrm{1}}{\mathrm{2}}{r}^{\mathrm{\text{'}}}{\mathit{M}}^{-\mathrm{1}}r)的平稳分布。

证 记H(\theta ,r)=U(\theta )+\frac{\mathrm{1}}{\mathrm{2}}{r}^{\mathrm{\text{'}}}{\mathit{M}}^{-\mathrm{1}}r，\mathit{G}=\left[\begin{array}{cc}\mathrm{0}& -\mathit{I}\\ \mathit{I}& \mathrm{0}\end{array}\right]，\mathit{D}=\left[\begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \mathrm{0}& B\end{array}\right]，\mathit{I}为单位矩阵，则（11）式可记为以下形式

\mathrm{d}\left[\begin{array}{l}\theta \\ r\end{array}\right]=-\left[\begin{array}{cc}\mathrm{0}& -\mathit{I}\\ \mathit{I}& B\end{array}\right]\left[\begin{array}{l}\nabla U(\theta )\\ {\mathit{M}}^{-\mathrm{1}}r\end{array}\right]\mathrm{d}t+N(\mathrm{0,2}\mathit{D}\mathrm{d}t)=-[\mathit{D}+\mathit{G}]\nabla H(\theta ,r)\mathrm{d}t+N(\mathrm{0,2}\mathit{D}\mathrm{d}t)（12）

由Fokker-Planck方程易知，在t时刻(\theta ,r)的联合概率密度p_t(\theta ,r)有如下等式

\frac{\partial p_t(\theta ,r)}{\partial t}={\nabla }^{\mathrm{T}}\{[\mathit{D}+\mathit{G}][\nabla H(\theta ,r)p_t(\theta ,r)]\}+{\nabla }^{\mathrm{T}}[\mathit{D}\nabla p_t(\theta ,r)]（13）

而由定理中的条件，可以得到\frac{{\partial }^{\mathrm{2}}{p}_t(\theta ,r)}{\partial \theta \partial r}=\frac{{\partial }^{\mathrm{2}}{p}_t(\theta ,r)}{\partial r\partial \theta }，从而有{\nabla }^{\mathrm{T}}[G\nabla p_t(\theta ,r)]=-\frac{{\partial }^{\mathrm{2}}{p}_t(\theta ,r)}{\partial \theta \partial r}+\frac{{\partial }^{\mathrm{2}}{p}_t(\theta ,r)}{\partial r\partial \theta }=\mathrm{0}，因此（13）式可化为

\frac{\partial p_t(\theta ,r)}{\partial t}={\nabla }^{\mathrm{T}}\{[\mathit{D}+\mathit{G}][\nabla H(\theta ,r)p_t(\theta ,r)+\nabla p_t(\theta ,r)]\}（14）

当p_t(\theta ,r)=\pi (\theta ,r)时，（14）式右边有\nabla H(\theta ,r){\mathrm{e}}^{-H(\theta ,r)}+\nabla {\mathrm{e}}^{-H(\theta ,r)}=\mathrm{0}，因此有\frac{\partial p_t(\theta ,r)}{\partial t}\left|{}_{p_t=\pi }\right.=\mathrm{0}，这表明\pi (\theta ,r)和时间t无关，因而\pi (\theta ,r)就是该系统下的平稳分布，证明完毕。

在实际应用中，通常不知道\mathit{B}的具体形式，因此要先利用样本对\mathit{B}做一个估计，记为\widehat{\mathit{B}}，然后人为指定一项\mathit{C}\ge \widehat{\mathit{B}}，从而（11）式的第二个微分方程变为如下形式

\begin{array}{l}\mathrm{d}r=-\nabla U(\theta )\mathrm{d}t-\mathit{C}{\mathit{M}}^{-\mathrm{1}}r\mathrm{d}t+N(\mathrm{0,2}(\mathit{C}-\widehat{\mathit{B}})\mathrm{d}t)+N(\mathrm{0,2}\mathit{B}\mathrm{d}t)=\\ -\nabla \tilde{U}(\theta )\mathrm{d}t-\mathit{C}{M}^{-\mathrm{1}}r\mathrm{d}t+N(\mathrm{0,2}(\mathit{C}-\widehat{\mathit{B}})\mathrm{d}t)\end{array}

为了能尽可能多地利用观察样本的信息，在每构建一个哈密尔顿系统前随机抽取一个样本子集进行随机梯度的计算，这样就克服了仅使用单一样本子集信息量不足的缺陷，具体算法如下：

步骤1 给定初值{\vec{\theta }}^{(\mathrm{1})}，要采集的样本数n，内循环次数m，步长\delta，矩阵\mathit{C}和\widehat{\mathit{B}}

步骤2 从观察样本全集\Omega中抽取一个样本子集\tilde{\Omega }。后面所用到的样本\vec{y}均为\tilde{\Omega }内的样本

步骤3 抽取样本r~N(\mathrm{0},M)，取{\theta }_{\mathrm{0}}为前一次采集得到的样本值，r_{\mathrm{0}}=r

步骤4{\vec{\theta }}_i={\vec{\theta }}_{i-\mathrm{1}}+\delta {\mathit{M}}^{-\mathrm{1}}{r}_{i-\mathrm{1}}

步骤5 抽取\vec{z}|{\vec{\theta }}_i,\vec{y}

步骤6r_i=r_{i-\mathrm{1}}-\nabla \tilde{U}({\vec{\theta }}_i)\delta -\mathit{C}{\mathit{M}}^{-\mathrm{1}}{r}_{i-\mathrm{1}}\delta +N(\mathrm{0,2}(\mathit{C}-\widehat{\mathit{B}})\delta )

步骤7 重复步骤4~6 m次，得到{\vec{\theta }}_m，将其作为该次采样过程中取得的样本。

重复步骤2~7 n次，即取得样本{\vec{\theta }}^{(\mathrm{1})},{\vec{\theta }}^{(\mathrm{2})},\cdots ,{\vec{\theta }}^{(n)}。

3　 实例分析

在实例分析中采用中证医药2019年3月13日—2020年11月3日的指数收盘价数据的变化率，数据集总计有398个数据，将其拆分为2019年3月13日—2020年1月2日199个数据以及2020年1月3日—2020年11月3日199个数据分别利用上述模型进行分析。这两部分的数据集的时间序列分别如图1和图2。

从图1和图2中可以发现这两部分数据的波动幅度随时间变化，因此适合用异方差模型进行估计，进一步绘制了两组数据的偏自相关函数（autocorrelation function，ACF）图，如图3和图4。

可以发现两组数据的偏自相关函数都是1阶截尾的，因此选用AR（1）-GJR-GARCH（1，1）模型比较合适。利用MATLAB软件使用SGHMC算法来对该模型的参数进行贝叶斯估计，每次内循环开始前从数据集中随机选择99个数据作为子集\tilde{\Omega }，并设置y_{\mathrm{1}}=\mathrm{0}。利用该算法获取1 000个样本，分别得到这两组数据的参数估计结果如表1和表2。

选择样本均值作为参数估计值。表1和表2的参数，这说明模型基本是稳定的。表1中\mu的均值为\mathrm{0.026}\mathrm{127}，表2中\mu的均值为\mathrm{0.026}\mathrm{237}，说明两组数据的收盘价均有一个缓慢上涨的趋势；二者数值差别不大，说明在这前后两个时间段该指数的收盘行情基本相差不大。表1中{\alpha }_{\mathrm{11}}和{\alpha }_{\mathrm{21}}的均值相差仅为0.000 205，变化率仅为0.08%，说明这段时期该指数的杠杆效应并不明显，好消息和坏消息对该指数产生的影响可以大致抵消；但在表2中，可以看到{\alpha }_{\mathrm{11}}的均值明显高于{\alpha }_{\mathrm{21}}的均值，说明在这段时间好消息对该指数的影响超过了坏消息。结合该时期分析，这段时间我国爆发了新冠疫情，疫情的出现使得原本相对平稳的医药生物产业走进了人们的视野，一时间导致大量投资者涌入了该产业，并且随着疫情转为长期存在的情形，许多投资者开始长期在该板块中投入。

综上，我们的模型一定程度上反映了该指数在一段时期内的背景信息，利用该模型可以对市场信息作一定的解释。

4　 结 语

本文利用随机梯度哈密尔顿蒙特卡洛方法在后验密度解析形式未知的情况下得到了含混合高斯项的AR-GJR-GARCH模型各参数的贝叶斯估计。该方法相比于前人的Griddy-Gibbs方法优势在于无需在每次迭代前再对一个复杂积分项进行数值估计，降低了计算复杂度。实例分析给出了AR（1）-GJR-GARCH（1，1）模型在中证医药2019年3月13日—2020年11月3日指数收盘价数据变化率下的参数贝叶斯估计值，并根据得到的估计值推算市场背后信息，发现其较好地反映了该指数的市场波动情况与好坏消息对该指数的影响，验证了模型和方法的合理性。

参考文献

[1] ENGLE R F. Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation ［J］. Econometrica， 1982， 50（4）： 987-1008. DOI：10.2307/1912773 .

[2] BOLLERSLEV T. Generalized autoregressive conditional heteroscedasticity ［J］. Journal of Econometrics，1986，31（3）：307-327. DOI: 10.1016/0304-4076(86)90063-1 .

[3] NELSON D B. Conditional heteroskedasticity in asset returns： A new approach ［J］. Econometrica， 1991， 59（2）： 347. DOI：10.2307/2938260 .

[4] GLOSTEN L R， JAGANNATHAN R， RUNKLE D E. On the relation between the expected value and the volatility of the nominal excess return on stocks［J］. The Journal of Finance， 1993， 48（5）： 1779-1801. DOI：10.1111/j.1540-6261.1993.tb05128.x .

[5] TAKAISHI T. Markov chain Monte Carlo versus importance sampling in Bayesian inference of the GARCH model ［J］. Procedia Computer Science， 2013， 22： 1056-1064. DOI：10.1016/j.procs.2013.09.191 .

[6] ALEXANDER C， LAZAR E. Normal mixture GARCH（1，1）：Applications to exchange rate modelling ［J］. Journal of Applied Econometrics， 2006， 21（3）： 307-336. DOI：10.1002/jae.849 .

[7] CHEN T Q， FOX E B， GUESTRIN C. Stochastic gradient Hamiltonian Monte Carlo ［EB/OL］. ［2020-11-12］.

[8] TALAY D. Stochastic Hamiltonian systems： Exponential convergence to the invariant measure， and discretization by the implicit Euler scheme ［J］. Markov Processes and Related Fields，2002，8：1-36.

[9] KIM M， LEE J. Hamiltonian Markov chain Monte Carlo for partitioned sample spaces with application to Bayesian deep neural nets ［J］. Journal of the Korean Statistical Society， 2020， 49（1）： 139-160. DOI：10.1007/s42952-019-00001-3 .

[10] 张新星， 唐亚勇. 基于Griddy-Gibbs抽样的混合高斯AR-GJR-GARCH模型的贝叶斯估计［J］. 四川大学学报（自然科学版）， 2016， 53（5）： 957-962. DOI: 103969/j.issn.0490-6756.2016.09.001 .

[11] ZHANG X X， TANG Y Y. Bayesian estimation of the Gaussian mixture AR-GJR-GARCH model with Griddy-Gibbs sampler［J］. Journal of Sichuan University（Natural Science Edition）， 2016， 53（5）： 957-962 （Ch）. DOI: 103969/j.issn.0490-6756.2016.09.001 .

基金资助

国家自然科学基金(11731012)