蒙特卡罗法模拟

2024-11-21

蒙特卡罗法模拟（精选7篇）

蒙特卡罗法模拟篇1

一、引言

作为资产评估的基本方法之一,收益法是依据资产未来预期收益经折现或资本化处理来估测资产价值的,其多数参数都须面向未来进行预测,所以对资产未来预期收益的合理预测成为收益法的核心内容之一。在对资产未来预期收益的预测过程中,评估师对评估对象的价值描述不仅要以一些基本假设为前提,而且还必须考虑评估对象所处未来环境的复杂性和不确定性,这就不可避免的会导致评估结果的波动性与不确定性。然而当前评估实务尤其是周期性公司估值实务中,在这一环节的做法差别较大,有时还要依赖于评估师一定的主观与经验判断,得到的评估结果也趋于单一化和绝对化,可信度往往会遭受质疑。由此可见,对资产评估收益法不确定性预测的处理方法研究迫在眉睫。蒙特卡罗模拟,又称随机模拟法,是以概率和统计理论为基础的一种计算方法。该方法将所求解的问题与某个概率模型联系在一起,在计算机上进行随机模拟,以获得问题的近似解及其分布情况。这是一种先进的数字仿真技术,其实质是利用服从某种分布的随机数来模拟现实系统中可能出现的随机现象。通常,对未来的情况都是不能确定的,但如果知道每一个输入变量的概率分布情况,就可以通过运用一个随机数发生器来产生具有相同概率分布的数值,重复次地给每个输入变量赋值,从而每次都会对应实际上可能发生的一种情况。通过大量次数的模拟,就可以得到结果的一个概率分布情况。可见,这种方法运用的是多次估值来表示结果,而不是用一个单一的点估计值来表示。蒙特卡罗法的这种特性,使其恰好可以弥补资产评估收益法中对不确定性预测的需求。鉴于此,本文以中国某周期性上市公司(以下简称“EC公司”)在评估基准日2009年12月31日的重大资产重组评估为例,详细说明在收益法评估过程中,通过蒙特卡罗模拟对未来数值预测的不确定性的处理。

二、蒙特卡罗模周期性公司收益法估值及其不确定性因素预测分析

(一)案例公司简介

EC公司是中国主营煤炭生产与经营的典型周期性上市公司,但由于公司整体规模较小。因此,决定将截至基准日除货币资金外的全部资产及负债出售给某投资公司,同时通过发行股份购买某煤炭集团持有的优质煤炭资产,以扩大上市公司资产规模,提升上市公司盈利能力。经评估机构对拟出售资产实施清查核实、实地查勘、市场调查和询证、评定估算等评估程序,采用收益法中常用的现金流折现方法(DCF)对拟出售资产进行评估。评估基准日2009年12月31日拟出售资产账面值为6,922.71万元,评估后的价值为14,158.12万元,评估增值7,235.41万元,增值率为104.52%。根据EC公司2005年至2009年的资产负债表、损益表和内部管理报表,综合考虑未来5年以及永续期各种相关因素的影响。本次评估中,加权平均资本成本WACC在2010年至2011年采用r=14.70%,2012年及以后采用r=14.63%。并且,EC公司拟出售净资产未来年度的现金流量预测如表(1)所示。

单位:万元

文中有关EC公司的相关数据信息均引自其上市公告年报信息。

(二)蒙特卡罗模拟

从本评估报告的收益法未来现金流预测中可以看出:(1)评估的详细预测期为5年,即2010年至2014年;(2)评估假设该公司可保持长时间的运行,故评估收益期按永续确定。且2014年后的永续收益趋于稳定,假定与2014年相同,即永续增长率为0%。显然,本次评估对EC公司2014年以后未来预期收益的处理过于简单。尽管煤炭行业相对其他行业来说收益较为稳定。但随着时间的长期推移,无论是国家出台的一系列对能源资源开采的鼓励抑或是限制性的宏观政策,还是公司内部进行资产结构调整以适应其更好发展的战略方针,都必然会使其经历或多或少的收益和成本的波动。不产生任何波动的未来收益预测,一定程度上是不具有说服力的。因此,对于EC公司在2014年后永续期的净现金流量现值的预测,十分有必要考虑其波动情况加入不确定性分析,并通过蒙特卡罗模拟重新调整其预测值。其中,对于详细预测期2010年至2014年所列的评估结果,本文仍采用原评估报告中的评估数据,不再另作分析。在进行蒙特卡罗模拟前,首先应当先对评估中所涉及的未来的不确定因素一一进行波动情况分析和确定。该实例中由于是采用收益法中常用的现金流折现方法(DCF)进行评估,因此涉及的未来不确定性因素主要有:主营业务收入、主营业务成本、销售税金及附加、营业费用、管理费用、财务费用、折旧、摊销、扣税后利息、资本性支出等。

单位:元

(三)EC公司主营业务收入和主营业务成本的波动性分析

主营业务收入和主营业务成本通常是对企业净利润额影响最大的两部分。因此,本文对这两个部分进行重点分析。但由于EC公司是煤炭行业的新进入企业,缺乏足够的历史数据支持,难以分析其在该行业中的未来营业收入和营业成本的波动情况。因此需要选取与其资产规模和营利能力基本相似的可比公司来进行收入和成本的波动性分析,并将分析结论作为分析EC公司时的参照。本文根据EC公司的资产规模和营利能力,选取与其在重大资产重组后行业类型相一致,且总资产和净利率相当的上市公司的历史利润表数据进行分析,以推测出EC公司未来利润表各个项目数值的概率分布。通过对煤炭行业35家上市公司的基本面情况的比较(特别是总资产规模及净利率的数值上),这里初步选定了中国另一家煤炭业上市公司(以下简称“SC公司”)作为EC公司的可比公司。本文通过SPSS软件对SC公司主营业务收入和主营业务成本这两部分的数据进行概率分布统计,分析其金额的大致分布特征,从而可推断出目标公司DY公司相应项目金额的未来概率分布。SC公司1998年至2009年的主营业务收入及主营业务成本情况如表(2)和表(3)所示。本文首先对SC公司的主营业务收入通过运用SPSS进行K-S单样本检验。通过该检验研究样本观察值的分布和指定的理论分布是否吻合,即利用样本数据推断其是否来自指定分布的总体。在SPSS软件中一共给出了4种指定分布,分别为正态分布、均匀分布、指数分布、泊松分布。在对上述营业收入进行检验时,选定的单侧显著性水平为0.05,且原假设和备择假设分别为:

H0:SC公司1998年至2009年的主营业务收入服从正态分布;

H1:SC公司1998年至2009年的主营业务收入不服从正态分布

检验结果如表(4)所示。显然,从表(4)的结果可以看出,均值为3,236,100,000,标准差为2,388,650,000,双侧渐近显著性水平为0.345;由于这里所选定的单侧显著性水平为0.05,且0.345>0.1,进而可得结论:检验不显著,无理由拒绝原假设,即认为SC公司1998年至2009年的营业收入和正态分布没有显著差异。由此可知,SC公司1998年至2009年的主营业务收入来自正态总体N(3,236,100,000,2,388,650,000^2)。同理,对SC公司的主营业务成本进行K-S检验,结果如表(5)所示。由此可知,SC公司1998年至2009年的主营业务成本同样服从正态分布。根据上述分析,由于SC公司与EC公司存在较好的可比性,假定EC公司未来可持续状态的主营业务收入及主营业务成本的金额也符合相应的正态分布,而并非是原评估报告中的分别保持在2014年的预测值26457.18万元和14602.47万元。为更好地表示出这种正态概率的波动性,这里假设EC公司在2014年之后每一年的主营业务收入和主营业务成本分别服从期望值为26457.18万元和14602.47万元的正态分布。即,在原评估值的基础上,赋予其一定程度的随机波动概率。

(四)EC公司其它不确定性项目的波动性分析

相对于主营业务收入与主营业务成本对净利润影响程度的显著性,对于处于煤炭行业中的EC公司而言,营业费用、管理费用、销售税金及附加数额则相对较小,且基本保持稳定。虽也会有波动,但波动范围不大,总体上来说概率分布均匀。因此,在原报告评估值的基础上,赋予其一定范围内的均匀分析概率,且具体范围以原评估报告值确定,即假设EC公司的营业费用、管理费用、销售税金及附加在2014年以后服从一定范围内的均匀分布,且该范围由原报告中预测值的最大值和最小值决定。于是,2014年之后每一年的营业费用在[998.95,1104.4]范围内服从均匀分布;每一年的管理费用在[6307.53,6525.42]范围内服从均匀分布;每一年的销售税金及附加在[929.26,977.71]范围内服从均匀分布。而财务费用、折旧率、摊销,扣税后利息、资本性支出的金额相对较少,对净利润的影响并不显著,因此,这里仍采用原报告中2014年的数值,未来保持不变。最后,对于EC公司未来永续增长率和加权平均资本成本的概率分布,假设两者在2014年以后,均服从三角形概率分布特征。三角形概率分布是一种简单的分布形式,适合于数据缺乏,但能得到变量的最高、最低和最可能值的情况,也是不确定性分析中常用的一种分布形式,尤其当变量的分布形式相当集中,分析者可以估计变量范围的极值、而极值的概率又很低时,这种分布更能确切地反映变量的分布。对于EC公司的永续增长率,假设其永续增长率介于[-1%,+1%]之间,且最可能值为原评估报告中的0%。而对于EC公司的加权平均资本成本,假设其在2014年后的最大值为16%,最小值为13%,最可能值即为原评估报告中的14.36%。综上,上述关于EC公司蒙特卡罗模拟前的收益波动性下不确定性因素的预测可小结如表(6)所示。

三、蒙特卡罗模拟运算与结果比较

(一)蒙特卡罗模拟分析

