基于经验似然贝叶斯计算的稳定分布参数估计

钱 瑾，钱夕元

（1.携程旅游网络技术有限公司；2.华东理工大学 理学院，上海 200237）

0 引言

金融数据分布不仅仅呈现简单的正态分布，大多具有“高峰厚尾”的特征，这是由于金融分布出现异常值的可能性更大，而单纯删除这些异常值是不科学的，会恰恰删除了数据中可以反映金融风险的部分。稳定分布当参数α＜2时方差无限，可以拟合方差较大甚至无限的数据，恰好符合金融数据的“高峰厚尾”特征。稳定分布已经被证明对于金融数据具有良好的拟合性：1960年稳定帕累托分布被Mandelbrot[1]用于拟合金融市场收益率截尾数据，且被证明较正态分布可以更好地拟合波动价格分布。Mantegna和Stanley（1995）[2]提出，方差无限的稳定列维分布可以拟合纽约股票交易所S＆P500指数高频数据的收益分布，且这些数据是尾部为近似指数下降的，并且位于中心（数据均值）6个方差范围内。

稳定分布又被称为列维稳定分布（Lévyalpha-stable distribution），它由列维于1920年提出，是一类没有显式密度函数表达式的连续概率分布。常见的一些分布如正态分布、柯西分布均是它的特例。广义中心极限定理指出，若大量标准化后的独立同分布随机变量之和的极限分布存在，则为稳定分布族中的一种，即稳定分布是唯一吸引场的分布。稳定分布的这一特殊性使得它被应用在各类实证分析中。稳定分布的特征函数有4个参数，可以刻画不同尾端厚薄、偏度、峰度特征的分布，具有较强的灵活性。但是稳定分布没有显式表达式，其分布参数估计存在一定困难。早期一些稳定分布参数估计方法包括：Du-Mouchel提出的极大似然法[3]、Nolan的直接积分法极大似然[4]、Zolotarev的特征变换法[5]以及McCulloch的分位数法[6]。可这些方法通常具有较大的局限性。顾娟和茆诗松[7]类似模拟矩法(SMM)的思想，提出并验证了一种可以得到强相合结果的参数估计方法。武东和汤银才[8]详细介绍了稳定分布的一些性质，并证明了稳定分布在金融风险度量中应用的有效性。Peters[9]利用近似贝叶斯计算（Approximate Bayesian Computation，ABC）给出一元及多元稳定分布的参数估计方法。

本文采用基于经验似然的贝叶斯计算（Bayesian Computation with empiricallikelihood，BCel）[10]对稳定分布进行参数估计：不仅通过模拟数据实验验证方法的有效性，还对2004年1月至2014年12月上证指数对数收益率数据对比正态分布进行实例研究，并利用稳定化的PP图方法进行拟合优度的检验。同时，参考Pradhan[11]中基于Gibbs抽样的贝叶斯预测的思路，本文提出基于经验似然的贝叶斯预测，并利用该方法预测2015年1月至2016年8月上证指数收益率特征分布情况。

1 稳定分布及基于经验似然的贝叶斯计算

稳定分布的定义分统计定义与参数定义两种，稳定分布的参数定义是指稳定分布特征函数的参数表示，具体如下：

稳定分布一般需要四个参数α，β，γ，δ来描述。α是特征参数，也称稳定化指数，是4个参数中最重要的参数，它可以刻画尾端的厚薄程度，α越小峰越高，尾越厚。β为偏度参数，反映了分布的峰偏离分布均值的方向与程度，β＞0则分布右偏，β＞0则分布左偏。γ为刻度参数，而δ为位置参数。当α＞1时，稳定分布均值才存在，并且此时δ就是分布的均值。当α=2，β=0时稳定分布就是正态分布；当α=1，β=0时稳定分布则为柯西分布。若X～S(α，β，γ，δ)，则 [(X-δ)/γ]～S(α，β，1，0)。Chambers等[12]给出的基于S(α，β，0，0)稳定分布随机数产生方法便是基于此种变换。

