蒙特卡罗方法(Monte Carlo method)
目录
|
蒙特卡罗方法又称统计模拟法、随机抽样技术,是一种随机模拟方法,以概率和统计理论方法为基础的一种计算方法,是使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系,用电子计算机实现统计模拟或抽样,以获得问题的近似解。为象征性地表明这一方法的概率统计特征,故借用赌城蒙特卡罗命名。
蒙特卡罗方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法,为它蒙上了一层神秘色彩。在这之前,蒙特卡罗方法就已经存在。1777年,法国Buffon提出用投针实验的方法求圆周率∏。这被认为是蒙特卡罗方法的起源。
Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪,人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定圆周率π。本世纪40年代电子计算机的出现,特别是近年来高速电子计算机的出现,使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。
考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”,如何求出这个“图形”的面积呢?Monte Carlo方法是这样一种“随机化”的方法:向该正方形“随机地”投掷N个点,有M个点落于“图形”内,则该“图形”的面积近似为M/N。 可用民意测验来作一个不严格的比喻。民意测验的人不是征询每一个登记选民的意见,而是通过对选民进行小规模的抽样调查来确定可能的优胜者。其基本思想是一样的。
科技计算中的问题比这要复杂得多。比如金融衍生产品(期权、期货、掉期等)的定价及交易风险估算,问题的维数(即变量的个数)可能高达数百甚至数千。对这类问题,难度随维数的增加呈指数增长,这就是所谓的“维数的灾难”(Curse of Dimensionality),传统的数值方法难以对付(即使使用速度最快的计算机)。Monte Carlo方法能很好地用来对付维数的灾难,因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率,科学家们提出了许多所谓的“方差缩减”技巧。
另一类形式与Monte Carlo方法相似,但理论基础不同的方法—“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)—近年来也获得迅速发展。我国数学家华罗庚、王元提出的“华—王”方法即是其中的一例。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般可比Monte Carlo方法提出高数百倍,并可计算精确度。
由概率定义知,某事件的概率可以用大量试验中该事件发生的频率来估算,当样本容量足够大时,可以认为该事件的发生频率即为其概率。因此,可以先对影响其可靠度的随机变量进行大量的随机抽样,然后把这些抽样值一组一组地代入功能函数式,确定结构是否失效,最后从中求得结构的失效概率。蒙特卡罗法正是基于此思路进行分析的。
设有统计独立的随机变量Xi(i=1,2,3,…,k),其对应的概率密度函数分别为fx1,fx2,…,fxk,功能函数式为Z=g(x1,x2,…,xk)。
首先根据各随机变量的相应分布,产生N组随机数x1,x2,…,xk值,计算功能函数值 Zi=g(x1,x2,…,xk)(i=1,2,…,N),若其中有L组随机数对应的功能函数值Zi≤0,则当N→∞时,根据伯努利大数定理及正态随机变量的特性有:结构失效概率,可靠指标。
从蒙特卡罗方法的思路可看出,该方法回避了结构可靠度分析中的数学困难,不管状态函数是否非线性、随机变量是否非正态,只要模拟的次数足够多,就可得到一个比较精确的失效概率和可靠度指标。特别在岩土体分析中,变异系数往往较大,与JC法计算的可靠指标相比,结果更为精确,并且由于思路简单易于编制程序。
通常蒙特·卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题,蒙特·卡罗方法是一种有效的求出数值解的方法。一般蒙特·卡罗方法在数学中最常见的应用就是蒙特·卡罗积分。
蒙特卡罗方法在金融工程学,宏观经济学,生物医学,计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。
在解决实际问题的时候应用蒙特·卡罗方法主要有两部分工作:
1. 用蒙特·卡罗方法模拟某一过程时,需要产生各种概率分布的随机变量。
2. 用统计方法把模型的数字特征估计出来,从而得到实际问题的数值解。
使用蒙特·卡罗方法进行分子模拟计算是按照以下步骤进行的:
1. 使用随机数发生器产生一个随机的分子构型。
2. 对此分子构型的其中粒子坐标做无规则的改变,产生一个新的分子构型。
3. 计算新的分子构型的能量。
4. 比较新的分子构型于改变前的分子构型的能量变化,判断是否接受该构型。
5. 如此进行迭代计算,直至最后搜索出低于所给能量条件的分子构型结束。
蒙特卡罗方法是以概率统计原理为基础,模拟事物的形成过程,以达到认识事物特征及其变化规律的方法。这种方法的前提假设是不确定性参数可以用概率分布来描述。蒙特卡罗方法实施步骤:
1、通过敏感性分析,确定随机变量;
2、构造随机变量的概率分布模型;
3、为各输入随机变量抽取随机数;
4、将抽得的随机数转化为各输入随机变量的抽样值;
5、将抽样值组成一组项目评价基础数据;
6、根据基础数据计算出评价指标值;
7、整理模拟结果所得评价指标的期望值、方差、标准差和它的概率分布及累计概率,绘制累计概率分布图,计算项目可行或不可行的概率。
从理论上来说,蒙特卡罗方法需要大量的实验。实验次数越多,所得到的结果才越精确。以上Buffon的投针实验为例、历史上的记录如下表1。
从表中数据可以看到,一直到公元20世纪初期,尽管实验次数数以千计,利用蒙特卡罗方法所得到的圆周率∏值,还是达不到公元5世纪祖冲之的推算精度。这可能是传统蒙特卡罗方法长期得不到推广的主要原因。
计算机技术的发展,使得蒙特卡罗方法在最近10年得到快速的普及。现代的蒙特卡罗方法,已经不必亲自动手做实验,而是借助计算机的高速运转能力,使得原本费时费力的实验过程,变成了快速和轻而易举的事情。它不但用于解决许多复杂的科学方面的问题,也被项目管理人员经常使用。
借助计算机技术,蒙特卡罗方法实现了两大优点:
一是简单,省却了繁复的数学报导和演算过程,使得一般人也能够理解和掌握;
二是快速。简单和快速,是蒙特卡罗方法在现代项目管理中获得应用的技术基础。
蒙特卡罗方法有很强的适应性,问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛,因此问题维数的增加不会影响它的收敛速度,而且存贮单元也很省,这些是用该方法处理大型复杂问题时的优势。因此,随着电子计算机的发展和科学技术问题的日趋复杂,蒙特卡罗方法的应用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题,而且在统计物理、核物理、真空技术、系统科学 、信息科学 、公用事业、地质、医学,可靠性及计算机科学等广泛的领域都得到成功的应用。
项目管理中蒙特卡罗模拟方法的一般步骤是:
1、对每一项活动,输入最小、最大和最可能估计数据,并为其选择一种合适的先验分布模型;
2、计算机根据上述输入,利用给定的某种规则,快速实施充分大量的随机抽样;
3、对随机抽样的数据进行必要的数学计算,求出结果;
4、对求出的结果进行统计学处理,求出最小值、最大值以及数学期望值和单位标准偏差;
5、根据求出的统计学处理数据,让计算机自动生成概率分布曲线和累积概率曲线(通常是基于正态分布的概率累积S曲线);
非权重蒙特卡罗积分,也称确定性抽样,是对被积函数变量区间进行随机均匀抽样,然后对被抽样点的函数值求平均,从而可以得到函数积分的近似值。此种方法的正确性是基于概率论的中心极限定理。当抽样点数为m时,使用此种方法所得近似解的统计误差恒为 1除于根号M,不随积分维数的改变而改变。因此当积分维度较高时,蒙特卡罗方法相对于其他数值解法更优。
一、问题的提出
随着社会主义市场经济体制的逐步完善、经济水平的逐步提高,我国社会经济活动日趋复杂,越来越多变,其影响越来越广泛,越来越深远,不确定性逐渐成为企业决策时所面临的主要难题。因此,如何在不确定条件下做出投资决策,就成为目前理论和实践工作者们广泛关注的一个核心课题。
传统的投资评价理论——以净现值法(NPV) 为代表的投资决策分析方法,其根本缺陷在于它们是事先对未来的现金流量做出估计,并假设其为不变或静态的状况,无法衡量不确定因素的影响,不能体现递延决策以应对所带来的管理弹性。所以,在不确定环境下的投资,用净现值法评估项目不能体现柔性投资安排决策所体现的价值,无助于项目在决策中回避风险。在多变的市场环境中,不确定性与竞争者的反应使实际收入与预期收入有所出入, 所以净现值法(NPV) 适用于常规项目,未来不确定性比较小的项目。
为此理论界对未来投资环境不确定性大的项目提出了实物期权法,但在实践中应用的还是比较少。实物期权法的应用对企业决策者的综合素质要求比较高,对企业资源能力要求也比较高。但是实物期权法改变了我国管理者对战略投资的思维方式。
基于以上的分析,我们得出这样的结论:传统的投资决策方法对风险项目和不确定性项目的评价有较多不完善之处,有必要对其改进;实物期权法理论上解决了传统决策方法对不确定性项目评价的不足,但其应用尚处于体系不成熟阶段,在实践中应用并不广泛。至此,引入蒙特卡罗模型的理论和其分析方法,此方法特别适用于参数波动性大,且服从某一概率分布的项目,例如地质勘察、气田开发等项目。
蒙特卡罗模型是利用计算机进行数值计算的一类特殊风格的方法, 它是把某一现实或抽象系统的某种特征或部分状态, 用模拟模型的系统来代替或模仿, 使所求问题的解正好是模拟模型的参数或特征量, 再通过统计实验, 求出模型参数或特征量的估计值, 得出所求问题的近似解。目前评价不确定和风险项目多用敏感性分析和概率分析,但计算上较为复杂,尤其各因素变化可能出现概率的确定比较困难。蒙特卡罗模型解决了这方面的问题,各种因素出现的概率全部由软件自动给出,通过多次模拟,得出项目是否应该投资。该方法应用面广, 适应性强。
惠斯通(Weston) 对美国1 000 家大公司所作的统计表明: 在公司管理决策中, 采用随机模拟方法的频率占29 % 以上, 远大于其他数学方法的使用频率 。特别, 该方法算法简单, 但计算量大, 在模拟实际问题时, 要求所建模型必须反复验证,这就离不开计算机技术的帮助, 自然可利用任何一门高级语言来实现这种方法。通过一案例具体实现了基于Excel 的Monte Carlo 模拟系统, 由于Microsof tExcel 电子表格软件强大的数据分析功能和友好的界面设计能力, 使系统实现起来颇感轻松自如。
二、理论和方法
蒙特卡洛模拟早在四十年前就用于求解核物理方面的问题。当管理问题更为复杂时,传统的数学方法就难以进行了。模拟是将一个真实事物模型化,然后对该模型做各种实验,模拟也是一个通过实验和纠正误差来寻求最佳选择的数值性求解的过程。模拟作为一种有效的数值处理方法, 计算量大。以前只是停留在理论探讨上, 手工是无法完成的。在管理领域由于规律复杂随机因素多, 很多问题难以用线性数学公式分析和解决, 用模拟则有效得多。在新式的计算机普及后, 用模拟技术来求解管理问题已成为可能。
计算机模拟技术和其它方法相比有以下优点:
1) 成本低、风险小, 在产品未投产, 实际生产未形成就可以对市场进行分析模拟, 极大地减少费用和风险。
2) 环境条件要求低, 工作人员不需要高深的数学能力, 完全依靠计算机进行, 在硬件和软件日益降价的情况下, 可以成为现实。
3) 可信度高, 常用的统计推理方法需要大量历史数据(如平均数法、最小二乘法) , 对无历史资料的场合就无能为力(如新产品) , 而且精度低。
模拟的最大特点是借助一个随机数来模仿真实的现实, 随机数的产生则由计算机来产生。称为伪随机数。即:
Rn = F (r - 1 , r - 2 ,……r - k)
在以对象为中心的软件中, EXCEL 有一个RANE()函数实现伪随机数功能。RANE( )实际上是一个会自动产生伪随机数的子程序。用产生的伪随机数模拟市场购买行为, 得出产品销售量, 在生产成本相对固定时进而推测出产品的利润。此方法不用编制复杂的程序, 思路假设为, 作为系统内部是可以控制的, 即企业内部生产成本可以人为控制, 但系统外部因素是不可控制的(消费心理导致的消费行为) , 则生产与销售就会产生矛盾。生产量小于销售量, 造成开工不足资源浪费;生产量大于销售量, 造成产品积压, 资金占用, 同样形成资源的浪费。最好生产量等于销售量, 则资源浪费最小, 自然经济效益就最高, 实际就是利润最大化。如果能科学地测算出在什么情况下利润最大, 则这时的产量就是最佳产量, 成本也就最低。这就是市场作为导向, 以销定产的公认市场经济的准则。实际工作中, 很多产品的消费是具有随机性的, 主要是一些需求弹性大、价格弹性大、价格低、与日常生活有关的中、小商品, 如副食品、日用消费品、玩具、轻工业产品。对企业而言利润较高的产品。
从以上分析可以看出, 蒙特卡洛模拟可以动态实现对产品利润的预测, 从而对产品产量科学控制,实现资源优化, 是一种较好的决策支持方法。
三、蒙特卡罗模型在Excel 表中的应用
某气田投资项目期投资、寿命期、残值以及各年的收入、支出,以及应付税金的税率、项目的资本成本等都是独立的随机变量,他们的概率密度函数如表1所示。
表 各变量对应概率密度函数表
A | B | C | D | |
2 | 概率 | 对应的随机数 | 可能值 | |
3 | 投资Yo | 0.2 | 0 | 450 |
4 | 0.5 | 20 | 500 | |
5 | 0.3 | 70 | 550 | |
6 | 寿命N | 0.5 | 0 | 6 |
7 | 0.3 | 50 | 7 | |
8 | 0.2 | 80 | 8 | |
9 | 残值F | 0.25 | 0 | 40 |
10 | 0.5 | 25 | 50 | |
11 | 0.25 | 75 | 60 | |
12 | 税率Te | 0.2 | 0 | 45 |
13 | 0.5 | 20 | 48 | |
14 | 0.3 | 70 | 51 | |
15 | 年收入R | 0.15 | 0 | 700 |
16 | 0.3 | 15 | 750 | |
17 | 0.4 | 45 | 800 | |
18 | 0.15 | 85 | 850 | |
19 | 年支出C | 0.2 | 0 | 100 |
20 | 0.4 | 20 | 150 | |
21 | 0.3 | 60 | 200 | |
22 | 0.1 | 90 | 250 | |
23 | 资本成本i | 0.1 | 0 | 10 |
24 | 0.2 | 10 | 12 | |
25 | 0.4 | 30 | 14 | |
26 | 0.2 | 70 | 16 | |
27 | 0.1 | 90 | 18 |
本案例用windowsXP 中的Excel2003 对该项目进行模拟如下:
1) 在A32 单元格(投资Yo 模拟:随机数) 输入:= RANDBETWEEN (0 ,99) ;在B32 单元格(投资Yo模拟:投资) 输入: = VLOO KUP (A32 , $C $3 : $D$5 ,2) ;
2) 在C32 单元格(寿命N 模拟:随机数) 输入: =RANDBETWEEN (0 ,99) ;在D32 单元格(寿命N 模拟: 寿命) 输入: = VLOO KUP ( C32 , $C $6 : $D$8 ,2) ;
3) E32 ,G32 , I32 , K32 ,M32 单元格分别输入: =RANDBETWEEN (0 , 99) ; F32 = VLOOPUP ( E32 ,$C $9 : $D $11 , 2) , H32 = VLOOPUP ( G32 , $C$12 : $D $14 ,2) ,J 32 = VLOO KUP ( I32 , $C $15 :$D $18 ,2) ,L32 = VLOO KUP ( K32 , $C $19 : $D$22 ,2) ,
N32 = VLOO KUP(M32 , $C $23 : $D $27 ,2)
4) O32 = (B32 - F32) / D32 , P32 = (J 32 - L32 -O32) * (1 - H32/ 100) + O32 ,Q32 = PV (N32/ 100 ,D32 , - P32) - B32 ;
5) H3 = AVERA GE ( Q32 , Q5031 ) , H4 =STDEV (Q32 ,Q5031) ,H5 = MAX ( Q32 , Q5031 ) , H6 = MIN ( Q32 ,Q5031) ,H7 = H4/ H3 ,H8 = COUN TIF (Q32 :Q5031 ,“ < 0”) / COUN T(Q32 ,Q5031) 。
在Excel 工具表中模拟5000次,结果输出见下表 :
表 结果输出表(1)
A | B | C | D | E | F | G | H | |
投资Yo模拟 | 寿命N模拟 | 残值F模拟 | 税率Te模拟 | |||||
随机数 | 投资 | 随机数 | 寿命 | 随机数 | 残值 | 随机数 | 税率 | |
32 | 17 | 450 | 78 | 7 | 51 | 50 | 2 | 45 |
33 | 31 | 500 | 84 | 8 | 87 | 60 | 67 | 48 |
34 | 22 | 500 | 63 | 7 | 97 | 60 | 88 | 51 |
35 | 95 | 550 | 70 | 7 | 40 | 50 | 81 | 51 |
36 | 31 | 500 | 96 | 8 | 20 | 40 | 12 | 45 |
37 | 16 | 450 | 1 | 6 | 41 | 50 | 66 | 48 |
38 | 79 | 550 | 33 | 6 | 87 | 60 | 51 | 48 |
39 | 0 | 450 | 97 | 8 | 78 | 60 | 17 | 45 |
40 | 35 | 500 | 43 | 6 | 22 | 40 | 5 | 45 |
41 | 3 | 450 | 70 | 7 | 52 | 50 | 87 | 51 |
42 | 78 | 550 | 39 | 6 | 69 | 50 | 30 | 48 |
43 | 20 | 500 | 36 | 6 | 90 | 60 | 2 | 45 |
44 | 96 | 550 | 5 | 6 | 92 | 60 | 40 | 48 |
45 | 51 | 500 | 36 | 6 | 90 | 60 | 20 | 48 |
46 | 58 | 500 | 39 | 6 | 1 | 40 | 11 | 45 |
47 | 4 | 450 | 79 | 7 | 22 | 40 | 29 | 48 |
48 | 83 | 550 | 36 | 6 | 40 | 50 | 62 | 48 |
… | … | … | … | … | … | … | … | … |
表 结果输出表(2)
I | J | K | L | M | N | |
年收入R模拟 | 年支出C模拟 | 资本成本i模拟 | ||||
随机数 | 年收入 | 随机数 | 年支出 | 随机数 | 资本成本 | |
32 | 12 | 700 | 88 | 200 | 4 | 10 |
33 | 11 | 700 | 88 | 200 | 59 | 14 |
34 | 3 | 700 | 79 | 200 | 7 | 10 |
35 | 68 | 800 | 20 | 150 | 77 | 16 |
36 | 23 | 750 | 21 | 150 | 53 | 14 |
37 | 98 | 850 | 73 | 200 | 40 | 14 |
38 | 37 | 750 | 23 | 150 | 99 | 18 |
39 | 72 | 800 | 92 | 250 | 16 | 12 |
40 | 81 | 800 | 96 | 250 | 46 | 14 |
41 | 32 | 750 | 17 | 100 | 74 | 16 |
42 | 70 | 800 | 73 | 200 | 17 | 12 |
43 | 39 | 750 | 78 | 200 | 68 | 14 |
44 | 12 | 700 | 46 | 150 | 92 | 18 |
45 | 79 | 800 | 75 | 200 | 15 | 12 |
46 | 10 | 700 | 52 | 150 | 54 | 14 |
47 | 45 | 800 | 1 | 100 | 87 | 16 |
48 | 75 | 800 | 47 | 150 | 4 | 10 |
… | … | … | … | … | … | … |
表 结果输出表(3)
O | P | Q | |
折旧Dt | 各年现金流量Yt | NPV | |
32 | 75 | 307175 | 840.3314803 |
33 | 76166667 | 348.18 | 856.3672298 |
34 | 64128571 | 342.8571 | 1064.716528 |
35 | 75 | 374 | 878.0912297 |
36 | 55.71429 | 364.7429 | 1114.128559 |
37 | 62.185714 | 368.1714 | 986.8844068 |
38 | 73.33333 | 355.9 | 883.9767691 |
39 | 58.57143 | 299.3714 | 1007.465496 |
40 | 83.33333 | 336.5 | 689.9136332 |
41 | 83.33333 | 326 | 717.7056104 |
42 | 66.66667 | 344 | 964.3241193 |
43 | 85 | 35218 | 749.9748285 |
44 | 57.14286 | 347.6429 | 1040.798547 |
45 | 64.28571 | 276.4286 | 761.5527004 |
46 | 76.166667 | 348.8 | 785.2358848 |
47 | 48.75 | 335.4 | 1105.87495 |
48 | 57.5 | 391.6 | 1200.950194 |
… | … | … | … |
所得结果如下:
表 净现值模拟计算结果表
F | G | H | |
2 | 净现值模拟计算结果 | ||
3 | 净现值期望值 | 952.13017 | |
4 | 净现值标准差 | 198.90501 | |
5 | 净现值最大值 | 1726.9833 | |
6 | 净现值最小值 | 405.54502 | |
7 | 变异系数 | 0.12089053 | |
8 | 净现值为负的概率 | 0 |
表 净现值概率分布统计表
净现值概率分布统计 | |||
系统分组 | 分布区间 | 概率 | 累计概率 |
300 | 3 | 以下0 | 0 |
400 | 3~4 | 0 | 0 |
500 | 4~5 | 0.0036 | 0.0036 |
600 | 5~6 | 0.0244 | 0.028 |
700 | 6~7 | 0.062 | 0.09 |
800 | 7~8 | 0.1322 | 0.2222 |
900 | 8~9 | 0.1898 | 0.412 |
1000 | 9~10 | 0.1992 | 0.6112 |
1100 | 10~11 | 0.1628 | 0.774 |
1200 | 11~12 | 0.1162 | 0.8902 |
1300 | 12~13 | 0.0548 | 0.945 |
1400 | 13~14 | 0.0338 | 0.9788 |
1500 | 14~15 | 0.0132 | 0.992 |
1600 | 15~16 | 0.005 | 0.997 |
1700 | 16~17 | 0.0018 | 0.9988 |
1800 | 17以上 | 0.0012 | 1 |
从分析结果得出,虽然此项目未来的不确定性很大,但由图可知,此气田开发项目服从正态分布,模拟5 000次的结果是净现值为负的概率为零,并且项目的期望净现值为952113 万元,说明项目值得开发。
由以上的案例分析可知,基于蒙特卡罗模拟的风险分析,对于工程实际应用具有较强的参考价值。随机模拟5 000 次,如果仅靠人的大脑进行计算,这在现实世界中是不可能的,但考虑到系统决策支持功能, 算法设计为由使用者自己设计方案, 采用人机交互, 这样可以发挥使用者的经验判断;系统实现模拟运算——系统对每一个设定的投资项目期投资、寿命期、残值以及各年的收入、支出,以及应付税金的税率、项目的资本成本等随机变量及他们的概率密度函数,通过蒙特卡罗模拟方法,得出了项目在不同概率发生的情况下净现值模拟计算结果。为人们解决不确定性项目的决策提供了简单的方法,节约了人们的工作量和时间。但是利用蒙特卡罗模型分析问题时,收集数据是非常关键的。