分别将表(6)中的12个项目设定为assumption,并分别设置好相应的概率分布情况;将2014年之后的净现金流量现值设定为forecast;将模拟次数设定为1000次,置信区间设定为95%,确定水平为100%;之后运行模型程序,输出结果如图(1)所示。图(1)(Frequency view)是对EC公司2014年后的净现金流量现值的预测图,共显示了997个模拟结果(997isplayed),结果中有3个异常值未列入模拟,即,模拟结果代表了对99.7%的数据的统计分析。100%的确定性水平说明997个模拟结果100%都落在了蓝色的区域范围之内。但图中显示的仅是结果数据概率分布的一个大致特征,基本是服从正态分布的。为了使结果更清晰地展现出来,可进一步分析图(2)中的数据输出结果(Statistics view)。根据图(2)所列数据,可以清楚的看出,模拟结果的平均值为8615.75万元,中位数为8674.38万元,标准差为10033.74万元,评估结果的波动范围在(-21555.19,47168.50)之间。因此,可以得出结论:EC公司2014年后的净现金流量现值最可能为8615.75万元,并且在95%置信水平下的价值区间-21555.19~47168.50充分反映了资产评估中不确定性的存在。

单位:万元

(二)估值预测结果比较分析

进一步将EC公司基于模特卡罗模拟的预测结果与原评估报告数据结果进行比较发现,原评估报告中对EC公司2014年后的净现金流量现值的预测值为9271.47万元,而蒙特卡罗模拟分析在考虑相关项目的波动性概率后,模拟的结果为8615.75万元,二者相差约655.72万元。这种差异即来自于未来相关因素的不确定性和收益、成本的随机波动,如表(7)所示。

四、结论与建议

当前评估实务尤其是周期性公司估值实务中的收益法预测,在某种程度上趋于单一化和绝对化,这或多或少会降低评估结果的可信度。而蒙特卡罗模拟则是通过确定未知参数恰当的波动范围,使预测值不会过于绝对化,并以此得出相应结果的波动范围和最可能值,从而提高了评估预测的合理性和评估结果的说服力。由此可见,蒙特卡罗模拟的不确定分析在很大程度上与资产评估收益法中的不确定性预测行为相一致,并从理论上较好满足了包括周期性公司估值在内的资产评估收益法预测的需求。但不可否认的是,在蒙特卡罗模拟中各个数值概率分布情况的假设方面,仍然需要更多的理论依据和数据支持,需要进行进一步的探索,以使其更合理的应用到资产评估收益法实践中来。

参考文献

[1]斯蒂芬A.罗斯:《公司理财》,机械工业出版社2009年版。

[2]戴维R.安德森:《商务与经济统计》,机械工业出版社2005年版。

[3]贾俊平、何晓群、金勇进:《统计学》,中国人民大学出版社2009年版。

蒙特卡罗法模拟篇2

质量稳定性是水泥的重要质量指标之一。许多用户对水泥质量稳定性的关注程度，甚至高于对质量指标水平的关注。水泥质量稳定性以均匀性试验得到的变异系数表示，GB12573—90《水泥取样方法》和JC/T578—1995《评定水泥强度匀质性试验方法》规定，均匀性试验采用系统抽样方法，样本容量为10。生产实践表明，该抽样方法的抽样误差较大，但是难以通过生产试验的方法确定抽样误差的数值。文献[1,2]应用蒙特卡罗方法对抽样误差进行数值模拟取得满意结果。文献[3]分析了抽样方差估计值的精度与样本数目的关系，指出抽样方差估计值的标准偏差与样本容量的平方根之积近似为一常数。

本文应用蒙特卡罗方法模拟了现行均匀性试验方法的抽样误差、扩大样本容量后的抽样误差。根据试验结果提出了对现行均匀性试验抽样方法的修改建议。

1 一般原理与模拟方法

1.1 一般原理

蒙特卡罗模拟的一般原理和方法已有介绍[4]。

在一个总数为N的检验批按系统抽样方法在其中抽取n个子样的方法是：

1)按下式计算抽样间隔k:

其中INT是取整函数。

2)确定抽样起点i:

其中RAN为随机数发生函数。式(2)随机产生均匀分布于(1,k)的自然数。

3)以i作为起点，依次在i,i+k,i+2k,...,i+(n-1)k的位置上抽取子样共n个。

1.2 模拟方法

设每个检验批(即一个编号的出厂水泥)为1 000t，按系统抽样方法在其中抽取n袋水泥。则总体个数N=20 000，样本容量(即每个出厂水泥编号的抽样数量)为n，抽样间隔k=20 000/n。设总体强度标准偏差(即一个编号内各袋之间的标准偏差，均匀性试验结果是总体标准偏差的估计值)为σ，总体强度平均值为μ。

产生N个服从标准正态分布的随机数R，即R～N(0,1)。均值为μ、方差为σ2的正态分布随机变量X可通过下列变换得到：

即X～N(μ,σ2)。

产生均匀分布于(1,20 000/n)的随机数i，取Xi,Xi+k,Xi+2k，…，Xi+(n-1)k，计算子样的标准偏差S，即完成一次模拟。连续模拟30次，计算30次模拟标准偏差估计值S的标准偏差。

2 模拟结果与讨论

按μ=50MPa(该数值大小与模拟结果无关)，σ=0.5MPa、1.0MPa、1.5MPa、2.0MPa,n=10次、20次、30次、40次、50次进行全部组合模拟，以观察不同标准偏差和不同样本容量对抽样误差的影响。σ=1.0MPa,n=10次的模拟结果见表1。

不同总体标准偏差σ、样本容量n的模拟结果的统计值见表2。

取置信概率为0.95，根据表2结果计算的不同总体标准偏差σ、样本容量n下，以标准偏差表示的抽样误差见表3。

MPa

在0.95的置信概率下，以标准偏差表示的抽样误差与样本容量的关系见图1。

表3和图1表明，当样本容量为10，总体标准偏差分别为2.0MPa、1.5MPa、1.0MPa、0.5MPa时，标准偏差表示的抽样误差分别为0.813MPa、0.676MPa、0.437MPa、0.220MPa，相对误差分别为40.7%、45.1%、43.7%、39.7%。如此大的抽样误差是无法接受的。表明现行均匀性试验方法的抽样数量偏少。

样本容量增加，抽样误差随之减小。当总体标准偏差为2.0MPa时，在样本容量10～50区间,随着样本容量的增加，抽样误差减小的速率基本保持恒定。当总体标准偏差为1.5MPa、1.0MPa、0.5MPa时，在样本容量10～50区间，随着样本容量的增加，抽样误差开始时下降较快，在样本容量为30以后抽样误差减小的速率变慢。这意味着样本容量小于30时，增加样本容量可以比较明显减小抽样误差;样本容量大于30时，增加样本容量对减小抽样误差的效果不明显。

总体标准偏差减小，抽样误差随样本容量增加而减小的速率减小。这意味着在总体标准偏差较大时，增加样本容量对减小抽样误差具有更明显的效果。

总体标准偏差增加，抽样误差显著增加。图2更加清楚地显示了这种关系。

图2表明，随着总体标准偏差的增加，抽样误差近乎呈线性增加。样本容量增加，抽样误差随着总体标准偏差增加的速率减小。这意味着当样本容量较小时，随着总体标准偏差的增加，抽样误差以更快的速度增加;当样本容量较大时，随着总体标准偏差的增加，抽样误差以较慢的速度增加。对于实际操作的意义是，在总体标准偏差较大时，增加样本容量对减小抽样误差的效果更加明显。

有关文件规定以变异系数表示均匀性试验结果。现行标准GB12573—90规定均匀性试验的样本容量为10，此时不同28d抗压强度下，以变异系数Cv表示的抽样误差与总体标准偏差的关系见图3。

当总体标准偏差为1.5MPa———这是一般水泥厂均匀性试验结果的最大值，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为1.7%、1.4%、1.1%。相对于《水泥企业质量管理规程》规定均匀性试验的变异系数不大于3.0%的要求，这个误差显然是不能接受的。再次证明现行均匀性试验方法的抽样数量偏少。

在总体标准偏差为1.5MPa和1.0MPa时，以变异系数Cv表示的抽样误差与样本容量的关系分别见图4和图5。

图4显示，当总体标准偏差为1.5MPa，样本容量为30，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.9%、0.7%、0.6%。

图5显示，当总体标准偏差为1.0MPa，样本容量为20，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.8%、0.6%、0.5%。

3 均匀性试验的适宜样本容量

总体标准偏差为1.5MPa，样本容量为30时，以变异系数Cv表示的抽样误差接近1%，这个误差依然有些偏大，但考虑到抽样成本和检验成本，以及在样本容量大于30以后继续增大对减小抽样误差的效果不显著，不宜继续增加样本容量。总体标准偏差为1.0MPa，样本容量为20时，以变异系数Cv表示的抽样误差也接近1%。

一般通用水泥28d抗压强度为40～60MPa，平均值约50MPa。当总体标准偏差为1.0MPa、1.5MPa时，变异系数分别为2%、3%。变异系数为2%对应着目前水泥厂的一般水平，变异系数为3%对应着目前水泥厂的较差水平。

建议对均匀性试验的抽样数量做如下规定：在连续生产的情况下，前次均匀性试验的变异系数不大于2%时，每次试验抽样数量为20个;当非连续生产或者前次均匀性试验的变异系数大于2%时，每次试验抽样数量为30个。

4 结论

1)当总体标准偏差为1.5MPa，样本容量为10，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为1.7%、1.4%、1.1%。证明现行均匀性试验抽样方法存在很大误差。

2)当总体标准偏差为1.5MPa，样本容量为30，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.9%、0.7%、0.6%。当总体标准偏差为1.0MPa样本容量为20，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.8%、0.6%、0.5%。

3)根据试验结果对均匀性试验的抽样数量作出了建议。

参考文献