本文提出利用BCel对稳定分布进行参数估计的方法。BCel是基于经验似然的一种贝叶斯计算方法，它的基础是近似贝叶斯计算方法。Owen在1988年提出经验似然，它是一类非参数似然，并且具有很多统计学的优势。首先经验似然构造置信区间有变换不变性、域保持性，而且置信域的形状也是由数据决定。经验似然的构建过程中，首先需要定义一个数据分布相关的功能函数θ：如分布的矩。具体来说，设独立同分布）为服从f(x)的随机变量，功能函数θ满足限制条件EF[h(Y，θ)]=0 ，则经验似然为p是集合中的值，而h的维数是θ的限制条件的个数。假设功能函数θ=Ef[Y]，在限制p1y1...pnyn=θ的条件下，经验似然即为p1...pn的乘积的最大值。经验似然已经被证明具有良好的收敛性，BCel是一种利用经验似然改进近似贝叶斯计算得到的一种算法。具体的算法逻辑如下：

for i=1:M

（1）由先验分布 π(·)抽样，得到参数θi；

（2）由θ模拟抽样得到对应功能函数；

i

（3）计算经验似然Lel()，并作为权重wi赋给θi。BCel的计算结果是M个权重为wi的参数θi，g(θ)的估计值为当i的值越大，越接近于g(θ)的真实值。BCel方法可以在参数估计中充分考虑参数的先验信息，而且计算较快较简便。利用BCel对于稳定分布进行参数估计时，类似McCulloch[6]的利用抽样的分位数来做估计的方法，本文将功能函数设为分位数，去对稳定分布进行参数估计。

2 基于经验似然的贝叶斯预测方法

Pradhan[11]利用Gibbs抽样与Monte Carlo模拟法对符合二参数的gamma分布的数据的未来观测值进行了预测。借鉴Pradhan[11]的预测思路，本文提出一种基于经验似然的贝叶斯预测方法。

设历史数据y与未来观测值y*，y*的后验预测密度函数为：p(y*|y)= ∫p(y*|y，θ).p(θ|y)dθ，p(θ|y)在历史数据y和先验条件下的θ的后验密度函数。因为y*，y独立同分布，所以设（m为数据长度），则m个顺序统计量为设次序统计量y*(r)（第r个）在已知数据y下的后验密度为则对于稳定分布来说是对应稳定分布的概率密度，而很难解析表示，故结合Monte Carlo模拟利用抽样得到的估计。

以上已经介绍过，以经验似然为似然的后验密度是p(θ|y)，而BCel的结果是M个有权重wi的参数θi，加权抽样得到的参数θ值即为估计值。这一过程又叫做Monte Carlo模拟，同样所 以的估计可以由以下抽样过程获得：

（1）利用参数先验获取参数 (θ1，θ2，...θn) ，并利用Chambers等[12]的稳定分布随机数产生方法产生长度为m的n组随机数

（2）计算n组随机数的r个次序统计量

对于第r个顺序统计量，其双侧100β%置信区间为下式，其中L代表上界，U代表下界：

得到估计。

3 针对稳定分布的BCel参数估计及模拟实验结果

因为α，β，γ，δ是严格限制定义域的，所以需要通过变换：ξ1=qnorm①qnorm:这里指R中stats包中函数，qnorm(x)代表取概率x在正态N(0,1)中对应的分位数；(α-1)②α的定义域为（1,2）；，ξ2=qnorm((β+1)/2)，ξ3=ln③ln(x)代表取底数为e的x的对数。将这4个参数转变为定义域为(- ∞，+∞ ) 的1×4的参数矩阵ξ。因为希望利用Chambers等[12]给出的基于S(α，β，0，0)稳定分布随机数产生方法产生稳定分布随机数据，所以需要利用X～S(α，β，γ，δ)，则 [(X-δ)/γ]

