第24卷第4期 统计与信息论坛 2009年4月 Vol.24 No.4 Statistics Information Forum Apr.,2009 【统计理论与方法】 基于MCMC方法的GARCH模型参数估计 潘海涛1，温小霓2 (1.西安财经学院统计学院，陕西西安710100：2.西安电子科技大学经济管理学院，陕西西安710061) 摘要：参数GARCH模型是最常用的度量金融市场波动性的模型。运用马尔科夫蒙特卡罗(MCMC)方法 对残差基于正态分布的GARCH(1,1)的参数进行估计，由沪市股指收益率数据的实证分析结果表明：基于马 尔科夫蒙特卡罗(MCMC)方法估计的GARCH模型比基于极大似然估计(ML)方法估计的GARCH模型具有 更好的拟合效果和预测能力。 关键词：GARCH(l,l)模型；MCMC算法；Gibbs抽样；Metropolis-Hasting算法；波动性；预测 中图分类号：0212.8文献标志码：A文章编号：1007-3116(2009)04-0012-05 在大多数研究中，对于参数GARCH模型系数 假定我们希望从分布f(日1X)中抽取一个随 的估计采用的都是经典的方法，例如极大似然法 机样本，它包含一个复杂的标准化常数，直接的抽取 (ME),当我们用基于梯度的方法来进行数值最优 要么太浪费时间，要么不可行，但若存在一个近似分 化求解过程中，目标函数应当是凸函数，否则，经典 布，且可以很容易的得到随机样本的话，那么则可以 的方法很难应用。此外，对于GARCH模型的系数 应用Metropolis算法。Metropolis算法就是从近似分 所进行的限制，也使得进一步对系数进行统计推断 布中产生一系列的随机抽取，使其分布函数收敛到 变得困雄。为了使得GARCH模型成为平稳过程和 f(01X)此算法如下： 满足其他的一些性质，从本质上要将GARCH模型 步骤1：抽取一个随机的初始值0o,满足f(001 表示为一些复杂的不等式组，这使得经典的方法很 X)>0; 难对GARCH的系数进行统计推断[1-3)。 步骤2：对t=1,2,… 本文引人基于Bayes计算的MCMC方法来对 (a)第t次迭代时，在给定前面的抽取0，-1下， GARCH模型的参数进行估计，通过构造GARCH 从已知分布中抽取一个候选样本日，。用J(61 模型系数的算法并结合沪市股指收益率实际数据进 0-1)表示已知分布，在Gelman中称为跳跃分布。这 行分析，证实了该方法的有效性。 个跳跃分布要求是对称的，即对于所有的日：，日，t, 一、Gibbs抽样和M-H算法 有J(018)=J(81a) (6)计算比率r= f(9.1X) Geman和Geman、Gelfand和Smith的Gibbs抽 f(0-11X) 样是应用最广的MCMC方法。Gibbs抽样具有将 8.;以概率min(r,l) 高维的估计问题利用所有参数的条件分布分解为几 (c)设定0.=0-1，其他 个低维问题的优点，具体抽样的思想和步骤可参看 在一些正则条件下，序列{日}依分布收敛到 上述学者的文章著作。下面仅简单介绍Metropolis f(01X),具体证明见Gelman等。 算法和Metropolis一Hasting算法。 算法的实施要求对所有的日，和日，-1计算比率 收稿日期：2009-02-19 基金项目：西安统计研究院课题《非参数回归的稳健Bot2rp在经济数据中的应用》(07JD13) 作者简介：潘海涛(1978-)，男，浙江奉化人，讲师，研究方向：统计计算，非参数方法； 沮小克(1955-)，女，山西孝义人，数授，研究方向：多元统计方法，经济时间序列分析。 12 万方数据
第24卷第4期 Vd．24 Nb．4 统计与信息论坛 Statistics&In{bnllation Fbrum 2009年4月 Apr-，2009 【统计理论与方法】 基于MCMC方法的GARCH模型参数估计 · 潘海涛1，温小霓2 (1．西安财经学院统计学院，陕西西安710100；2．西安电子科技大学经济管理学院，陕西西安710061) 摘要：参数GARCH模型是最常用的度量金融市场波动性的模型。运用马尔科夫蒙特卡罗(McMC)方法 对残差基于正态分布的GARCH(1，1)的参数进行估计，由沪市股指收益率数据的实证分析结果表明：基于马 尔科夫蒙特卡罗(McMc)方法估计的GARCH模型比基于极大似然估计(ML)方法估计的GARCH模型具有 更好的拟合效果和预测能力。 关键词：GARCH(1，1)模型；MCMC算法；GibbS抽样；Met删dis—HaSting算法；波动性；预测 中图分类号：0212．8 文献标志码：A 文章编号：1007—3116(2009)04—0012—05 在大多数研究中，对于参数GARCH模型系数 的估计采用的都是经典的方法，例如极大似然法 (MLE)，当我们用基于梯度的方法来进行数值最优 化求解过程中，目标函数应当是凸函数，否则，经典 的方法很难应用。此外，对于GARCH模型的系数 所进行的限制，也使得进一步对系数进行统计推断 变得困难。为了使得GARCH模型成为平稳过程和 满足其他的一些性质，从本质上要将Q娘CH模型 表示为一些复杂的不等式组，这使得经典的方法很 难对G根CH的系数进行统计推断[卜31。 本文引入基于Bayes计算的MCMC方法来对 GARCH模型的参数进行估计，通过构造Q啵CH 模型系数的算法并结合沪市股指收益率实际数据进 行分析，证实了该方法的有效性。 一、Gibbs抽样和M—H算法 G锄ail和G电man、G毛l缸d和Smith的GibbS抽 样是应用最广的MCMC方法。GibbS抽样具有将 高维的估计问题利用所有参数的条件分布分解为几 个低维问题的优点，具体抽样的思想和步骤可参看 上述学者的文章著作。下面仅简单介绍Metroldis 算法和Metropolis—HaSting算法。 假定我们希望从分布厂(口l X)中抽取一个随 机样本，它包含一个复杂的标准化常数，直接的抽取 要么太浪费时间，要么不可行，但若存在一个近似分 布，且可以很容易的得到随机样本的话，那么则可以 应用Metropolis算法。Men即lis算法就是从近似分 布中产生一系列的随机抽取，使其分布函数收敛到 厂(口I X)。此算法如下： 步骤1：抽取一个随机的初始值日o，满足厂(臼o X)>0； 步骤2：对￡=1，2，… (a)第￡次迭代时，在给定前面的抽取f9f一1下， 从已知分布中抽取一个候选样本口。。用^(岛I 吼一1)表示已知分布，在叫man中称为跳跃分布。这 个跳跃分布要求是对称的，即对于所有的皖，乱，￡， 有^(岛I谚)=^(岛I巩) (b)计算比率r=舷隅 (c)设定馥=钯雾嘶出’1’ 在一些正则条件下，序列{以}依分布收敛到 ．厂(口l X)，具体证明见蹦man等。 算法的实施要求对所有的a。和以一1计算比率 收稿日期：2009一02—19 基金项目：西安统计研究院课题《非参数回归的稳健Boo晓舰p在经济数据中的应用》(07JDl3) 作者简介：潘海涛(1978一)，男，浙江奉化人，讲师，研究方向：统计计算，非参数方法； 温小霓(1955一)，女，山西孝义人，教授，研究方向：多元统计方法，经济时间序列分析。 万方数据
潘海游，温小宽：基于MCMC方法的GARCH模型参数估计 r,以便从跳跃分布中抽取日，，并从均匀分布中抽 的，这样就无法找到参数的共轭先验分布，因此，借 取一个随机实现，决定接受或者拒绝0，。不需要 助于Metropolis-Hasting(M-H)算法。具体方法 f(01X)的标准化常数，因为这里只利用比率。 如下： 此算法的接受和拒绝准则思想如下：1.如果从 参数的初值为⊙o】=(ao,po),因为，无法 日，-1到0，的跳跃增加了条件后验密度，则接受0。 得到条件分布的显示表示，从参数α，B各自的分布 作为日，；2.如果这个跳跃降低了后验密度，则以等 进行随机数的抽取，运用(M-H)算法J次，可以得 于密度比r的概率设定8，=0，，否则设定日= 到第J次参数的迭代结果为： 9-1这个思想显然是合理的。 alil ~p(a l pli-n,y) Hasting以两种方式推广了Metropolis算法。首 g1-p(B1a,y）) 先，跳跃分布无需对称；其次，跳跃准则修正为： 其次，注意到GARCH(1,1)模型可以写为如下 r=fa.1X00.(0、102 的ARMA(1,1)模型： f0-11X)/J(0-116.) h=a0+a1y2-1+:-l 这个修正的算法称为Metropolis-Hasting算法。 台y2=ao+(a1+)y经-1--1+w 下面应用MCMC算法的思想来构造GARCH (1,1)模型中参数的估计方法。 其中m,=好-h=(学-1Dh h 二、MCMC算法在GARCH(1,1) =（x7-1)h (1) 由{}的构造可知，序列{w}是一个鞅差过 参数估计中的应用 程(Martingale Difference Process),可以得到其方差 新息(innovations)为正态分布的GARCH(1,l) 为2h。据Nakatsuma知，很难直接从式(1)中生成 模型如下： 参数⊙的随机数。因此，用一个期望为0，方差为2h2 =ehn,t=1,…,T 的正态随机变量：来近似w。据此，可以构造如下 e:~idN(0,1) 的辅助方程： h:=a0+a1yi-1-,-1 y2=a0+(a1+)y2-1--1+4 这里，a0>0,a1≥0，B≥0。条件方差h,是滞 注意到，名：和h:都是⊙的函数，即有 后一阶平方观察值和滞后一阶方差的线性函数。参 z(8)=y2-a0-(a1+)y2-1+-(⊙) 数的正定性的限制保证了方差为正定的。为了写出 h,(⊙)=a0+a1y2-1+hr-1(⊙) (2) 似然函数，定义以下的向量：y=（y1,…,yT）',a= 再记(T×T)的对角矩阵A=A(⊙)= (ao,a1)',然后将所有参数记为：⊙=(a,B),除此 diag({2h(⊙)T1)和z=（z1,…,zr)',我们得到 之外，再定义一个(T×T)的对角矩阵∑= 了似然函数的近似为： ∑(o)=diagh1(8)证1，这里：h,(8)=co+ (detA)xp](3) a1y2-1-伍-1（⊙），由此，可以写出似然函数如下： 下面会看到，抽取a,B的随机数的分布就基于 l(e1y)o(det∑)nepl-2y∑'y】 上述的似然函数。 为了方便起见，将第一个观察值作为初值，方差 (一)ARCH系数的生成 的初值固定为a0。给予参数a,B如下的先验： 由Chib&Greenberg迭代变换可知，式(2)中的 p(a)cc Na(a l,)I 2,(⊙)可以写为a的线性函数。为了方便，记v:= y2,迭代变换为：l:=1+仪-1，v=-1+1, p(B)cN(BI9,∑BILB>0I 这里0，vg的初值均设为0。记=（v1,…,r）', 这里，4，∑.是超参数，是示性函数，N4是 c:=(l:,v,),构造(T×2)的矩阵C,其第t行是 d维正态分布。此外，假设参数a,B是独立的，即 c。可以看出z=v一Ca,这样可以将α的似然函数 (⊙)=(a)p(B),这样，由贝叶斯法则，可以得到 记为： 参数的联合后验分布为： p(⊙1y)cl(日Iy)p(⊙)。 I(aB)c(detA)-1exp- GARCH(1,1)模型中的方差方程是递推形式 Ca)'A-1(v-Ca)] 13 万方数据
潘海涛，温小霓：基于MCMC方法的GARCH模型参数估计 r，以便从跳跃分布中抽取口。，并从均匀分布中抽 取一个随机实现，决定接受或者拒绝口，。不需要 厂(口I X)的标准化常数，因为这里只利用比率。 此算法的接受和拒绝准则思想如下：1．如果从 B—l到口。的跳跃增加了条件后验密度，则接受p。 作为巩；2．如果这个跳跃降低了后验密度，则以等 于密度比r的概率设定以=口。，否则设定B= 幺一1。这个思想显然是合理的。 HaSting以两种方式推广了Met砌polis算法。首 先，跳跃分布无需对称；其次，跳跃准则修正为： 厂(口。I X)仉(口。l仇一1) ’一，(仇一1 X)∥。(良一l l口．) 这个修正的算法称为Met叫diS-Hasting算法。 下面应用MCMC算法的思想来构造G讯CH (1，1)模型中参数的估计方法。 二、MCMC算法在GARCH(1，1) 参数估计中的应用 新息(innovations)为正态分布的GARCH(1，1) 模型如下： M=e正：忍，￡=1，…，T e。～i以N(0，1) ^f=口o+口1 3，；一l一犀￡一l 这里，口o>0，口1≥0，p≥0。条件方差^。是滞 后一阶平方观察值和滞后一阶方差的线性函数。参 数的正定性的限制保证了方差为正定的。为了写出 似然函数，定义以下的向量：y=(y1，…，了T)7，a= (ao，口1)7，然后将所有参数记为：@=(口，卢)，除此 之外，再定义一个(T×T)的对角矩阵∑= ∑(@)=diag{^1(@)}二1'这里：^。(@)=口o+ 口ly；一l一摩f_1(@)，由此，可以写出似然函数如下： z(@I了)。C(det∑)1尼eXp[-告了’∑～y] 为了方便起见，将第一个观察值作为初值，方差 的初值固定为口o。给予参数口，卢如下的先验： 户(a)。C N2(a I心，∑。)叱>o] p(卢)。C N(J9 J邱，∑日)J[p>o】 这里，∥，∑．是超参数，，[．]是示性函数，N0是 d维正态分布。此外，假设参数口，卢是独立的，即 p(@)=p(a)户(p)。这样，由贝叶斯法则，可以得到 参数的联合后验分布为： 户(@I y)oc Z(@I了)夕(@)。 GARCH(1，1)模型中的方差方程是递推形式 的，这样就无法找到参数的共轭先验分布，因此，借 助于Me仃c∞lis—HaSting(M—H)算法。具体方法 如下： 参数的初值为@[o]=(∥01，∥01)，因为，无法 得到条件分布的显示表示，从参数口，卢各自的分布 进行随机数的抽取，运用(M—H)算法_r次，可以得 到第J次参数的迭代结果为： 口[j】～户(口I∥j一1|，y) ∥～p(卢I口[川，y) 其次，注意到GARCH(1，1)模型可以写为如下 的ARMA(1，1)模型： ^￡=口o+口1y；一l+陋￡一l {jy；=口o+(口l+J9)了；一l—J9础f—l+t‰ ．．2 其中姒=y；一^。=(半一1)^。 =(z}一1)^。 (1) 由{姒}的构造可知，序列{毗}是一个鞅差过 程(Maningale Diffe嘲撇P懒)，可以得到其方差 为2^；。据NakatSuma知，很难直接从式(1)中生成 参数@的随机数。因此，用一个期望为o，方差为2^； 的正态随机变量≈来近似砒。据此，可以构造如下 的辅助方程： y；=口o+(口1+J9)y；一l一盘卜．1+≈ 注意到，砀和^f都是@的函数，即有 ‰(@)=y；一口。一(口l+卢)了；一1+艮卜．1(@) ^￡(@)=口o+aly；一1+肛f—l(@) (2) 再记(T×T)的对角矩阵以=A(@)= diag({2^；(@)}五1)和z=(2l，…，zT)7，我们得到 了似然函数的近似为： z(@I，)。C(detA)-1／2eXp[一告z么．1z] (3) 下面会看到，抽取口，卢的随机数的分布就基于 上述的似然函数。 (一)ARCH系数的生成 由Chib&Q代nberg迭代变换可知，式(2)中的 z。(@)可以写为口的线性函数。为了方便，记仇= 了；，迭代变换为：z，=1+犀0l，口f=仇一1+加工l’ 这里z孑，口彳的初值均设为0。记口=(口1，…，钞丁)7， ct=(zf，口?)，构造(T×2)的矩阵C，其第￡行是 c。。可以看出z=口一巳，这样可以将a的似然函数 记为： z(口I卢，y)∞(det以)一1尼exp[一告(口一 ＆)么-1(口一已)] 13 万方数据
统计与信息论坛 结合此似然函数，由通常的贝叶斯方法即可抽 的结果进行比较。 取a的随机样本： 三、应用MCMC方法的实证分析 (a)c Na(a)It 取上证指数为研究对象，数据全部来自上海证 2.=cac+. 券交易所网站。收益率指数采用对数百分收益率，即 a.=2.(ci+∑h) r:=(lnp-lnp-i)×100,其中，p:和p-1分别是 第t日和第t-1日指数收盘价，数据已进行复权处 这里(T×T)的对角阵A为diag({2h1), 理。为了避免股市的暴涨暴跌对模型拟合造成的影 其中方=a0+a1y2-1+尻，-1，a指的是前一次M- 响，数据的起始点选在1997年以后，即股市实行涨 H抽得的样本。其中某一个a被接受的概论为： 跌停板制度以后，时间跨度为1997年2月日1至 p(a",B8Iy)qa(a",a) 2008年4月14日，上证指数共2717个数据。通过经 min p(a,BIy)qa(a,a") 典的ML方法[4，对于此组数据序列建立新息 (二)GARCH系数的生成 (innovations)基于正态的GARCH(1,1)模型如下： 由前面知道，式(2)中的x,(8)可以写为a的线 %=0.01829+e4 性函数，但是（⊙）不能表示为B的函数。为了克服 e =ai n,n~iidN(0,1) 这个困难，将，（β）在点B处一阶泰勒展开如下： o1=0.08046+0.14048e2-1+0.84044a-1 ()≈()+ .·(B-B) (4) d8 B-B 下面根据上文构造的MCMC算法来重新估计 这里的B指的是前一次的M~H抽样值，采用 上面的模型参数向量a=(ao,a1)'的先验分布中 如下记号： r()=,()+v, 初值取为4，=(0,0y和∑。=diag(1000, 1000),这里方差取为1000，是用了所谓的扩散先 可，=-那 t 验的思想(diffuse prior)。参数B的先验分布中的初 其中的梯度可由如下的迭代来计算： 值取为=0和∑g=1000,其中方差取为1000 V,=y2-1--1(8)+V:-1,初值为V0=0。 的理由同上。 然后，再记r=(r1,…,rT）}',V=（1,…,r)', 为了确保收敛性，模拟了两条马氏链，每条各迭 这样可以将式(3)中的指数近似为： 代1000次，下面是各次迭代完成后的踪迹图(trace exp[-2r-Ya(r-o)】 plot),左边为踪迹图，右边为各个参数的后验核密 度估计(kernel density)。此图可以告知迭代的次数 然后，同样地，β可由通常的贝叶斯方法结合它 是否足够。 自己的先验来生成： 由图1可以看到，在迭代600左右附近，两条链 gA,BCN(p1Rp,Ipn1 已较接近，但仍有波动。再计算由Gelman&Rubin 提出的方差比(variance ratio)法。若R>l,说明链 :=可1g+, 收敛的不好，若R=1,提示链达到了静止状态。计 =2,口ir+∑gm) 算结果为：R0=1.38,R。,=1.55,Rg=1.38,R =1.47,全部大于1，说明两条链收敛效果并不理 这里(T×T)的对角阵A为diag({2h引Ti), 想。 其中，=a0+a1y匠-1+五-1其中某一个B”被接 下面增大迭代次数，迭代3000次，计算的各个 受的概率为： 参数的方差比为：R,=1.00,R。,=1.02,Rg= J(B",a【y)g(β"，B) minpBaly)) 1.00,R.=1.07,可以看到，方差比的结果均接近 于1，说明在迭代3000次时，MCMC算法到达收敛 下面运用上文提出的算法对于基于数据为沪市 (耗时169.30秒，本机配置为Intel CPU1.73GHz, 股票收益率的序列构造的GARCH(1,1)(其新息为 1G内存)。 正态分布)模型的系数进行估计并同经典方法得到 其踪迹图如下： 14 万方数据
统计与信息论坛 结合此似然函数，由通常的贝叶斯方法即可抽 取口的随机样本： q口(二，口)。c N2(a愧，∑。)亿>o] 佥：1：c7三．1c+∑：1 舷=宝。(c7五一·口+∑■) 这里(丁×丁)的对角阵以为d诹({2五；}工1)， 其中元。=二o+二1y；一1+岳一，二指的是前一次M— H抽得的样本。其中某一个口。被接受的概论为： ． {p(口。，星I 3，)q。(口。，口)l 舢n1及i打万i磊万f (二)GARCH系数的生成 由前面知道，式(2)中的‰(@)可以写为口的线 性函数，但是‰(@)不能表示为卢的函数。为了克服 这个困难，将z。(p)在点声处一阶泰勒展开如下： “p)≈砀(声)+雾b邯一声) 这里的声指的是前一次的M—H抽样值，采用 如下记号： r。(卢)=z。(卢)+卢V。 。一 如￡I q一邵J卢：芦 其中的梯度可由如下的迭代来计算： V。=y；一l一≈一1(声)+声VH，初值为Vo=o。 然后，再记r=(rl，…，rT)7，V=(Vl，…，VT)7， 这样可以将式(3)中的指数近似为： e冲[一吉(r一户V)么(r—pV)] 然后，同样地，p可由通常的贝叶斯方法结合它 自己的先验来生成： 邻(声，p)。c N2(p 1邱，∑日)f[J9>o】 奎；1=v7五-1V+∑；1 邱=奎卢(v7芡一·r+∑；k) 这里(T x T)的对角阵A为d谊g({2五；}墨1)， 其中左。=二o+二ly；一l+屏H。其中某一个J9。被接 受的概率为： ．I丛p。，口I y)口巨(卢。，卢)l mlnli万万万i丽f 下面运用上文提出的算法对于基于数据为沪市 股票收益率的序列构造的GARCH(1，1)(其新息为 正态分布)模型的系数进行估计并同经典方法得到 】4 的结果进行比较。 三、应用MCMC方法的实证分析 取上证指数为研究对象，数据全部来自上海证 券交易所网站。收益率指数采用对数百分收益率，即 九=(1nA—lnA一1)×100，其中，A和A—1分别是 第￡日和第￡一1日指数收盘价，数据已进行复权处 理。为了避免股市的暴涨暴跌对模型拟合造成的影 响，数据的起始点选在1997年以后，即股市实行涨 跌停板制度以后，时间跨度为1997年2月日1至 2008年4月14日，上证指数共2 717个数据。通过经 典的ML方法[41，对于此组数据序列建立新息 (innovationS)基于正态的GARCH(1，1)模型如下： M=0．018 29+￡f 屯=盯；以仇，琅～i劫V(0，1) 吼=o．080 46+o．140 48￡；一1+o．840 44吼一l (4) 下面根据上文构造的MCMC算法来重新估计 上面的模型。参数向量口=(口o，口1)7的先验分布中 初值取为心=(o，o)7和∑．=d缸g(1 000， 1 000)，这里方差取为1 000，是用了所谓的扩散先 验的思想(diffuse prior)。参数p的先验分布中的初 值取为脚=o和∑。=1 000，其中方差取为1 000 的理由同上。 为了确保收敛性，模拟了两条马氏链，每条各迭 代1 000次。下面是各次迭代完成后的踪迹图(trace plot)，左边为踪迹图，右边为各个参数的后验核密 度估计(kemd denSity)。此图可以告知迭代的次数 是否足够。 由图1可以看到，在迭代600左右附近，两条链 已较接近，但仍有波动。再计算由Gelman&Rubin 提出的方差比(、枷ance ratio)法。若食>1，说明链 收敛的不好，若R=1，提示链达到了静止状态。计 算结果为：R。o=1．38，食。，=1．55，R卢=1．38，R产 =1．47，全部大于1，说明两条链收敛效果并不理 想。 下面增大迭代次数，迭代3 000次，计算的各个 参数的方差比为：R。。=1．00，R。，=1．02，郫= 1．00，R。=1．07，可以看到，方差比的结果均接近 于1，说明在迭代3 000次时，MCMC算法到达收敛 (耗时169．30秒，本机配置为Intel CPU 1．73GHz， 1G内存)。 其踪迹图如下： 万方数据
潘海涛，温小霓：基于MCMC方法的GARCH模型参数估计 Traee ofalpha0 Density ofalphao 为了更进一步与前面运用极大似然法估计的结 8学 果进行比较，介绍常用的六种预测评判指标。 02004006008001000 0.000.050.100.150200.25 Iterations N-1000 Bandwidth-0.005602 表1MCMC方法对基于正态分布的 Trace of alphal GARCH(1,1)模型系数的估计表 Demsity ofalphal 估计值 标准误 95%置信区间 02004006008001000 0.100.150.200.25 N-1000 Bandwidth-0.004942 0.09816 0.02176 (0.06075,0.1441) Iterations aj 0.14278 0.01961 (0.10768,0.1835) Trace ofbeta Density ofbata 0.836340.02052 (0.79426,0.8721) 02004006008001000 0.700.750800.850.90 N-1000 Bandwidth=0.006179 将股市在第t日收益率减去均值的平方(：一 )2作为股市在该日的实际波动，由模型所计算 Trace of nu Deasity ofnu 的条件方差？作为股市的波动预测，Buhlmann和 0200400600s001000 020 406080100 MeNeil以如下的六种预测误差度量指标，即平均平 N-1000 Bandwidth-0.1774 方误差(MSE1、MSE2)、平均绝对误差(MAE1、 图1MCMC模拟1000次的踪迹图 MAE2)、高斯准极大似然损失函数误差(QLIKE) Trace ofalphs0 Density of alphao 及对数损失函数误差(RLN)作为判断标准： 05001000 2000 3000 0.000.050.100.i50.20025 【ter5oet N-3000 Bandwidth-0.004465 MsE1=2(6-月， Trace of alphal 0500100020003000 0.100.130200.23 MSE2=元之（好-）2 Iterations N-3000 Bandwidth-0.00381 MAE1=121a,-,l, 0500100020003000 0.700.750.800.850.90 Iterationa N=3000 Bandwidth=0.004169 M2=容弟-1 Density of nu Trace ofau 0204060801600 02004006008001000 LE(() Iterations N-3000 Bandwidtb-0.004169 R'LN-1>[tn(ri)P (6) 图2MCMC模拟3000次的踪迹图 nt=l 其中MSE1,MSE2,MAE1,MAE2越小表示预测 利用最后的这个改变初值后的迭代结果，并且 精度越高，QLKE和RLN越大表示预测精度越 舍弃选代的前1500次(bum-in)样本，利用后1500 高。 迭代的样本对基于正态分布的GARCH(1,1)模型 然后运用前面这六种预测指标来评价这两个模 的MCMC方法的系数估计的结果，即方差方程为： 型，计算结果见表2。 a=0.09816+0.14278e-1+0.83634a-1(5) 由上面的计算可知，在各个指标上均显示出基 考察此模型参数的联合分布图（图略），发现参 于MCMC估计的GARCH模型要优于ML方法估 数α1和β的联合分布与椭圆稍有偏差，这预示我们 计的模型。 基于经典方法的推断可能存在问题[5)。 表2基于ML和MCMC方法估计GARCH(1,1)模型的样本期内拟合误差度量指标表 模型 MSEI MSE2 MAEI MAF2 QLIKE RLN ML 0.000769230 0.001347940 0.02369042 0.03353116 0.2980602 -0.0602726 MCMC 0.000565027 0.001049624 0.01607877 0.0249226 0.5393854 0.04185402 四、结论 残差基于正态分布的GARCH(1,1)进行估计，通过 构造一个收敛于待估计参数的后验分布而绕开了 本文运用马尔科夫蒙特卡罗(MCMC)方法对 GARCH模型中参数的约束使得一般优化方法难以 15 万方数据
潘海涛，温小霓：基于McMc方法的G娘cH模型参数估计 Tr_cc of-lDh-O §难辜罩驿醉回 O 2∞400 600 800 1000 I“nboⅡ- Trlct ofIlph^l o■。t’。。。___-____——-_______’一 苫非生■■笠阜蛐I O 200 400 600 8001000 Itmtions N-l000 BtndⅣidth—o．005602 Traccofbct· D删时ofb^忸 !年军=粤=马 O 200 400 600 800lOoo i[二；军=马 0．70 0．75 O．80 O．85 O．90 T．一“一。 N一10∞B-ndwidthlO．O嘶179 Dcmi舛of¨ 詈匹====] 0 20 40 60 80 lOn NzlO∞BandwidthlO．1774 图l MCMC模拟1 000次的踪迹图 詈讴草晒革函 Ttacc ofa】Dh8l O 5001000 2000 3000 Ilemion● 下racc ofheta 蓥犟覃霎霉翠 O 500l000 2000 3000 ItcfnIo¨ N-3000 B蛆dwidthtO．0038l N一3000 B柚dwidnll0．004169 图2 MCMC模拟3 000次的踪迹图 利用最后的这个改变初值后的迭代结果，并且 舍弃迭代的前1 500次(bum．in)样本，利用后1 500 迭代的样本对基于正态分布的GARCH(1，1)模型 的MQ以c方法的系数估计的结果，即方差方程为： 吼=o．098 16+o．142 78e；一1+o．836 34吼一1(5) 考察此模型参数的联合分布图(图略)，发现参 数口1和卢的联合分布与椭圆稍有偏差，这预示我们 基于经典方法的推断可能存在问题[51。 为了更进一步与前面运用极大似然法估计的结 果进行比较，介绍常用的六种预测评判指标。 表1 M口佗方法对基于正态分布的 GARCH(1。1)模型系数的估计表 将股市在第￡日收益率减去均值的平方(^一 ≠)2作为股市在该日的实际波动口；，由模型所计算 的条件方差孑；作为股市的波动预测，Buhlmann和 McNeil以如下的六种预测误差度量指标，即平均平 方误差(MSEl、MSE2)、平均绝对误差(MAEl、 №)、高斯准极大似然损失函数误差(QLIKE) 及对数损失函数误差(R2LⅣ)作为判断标准： MsEl=告∑(五一吼)2， ’’￡=l MSE2=号∑(孑；一d；)2 。’￡=l 脱1=号∑I况一吼I， 。’‘=l 脚2=昙∑I z一盯；l ’’￡=1 QLJⅫ=号∑(1n(蕾)+盯多A ’。t=l R2LN=吉∑[1n(孑搿)]2 (6) 其中MsEl，MSE2，^硝E1，MAE2越小表示预测 精度越高，QLJKE和R2LN越大表示预测精度越 高。 然后运用前面这六种预测指标来评价这两个模 型，计算结果见表2。 由上面的计算可知，在各个指标上均显示出基 于MCMC估计的GARCH模型要优于ML方法估 计的模型。 表2基于ML和M回肥方法估计GARCH(1，1)模型的样本期内拟合误差度量指标表 本文运用马尔科夫蒙特卡罗(MCMC)方法对GARCH模型中参数的约束使得一般优化方法难以 四、结论 蔫耋鬈轰器鬈裟燃霁莩15 n—o n—o一蕴 o∞o— n—o瓣 一 一 万方数据
统计与信息论坛 运用的困难。实证分析结果表明：对于此组数据，基 模型比基于极大似然估计(ML)方法估计GARCH 于马尔科夫蒙特卡罗(MCMC)方法估计的GARCH 模型更接近真值。 参考文献： [1]David Arida.Bayesain estimation of a markov-switching threshold asymmetric GARCH model with student-t innovations[J]. Econometrics Joumal,2008(1):1-22. [2]Teruo,Nakasuma.A Markov-Chain Sempling Algorithm for GARCH models[J].Studies in Nonlinear Dynamics& Econometrics,Berkeley Electronic Press,1998(2):46-53. [3]Teruo,Nakasuma.Bayesian analysis of ARMA-GARCH models:A Markov chain sampling approach[]].Joumal of Eoonometrics,2000(95):57-69. [4]陆懋祖.高等时间序列计量学[M].上海：上海人民出版社，1999：231-302. [5]Gelman A,Carlin J B,Stern HS,Rubin D B.Bayesian data analysis[M].CRC Press:London,1995:222-240. [6]Box GEP,Tiao GC.Bayesian inference in statistical analysis[M].Addison-Wesley:Reading,1973:340-367 [7]Carlin B P,LouisTA.Bayes and empirical bayes methods for data analysis[M].2nd ed.London:Chapman and Hall,2000: 231-302. [8]Hasting W K.Monte carlo sampling methods using markov chains and their applications(].Biometrica,1970(57):97-109. [9]Jonathan D.Cryer,Kung-Sik Chan.Time Series Analysis with Applications in R[M].New York:Springer,2001:331- 342. [10]Tiemey L.Markov chains for exploring posterior distributions[J].Annals of Statistics,1994(22):1701-1762. (贵任编辑：王南丰) The Estimation of GARCH Model Parameters Based on MCMC Algorithms PAN Hai-taol,WEN Xiao-ni2 (1.School of Statistics,Xi'an University of Finance and Eoonomics,Xi'an 710061,China; 2.School of Management,Xidian University,Xi'an 710061,China) Abstract:Traditional method is to use ML method to estimate the parameters,and ML method is essential- ly an optimization method.But GARCH model ordinarily has many constraints among parameters,which result in the failure of trust of MLE results.This paper use Markov Monte Carlo (MCMC)method to estimate the parameters of normal-based GARCH(1,1)model.The results based on MCMC are more reliable and we also show results based on MCMC are better than that of ML based by using real financial data. Key words:GARCH model;MCMC algorithm;Gibbs sampling;Metropolis -Hasting algorithm; volatility;forecasting 16 万方数据
统计与信息论坛 运用的困难。实证分析结果表明：对于此组数据，基 于马尔科夫蒙特卡罗(MCMC)方法估计的GARCH 参考文献： 模型比基于极大似然估计(ML)方法估计GARCH 模型更接近真值。 [1]&丽d Arida．Bayesain estirIlation of ama如、，一s诵tching thr∞hdd asymmenic GARCH r∞dd诵也stud朗t．t inrK舰tionS[J]． Ec删trics Joumal。2008(1)：1—22． [2] Ten的，I、iak姻lrlla．A Md唧一Q面n s啪pling～鲥thm妇GARCH m硇els【J]．Studies in N袖near功咖遗& Ec踟帅etri鸥，酬ey Elea旺H1ic Pre鹞，1998(2)：46—53． [3]T印时，Nal【a双叻a．Bay鹳ian analysis of川RMA—GARCH mDdds：A Marl鼢，chain s锄pling appfoach[J]．．『a鼢a1 0f Ec雠trics，2000(95)：57—69． [4]陆懋祖．高等时间序列计量学[M]．上海：上海人民出版社，1999：23l一302． [5]Gdm锄A，Carlin J B，st锄H S，Rub．m D B．Bay菌锄data锄al酒s[M]．CRC Pr髑：Lond∞。1995：222—240． [6]BOX GE P，Ti∞G C。BayeSi孤infef朗ce irI statiStical anal姆s[M]．Addis∞一wesley：Re8diIlg，1973：340一367． [7]Cdin B P，L∞如T A．Bayes and唧idcal bay∞methods for data analySis[M]．2nd ed．k叽don：(、hpnm and H“，2000： 231—302． [8】HaSting wK．M。nte跚do sampling methods LtsiIlg markov c鼢nS and their applicationS[J]．Biofne伍ca，1970(57)：97—109． [9] Jonathan D．Cfy口，Kung—sik Cha工1．Time S萌皓AnalysiS诵th Applications iIl R[M]．New Y甜【：S研ng盯，2001：331— 342． [10] m咖ey L．～Iarl州c陆ns for唧嫡I】g p06td恼如tributiollS[J]．Anna】s of statiSti心，1994(22)：170卜1762． (责任编辑：王南车) ，nIe Esti眦ti蚰of Q堰CH ModeI Pa憎mete璐Based蚰M例C A120rith瞄 PAN Hai—taol。WEN Xia肛nP (1．Sch∞l 0f stat敏i∞，xi’柚U11iversi够0f R咖ce and Eo咖ics，)(i’锄710061，(Xm； 2．sch001 of M阻ag锄∞t，ⅪdiaJl Urliversi够，xi’锄710061，C扭m) Abstmct：7l、raditiom】method is to use ML method t0 eStinlate the parameters，and M[L meth()d is eSSlential． 1y a11 optimization method．But GARCH model ordimdly h弱many∞nSt商ntS锄ong pa埘neterS，which reSult in the fajlure of trust of Mu三results．This paper uSe Markov Monte Caurlo(MCMC)method t0 eStimate the parameterS of nonTla】一based GARCH(1，1)mOdel．T11e results based on MCMC are rnore reliable and、他also show results baSed on MCMC are better than that of ML b够ed by uSing real finaJlcial da乜． Key wordS：GARCH mDdel；MCMC algonthm；GibbS sarnpling；Met印oIis—HaSting a190rithm； v。latility；fc∞ecasting 16 万方数据
