基于极值po模型的我国生猪疫病风险评估研究

基于极值po模型的我国生猪疫病风险评估研究

随着社会经济的快速发展，中国居民的动物消费显著增加，家禽行业也显著增加。201年,我国肉类产量达7925.8万t,比2005年增长14.2%,连续21年居世界第1位;禽蛋产量2762.7万t,增长13.4%,连续26年居世界第1位;奶类产量374万t,增长30.8%,居世界第3位。我国已发展成为世界畜禽产品生产与消费的第一大国,畜禽产业在改善国人食物结构、增加农民收入及保障食物供给安全方面发挥了越来越大的作用。与此同时,频繁发生的禽流感、蓝耳病、口蹄疫等动物疫病,也极大地影响着食物卫生与公共卫生安全,以及畜禽生产安全。因此,清楚认识我国重大动物疫病风险,并构建有效的动物疫病风险防控体系的任务显得尤为紧迫。总结近些年来国内外对动物疫情风险分析的研究文献,发现学者们对动物疫病的风险分析较多的集中在定性研究方面,而定量研究较少。其中比较具有代表性的是吴春燕等、李静等、王靖飞等利用专家打分法和层次分析法相结合的方法确定疫病的风险因素及权重,对我国疫病风险进行了评估,并且同实际的疫情发展情况进行对比,得出此研究方法具有一定的可行性。而李亮等则利用蒙特卡罗法进行了定量研究,以安徽省生猪为例,针对猪瘟、猪丹毒、猪肺疫等3种主要生猪疫病的风险分布进行拟合,利用蒙特卡洛模拟法得出了各自的风险分布类型及参数,得出的研究结论为Lognormal分布为3种疫病风险的最优分布,安徽省猪瘟、猪丹毒、猪肺疫的死亡风险并不服从于正态分布。相比较国内研究来说,国外对动物疫病风险分析的研究较为成熟,研究方法也相对丰富得多,其中Schlosser等提出利用以往的历史数据可以对疫病的爆发情况进行相关的预测,运用蒙特卡洛模拟法对Z病进行分析,得出Z疫病发生率符合Beta分布,并计算出了风险分布参数,对Z病成本收益进行了评估。Goldbach等通过蒙特卡洛模拟和成本收益相结合方法对荷兰猪肉生产中沙门氏菌的3种防控措施进行了评估,利用蒙特卡洛模拟得到3种措施的NPV分布,最后比较得出最优方案。目前对我国疫病的定量分析主要是利用蒙特卡洛法对疫病分布进行拟合,计算出分布参数,得出疫病风险的最优分布类型。这主要产生了两方面的不足:其一,动物疫病灾害的发生具有低频高损的特点,重大疫病的突发表现为发生率低,但是一旦发生则损失巨大,这也是我们防控疫病的重点所在。而仅仅对历年疫病情况进行简单拟合,往往会忽视重大疫病的发生情况,所得到的最优分布也不能完整地反映疫病巨灾的信息;其二,仅仅得出疫病风险的最优分布在构建简单明了的疫病防控体系方面还稍显不足。为此,本研究利用现代极值方法对我国重大生猪疫病风险进行分析评估,利用极值POT方法对生猪疫病损失尾部进行有效拟合,避免了传统估计方法在分析极端事件方面由于关注与均值而对尾部数据估计不足的缺陷;将风险价值(VaR)引入我国疫病风险评估体系中,实现我国疫病风险的数字化,并计算出不同程度大小疫病风险的95%置信区间,为构建我国疫病风险防控体系提供有益的参考。1数据来源与研究方法1.1基于蒙特卡洛模拟方法的样本数据的分析疫病风险对我国生猪养殖的影响主要表现在生猪感染疫病后所产生的死亡,因此本研究采用生猪感染疫病后产生的最终死亡数的波动形式及大小来衡量疫病的风险。其中最终死亡数的统计选取生猪感染疫病后的自然死亡数以及为防止疫病传染对生猪进行的捕杀数量之和作为最终指标,数据来源于农业部兽医局发布的《兽医公报》,选取的生猪疫病为公报上统计的所有生猪疫病,包括口蹄疫、猪瘟、猪丹毒、猪肺疫、猪水泡病等共9种疫病,样本取值区间为2000年1月至2012年12月,其中有些年份的数据缺失,最终得到了139个样本数据。鉴于资料所限,本研究统计数据仅有139个,样本数据过少,会产生在极大似然估计中误差增大问题,因此采用蒙特卡洛模拟法来增加我国生猪疫病死亡数的样本空间。蒙特卡洛模拟方法也称为随机模拟或随机抽样法,其基本思想是:首先建立一个概率模型或随机过程,使其参数等于问题的解;然后通过对这个模型进行随机抽样观察,即产生随机变量;最后给出问题的近似估计值。其中构建与描述概率过程是其核心。本研究假设我国生猪疫病损失的历史数据服从Beta分布、Gamma分布、Weibull分布、Lognormal分布、Normal分布、Logisitic分布以及Log-Logistic这7种在国内外研究中经常使用的分布类型,对其进行拟合,然后根据计量经济学中的拟合优度检验来判定所选分布对样本数据的拟合情况。最后根据检验结果确定最优分布,对其进行随机抽样生成大量模拟的随机变量,构建我国生猪疫病损失新的样本空间。计量经济学中的拟合优度检验包括AD检验、K-S检验与卡方检验3种检验方法。对分布类型选取的评判方法是:如果3种统计量的检验结构相同,就以该检验的结果为准;如果其中2种统计量的检验结果一致,则以数量多的结果为准;如果3种统计量检验结果均不一致,考虑到AD统计量是最为稳健的检验方法,则选取AD统计量的检验结果为准。1.2基于pot方法的模型分析极值理论是研究随机过程的极值分布及其特征的模型技术,对随机过程中厚尾现象具有突出的针对性,并可在总体分布未知情况下,依靠样本数据外推得到总体极值的变化性质,克服了传统统计方法不能超越样本数据进行分析的局限。目前极值理论中确定极值有两类主流方法:区间样本极大值法(BMM)和超越阈值方法(POT)。由于POT方法在有效性和操作性方面优于BMM方法,因此本研究选取POT方法来构建我国重大生猪疫病损失的概率分布模型。设u为我国生猪巨大疫病所导致的损失,依据现代极值理论的超越阈值模型,位于阈值之上的观测数据将服从广义Pareto分布(GPD)。(1)式为我国生猪重大疫病损失的分布函数,其中,y=x-u,σ>0是尺度函数,k为形状参数。GPD分布实质上是对超过阈值u以上的次序统计量的观察,合理的选取阈值u是POT模型的关键,阈值u值过高,则超过阈值u的样本数过少,这就可能会造成参数估计的方差增大。反之,如果选取的阈值u过低,虽然使得可观测样本数增大,但样本数目的增加将增强样本的中心分布特征,从而造成参数估计的走偏。目前学术界普遍采用平均残差寿命图来确定阈值,如果选取的阈值使得平均残差寿命图趋于线性时,这个值就可以作为样本的阈值。因此,本研究利用平均残差寿命图来确定阈值,并在此基础上,采用极大似然估计法来对广义Pareto分布的参数进行估计,从而获得我国重大生猪疫病损失的分布函数。1.3var的初始定义根据重大生猪疫病风险损失分布模型,运用当前金融市场风险度量的主流方法—风险价值法(简称VaR),对生猪重大疫病风险进行有效度量。VaR是ValueatRisk的缩写,可译为风险价值或在险价值,飞利浦·乔瑞(2005)将其定义为在市场正常波动下,某一金融资产或资产组合的最大可能损失。更为确切的是指在一定置信度水平下,某一金融资产或资产组合在未来特定时间段内的最大可能损失。可表示为:式中,ΔP为金融资产或资产组合在特定持有期内的损失,C为显著水平,VaR为置信水平在1-C下的风险价值。基于VaR的初始定义,在GPD分布中VaR的估计表示为:式中,u为阈值,k与σ分别是GPD中的估计参数,q为相应的置信水平,n为样本数,nu为大于阈值u的样本x的个数。结合(3)式,可计算出我国生猪养殖在遭遇10年一遇、20年一遇、50年一遇甚至100年一遇的重大疫病灾害时的损失大小,从而能够实现对重大生猪疫病灾害风险的有效分析和评估。2中国猪疫情风险的实证分析2.1我国冷却系统模型上的时间序列分析本研究运用生猪最终死亡数来反映我国的生猪疫病风险,根据上文所提到的数据选取方法,最后共得到了139个样本数据,编号为1~139(图1)。从图1可以看出,由于我国生猪疫病在2006年6、7月的疫情比较严重导致的死亡数较多,对这两个极端数据予以剔除,然后对剔除后的样本数据进行平稳性检验。本研究针对样本数据进行ADF检验和PP检验双重稳定性检验,结果见表1。根据检验结果可以看出,我国生猪疫病所导致的生猪死亡数序列在1%、5%、10%的显著水平下拒绝了存在单位根的原假设,由此可以判断改时间序列数据不存在时间趋势,是一组平稳数据,可以对上述时间序列进行下一步的分析。为了反映我国生猪疫病风险的总体特征,进一步对比样本时间序列做简单的统计分析,结果(表2)表明,近几年来我国生猪疫病总体相对较轻,最重的一场生猪疫病发生在2006年,生猪月最高死亡数达到33376头;我国疫病损失的时间序列在5%的显著水平下拒绝了正态分布的假设,这点可以由统计数据简明的看出,该时间序列的峰度明显大于3,呈现出典型的“尖峰厚尾”的分布特征,与重大灾害“低频高损”的特点相吻合,因此完全可以用极值理论的POT方法对样本数据进行分布拟合。由于本研究所得到的统计数据仅有139个,样本数据过小,因此利用上文提到的蒙特卡洛模拟方法来增加我国生猪疫病死亡数的样本空间。表3是我国生猪疫病风险的概率分布检验结果。从表3的拟合优度检验结果可以看出,Gamma分布是我国生猪疫病风险的“最优分布”形式。根据所得到的Gamma概率分布模型,对其进行随机抽样1000次,即进行模拟1000个月,据此产生的已知Gamma概率分布的随机变量并且建立了新的样本空间(图2)。2.2基于u3000提取通过R统计学软件,利用蒙特卡洛模拟法所取得的我国生猪疫病所导致的生猪死亡数新的样本空间集合,运用POT方法将其超过阈值的部分拟合成为广义Pareto分布(GPD),以此取得我国重大生猪疫病灾害的概率分布函数,图3是我国生猪疫病死亡数时间序列的平均残差寿命图。由图3可知,当u>20000时,图形产生剧烈的波动,则证明u值应该位于20000左右,当u=20000时,共有50个数据超过阈值用于判断,数据较少,当u=15000时,其右侧图形为直线,共有116个数据超过阈值,因此暂定u=15000。为了推断u=15000作为阈值是否合适,需要进行更细致的检验,即选择一系列阈值,对不同的阈值利用几大似然估计得到一系列参数值,如果参数估计值在所选阈值附近是稳定的,就说明此阈值合适。在阈值8000~24000均匀地选取100个值作为阈值,用GPD模型估计得到100组参数值(赞σ,赞k),见图4。从图4可以看出,选取u=15000作为阈值是合理的。在取得了我国生猪疫病损失阈值的基础上,采取极大似然估计法对我国生猪重大疫病风险分布函数的参数值进行估计,据此得出了相应的重大疫病灾害概率分布函数:为了对我国生猪重大疫病风险灾害损失概率分布模型拟合的准确性进行进一步的检验,本文将采用经验概率图(PP图)、分位数图(QQ图)、重现水平图和概率密度图对拟合的广义Pareto分布进行诊断比较。从图5可以看出,PP图的点分布基本上都位于对角线的附近,QQ图中右上方的点分布虽然有些偏离,但是由于重现水平图中显示估计的重现水平全部位于重现水平的置信区间内,且密度函数中分布函数的估计同概率密度图的拟合也较好。由此可以得出,本文采用POT方法拟合的我国生猪重大疫病风险损失的概率分布模型是合适的。2.3重大疫情灾害对我国养殖业的影响以及应对根据已经建立的我国生猪重大疫病损失的概率分布函数,运用VaR方法来计算我国生猪养殖在遭遇10年一遇、20年一遇、50年一遇甚至在100年一遇的重大疫病灾害事件下的损失大小,并计算出上述情况下的各自95%置信区间,所得结果见表4。从表4可以看出:(1)我国生猪重大疫病在正常年景的损失(即生猪死亡数)估计在16109头左右,表明我国生猪疫病对我国生猪养殖造成的影响还是比较严重的;而且在遭遇10年一遇和20年一遇的重大疫病灾害时,生猪死亡数就迅速升至3万头左右,如果在遭遇50年甚至100年一遇的疫病巨灾时,死亡数将攀升至4万头左右,这将对我国的畜牧业造成极其重大的影响。(2)不同等级的疫病风险应采取不同的应对措施,根据VaR值的置信区间体系可以构建更加有效的疫病防控体系。例如,我国生猪疫病损失的均值为16109.15,其95%置信区间为(16075.24,16288.62),由于疫病损失越低越好,不存在疫病损失下降的风险,故更关注置信区间上限;如果疫病超过16288.62的置信上限,则疫病损失已超过警戒,需引起高度关注;如果继续上升到27748.05即10年一遇疫病风险95%置信区间的下限,则表明我国已出现了10年一遇的疫病风险,应采取相应的措施。因此根据VaR值的置信区间体系可以很方便的构建疫病防控体系,这对我国畜牧业的稳定发展将产生巨大的作用。3建立冷却压力型养殖业猪马立克氏模型本研究探讨了我国生猪重大疫病灾害风险的评估方法,具体内容包括我国生猪