下面讨论利用BCel对稳定分布进行参数估计的具体方法及模拟实验结果。～S(α，β，1，0)这一变换处理参数。而利用以上2种技巧就可以利用模拟产生需要的随机数。

以上介绍过，X～S(α，β，γ，δ)，则 [(X-δ)/γ]～S(α，β，1，0)这一稳定分布的标准变换。又Famma（1971）年提出可以利用数据的分布特征直接得到参数γ，δ的估计值。所以利用BCel对于稳定分布进行参数估计的思路两种，一是利用稳定分布标准变换得到S(α，β，1，0)，利用BCel估计参数α，β的值，而γ，δ的估计值Famma（1971）的方法得到；二是利用BCel直接估计稳定分布四个参数α，β，γ，δ的值。下面先介绍第一种思路的具体算法过程，假设原始数据集y服从稳定分布S(1.7，0.9，10，10)，且y的数据长度为200，利用标准化变换z=[(y-δ)/γ]可得z，即z～S(1.7，0.9，1，0)，对于变换后得到的数据集z计算累计概率密度为 0.1，0.2，...,0.9 的分位数qz1，qz2，...,qz9。设9×100的矩阵当功能函数θ′=(qz1，qz2，...qz9)时的经验似然可以定义为：在限制p1Nz1...pnNzn=θ的条件下，p1...pn的乘积的最大值。由稳定分布标准变换与ξ1=qnorm①qnorm:这里指R中stats包中函数，qnorm(x)代表取概率x在正态N(0,1)中对应的分位数。(α-1)②这里α的定义域为(1 ，2] 。，ξ2=qnorm((β+1)/2)这一变换，知参数α，β，γ，δ(1.7,0.9,1,0)对应ξ(0.52，1.64，0，0)，取正态分布N(0，1)与N(1，1)作为参数ξ1与ξ2的先验。利用以下步骤：

for i=1:10000：

（1）从先验分布N(0，1)与N(1，1)中取，；

（2）由ξ1i，ξ2i利用变换得到αi，βi，再利用Chambers的随机数产生方法，产生长度为200的随机数mzi；

（3）计算9×100的矩阵Nmz为(t=1，..，9，j=1，..，200)；

（4）利用R包‘emplik’求出时对于Nmz的经验似然Lel(θ′|mz)；

（5）对于参数αi，βi，设其权重wi=exp(log(Lel(θ′|z)))-

图1的箱线图在利用以上算法进行50次重复实验后获得。图1箱线图表示此方法的估计效果是良好的：α与真值很接近，而β的估计基本吻合参数真值，并且50组模拟实验得到的参数α，β的估计值方差基本均小于0.1。若数据长度n=500，利用以上方法得到的α与β的估计值为1.72，0.78，具体的累计概率密度函数分布图见图2，可以看出此时拟合分布与真实分布更为接近。

图1 50次重复实验得到的稳定分布参数估计结果

图2 n=500时的累计概率密度分布拟合图

下面讨论不利用稳定分布标准化变换，直接对稳定分布利用BCel进行参数估计的方法。依上取稳定分布S(1.7，0.9，10，10)，原始数据集y长度为200，对于稳定分布直接利用BCel进行参数估计不同与之前的算法过程的几点主要在于：

（1）不需要将y标准化转为z，进而求z的分位数的值。而是直接求y的概率为0.1，0.2，...0.9的分位数qy1，qy2，...qy9，并在计算经验似然的过程中将其作为功能函数θ′；

（2）通过变换ξ1=qnorm(α-1)，ξ2=qnorm((β+1)/2)，ξ3=ln③ln(x)代表取底数为e的x的对数。(γ)，ξ4=δ，将参数α，β，γ，δ(1.7，0.9，10，10)转化为对应的ξ=(0.52，1.64，2.30，10)，并且对于 4个参数2，ξ3，ξ4分别取N(0，1)，N(1，1)，N(0，2)及N(0，10)作为先验。