[1]高志,何锡文,李一峻.分层性物质的组合取样误差与份样数目之间的关系及其Monte Carlo模拟[J].分析科学学报,1999,15(5):353-357.

[2]高志,何锡文,李一峻,等.组合取样的误差理论及其Monte Carlo模拟[J].高等学校化学学报,1999,20(12):1853-1857.

[3]高志,李一峻,何锡文,等.取样方差估计值的精度与样本数目之间的关系[J].分析化学,2001,29(2):171-174.

[4]张大康.基于蒙特卡罗方法的率值稳定性定量分析(Ⅰ)——随机检验误差对生料率值稳定性的影响[J].水泥,2007,(8):1-6.

蒙特卡罗法模拟篇3

Black-Scholes模型的假设条件过于完美, Papapantoleon对Levy过程进行了概述[1], 指出一类无限可分、左极限右连续的随机过程, 包括了Brown运动、泊松过程、指数过程及其它无穷跳跃过程, 一般没有确切的密度函数解析式, 而特征函数形式可根据Laplace分解得到具体形式, 即φ (u) =E (eiuX) =e-tφX (u) , 特征指数φX (u) 可以表达为Levy-Khinchine公式:

$φ_{X} (u) = - i μ u + \frac{σ^{2}}{2} u^{2} + \int_{R} (1 - e^{i u x} + i u x 1_{| x | < 1}) k (d x) (1)$

(μ, σ, k (dx) ) 称为Levy过程特征三项, 分别表示线性漂移率、连续布朗运动扩散、纯跳跃测度, 这三项测度能全面概括随机过程的运动特征。有限的跳跃到达率为: $Ι = \int_{R} k (d x) = λ < \infty$ , 无穷到达率为:I=∞.

随机过程的一次、二次变差表达式为:

$\begin{array}{l} J_{1} = \int_{- \infty}^{+ \infty} | x | k (x) d x \\ J_{2} = \int_{- \infty}^{+ \infty} x^{2} k (x) d x (2) \end{array}$

Madan, Carr和Chang[2]建立的VG (C, G, M) 模型可以表示为正Gamma (C, M) 和负Gamma (C, G) 两个不同过程的代数和, 是一类无穷活动率过程, 存在二次变差。模型能够描绘无穷小跳跃的活动过程, Carr, Geman, Madan和Yor[3]对VG模型的进行了扩展, 在VG测度中加入了参数Y, 以此刻画不同活动水平下的随机过程。即:

