基于gh分布的极值分布拟合新方法统计与决策第8期.docx
《基于gh分布的极值分布拟合新方法统计与决策第8期.docx》由会员分享，可在线阅读，更多相关《基于gh分布的极值分布拟合新方法统计与决策第8期.docx（15页珍藏版）》请在冰豆网上搜索。
基于gh分布的极值分布拟合新方法统计与决策第8期
一种新的极值分布拟合方法：
基于g-h分布的拟合法
吴建刚
摘要：
对诸如重大危机的极端事件的分析和预测中常常涉及到对极端值的统计分析，这些分析中最主要的部分是对极值分布进行拟合。
这次源于美国的金融危机更是提醒我们极值研究的重要性。
西方学术界对变量极端值的分布拟合的研究相当广泛，其中最主要的方法是被称为“极值理论”（EVT）的方法。
g-h分布是一种能对具有尖峰、厚尾、偏态特征的分布进行很好拟合的分布，还未有人将其专门用于分布尾部的局部拟合，同时g-h分布传统拟合方法是用分位数分别对其参数进行拟合的，这很难做到使四阶矩同时与目标分布一致。
本文提出了g-h分布拟合的蒙特卡罗模拟法并给出了算法，它克服了传统方法的缺点。
这一方法不仅能整个分布进行拟合，在对分布的部分进行拟合时也取得了非常好的效果。
本文首先提出了g-h分布的蒙特卡罗算法，然后利用它进行股票收益率的极端值进行拟合，最后和极值分布拟合方法进行了对比分析。
实证表明，g-h分布的蒙特卡罗算法比极值理论更加方便、灵活和准确。
关键字：
分布拟合；g-h分布；极值理论
中图分类号：
O212文献识别码：
A
Title:
ANewSimulationMethodofExtremeValueDistribution:
SimulationBasedong-hDistribution
Author:
WuJiangang,ManagementSchooloftheShanghaiUniversity
Abstract:
Themailpartofstatisticanalysisintheanalysisandpredictionofextremeevents,suchasvarietiesofgreatcrisis,isdistributionfittingofextremevalues.Thefinancialcrisishasagainremindedustheimportanceofextremevaluestudy.Westernacademiahashadawiderangeofresearchonthisissue,themainmethodiscalledextremevaluetheory（EVT）.Thoughtheg-hdistribution methodhaveitsadvantagesoftakingpeakandasymmetryintoconsideration,theexistingmethodofcalibration,whichusesquantilestocalibrate,donotfittheirparameterswell,whengivesthefourth-ordermomentsofadistribution,.Inordertosolvetheseproblems,thispapergivesaMonteCarlofittingmethod.Evidenceshowsthatthisapproachcanwellfitnotonlyasingledistributionbutalsoapartofadistribution.ThispaperfirstlygavetheMonteCarlofittingmethodofg-hdistribution,thenusedittofitthereturnsofstocks,andfinallycompareditwithEVTandfoundthatg-hmethodisabettermethodfordistributionfittingofextremevalues.
Keyword:
g-hDistribution;Distributionfitting;ExtremeDistribution
在诸如天气预测、危机管理、保险分析、金融风险管理中经常需要对极端值进行统计分析，这次源于美国的金融危机更是提醒我们对于极值的统计分析研究的重要性。
类似分析中往往涉及到对极值分布的拟合，即在难于对整个分布进行精确拟合时，我们常常采用对分布的尾部进行拟合的方法。
西方学术界对极值分布的拟合作过广泛而且持久的研究，主要成果体现在分位数回归和极值理论（EVT）等方法上，其中极值理论的研究和使用最为广泛。
g-h分布最早是由Tukey在1977年一篇工作论文中提出的。
这一分布的优势主要体现在它能对尖峰厚尾及偏态分布进行较好的拟合。
g-h分布传统拟合方法是用分位数分别对其参数进行拟合。
传统方法的明显缺陷是它很难做到使四阶矩同时与目标分布达成一致。
本文提出了g-h分布新的拟合方法，即基于蒙特卡罗模拟的拟合方法，这一方法不仅能对证券收益率的整个分布进行拟合还可以对分布的部分进行拟合。
本文利用尖峰厚尾的股票收益率数据对g-h分布的蒙特卡罗算法进行了验证，并且与极值理论进行了对比分析。
实证表明，g-h分布的蒙特卡罗算法比极值理论更加灵活和准确。
1.g-h分布拟合方法概述
1.1g-h分布的定义
g-h分布是在1977年由Tukey没有发表的论文最早提出，由Hoaglin（1983，1985）[1][2]加以发展的。
这一分布的统计特征也可以从Martinez和Iglewicz（1984）[3]、MacGillivray和Balanda（1988）[4]以及MacGillivray（1992）[5]中找到。
假设Z服从标准正态分布，Tukey通过变换服从标准正态分布变量Z得:
（1）
其中g和h分别代表影响偏度和峰度的参数。
方程右边第一个式子
被称为g分布，可以用
表示，它是影响g-h分布对称性的主要来源；当g趋向零时的极限分布
被称为h分布，用
表示，它是影响g-h分布厚尾性的主要来源。
g分布只允许不对称存在，h分布只允许厚尾性存在，而g-h分布允许两者并存。
由于上述方程没考虑位置参数和尺度参数，可以通过以下方程将其引入：
（2）
其中a、b分别表示位置参数和尺度参数。
1.2g-h分布的特征
通过调整g、h参数可以产生许多其它分布。
正如Martinez和Iglewicz（1984）[3]指出，g-h分布可以产生或拟合很多其它分布。
比如当g、h都取零时可以得到正态分布（通过罗毕塔定理可以求得
）；当g和h分别取1和0时，得到对数正态分布（即
）；当g和h分别取0和0.97时可以得到柯西分布；而不同的g、h取值还可以很好的近似均匀分布、t分布、指数分布、威布尔分布、logistic分布、双指数分布等等。
1.3g-h分布参数的传统估计方法及缺点
常用的估计方法是通过历史数据的分位数来估计g和h。
这一估计时只用了历史数据的少数分位数，这样估计是比较粗糙的；选择什么样的分位数计算并没有客观的标准，这使这一方法主观性比较强；估计结果的好坏没有客观标准，难以进行误差修正。
另外，这一方法要求先后估计a、g、h、b，这样估计出的分布的各阶矩与历史分布是有很大不同的。
这是因为g-h分布的各阶矩与至少g、h、b都有关系，由于每估计下一个参数前都不再对前一个参数进行修正，从而造成估计出来的分布并不能在矩上与被估计分布相符，从而极大了影响了分布的拟合优度。
2.基于蒙特卡罗的g-h极值分布拟合方法
2.1基于蒙特卡罗的g-h拟合法的原理
由于g-h分布是建立在正态分布之上的，只是其均值、方差、偏度和峰度经过位置参数a、尺度参数b、偏度参数g、峰度参数h的调整而使其在均值、方差、偏度和峰度上与正态分布有了不同的矩特征。
能否由蒙特卡罗模拟产生一组正态分布的随机数，然后利用g-h分布调整这个四个参数以便达到均值、方差、偏度和峰度与历史分布完全匹配呢？
匹配后的四个参数的值即为这一历史分布相拟合的g-h的值，只要误差达到一定的精度这样做出来的g-h分布将在均值、方差、偏度和峰度上与历史分布几乎完全匹配。
进一步，如果这一设想能实现，那么只要任意给出一组均值、方差、偏度和峰度的参数就可以计算出g-h分布的四个参数，并且利用g-h分布随机模拟出符合这一条件的历史分布。
2.2完整分布的g-h拟合法的估计步骤
经过实验，可以发现用以下方法是可以估计出历史分布的g-h分布的，而且精度可以控制，步骤如下：
第一步：
输入目标参数
输入历史分布的矩特征，假设均值为m，标准差为v，偏度为s，峰度为k。
第二步：
产生随机正态模拟数
用计算机程序产生一个k维标准正态分布的向量z。
这在许多计算程序中可以很容易做到。
k的大小适当就可以了，太小了估计精度降低，太大了会延长计算时间。
第三步：
预估计g
利用g分布，如下式：
（3）
将正态分布的随机向量z代入上式，并在g常规取值范围内用牛顿算法或用遍历的方法计算y向量，并计算y向量的偏度，直到这一分布的偏度为输入参数s时，这时g的取值为预估计的g值，记为g_temp。
由于历史分布可能为左偏或右偏，g的取值为[-x,x]内，其中x是可以调节的，通常情况下是在正负1之间。
第四步：
预估计h
利用h分布，如下式：
（4）
用与估计g同样的方法可以估计出h。
由于股票历史分布一般峰度都大于3，所以h取值为正。
第五步：
同时逼近a、b、g、h
假设计算均值、标准差、偏度和峰度的函数分别为mean、std、skewness、kurtosis，则第一次估计的a、b分别为：
（5）
（6）
这样就求出了a、b、g、h的初值，这样就可以求出拟合出的g-h分步的随机向量y。
（7）
接下来要用逼近算法，调整这四个值使模拟出的分布y的均值、标准差、偏度和峰度分别逼近目标参数m、v、s、k，直到达到所要求的精度为止。
在实际运算时采用的是万分之一的精度，这已经足够好了。
逼近时需要调整a、b、g、h的取值，先调整g、h的取值。
由于g-h分布里g和h分别控制偏度和峰度，这为逼近带来了方便。
我们用初次估计模拟出的随机变量y的偏度和峰度与目标偏度和峰度进行比较，用这些值来做逼近运算。
具体算法是，新的g′和h′分别取下式的值：
（8）
（9）
这样就可以用迭代的方法，将g′和h′分别代入（6）式以替换原来的g和h，求出新估计出的b′；再将b′、g′和h′分别代入（6）式以替换原来的b、g和h，求出a′；最后，将a′、b′、g′和h′代入（7）式，求出新的g-h分布的随机向量y′。
这样就完成一次迭代运算。
通过这样反复的迭代可以使g-h分布的随机向量的均值、标准差、偏度和峰度逼近目标参数。
2.3计算步骤的相关说明
g-h分布参数估计的常用方法是分位数法，它的明显缺陷是不能保证估算出的分布与原分布在前四阶矩上是一致的。
这里提出的方法是基于蒙特卡罗模拟的方法，虽然经过多次试验都能成功拟合原分布，但为了从理论上说明以上计算步骤的适用性，特作以下说明。
需要首先说明的是，这一计算过程经过多次实验是成功的，但是理论上严格证明将在后续研究中继续进行，这里只是对相关理论问题作一个简单说明：
（1）关于第三步与第四步预估g和h的方法的说明
由于采用蒙特卡罗模拟法是对最先产生的正态分布的向量进行修正以适合需要模拟的目标分布，这里主要采用迭代逼近的算法，所以为了减少迭代次数，先利用g分布和h分布对g和h进行预估计是必要的。
（2）关于为什么先预估g、h再预估a、b的说明
g-h分布中的a是对分布的均值的主要决定因素，b是分布的标准差主要决定因素，g和h分别是对偏度和峰度的主要决定因素，但是由于a、b、g、h的变化对所使用的向量y′的均值、标准差、偏度和峰度都同时有影响，所以，只能使a、b、g、h同时调整以达到前四阶矩与原分布一致的目的。
在计算时，由于通过对a和b的增大或缩小很容易调整该向量使其与目标均值m和标准差v相等而又对偏度和峰度影响不大，所以先预估g、h，然后再预估a、b的方法将更加容易计算。
计算步骤中体现了这一点。
（3）关于第五步是否能收敛的说明
试验设想是通过（8）式和（9）式每修正一次就使y′的偏度和峰度更接近目标值s和k，这一结果得到实验的证实。
（4）关于在利用g-h分布对蒙特卡罗产生的正态分布的向量进行修正过程中产生的a、b、g、h形成的g-h分布是否在前四阶矩上与原有分布相等的说明
由于a、b、g、h形成过程，即是通过调整正态分布并使之与目标分布在前四阶矩一致的过程，而且这一计算步骤是可逆的，即利用估计出的g-h分布的参数a、b、g、h可以产生一个其各阶矩与目标值相等的向量。
2.4基于蒙特卡罗的g-h尾部分布拟合方法计算步骤
以下步骤可以对尾部进行拟合：
1.选择用于模拟的尾部
从第一部分分析可知，由于收益率分布可能是多峰的，这就使单峰的g-h分布很难会处处拟合得准确，但这可以通过分段拟合来达到目的。
对于尾部的拟合，可以选择整个分布的一定比例进行拟合，只要所选择部分是单峰的就可以取得比较好的拟合效果。
2.用g-h分布产生模拟数
计算出所选部分的均值、标准差、偏度和峰度，并使g-h分布与之匹配，通过g-h分布模拟出随机数。
3.计算分位数和下尾期望
由于这里计算的分位数是整个分布的分位数，所以要把尾部拟合算出的分位数与整个分布的分位数进行转换。
假设尾部占整个分布的1/n，那么尾部分布的q×n分位数即为整个分布的q分位数。
同时，由于只使用了整个分布的1/n进行拟合，所以模拟出的g-h分布只能估计置信度小于1/n的分位数，即q×n必须小于等于1。
计算出分位数后，这时就可以计算所有小于分位数的数的平均数，假设与分位数相等的模拟数不超过两个（这是通常的情况），这个平均数即为这一置信度下的下尾期望（ExpectedShortfall，缩写为ES，下同）。
3.极值理论方法介绍
3.1极值理论的提出
研究极值分布的理论称为极值理论（ExtremeValueTheory，EVT）。
经典的极值理论起源于1928年Fisher和Tippett的研究。
近些年来，极值理论有了迅速发展，并引起重视。
在可靠性工程领域及人寿保险中，用来计算产品及人员的寿命；在环境监测中，用来估计降雨量，洪峰流量，水库容量等；在金融领域中，用来做风险分析。
一般来说，极值的精确分布难以估计，主要研究的是其渐进分布。
极值理论主要包括两部分：
经典极值理论和近几十年来发展起来地POT理论。
经典极值理论又称BMM模型（BlockedMaximaMethod）用于对大量独立同分布的样本分组后的极大值建模；POT理论（PeaksOverThreshold）用于对超过给定阈值（threshold）的样本数据建模，从而描述分布尾部的特征。
由于BMM方法在使用数据时把数据分成多段，在数据不足的情况下数据使用效率不高，这里采用POT方法。
3.2POT方法介绍
1.定义
假设序列{zt}的分布函数为F（x），定义Fu（y）为随机变量Z超过阈值u的条件分布函数，可以表示为：
（10）
根据条件概率公式我们可以得到：
（11）
定理：
（Pickand（1975）[6]，Balkema和Haan（1974）[7]）对于一大类分布F（几乎包括所有的常用分布）条件超量分布函数Fu（y），存在一个
使得：
（12）
当ξ≥0时，y∈[0，∞）；当ξ<0时，y∈[0，-σ/ξ）。
图1为广义Pareto累积分布在σ=1，ξ取0.3，0，-0.3的图形。
图1广义Pareto累积分布（横轴：
y；纵轴：
概率值）
从图1可以看出，ξ可以看作形状参数（σ这里是尺度参数）。
当ξ>0时，就是通常的分布Pareto分布，这种分布具有厚尾的特点，在金融时间序列分析中备受关注；当时ξ=0时是广义指数分布，它具有正常的尾部。
当ξ<0时，称为ParetoII型分布，具有较薄的尾部。
实际上广义分布的尾部具有与广义极值分布相同的性质，因此可以利用超过量来估计极值分布的尾部。
分布函数
被称作广义的Pareto分布。
根据定理我们可以得到广义的Pareto分布的概率密度函数
：
（13）
2.分布函数
在确定阈值后，利用{zt}的观测值，同时可以得到比阈值u大的个数，记为N，根据前面的条件概率公式，用频率代替F（u）的值，可以得到F（z）在z>u时的表达式：
（14）
3.分位数的计算
对于给定某个置信水平p，可以由分布函数F（z）得到分位数：
（15）
3.3阈值的确定
在利用最大似然法估计参数，但在估计之前，先要确定定理中的阈值u。
它是正确估计参数σ和ξ的前提。
如果阈值u选取的过高，会导致超额数据量太少，使估计出来的参数方差很大；如果阈值选取的过低，则不能保证超量分布的收敛性，使估计产生较大偏差。
有学者（Danielsson和Vries（1997）[8]以及Dupuis（1998）[9]）给出了对阈值u的估计方法，有两种：
一是根据Hill图。
令
表示独立同分布的顺序统计量。
尾部指数的Hill统计量定义为：
（16）
Hill图定义为点
构成的曲线，选取Hill图形中尾部指数的稳定区域的起始点的横坐标K所对应的数据Xk作为阈值u。
计算出阈值后可以利用Hill给出的一阶矩估计形式对ξ进行估计，方程为：
（17）
其中k表示超过阈值的样本数。
这一估计的缺点是对阈值有依赖性。
二是根据超额限望图。
令
表示独立同分布的顺序统计量。
样本的超限期望函数定义为：
（18）
超限期望图为点（u,e（u））构成的曲线，选取充分大u的作为阈值，它使得当x≥u时e（x）为近似线性函数。
4.实证分析
由于股票收益率大起大落，往往有尖峰厚尾及偏态的特征，用股票收益率数据对以上理论作实证检验是合适的。
以下，首先对数据进行说明，其次使用极值理论进行了拟后，然后利用了g-h分布进行了拟后，最后对各种方法的效果进行了对比。
4.1数据相关说明
1.选样说明
数据采用了中国股票市场交易数据库（CSMAR）股票行情数据，数据期间为1996年1月1日到2008年6月30日之间的日、周、月的股票收盘价格，所有的价格都采用了向前复权复息处理。
收益率用价格的对数差计算而得。
考察的数据应该越多越好，特别是考虑周和月的收益率时，数据量就会减少很多，而采用1996年开始的数据的原因是，1996年后的中国股票市场相对来说从初创走上正轨，用这期间的数据进行计算比较有实际意义。
2.收益率计算方法说明
计算当日收益率的一般方法是rt=ln（Pt）-ln（Pt-1），其中r表示收益率，P表示收盘价，下标表示时间。
但在实证中一般采用当日收盘价的自然对数减去前一天收盘价自然对数的方法，即Rt=（Pt-Pt-1）/Pt-1。
这是由于两者在数量上差别不大，而且也不存在系统性的偏差，用对数收益率时间序列就会显得方便和合逻辑一些。
在方便性上的考虑主要是基于对数收益率在时间上可加性，即多期收益率可以由单期收益率直接加总取得的，而简单收益是需要用复利计算公式来取得，即其计算要采用连乘的方法。
合逻辑是因为时间上可加的收益率序列在时间序列处理上更加合理一些。
3.静态矩分析
计算得，均值为0.092225，标准差为1.7115，偏度为-0.058355，峰度为8.2675。
可见样本是严重左偏厚尾的，这一分布是极为偏离正态的。
4.2利用极值理论进行尾部拟合
1.估计阈值
做出上证指数Hill图，如图2所示。
图2上证指数Hill图（横轴：
由大到小排列的收益率的序号；纵轴：
hill统计量）
从上图可以看出，最小的105个收益率为阈值，因为它是稳定的起始点。
这一点对应的收益率为-2.9226%。
做出超限图，如下图所示：
图3上证指数超限期望图（横轴：
阈值；纵轴：
超限期望函数值）
由上图可知，取充分大的u（为x轴为上取值为0.028，即阈值为收益率为-2.8%时）使u和e趋向于线性关系。
这一方法选择的阈值与Hill图基本一致。
2.估计广义帕累托分布的尺度参数和形状参数
当u确定以后,利用{Xt}的观测值，根据最大似然估计得到ξ和σ。
最大似然函数
如下：
（19）
其中yi表示极值与阈值的差。
由于计算复杂，这里先利用Hill对ξ的一阶估计计算出ξ，再估计σ。
Hill还给出了方程对ξ的估计再通过最大似然函数估计。
通过计算可知ξ=0.41883，σ=0.0138。
4.3利用g-h分布模拟随机向量
1.全部拟合
将全部2900个数据进行拟合，模拟出5000个g-h分布的随机数。
用g-h分布四参数法，拟合结果是：
a=0.0017469，b=0.73529，g=-0.0046535，h=0.054883。
由于每一次模拟所使用的标准正态分布模拟数的数值是不一样的，所以迭代次数不一样，大约在30次左右。
2.尾部1/4拟合
这里选择分布的损失的前1/4的数据作为模拟对象，模拟出的数据是2000个。
通过计算可知，其均值为-1.8772，标准差为1.348，偏度为-2.6511，峰度为11.92，这一分布也是极为偏离正态的。
用g-h分布四参数法，拟合结果是：
a=-1.4213，b=0.82851，g=-0.91027，h=0。
大约迭代15次可以得到结果。
4.4g-h极值分布拟合法与极值理论方法的对比分析
有了分布函数就可以对相关指标进行计算了。
由于计算出多个分位数并观察与实际数据的离异情况就可以看出由其计算出的ES是否准确，所以这里只计算分位数，并将其与正态假设下以及实际的分布的分位数一起对比。
如下表所示：
表1上证指数三种分位数方法比较
实际数值
正态分布
极值分布
全拟合g-h
1/4拟合g-h
95.00%
2.66%
2.82%
2.32%
2.53%
2.52%
95.50%
2.77%
2.90%
2.44%
2.71%
2.66%
96.00%
2.92%
3.00%
2.58%
2.85%
2.82%
96.50%
3.07%
3.10%
2.75%
3.03%
2.98%
97.00%
3.33%
3.22%
2.96%
3.19%
3.17%
97.50%
3.47%
3.36%
3.22%
3.35%
3.46%
98.00%
3.77%
3.52%
3.57%
3.69%
3.75%
98.50%
4.12%
3.71%
4.08%
3.95%
4.18%
99.00%
4.78%
3.98%
4.90%
4.42%
4.73%
99.50%
6.26%
4.41%
6.68%
5.71%
6.07%
从表可以看出，极值分布比正态分布好，而g-h分布又比极值分布好。
以下通过图来比较直观的展示三者的区别。
如下图所示：
图4上证指数三种分位数方法对比图（横轴：
百分比；纵轴：
分位数）
图中最上面的直线为真实值，“短线”虚线为正态分布的值，“短线加点”虚线为极值分布的值，点为全部拟合的g-h分布的值，星号为部分拟合g-h分布的值。
从上图中可知，g-h分布有极出色的拟合效果，而部分拟合法又
升级会员 