除此之外过程大体与只估计2参数的相似，表1表示各个不同的方法得到的四参数稳定分布的参数估计值，其中第一列至第四列分别表示Buckle（1995）[13]利用Gibbs抽样，McCulloch（1986）[6]利用分位数法，Lombardi（2007）[14]利用反向MCMC法，本文利用BCel法，在10次重复实验的基础上得到的结果。

表1 四参数稳定分布的估计结果

表1数据表明：相较与其他三种方法，BCel对于参数α，γ，δ的估计较为精准，尤其对于参数α，δ，方法BCel的估计精度是最高的。可是BCel对于参数β的估计有些偏差，只能做到与其他方法持平。但总体来说BCel的估计结果的总体标准差更小，估计效果更加稳定。

4 实证

本文对于2004年1月2日至2014年12月31日的上证综合指数日收盘价数据（共2640个）利用BCel方法进行了分析。假设时间序列{Pt}对应这2640个数据，应用对数收益率变换：yt=100·(lnPt-lnPt-1)，可得序列{yt}。{yt}的均值为0.029，方差为2.923，序列的偏度与峰度为0.877与22.376，满足上文阐述的金融数据的“高峰厚尾”的特征。利用上文说明的2种BCel估计方法对于序列进行分布拟合。首先说明第一种方法的应用过程：采用数据均值，Famma（1971）的估计方法得到参数γ，δ的估计，利用BCel估计参数α，β，且这一过程中采用的先验分布，模拟次数均与之前的模拟实验相同。这一方法得到的参数估计值为(1.88,0.12,0.84,0.029)。而采用直接利用BCel估计4个参数的方法得到的参数估计结果为(1.56,0.122,0.837,-0.12)。同时，采用正态分布对比进行拟合，得到的正态分布为N(0 .029，2.913) 。

PP图是检验分布拟合数据优度的一种重要的方法，可是传统的PP图存在很多问题：采用的均匀分布的次序统计量方差不稳定，且波动较大；曲线尾部无法准确反映对数据的拟合优度等。Michael针对PP图的此类缺陷提出了稳定化的PP图，稳定化的PP图采用均匀分布U经变换后得到的分布的次序统计量，这些次序统计量的渐进方差相同，且上下波动较小，故稳定化的PP图可以更好地进行拟合优度的检验。图3是利用稳定化的PP图对于上段得到的估计结果：稳定分布S(1.88，0.12，0.84，0.029)（标准化后稳定分布BCel拟合结果），S(1.56，0.122，0.837，-0.12)（四参数稳定分布BCel拟合结果），正态分布N(0 .029，2.913)进行拟合优度检验后的结果。图3表明正态分布的拟合优度远不及其他2类分布。而稳定分布标准化后的估计结果S(1.88，0.12，0.84，0.029)较稳定分布直接估计的结果S(1.56，0.122，0.837，-0.12)对于数据的拟合优度更好。

图3 PP图对于3类分布拟合检验结果

为了进一步进行数据分布拟合检验，图4绘制了数据的频率直方图及估计得到的正态分布N(0 .029，2.913) ，稳定分布S(1.88，0.12，0.84，0.029)的密度函数图。图 4表明稳定分布S(1.88，0.12，0.84，0.029)对于数据的拟合程度与数据特征的表现程度远远优于正态分布，这不仅说明了稳定分布可以更好地拟合金融数据“高峰厚尾”的特征，也说明BCel方法可以较好地估计稳定分布参数的值。基于稳定分布S(1.88，0.12，0.84，0.029)对数据的较优的拟合性，可以得到数据本身的一些特征：得到的稳定分布参数α估计值小于2说明上证指数对数收益率存在相关性，并不满足随机游走；而得到的参数估计值β大于0说明上证指数对数收益率呈现右偏的厚尾分布。

图4 上证指数收益率的频率直方图及正态分布、稳定分布密度函数图