$k_{C G Μ Y} (x) = {\begin{matrix} \frac{C e^{G x}}{| x |^{Y + 1}}, & x < 0 \\ \frac{C e^{- Μ x}}{x^{Y + 1}}, & x > 0 \end{matrix} (3)$

VG模型是Y=0的特殊情形。当Y<0时, 属于类似于复合泊松分布的有限活动率过程, 0<Y<1为无限活动率过程, 1<Y<2为无限活动率且不存在有限变差, Y>2为无限活动率不存在二次变差。CGMY模型能有较好的兼容效果, 能通过参数的设定得到不同活动率形式下的随机过程, 研究这个过程下的期权定价结果具有重要意义。测度的复杂性不仅给参数估计带来了一定的麻烦, 也使得计算机模拟变得非常不便。 Carr和Wu[4]提出了CGMY过程的模拟步骤, 便是基于Levy测度下的随机数生成规则, 通过差分N阶变换, 得到Kummer等式[5]的解, 属于合流超几何函数 (confluent hypergeometric function, CHF) 。Zhang[6]提出了两种模拟方案, 一是复合泊松过程模拟方法, 二是随机时变布朗运动 (Time-Changed Brownian Motion, TCBM) 。两者依赖CHF计算, 即Francesco[7]得出的解析式。Ji[8]对CGMY过程的模拟采用的Madan和Yor (2006) [9] 的方法。Zhang[6]和Ji[8]两种判断随机数淘汰的条件都是通过Kummer等式解的变形而来, 都包含CFH, 各有优势和不足, 前者极端值偏离, 后者尾部较重。

看涨期权而言, 通常以现货价值的上下浮动5%为界区分价外、价内和平价。由于CGMY测度复杂, 模拟中随机数的淘汰涉及CHF计算, 对参数的依赖程度高。Carr和Madan[10], Carr和Wu[11], Ji[8]快速傅立叶变换 (FFT) 数值方法便利而快捷, 但在虚值期权定价时失效。

基于上述现状, 本文试图解决虚值期权定价精度问题。CGMY模型特别是参数Y比较小的情形下, 探索提高模拟效率的方法, 引入复数合流超几何函数 (CCHF) , 改进计算效率, 克服模拟困难。

2 模型与参数

2.1 VG模型

基于两个相反Gamma过程的代数和建立的无穷纯跳跃VG模型, 以参数 (C, G, M) 表示。

其中Gamma (a, b) 分布的密度函数:

$f_{G} (x; a, b) = \frac{b^{a}}{Γ (a)} x^{a - 1} e^{- b x} (4)$

若G (a, b) 表示Gamma过程, 价格增量演变过程可以用等式描绘:

$Δ L_{Τ} = G_{Τ}^{1} (C, Μ) - G_{Τ}^{2} (C, G) (5)$

2.2 CGMY模型

作为无穷活动率下纯跳跃VG模型的拓展, 仅考虑纯跳跃来刻画资产价格演变过程, 那么它的特征三项中连续布朗运动扩散可以为0, 即 (w, 0, k (x) ) 。跳跃的活动到达率水平取决于Y的大小。Y在 (0, 1) 区间就能体现无穷小跳跃、无穷到达率水平。w为线性漂移率, 非随机趋势, 根据式 (3) Levy测度k (x) , 可以得到CGMY过程的跳跃频率和强度。

由于CGMY模型没有确切的密度函数解析式, 而特征函数表达式比较简洁, 通过Levy-Khinchine公式, 模型的特征函数为:

$φ_{C G Μ Y} (u) = e^{C Γ (- Y) [(Μ - i u)^{Y} + (G + i u)^{Y} - Μ^{Y} - G^{Y}]} (6)$

其中伽马函数:

$Γ (a) = \int_{0}^{\infty} x^{a - 1} e^{- x} d x (7)$

关于参数Θ下随机过程Xt的Laplace变换公式是:

$φ_{X_{t}} (u; Θ) = E (e^{i u X_{t}}) = e^{- t φ_{X} (u)} = l_{t}^{Θ} (φ_{X} (u)) (9)$

指数Laplace形式是自然底数e的指数形式, 往往称为特征指数、累积特征指数, 转化为资产价格对数收益率特别方便, 利率期限结果模型中用于计算零息债券的价格公式 (假设有价证券面值为一个单位) 即为E (e-tRt) =e-tR (λ) =lΘt (λ) , 是关于利率仿射指数变换的函数表达式。

另外, 令u=-i, 有φXt (-i;Θ) =E (eXt) =e-tφX (-i) , φX (-i) 是随机过程的指数线性漂移率、单位时间的预期对数收益率, 用于建立风险中性资产价格模型和鞅过程。

由CGMY模型的特征指数得到CGMY过程分布的累积函数:

$\begin{array}{l} C_{n} = \frac{1}{i^{n}} \frac{d^{n} \log (Φ (u))}{d u^{n}} |_{u = 0} \\ = C Γ (n - Y) [Μ^{Y - n} + (- 1)^{n} (G^{Y - n})] (10) \end{array}$

依据特征指数可以计算风险中性漂移率修正参数:

$\begin{array}{l} w = φ (- i) \\ = C Γ (- Y) [(Μ - 1)^{Y} + (G + 1)^{Y} - Μ^{Y} - G^{Y}] (11) \end{array}$

在非单位时间段t内, CGMY模型的参数变换关系:

$C_{t} = t C; G_{t} = G; Μ_{t} = Μ; Y_{t} = Y (12)$

为便于计算和模拟, 设立一个正的极小值ξ→0, 作为无穷小跳跃的门槛, 高于小跳门槛的以CGMY的Levy测度表示, 低于该门槛则以无穷小波动的连续布朗运动代替。-ξ<x<ξ称为小跳跃, |x|>ξ为CGMY的跳跃, 通过设定和排除趋近于正态分布的无穷小跳跃, 能更加方便地模拟随机数的分步以及进行数值计算。

2.3 特征函数的矩估计

有了矩条件表达形式, 参数的矩估计采用Y.Miyahara的方法[18]。

设z分布服从CGMY过程的增量分布, k阶原点矩表达式:

$\begin{array}{l} \hat{m}_{k} = E (z^{k}) \\ \hat{A} (u) = E (e^{i u z}) \\ \hat{A}^{(k)} (u) = d^{k} E (e^{i u z}) / d u^{k} \\ \hat{m}_{k} = \hat{A}^{(k)} (u) / i^{k} |_{u = 0} (13) \end{array}$

根据式 (7) 的矩条件: $\hat{C}_{n} = C Γ (n - Y) [Μ^{Y - n} + (- 1)^{n} (G^{Y - n})], Ν$ 取一个比较大的整数解四阶方程组, 参数结果表达式:

$\begin{array}{l} \hat{Y} = \frac{\ln (Re (\ln (\hat{A} (Ν))))}{\ln Ν} \\ \hat{A} (u) = E (e^{i u x}) \\ \hat{C} = \frac{Re (\ln (\hat{A} (Ν)))}{2 Ν^{\hat{Τ}} Γ (- \hat{Y}) \cos (\frac{σ^{2} \cdot \hat{Y}}{2})} (14) \end{array}$

代入下列方程组求解:

$\begin{array}{l} \hat{C} Γ (- \hat{Y}) (\hat{Y} - 1) \hat{Y} [G^{Y - 2} + Μ^{Y - 2}] = C_{2} - C_{1}^{2} \\ \hat{C} Γ (- \hat{Y}) \hat{Y} (\hat{Y} - 1) (\hat{Y} - 2) [G^{Y - 2} + Μ^{Y - 2}] \\ = C_{3} - 3 C_{1} C_{2} + 2 C_{1}^{3} (15) \end{array}$

方程 (15) 由于指数复杂, 普通计算很难有较好的准确性。可在Matlab中采用函数最小值法求参数。X表示参数集, f1 (X) 、f2 (X) 分别表示模型矩条件, 参数估计原理为

$X = a r g m i n (f_{1} (X) - Μ_{1})^{2} + (f_{2} (X) - Μ_{2})^{2} (16)$

通过上述低阶矩估计方法获取参数后, 进行期权的数值定价或蒙特卡罗模拟比较在CGMY模型下的实值期权与虚值期权之间的差异。

3 分布模拟原理

3.1 传统模拟计算技术

根据Madan和Yor[9]的模拟方法, 在随机泊松跳跃时间内先产生一个逆CGMY随机变量, 可称之为levy测度下的随机跳跃y, 再通过计算CHF (y) 作为判断条件进行拒绝-接受检验, 可以接受的随机数y作为随机时间增量生成漂移率为A的布朗运动, 步骤:

①计算以下参数:

A= (G-M) /2

B= (G+M) /2

$k = C \sqrt{π} / [(Γ ((Y + 1) / 2) \cdot 2^{Y / 2}]$

λ=2k/ (YεY/2)

d=kε1-Y/2/ (1-Y/2)

ε→0

②生成以λ为参数的指数分布a≤ΔT, 以此作为判断跳跃的时间, 反之退出模拟

③生成两个独立的均匀分布 :u1, u2～U (0, 1)

④生成列维逆正态分布, 即随机跳跃的幅度: $y = \frac{ε}{u_{1}^{\frac{2}{Y}}}$

⑤根据随机数生成接受-拒绝的判断条件 (以下称U法则) :

$u = e^{\frac{- A^{2} + B^{2}}{2}} \frac{Γ (\frac{Y + 1}{2})}{\sqrt{π}} U (\frac{Y}{2}, \frac{1}{2}, \frac{B^{2}}{2} y) > u_{2}$

⑥若满足上述判断, 以随机数作为跳跃时间幅度的累积:sum=sum+y

⑦把跳跃作为时变布朗运动的随机时间, 生成CGMY过程的增量: $Ζ = A \cdot s u m + ξ \sqrt{s u m}$ , 这样价格过程就可以表示成:

$S_{t} = S_{0} e^{(r - q) t + Ζ (t)} (17)$

U (a, b, x) 是一个合流超几何函数。

另外, Zhang[6]运用了类似的淘汰规则, 是基于Kummer方程的第二个解形式 (下简称H法则) :

sum=sum+yi1h (yi) >ui

$\begin{array}{l} h (y) = 2^{Y} \exp {\frac{(B^{2} - A^{2}) y}{2}} [F (\frac{Y}{2}, \frac{1}{2}, \frac{B^{2}}{2} y) \\ - \frac{\sqrt{B^{2} y} \cdot π \cdot Γ (\frac{Y + 1}{2})}{\sqrt{2} Γ (\frac{Y}{2}) Γ (\frac{3}{2})} F (\frac{Y + 1}{2}, \frac{3}{2}, \frac{B^{2}}{2} y)] (18) \\ F (a, b, x) = U (a, b, x) \end{array}$

3.2 复数合流超几何函数计算

特别地广义超几何函数表达式为:

$F (n, d; x) = \sum_{k = 0}^{\infty} \frac{\prod_{i = 1}^{j} \frac{Γ (n_{i} + k)}{Γ (n_{i})}}{\prod_{i = 1}^{m} \frac{Γ (d_{i} + k)}{Γ (d_{i})}} \cdot \frac{x^{k}}{k!} (19)$

Y越小产生的随机数y可能越大, 函数计算复杂度增强, 出现的极端值非常不稳定。考虑到计算过程受到参数影响, 对参数进行复变化可以收敛实数分布域, 本文采用复数变换合流超几何函数 (Kummer Complex函数) 代替传统CHF, 是Stepan[23]的一种广义傅立叶变换算法, 即:

$C (z; Θ) \equiv \int_{- \infty}^{\infty} e^{i z v} \cdot F (x; Θ) d v (20)$

z可以为一切复数, 通过变换, 可以简化相关计算。特征函数就是这类转换的一种特例。采用了傅立叶变换方式, 该函数可以扩展到复数领域参数计算, 反之也可以将虚数转换到实数空间, 是傅立叶逆转换, 表达式为:

$F (x; Θ) \equiv \frac{1}{2 π} \int_{i z - \infty}^{i z + \infty} e^{- i z k} \cdot C (z; Θ) d z (21)$

虽然需要通过正反两次变换, 但无论参数多大的实数督可以变成二元复数, 大大减少了CHF的计算时间和复杂度。

Stepan[23]指出通过函数计算对比F函数和U函数的值误差不超过10-6, 计算速度却快了很多。同样也解决了y非常大的时候广义超几何函数无结果的问题, 计算效率得到较大程度的提高。

为进行比较, 针对不同Y模拟生成1000个随机数, 并进行CGMY分布的模拟, 比较两种计算方法的平均时间。当Y越小时, 需要的时间越长, 随机数出现超极端值的概檬就越高, CHF函数的计算就越复杂, 可能导致死循环, 通过无效值比例可以看出模拟1000个随机数时, 就可能导致一次死循环。两种函数计算时间比较参见表1[10]引入转换因子, 将期权价格进行广义傅立叶变换, 期权价格的傅立叶形式为 (α为转换因子也称修正因子) :

$F (v) = φ (v) = \int_{- \infty}^{\infty} e^{i v k} e^{α k} C_{Τ} (k) d k (24)$

式中, k为执行价 (对数形式) , CT (k) 表示到期日的期权价格, 关于执行价序列的期权价格向量转换成关于v向量的傅立叶序列, 该序列有许多优点, 代入欧式期权定价公式并进行积分秩序变换, 密度函数积分成特征函数得到以下结果 (小写的价格和执行价均表示对数形式) :

F (v) =∫∞-∞eivkeαkCT (k) dk

=∫∞-∞eivkeαk∫∞ke-rt (es-ek) f (s) dsdk

通过推导得到了特征函数形式下的期权傅立叶形式的数值解:

$φ (v) = \frac{e^{- r t} Φ (v - (α + 1) i)}{α^{2} + α - v^{2} + i (2 α + 1) v} (25)$

只需要进行傅立叶逆变换就可以得到期权的实数解:

$C_{Τ} (k) = \frac{e^{- α k}}{2 π} \int_{- \infty}^{\infty} e^{- i v k} φ (v) d v (26)$

解析式代表了在特征函数下的价格封闭解形式, 连续函数给电脑计算却带来了不便, 而且现实中的期权执行价序列也是离散的多个有限点, 利用Trapeezoidal计算规则进行离散化可以得到:

$\int_{- \infty}^{\infty} e^{- i v k} φ (v) d v \approx \sum_{j = 0}^{Ν - 1} e^{- i v_{j} k} φ^{'} (v_{j}) η (27)$

特别注意到积分结果依赖于积分间隔η和积分上限a=Nη, 通过设立k、v和N的关系, 可以得到快速傅立叶形式:

$Y = \sum_{j = 0}^{Ν - 1} e^{- i \frac{2 π}{Ν} j} X (28)$

通过特征函数求期权价格的傅立叶形式, 并计算与执行价无关的序列X, 即 $h_{j} = φ^{'} (v_{j}) η e^{i \frac{Ν λ}{2} η j}$ , 再进行快速傅立叶逆变换, 是FRFT数值方法 (Fractional FFT) : $Y = \sum_{j = 0}^{Ν - 1} e^{- i \cdot 2 π j ε} X$ 的特殊情形, 令 $ε = \frac{1}{Ν}$ , 得到新的序列Y就是关于执行价的期权价值。通过CGMY的参数、不同的执行价、恒生指数、无风险利率、时间参数设定后, 历时不到2s, 数值计算结果如图1所示。

图1显示, 执行价在20000点以内, 实值期权 (in the money) 定价结果比较精确, 误差较小。以20000点为界的价外期权 (out of the money, 虚值期权) 深度虚值甚至出现负数, 不符合期权的实际价值, 需要进一步研究虚值范围内蒙特卡罗模拟定价效果, 假设K>S即为虚值期权。

4.3 蒙特卡罗模拟

通过第3节的模拟比较, 为了提高模拟效率, 一例以Kummer Complex函数计算合流超几何函数值。分别采用H和U淘汰法则来模拟生成CGMY分布如图2所示。

数据拟合统计检验采用MATLAB中t-test方法的统计检验见表3。

H0=1接受原假设, 通过t-test二维拟合检验。

现若进行期权定价只需要对C进行参数变换:

$C_{Τ} = C \cdot Τ = 254 C (29)$

风险中性修正的漂移率为:

CGMY过程的风险中性下模拟路径:

$S_{t} = S_{0} e^{(r - w) t + Ζ (t)} (31)$

欧式看涨期权蒙特卡罗模拟定价关系式:

$c (Κ) = E (\max (S_{t} - Κ, 0)) = \frac{1}{n} \sum_{i = 1}^{n} m a x ((S_{t} - Κ, 0)) (32)$

通过MATLAB编程实现模拟路径及期权价格计算, 设定路径数np=10000, 每条路径步数ns=100, 共需产生CGMY随机数1000000个, 历时小于30s. 而采用普通广义超几何函数进行计算, 因一周内产生死循环重新运行MATLAB数次, 故不做时间统计。通过模拟计算得到的期权价格参见图3。

如图3, 通过模拟比较, 无论哪种法则, 相对数值方法误差更小, 深度虚值领域也存在正的价值, 这比较符合期权理论和现实。相对数值结果, 模拟定价的效果要精确很多, 证明该方法的有效性。相比起数值方法, 这种模拟方法在虚值期权定价上更具先进性。同时, 两种法则中H法则下的模拟结果较为发散, 在恒生指数期权定价方面, Kummer函数显然比Francesco函数好。

5 结论

通过CGMY模型下数值和模拟两种常用定价技术的比较, 进一步提出改进该模型下蒙特卡罗模拟方法, 引入复数合流超几何函数计算, 以提高效率, 结果显示:①运用于恒生指数期权定价, 数值计算期权价格十分快捷和便利, CGMY模型下实值期权内定价精确有效, 但不适合虚值期权定价。②CGMY模型下蒙特卡罗模拟虚值期权定价比较有效, 深度虚值期权也不例外。相比起数值计算有很大的改进。误差较小, 速度较慢。③CGMY模拟中采用广义超几何函数或实数域内的合流超几何函数在Y<0.5的计算都是低效率的, 通过合流超几何函数的傅立叶变换使得计算变得快捷而且有效, 在Y越小时越明显, 采用复数域的合流超几何函数代替合流超几何函数计算可以提高计算速度。④Kummer微分方程的两个非线性相关解都属于随机数淘汰条件, 复数代替广义合流超几何函数, 在期权定价模拟采用Kummer方程第一类解, 相对第二类解稳定和收敛一些。

根据以上结论, 提出如下建议:①新函数计算速度不依赖于模型参数变化, 不因参数无穷小而影响计算速度, 极大提高随机数生成的效率, 从而避免随机数生成死循环。CGMY模型下, Y接近于零或者为负数的时候最为适用。②如果只从单向期权进行定价, 实值期权采用数值计算而虚值期权采用模拟的方法比较快捷和稳健。

蒙特卡罗法模拟篇4

房地产开发是一种投资巨大、开发周期长的项目, 这就导致房地产成为一个有较高风险的行业, 所以对房地产开发项目进行风险分析有很强的现实意义。房地产项目风险是指在房地产开发活动过程中存在影响开发收益的多种因素, 这些风险使得企业实际开发利润可能与预计的利润发生偏差, 从而使企业有产生经济损失的可能。按风险的成因可将房地产项目的风险分为:政治风险、经济风险、自然风险、技术风险、管理风险、社会风险、国际因素风险等。房地产项目的风险有如下特点:

(1) 风险的综合性。由于和房地产行业相关的产业众多, 包括建筑业、建材业、机械制造业、钢铁业、交通运输业等, 这些行业原材料价格的变动、供需的变动都会对房地产行业产生影响。

(2) 资金规模大。房地产项目开发一般是投资规模较大的大型工程。从购买土地到房产建成销售, 涉及很多环节, 需要投入巨额的资金。开发一个房地产项目有时需要资金多达上亿元甚至几十亿到上百亿元。巨大的资金需求使房地产开发商面临利率风险。同时由于房地产开发项目周期较长, 开发商也面临通货膨胀风险。

(3) 房地产售价变动的风险。房地产开发商从选地、购地、建筑设计、房地产开发到销售, 往往需要几年时间, 一般来说, 房地产开发周期越长, 房地产的售价越难以预测, 房地产开发商面临的房地产售价变动风险越大, 同时由于房地产开发项目周期长, 房地产项目后续建设和建筑安装费用不确定性的风险也增大。

2 蒙特卡罗模拟的Excel实现及敏感性分析

蒙特卡罗模拟也称统计模拟方法, 它是一种以概率统计理论为指导的一类非常重要的数值计算方法。蒙特卡罗方法在数学、金融工程学、宏观经济学、计算物理学等领域有着广泛应用。应用蒙特卡罗模拟解决问题一般包括以下几个步骤:

(1) 分析问题, 建立实际问题的模型。一般是先分析实际问题, 然后构造随机变量的统计量, 比如Y=f (X1, X2, …, Xn) , 使这个统计量Y的某个统计参数等于要求解的问题的值。

(2) 当不知道随机变量的分布类型和参数时, 用已知的样本数据拟合随机变量分布和参数。

(3) 重复产生大量随机向量X= (X1, X2, …, Xn) , 将其代入Y=f (X1, X2, …, Xn) 中即可得到Y的值。然后计算Y的相关统计参数就可以得到求解的问题的值。

(4) 利用模拟中产生的数据Y和X进行其他分析, 比如敏感性分析等。

Excel是office办公套装软件的一个重要组成部分, 它可以进行各种数据的处理、统计分析和辅助决策操作, 广泛地应用于管理、统计财经、金融等众多领域。Excel软件特别适合数据分析与建模, 同时Excel还具有强大的绘图功能, 方便将数据以图表的方式展现出来。同时Excel还支持VBA编程功能, 如果需要一些Excel没有的功能和函数, 可以使用VBA编程来增强Excel的功能。所以使用Excel进行蒙特卡罗模拟是很方便和快捷的。

现在通过一个实例来说明蒙特卡罗模拟在房地产项目风险分析中的应用。

某房地产公司准备开发一个新的商品房项目。该房地产公司先支付土地出让金及前期费用后进行商品房开发, 商品房分三期, 预计每年建设一期, 销售一期。房地产公司经过分析认为前期费用与土地出让金服从3 500万~4 500万元之间的离散均匀分布, 期初的建筑与安装费用服从800万~1 300万元之间的离散均匀分布, 第一年的销售费用服从80万~120万元的离散均匀分布, 销售收入年增长率服从均值为10%, 标准差为2%的正态分布, 建筑安装费用年增长率服从均值为15%, 标准差为3%的正态分布, 销售费用年增长率均值为10%, 标准差为2%的正态分布。在Excel中后建模如表1。

然后利用Excel的模拟运算表模拟1 000次后得到NPV的直方图。

该房地产项目NPV的相关统计参数和表2。

由于该房地产项目NPV大于等于0的概率为84%, 平均值为610.52万元, 可见该房地产是一个挺好的项目, 但同时也存在一定的风险。开发商准备控制风险, 所以要分析该项目各个风险因素影响的程度, 从而得到影响NPV的主要风险因素, 这就是房地产项目敏感性分析。

本文提出基于Spearman秩相关系数的敏感性系数的构建方法。假设风险因素有k个, 分别用x1, x2, …, xk表示。指标Y表示投资项目优劣判断指标, 比如NPV, IRR等, Y=f (x1, x2, …, xk) 。计算Spearman秩相关系数时, 只要x和Y具有单调的函数关系, 那么x和Y就是完全Spearman相关的, 即其绝对值为1。

Spearman秩相关系数通常被认为是排列后的变量之间的Pearson线性相关系数,

式中, ρs表示xi和Y的Spearman秩相关系数;xi′, Yi′, 表示原xi, Yi在排列后数据所在的位置;xi′, Yi′称为变量xi, Yi的秩次。

Spearman秩相关系数的符号表示x和Y之间联系的方向, 大小表示联系的紧密程度。当x和Y有严格单调增加的关系时, 它们之间的Spearman秩相关系数为1;反之, 在x和Y有严格单调减少的关系时, Spearman秩相关系数为-1。所以用Spearman秩相关系数作为多因素敏感性系数指标是可行的。当分别计算出敏感因素x1, x2, …, xk对输出Y的Spearman秩相关系数ρ1, ρ2, …, ρk后, 敏感因素x1, x2, …, xk对输出Y的敏感性系数就是其对输出Y的Spearman秩相关系数ρ1, ρ2, …, ρk。

求出用Spearman秩相关系数表示的敏感性系数之后, 也就得到了输入xi对输出Y的影响程度。在实际管理中, 管理者还希望知道输入xi对输出Y影响程度的百分比及方向。下面介绍如何得到敏感性系数xi对输出Y影响程度的百分比及方向。

假设输入x1, x2, …, xk对输出Y的敏感性系数为ρ1, ρ2, …, ρk。设输入xi的敏感性系数对平方和贡献的百分比用Mi表示, 令

Mi的绝对值大小表示了输入xi对输出Y影响程度的百分比, 其符号表示输入xi对输出Y的影响方向。

将进行蒙特卡罗模拟时生成的数据带入式 (1) 和 (2) 可以得到表3和图2。

从图2中可以清楚地看出对NPV最敏感的是期初投入的建筑安装费用, 为-64.08%, 其次为前期费用与土地出让金, 占-23.69%。所以为了提高该房地产项目的NPV, 管理者应首先加强建筑安装的原材料采购、施工等环节的管理, 降低建筑安装费用。其次应分析前期费用与土地出让金的使用情况, 对资金使用流向加强监管, 尽量减低前期费用与土地出让金, 从而使该房地产项目的NPV达到最大。

3 小结

房地产行业是目前我国经济的支柱产业之一。由于房地产项目周期长、占用资金额巨大, 所以对房地产项目进行风险分析很重要。本文用Excel软件对房地产项目进行了风险分析。通过对房地产项目的风险分析得到NPV的直方图。通过敏感性分析, 找到影响房地产项目评级指标NPV的关键因素, 从而为决策者提供参考, 使项目获得最大收益。

摘要：由于房地产项目周期长、占用资金额巨大, 所以对房地产项目进行风险分析很重要。本文用Excel软件对房地产项目进行了风险分析。通过对房地产项目的风险分析得了NPV的直方图。通过敏感性分析, 找到影响房地产项目评价指标NPV的关键因素, 从而为决策者提供参考, 进而使项目获得最大收益。

关键词：蒙特卡罗模拟,房地产项目,风险分析

参考文献

[1]王学强, 庄宇.基于蒙特卡罗模拟模型的投资项目风险分析[J].工业工程, 2007, 10 (5) :93-96.

[2]王跃虹.多因素敏感性分析的积分法[J].云南工业大学学报, 1999, 15 (2) :75-78.

[3]蔡毅.敏感性分析综述[J].北京师范大学学报:自然科学版, 2008, 44 (1) :9-14.

[4]Burmaster D E, P D Anderson.Principles of Good Practice for the Use ofMonte Carlo Techniques in Human Health and Ecological Risk Assessments[J].Risk Analysis, 1994, 14 (4) :477-481.

[5]Saltelli A, Marivoet J.Non Parametric Statistics in Sensitivity Analysis for Model Output:A Comparison of Selected Techniques[J].Reliability Engineering and System Safety, 1990, 28 (2) :229-253.

[6]Hammersley J M, D C Handscomb.Monte Carlo Methods[M].New York, NY:Chapman and Hall, 1964.

[7]Kalos Malvin H, Paula A Whitlock.Monte Carlo Methods[M].New York, NY:John Wiley&Sons, 1986.

蒙特卡罗法模拟篇5

关键词：大气氡,蒙特卡罗方法,伽马能谱测量

0 引言

天然放射性伽马能谱分析法是分析天然放射性核素的比较主流的方法之一, 它能比较直观的分析出岩矿石或是环境中天然放射性核素U、Th、K的各项参数, 并且通过观察谱线可以很好的解决问题。

我们生活的大气中也充满了各种各样的放射性元素。譬如大气中的氡及其子体, 222Rn是238U系中核素226Ra的衰变子体, 220Rn是232Th系中核素224Ra的衰变子体, 还有214Bi和214Pb, 虽然它们非常微量, 却都是天然放射性核算U、Th系的衰变子体, 而且它们的浓度会随着周围环境和气象的不同而有明显的变化, 这些对伽马能谱测量分析会造成很大的影响。因此, 了解大气氡的浓度对伽马能谱测量影响就显得非常重要了。

蒙特卡罗方法是一种使用随机数 (或更常见的伪随机数) 来解决很多计算问题的方法, 它应用领域广泛, 可以模拟各种粒子运输的途径。蒙特卡罗方法的优势在于, 有些伽马能谱数据分析工程量巨大, 需要较长的时间, 人为计算也会造成很大的误差, 理论上很难推算出来。那么通过蒙特卡罗模拟程序可以很好的进行模拟分析, 可以较为直观的进行模拟, 减轻工作量。

1 蒙特卡罗方法

蒙特·卡罗方法 (Monte Carlo method) , 也称统计模拟方法, 其基本思想是当所求解问题是某种随机事件出现的概率时, 或者是某个随机变量的期望值时, 通过某种“实验”的方法, 以这种事件出现的频率估计这一随机事件的概率, 或者得到这个随机变量的某些数字特征, 并将其作为问题的解。

蒙特卡罗方法的优点在于其模拟数据曲线与真实曲线还是比较相同的;由于是用参数定义各项资料, 因此模拟过程受到的几何条件限制比较小, 很多复杂的地形或是空间可以很好的进行描述;收敛速度与问题的维数无关;误差容易确定, 在模拟结束的窗口里面可以直观的看出误差范围和其他可能存在的错误;程序结构简单, 容易上手, 易于实现。同时, 它也有弊端, 收敛速度慢, 运行时间比较长, 若是减少运行时间, 则会增大误差;误差具有概率性;模拟结果与真实值是有出入的, 若是需要精确值还是需要实际测量。所以使用蒙特卡罗方法的时候需要考虑其优缺点, 选择适当的运行时间, 发挥其特点, 为实际测量结果提供参考。

本论文中, 将运用在伽马能谱数据模拟, 通过蒙卡模拟的数据来解决问题。

2 放射性伽马能谱测量基本知识

伽马能谱测量可以在野外分别测定矿石 (岩石) 中铀 (镭) 、钍、钾三种放射性元素的含量。因而可以为放射性矿床的普查和勘探工作直接提供更多的地质信息。

2.1 铀、钍矿石的伽玛射线仪器谱

天然放射性核素的伽马射线大多是复杂谱, 其伽马射线仪器谱就更为复杂:常见的放射性铀、钍矿石就是具有多组能量的复杂伽马源, 因此用伽马能谱仪测得的仪器谱也是非常复杂的。图一、图二分别是用高分辨率的锗 (锂) [Ge (Li) ]半导体探测器在800道脉冲幅度分析器上测得的平衡铀矿石的伽马射线仪器谱与小块钍矿石的伽马射线仪器谱, 图中数字是能量峰值。

目前野外伽马能谱测量中, 一般常用Na I (TI) 晶体探测器, 其能量分辨率ε较之Ge (Li) 半导体探测器要差些。所谓能量分辨率ε就是用闪烁探测器测定35137Cs光电峰 (0.661Me V) 时, 其峰值半宽度△V与最大能量V0之比:

由式 (1) 可知, ε越大, 曲线越平缓, 则分辨率越低;ε越小, 曲线越尖, 则分辨率越高。因此我们总是希望ε小些, 即所测得能量分布曲线窄些。这样就有可能将两种能量相近的谱线分辨开。野外伽马能谱仪的分辨率一般在10%~15%左右。

根据四道伽马能谱仪在U、Th、K矿石模型表面中心测得铀、钍、钾微分谱线, 可得如下结论: (1) 康-吴散射的散射伽马射线为连续谱。该连续谱虽然掩盖了很多组分的全能峰, 但不论是铀道谱线还是钍道谱线, 铀系、钍系、钾在高能区的几个主要伽马特征峰仍清晰可辨。这几条谱线是:铀系, 1.12Me V, 1.38Me V, 1.76Me V, 2.20Me V;钍系, 2.62Me V以及单逃逸峰2.1Me V、双逃逸峰1.60Me V;钾, 1.46Me V。这就为铀钍钾道的道中心的选择提供了依据。 (2) 铀、钍微分谱在能量为几千电子伏是8 Me V的相当宽的范同内连续分布, 而钾的微分谱能量分布不超过1.46Me V。因此, 在1.0~1.60Me V能谱段内, 钾对铀、钍的干扰不容忽视。 (3) 对比铀、钍、钾的微分谱线, 往往不同能谱段内它们相互间的影响是有差别的, 或者说, 在选定谱段内各谱线与横轴所夹的面积相差很大。这就为野外γ能谱测量谱段的选择提供了依据。

2.2 伽马能谱测量的影响因素

由于伽马能谱测量并不是贴地测量, 受到的影响因素较多, 特别是在航空伽马能谱测量中, 采集的数据需要作相应的校正, 主要如下: (1) 周围环境和宇宙射线产生的本底; (2) 康普顿散射的影响; (3) 大气氡的影响; (4) 探测高度变化的影响。

本文的蒙特卡罗模拟也就是主要为了说明大气氡的影响。

3 大气氡对伽马能谱测量影响的蒙特卡罗模拟

根据蒙特卡罗方法, 分别模拟Na I (Tl) 伽马能谱仪对地面天然放射性核素源、大气氡源响应谱线的技术路线, 如下所示:

探测器栅元卡片填写→建立模型系统→源定义→输入卡片→模拟输出结果

3.1 蒙特卡罗模拟

(1) 模型系统的建立

假设地面平坦并且没有地形起伏, 考虑到实测以及蒙卡模拟的各种限制, 因此, 将体源模型设置为圆柱体结构, 圆柱体下部为土壤, 上方为大气, 建模结构示意图如图三所示。天然放射性伽马粒子均匀分布于土壤中。土壤的厚度设为100cm, 大气的高度设为100m, 探测半径R设为100cm, 将探测器置于模型中间, 并以地面作为坐标轴Z=0的面。

所以根据上述条件建立的模型系统, 填写信息卡:

(2) 源定义

模拟土壤作为天然放射性核素源的时候, 分别对铀系、钍系及40K进行模拟, 再将大气氡作为源进行模拟, 最后将铀钍钾氡整合处理。

(3) 输入卡片

单U系、单Th系、单K以及单Rn系模拟程序工作过程以单Rn系为例。

(4) 模拟输出结果

本次模拟共耗时2880分钟, 获得了4份模拟数据, 分别是单U系、单Th系、单40K以及单Rn系的模拟数据。能量变化范围是0-3Me V, 总共有1216个数据, 只含U、Th、K时总计数率为0.000714, 在此基础上加入Rn时总计数率为0.000735, 扩大了2.94%。模拟结果如图四至图九所示。

从模拟数据可以看出, K和Rn的模拟数据效果比较好, 而U和Th的模拟数据中有许多0值, 在成图的时候人为的进行了整理, 因而图中谱线可能存在误差, 并非是精确值。尽管存在误差, 但是U、Th、K的特征能量峰还是非常明显的, 譬如U的1.76Me V时的能量峰 (见图四) , Th的2.62Me V时的能量峰 (见图六) , 以及K的1.46Me V时的能量峰 (见图五) , 说明模拟曲线还是与实际曲线比较吻合的。最后的不含氡与含氡模拟数据对比图 (见图九) 可以看出, 红线为含氡时模拟图, 也就是加入了Rn的数据, 而蓝线为不含氡模拟图, 两者的形状几乎相同, 能量在1.5Me V~2.0Me V的时候有明显的区别, 而能量在1.5Me V之前差别不是非常的明显, 而能量在2.5Me V之后就几乎没有区别了。

3.2分析结果

对上述模拟数据进行归一化处理[11], 分别处理铀系、钍系、40K、氡及其子体的放射性核素系列的贡献。归一化的公式如下:

式 (2) 中, SMC (E) 是MCNP程序模拟的输出结果, Snorm (E) 为归一后的所求数据, ρs为源的密度, Vs为源的总体积, ∏是与源的产生率有关的扩展因子, 扩展因子通常通过与实验数据直接对比获得, ACi是天然放射性核素i的比活度 (单位:Bq/Kg) , pi为放射性核素i的总伽马射线发射概率, ε为探测器的探测效率, 求ε公式如下:

将所得数据进行归一化处理, 处理中所用到的各种放射性核素的比活度如表一所示, 数据来源为中国地质大学 (北京) 校园土壤[12]。

通过计算得出各核素的数据, 归一化后, 只含U、Th、K时的总计数率为:0.659, 在此基础上加入Rn时的总计数率为:0.659562, 扩大了0.085%。将U、Th、K归一化数据曲线进行整合, 得到U、Th、K综合影响数据图 (见图十) 。

在此基础上将Rn曲线图进行整合, 得到U、Th、K、Rn综合影响曲线图, 并与之前不含Rn时的综合影响曲线图对比 (见图十一) 。

从图中可以看出, 归一化后的不含氡与含氡曲线的差别不大, 同样是能量在1.5Me V~2.0Me V的时候能看出区别, 当能量大于2.0Me V以后, 几乎就没区别了。

为了反映氡浓度的影响, 因此将大气氡的浓度扩大10倍、20倍, 扩大10倍时总计数率为0.665, 较之原来的总计数率扩大了0.76%。而扩大20倍时总计数率为0.67, 较之原来的总计数率扩大了1.51%。对比数据如图十二所示。

图中蓝色红色和绿色曲线分别对应原本浓度的氡, 以及扩大10倍和20倍氡浓度时候的含氡曲线图, 紫色为不含氡的曲线。当能量小于1.5Me V时, 氡浓度的变化对曲线几乎没影响, 而当能量在1.5Me V和2.5Me V之间的时候, 浓度的变化对曲线的影响是非常巨大的, 当能量大于2.5Me V时则又几乎没什么影响了 (图中紫色曲线有许多断点, 这些点是0值, 而对应含氡曲线的点是非0值, 这个值即为大气氡的影响值) 。由此可知, 随着大气氡的浓度增加, 对伽马能谱曲线的影响也越大, 尤其当能量在1.5Me V~2.5Me V之间, 影响最为明显。

4结束语

本文利用蒙特卡罗模拟方法中的MCNP程序, 针对大气氡对现场伽马能谱数据的影响进行了模拟, 利用模型理论来尝试说明大气氡的影响。从模拟所得数据, 以及数据成图的谱线中可以看出, 含氡的数据 (计数率) 及曲线要比不含氡大, 原因就在于氡的存在, 所以一般在进行现场测量的时候都要对大气氡进行校正。

蒙特卡罗法模拟篇6

一、投资项目敏感性分析的方法

敏感性分析能够得出风险因素对投资项目目标影响的程度, 因而是投资项目风险分析中一个至关重要的环节。假设风险因素有k个, 分别用x1, x2, …, xk表示。指标Y表示投资项目优劣判断财务指标, 比如NPV、IRR等。Y是这些风险因素的函数, 即:Y=f (x1, x2, …, xk) , 其中xi表示模型的第i个风险因素。敏感性分析就是研究和预测这些输入对输出值的影响程度。影响程度的大小可以称为该输入的敏感性系数。显然, 敏感性系数的绝对值越大就说明该风险因素对投资项目的影响程度越大。敏感性分析的关键就是找到各个风险因素的敏感性系数, 进而得到各个风险因素的重要程度排序。

根据敏感性分析的作用范围, 我们可以将其分为单因素敏感性分析和多因素敏感性分析。单因素敏感性分析是在假设其他输入不变的情况下分析其中单个风险因素变动对模型的输出的影响程度, 而多因素敏感性分析是分析多个风险因素共同作用对模型输出的影响程度。单因素敏感性分析的方法简单, 但其不足之处在于忽略了风险因素之间的相关性。实际上, 一个因素的变动往往也伴随着其他因素的变动, 多因素敏感性分析考虑了这种相关性, 因而能反映几个风险因素同时变动对项目产生的综合影响, 弥补了单因素分析的局限性, 可以更全面地解释输出Y变动的原因。本文讨论的是投资项目的多因素敏感性分析。

1. 基于Pearson相关系数的敏感性系数。

Saltelli和Marivoet在1990年就详细地提出了多因素敏感性分析的方法。该方法首先对数据建立多元线性回归模型, 再利用下面公式分别求出每个输入的敏感性系数:

其中:PEAR (xi) 表示xi与Y的总体相关系数, σi和σj分别表示xi和Y的标准方差, Cov (Y, xi) 表示Y和xi的协方差, bi则表示Y关于xi的回归系数。

若用样本计算的协方差和标准差代替总体的协方差和标准差, 则为样本相关系数, 一般用r表示:

该方法的缺点是其只能用于线性模型, 当Y与x1, x2, …, xk之间的关系呈非线性形态时, 用r表示Y与xi的敏感性系数将与实际情况相差很大。

2. 基于Spearman秩相关系数的敏感性系数模型。

由于用Pearson相关系数表示的敏感性系数只能用于线性模型, 本文提出基于Spearman秩相关系数的敏感性系数的方法。计算Spearman秩相关系数时, 只要x和Y具有单调的函数关系的关系, 那么x和Y就是完全Spearman相关的, 即其绝对值为1。例如:Y=x10, 在区间[1, 2]内生成步长为0.01等差向量x, 然后计算Y=x10, 得到向量Y。可以计算出向量与向量的Pearson相关系数为0.840 8, 而向量x与向量Y的Spearman秩相关系数为1。

Spearman秩相关系数通常被认为是排列后的变量之间的Pearson线性相关系数:

其中:ρS表示xi和Y的Spearman秩相关系数。xi', Yi'表示原xi, Yi在排列后数据所在的位置, xi', Yi'称为变量xi, Yi的秩次。

Spearman秩相关系数的符号表示x和Y之间联系的方向, 其大小表示联系的紧密程度。当x和Y有严格单调增加的关系时, 它们之间的Spearman秩相关系数为1。反之, 在x和Y有严格单调减少的关系时, Spearman秩相关系数为-1。

所以用Spearman秩相关系数作为多因素敏感性系数指标是可行的, 当分别计算出敏感因素x1, x2, …, xk对输出Y的Spearman秩相关系数ρ1, ρ2, …, ρk后, 敏感因素x1, x2, …, xk对输出Y的秩相关系数ρ1, ρ2, …, ρk就是其对输出Y的Spearman敏感性系数。

3. 敏感性系数的标准化模型。

求出用Spearman秩相关系数表示的敏感性系数之后, 也就得到了输入xi对输出Y的影响程度。在实际管理中, 管理者还希望知道输入xi对输出Y影响程度的百分比及方向。本文提出两种方法:

第一种, 敏感性系数xi对绝对值之和贡献的百分比。

假设输入x1, x2, …, xk对输出Y的敏感性系数为ρ1, ρ2, …, ρk。

设输入xi的敏感性系数对绝对值之和贡献的百分比用mi表示, mi的求解公式如下:

mi的绝对值大小表示了输入xi对输出Y影响程度的百分比, 其符号表示输入xi对输出Y的影响方向。

第二种, 敏感性系数对平方和贡献的百分比。

设输入xi的敏感性系数对平方和贡献的百分比用Mi表示, 令:

Mi的求解公式如下:

类似的, Mi绝对值大小表示敏感因素xi对输出Y影响程度的百分比, 其符号表示敏感因素xi对输出Y的影响方向。

二、投资项目敏感性分析实例

假设某投资项目初始投资金额服从均值为1 000万元, 标准差为50万元的正态分布;项目使用期为5年, 期末残值为0元。期初投入的营运资本服从均值为50万元, 标准差为10万元的正态分布;项目贴现率为10%, 所得税税率为25%, 每年固定经营成本服从均值为150万元, 标准差为10万元的正态分布;每年销售量服从均值为200万元, 标准差为20万元的正态分布;销售价格服从均值为6元, 标准差为1元的正态分布;单位变动成本服从3元到4元之间的均匀分布, 销售量年增长率为8%。

首先建立这个投资项目的电子表格模型, 然后在Excel中模拟1 000次得到表1。

根据表1中NPV的数据可以计算出NPV的统计参数 (见表2) 和直方图。

根据表1的数据可以计算出其Spearman秩, 结果如下表:

根据公式 (1) 和公式 (3) 可以计算出投资项目的相关参数的敏感性系数表, 即表4。敏感性系数按照绝对值从大到小排序, 可以得到投资项目的敏感性系数的条形图。

从图3中可以清楚地看出销售价格对NPV的影响程度最大, 达到86.12%, 然后是单位变动成本为-7.96%。所以公司管理者应该确保项目建设工期如期或提前完工, 加大市场营销力度, 在竞争者产品进入市场前让产品以较高的利润销售出去。所以该投资项目产品的快速撇脂定价策略能否成功是该项目能否成功的关键。

下一步, 分析产品工艺, 降低产品单位变动成本, 从而使企业项目投资获得较好的回报。

三、结束语

在市场经济条件下, 企业投资项目面临诸多风险因素。为了确定风险因素的大小, 将Spearman秩相关系数引入到多因素敏感性系数确定中, 并在此基础上建立了敏感性系数的标准化模型。然后结合实例进行了敏感性分析。

由于敏感性分析在投资项目风险分析中用着广泛的应用, 所以基于Spearman秩相关系数建立的多因素敏感性系数方法有很强的现实意义。

参考文献

[1].王学强, 庄宇.基于蒙特卡罗模拟模型的投资项目风险分析.工业工程, 2007;10

[2].王跃虹.多因素敏感性分析的积分法.云南工业大学学报, 1999;2

[3].Burmaster, D.E., and P.D.Anderson.Principles ofGood Practice for the Use Of Monte Carlo Techniques inHuman Health and Ecological Risk Assessments.Risk Analysis, 1994

蒙特卡罗法模拟篇7

关键词：城市交通,蒙特卡洛,排队论

引言

南京市将于2014年举办青奥会, 随之大量游客将涌入该市。为缓解交通压力, 届时贯穿南京市中心区域的地铁3号线即将建成运营。受地铁3号线的影响, 公交系统应做相应调整, 以便于充分利用公交资源。而其中公交客流量的变化和乘客排队长度是影响公交线网布局的一个重要依据。高自友[1]利用拟动态均衡配流模型估计了不同时段内拥挤条件下公交客流量和乘客排队长度变化情况。四兵锋[2]则根据影响城市交通网络客流量分配的主要因素, 构造广义费用函数来分析乘客的路径选择行为, 并应用Logit模型, 通过搜索网络有效路径来解决客流的分配问题。本文在对南京公交及地铁运营现状分析的基础上, 基于蒙特卡洛模拟算法以及排队论仿真方法, 模拟地铁3号线建成开通后沿线公交站点的客流量和乘客平均候车时间等变化情况, 分析地铁开通后对公交客流的分流影响, 为改善公交布局以及地铁运营调整提供参考依据。

一、对当前公共交通系统的评价

为了对比分析地铁建设前后交通客流量变化, 这里选取3号线预设站点附近300m的122条线路、1 105个公交站作为研究对象, 并按照地铁站点所在地域将其分为交通汇聚地、商业区和名胜风景区。普通站点 (即不具备前两种特点的站点) 三类站。采用分层抽样方式选取了3个代表站点:南京火车站、夫子庙站、火炬南路站, 而其相邻公交站点:南京车站、夫子庙站、桃园站, 前两个站点基本与地铁站点位置吻合, 而火炬南路站为新开辟站点, 这里选取了最相近的且覆盖范围最大的桃园站作为分析对象。这里通过分时段抽样对上述3个站点进行实地调研, 获得了3个代表站客流量和平均等待时间、等待队长等资料。

(一) 地铁建成前公交站台客流量的分析

图1给出了3个代表站的客流量日逐小时变化图, 可以看到这3个代表性站点由于地理位置及周围环境的不同, 客流随时间的分布不尽相同。夫子庙站的客流量最大, 下午和傍晚时间段客流量明显高于上午, 这与其为南京市主要的商业区和风景区有关, 其客流变化特征与其他两站存在明显差异。南京车站公交平均客流量较少, 这是因为该站周围的公交站台多且密集, 每个站台所需服务的交通范围面积狭小, 另外已经开通的地铁1号线也分流了不少乘客。桃园站和南京车站的客流总体趋势基本一致, 6:00-9:00属于上班高峰期, 之后有所回落, 保持一段时间的稳定后, 晚间7:00后又再次进入高峰期。因桃园位于江北新区, 居民多需要乘车前往城区上班, 所以该区域存在客流量大、平均乘车时间长, 且上午高峰期要早于南京车站, 下午高峰期迟于南京车站的特点。

(二) 地铁建成前平均等待时间

平均等待时间为所有乘客从到达站台至上车所用时间的平均值, 考虑乘客心理承受极限, 这里只统计等待时间在20分钟以内的乘客数。表1给出了3个公交站点乘客平均等待时间不同时长的分布, 可以看出:南京车站、夫子庙站的乘客等待时间大多在4分钟内, 而桃园站的等待时间两极分化严重;3个站点的乘客平均等待时间分别为2.83分钟、2.36分钟、4.17分钟;南京车站为交通汇聚地, 客流虽大但车次较多, 乘客平均等待时间较短, 而以夫子庙站为代表的商业区和风景区客流大, 但站台分布密集, 故平均等待时间反而最短, 而桃园站则由于车次少, 人流量大平均等待时间最长。

二、地铁3号线建成后的预测与分析

(一) 乘车选择优化模型与蒙特卡洛仿真

地铁3号线建成后, 乘客可选择的交通方式将会发生改变, 这里用G (v, e) 表示包含地铁和公交两种乘车方式的南京市公共交通网络, 其中v, e分别为南京公交网络所有站点、线路的集合。现在假设某市民乘车起点为i, 终点为j, xij表示以起点为i, 终点为j的路线, 而T (xij) 为所对应的乘车时间, n为乘客所能容忍的最大换乘次数。在实际生活中, 当换乘次数超过一定数量时就会与乘客方便出行的心理相违背, 因此, 这里仅考虑换乘次数最多2次, 若超过2次, 则视为不可达。由此, 建立乘客乘车选择优化模型[3]为:

建立上述乘车选择优化模型后, 这里对南京地铁3号线建成后进行蒙特卡洛仿真[4], , 算法步骤如下:

(1) 确定起始站i, 随机模拟一个目的站j, 将含有站的路线存储在中x0, 含有j的路线存储在x1中;

(2) 除i站外, 搜寻x0与x1中是否有重复站点k。若k存在, 且k=j, 则i到j直达;若k≠j, 则视为一次换乘, 并将该路线存入r后转4) 。若k不存在, 则视i到j为两次换乘, 转3) ;

(3) 将经过i后面各个站点i1, i2, …, im的路线存入x3, 经过j前面的各个站点j1, j2, …, jn的路线存入x4, 依次将i1, i2, …, im作为起点站, j1, j2, …, jn作为终点站, 转2) ;

(4) 从r中选择换乘次数最小且时间短的路线为最佳路线;

(5) 最佳路线中, 若i站的下一站是地铁站, 选坐地铁, 否则选公交, 依此计算乘坐公交车的概率。

(二) 地铁建成后客流量的分析

运用上述算法, 分别以南京车站、夫子庙站、桃园站为起点模拟2 000次, 图2给出了建立地铁后公交站客流量仿真的结果, 客流量分布特征与建成前基本一致。地铁3号线建成后, 南京车站在7:00—8:00的客流量由地铁3号线建成前的240、245人次减到200人次左右, 而其余时段减少的客流量更是微乎其微, 这是因为南京车站是公交、地铁、出租等的交通枢纽区, 公交系统本身比较完善, 地铁线路的增加对公交的分流并不大。作为商业区和风景区的夫子庙站, 客流较大。地铁3号线建成后, 夫子庙站在上班高峰期8:00—9:00客流量相对于地铁建成前减少500人次, 同样处于客流高峰期的15:00—16:00大约减少100人次, 其余时间段公交客流也相应减少, 说明地铁3号线建成对夫子庙站客流分流较明显, 公交系统的压力明显减小。而远离市区且附近没有轨道交通的桃园站总客流量明显下降, 其中6:00—7:00的客流量由500人次减少为近300人次。这说明地铁3号线的建成对桃园站客流量的分流作用明显, 很大程度上缓解公交压力。总体上, 3个公交站点的客流量相对减少了16.20%, 18.91%, 31.30%。

(三) 地铁建成后平均等待时间及排队长

为进一步分析地铁建成对公交系统的影响, 本文根据公交的服务特点, 利用排队论[5]的相关知识选取乘客平均等待时间和排队长为影响公交系统的指标。由于经过同一站台的a条不同的线路相互独立, 则顾客的排队可看做是a个M/M/1系统。平均等待时间Wq是从乘客到达站台至上车所耗时间, 平均队长Lq是排队系统中未接受服务的人数。假设公交排队系统进入统计平衡状态, 系统中乘客的平均到达率λ以及乘客的平均服务率μ由调研可得 (见表2) 。

地铁3号线建成后, 南京车站、夫子庙站、桃园站将分别有16.2%、18.91%、31.30%的客流选择乘坐地铁, 故地铁建成后这3个公交站点顾客平均到达率变为112.579、300.025、144.257。借用R软件[6]对公交排队系统中乘客的平均等待时间进行客流量模拟, 得到地铁建成后的平均等待时间 (见表3) , 可以看出3个站点的平均等待时间均减少, 夫子庙站和桃园站尤为明显, 乘客平均等待1分钟即可上车。

根据Little公式 (Lq=λWq) , 计算得到各站点平均排队长 (见表4) 。地铁建成后, 南京车站平均排队长的减少幅度最小, 而桃园站减少幅度最大, 平均排队长减少率与客流量的分流有相似之处, 南京车站是交通枢纽, 增加一条

三、结论

通过上述分析, 不难看出, 地铁3号线的开通后公交站台客流量和平均等待时间和等待队长将大幅降低, 能够有效减轻和分流地铁沿线公交客流;同时可以看出, 地铁3号线的开通对各个区域公交站点的影响也存在着明显的区域差异, 公交线路调整和分流应该视所在区域区别对待。如南京火车站由于位于交通汇聚地, 虽客流量大但交通方便, 且附地铁对客流量的影响并不大, 自然乘客平均排队长的改变相对较小。相反, 桃园远离市区, 附近没有轨道交通, 增加一条运送客流量大的交通路线, 乘客的排队长会得到明显的减少。近公交站台密集, 地铁3号线的建成对其影响不大, 附近公交线路无需进行改动;而以夫子庙站为代表的商业区和名胜风景区特点是客流大相邻站台间隔长, 地铁建成后其公交客流明显减少, 于是可适当延长公交发车时间或减少部分公交班次。对于既不是交通汇聚地又不是商业区和名胜风景区的普通站点, 如火炬南路, 地铁的建成使其附近公交站客流急剧下降, 乘客排队长有效减少, 应采取增大相邻站台间隔、延长发车时间甚至取消离地铁站最近的公交站等措施。

参考文献

[1]高自友.拥挤条件下公交系统的拟动态均衡配流模型[J].交通运输系统工程与信息, 2002, (2) :38-45.

[2]四兵锋.无缝换乘条件下城市轨道交通网络客流分配模型及算法[J].铁道学报, 2007, (6) :12-18.

[3]周媛.基于蒙特卡洛模拟的地铁对公交线网布局影响分析[J].湖南交通科技, 2012, (1) :151-154.

[4]谢赤.大型水电工程造价风险评估及其关键因素识别[J].水利发电学报, 2010, (29) :63-68.

[5]胡运权, 郭耀煌.运筹学教程[M].北京:清华大学出版社, 2007:306-319.

【蒙特卡罗法模拟】推荐阅读：

蒙特卡罗模拟08-16