由于时间序列数据受周期性的时间因素影响较大，为了对2015年1月至2016年8月的上证指数对数收益率进行预测，历史数据选取2012年1月至2014年8月的上证指数对数收益率数据。利用前文给出的基于经验似然的贝叶斯预测方法与4参数的稳定分布的BCel方法，利用历史数据预测了未来数据的一些特征的分布情况。图5是预测的未来数据的最小值、中位数、最大值的概率密度图，预测的最小值、中位数、最大值的95%的置信区间为(-46.01，-6.482 )，(- 0.53，0.51) ，(6 .31，46.74 )。而 2015中位数、最大值的真值(- 8.87，0.15，5.60 )用虚线标出，可以看出最小值与中位数的真值落在预测的最小值，中位数的95%的置信区间内，预测结果较好。而最大值真值略小于预测的区间的最小值，应该是因为2012年1月至2014年8月的数据的“高峰”特征更明显，所以预测结果存在偏差。年1月至2016年8月的上证指数对数收益率的最小值、

图5 预测数据的最小值、中位数、最大值的概率密度图

5 结论

本文介绍了基于经验似然的贝叶斯计算方法，并利用该方法对稳定分布进行参数估计并且得到良好的模拟实验参数估计效果。同时本文借鉴Pradhan的基于Gibbs抽样的贝叶斯预测方法的思路，提出一种基于经验似然的贝叶斯预测方法，并从理论层面解释了该方法。在上证指数对数收益率的实证研究中，先是利用BCel方法结合稳定分布拟合数据，其拟合结果不仅证实了稳定分布较正态分布更适合拟合“高峰厚尾”特征的数据，也证实参数估计方法的准确性；之后又通过历史数据预测未来数据的特征情况，得到了符合预期的结果。

参考文献：

[1]Mandelbrot B.The Pareto-Levey Law and the Distribution of Income[J].International Economic Review,1960，（1）.

[2]Mantega R M,Stanley H E.Scaling Behaviour in the Dynamics of an Economic Index[J].Nature,1995，（367）.

[3]DuMouchel W.Stable Distributions in Statistical Inference:Informa⁃tion from Stably Distributed Samples[J].Journal of the American Sta⁃tistical Association,1975，（70）.

[4]Nolan J P.Numerical Computation of Stable Densities and Distribu⁃tions[J].Comm.Stat.Stochastic Models,1997，（13）.

[5]Zolotarev V M.One-dimensional Stable Distributions[J].Translations of Mathematical Monographs,1986，（65）.

[6]McCulloch J H.Simple Consistent Estimators of Stable Distribution Parameters[J].Comm.Stat.Simulation and Computation,1986，（15）.

[7]顾娟,茆诗松.稳定分布的参数估计[J].应用概率统计,2002，（18）.

[8]武东，汤银才.稳定分布及其在金融中的应用[J].应用概率统计,2007，（23）.

[9]Peters G W,Sisson S A,FanY.Likelihood-free Bayesian Inference Forα-Stable Models[J].Computational Statistics and Data Analysis,2012，（56）.

[10]Mengersen K L,Pudlo P,Robert C P.Bayesian Computation via Em⁃pirical Likelihood[J].Proceedings of the National Academy of Sci⁃ences,2012，（110）.

[11]Pradhan B,Kundu D.Bayes Estimation and Prediction of the two-parameter Gamma Distribution[J].Journal of Statistical Com⁃putation and Simulation,2011，（81）.

[12]Chambers J,Mallows C,Stuck B A.Method for Simulating Stable Random Variables[J].Journal Am.Statist.Assoc.,1976，（71）.

[13]Buckle D J.Bayesian Inference for Stable Distributions[J].Journal of the American Statistical Association,1995，（90）.

[14]Lombardi M J.Bayesian Inference for Alpha Stable Distributions:A Random Walk MCMC Approach[J].Computational Statistics&Da⁃ta Analysis,2007，（51）.