1、模拟退火（Simulated Annealing）、进化策略（Evolution strategy）等等。一个最简单的例子是，已知某函数，我们要求此函数的最大值，那么我们可以不断地在该函数定义域上随机取点，然后用得到的最大的点作为此函数的最大值。这个例子实质也是随机数值积分，它等价于求此函数的无穷阶范数（-Norm）在定义域上的积分。由于在金融产品定价中，这部分内容用的相对较不常见，所以此课程就不介绍随机最优化方法了。十二、十二、MontMonte Carloe Carlo 形式与一般步骤形式与一般步骤（一）积分形式（一）积分形式做 Monte Carlo 时，求解积分的一般形式是：X 为自变
2、量，它应该是随机的，定义域为（x0,x1），f（x）为被积函数，（x）是 x 的概率密度。在计算欧式期权例子中，x 为期权到期日股票价格，由于我们计算期权价格的时候该期权还没有到期，所以此时 x 是不确定的（是一随机变量），我们按照相应的理论，假设 x的概率密度为（x）、最高可能股价为 x1（可以是正无穷）、最低可能股价为 x0（可以是 0），另外，期权收益是到期日股票价格 x 和期权行权价格的函数，我们用 f（x）来表示期权收益。（二）一般步骤（二）一般步骤我将 Monte Carlo 分为三加一个步骤：1依据概率分布（x）不断生成随机数 x,并计算 f（x）由于随机数性质，每次生成的 x
3、的值都是不确定的，为区分起见，我们可以给生成的 x赋予下标。如 xi表示生成的第 i 个 x。生成了多少个 x，就可以计算出多少个 f（x）的值2将这些 f（x）的值累加，并求平均值例如我们共生成了 N 个 x，这个步骤用数学式子表达就是3到达停止条件后退出常用的停止条件有两种，一种是设定最多生成 N 个 x，数量达到后即退出，另一种是检测计算结果与真实结果之间的误差，当这一误差小到某个范围之内时退出。有趣的类比：积分表达式中的积分符合类比为上式中累加符号，dx 类比为 1/N（数学知识告诉我们积分实质是极限意义下的累加；f（x）还是它自己，积分中的（x）可类比为依据（x）生成随机数4误差分析
4、Monte Carlo 方法得到的结果是随机变量，因此，在给出点估计后，还需要给出此估计值的波动程度及区间估计。严格的误差分析首先要从证明收敛性出发，再计算理论方差，最后用样本方差来替代理论方差。在本课程中我们假定此方法收敛，同时得到的结果服从正态分布，因此可以直接用样本方差作区间估计。详细过程在例子中解释。这个步骤的理论意义很重要，但在实际应用中，它的重要性有所淡化，倘若你的老板不太懂这些知识，你报告计算结果时可以只告诉他点估计即可。注意，前两大步骤还可以继续细分，例如某些教科书上的五大步骤就是将此处的前两步细分成四步。十三、最简单的例子十三、最简单的例子举个例子:计算从函数从 0 到 2
5、的定积分值。数学方法：我们已知的原函数是，那么定积分值就是：=6.38905609893065。计算这个数值可以在 Matlab 中输入代码：exp（2）-exp（0）上面得到的值是此不定积分的真实值。常规数值积分：在区间内取 N 个点，计算各个点上的函数值，然后用函数值乘以每个区间宽度，最后相加。Matlab 代码：N=100;x=linspace（0,2,N）;sum（exp（x）.*2/N）试着调大 N 的值，你会发现，最后的结果将更接近真实值。Monte Carlo 数值积分法：在内随机取 N 个点，计算各个点上的函数值，最后求这些函数值的平均值再乘以 2（为何要乘以 2 在后面小节详
6、细讲）。看 Matlab 代码：x=unifrnd（0,2,N,1）;mean（2*exp（x）同样的，通过增大 N，这种方法得到的结果也将越来越接近真实值。解释解释这个例子要求的积分形式是：，还不完全是形式，我们先做变换，这里是 f（x）；1/2 是（x），它表示，在取值范围（0,2）区间内，x 服从均匀分布。前一例子共三条语句，逐句解释如下：设定停止条件，共做 N 次 Monte Carlo 模拟。按照（0,2）区间均匀分布概率密度对 x 随机抽样，共抽取 N 个 xi。此句相当于第一个步骤中的前半部分。mean（2*exp（x）2*exp（x）作用是对每个 xi计算 f（xi）的值,共可
7、得到 N 个值，这个相当于第一个步骤后半部分；Mean（）函数的作用是将所有的 f（xi）加起来取平均值，相当于第二个步骤。这段代码中的停止条件隐含于 N 值设定中，它一次性生成 N 个 x 值，完成此次计算后整个程序就结束了。十四、十四、Monte CarlMonte Carlo o 方法的优点方法的优点对比前面常规数值积分和 Monte Carlo 数值积分代码，同样数量的 N 值也就意味这几乎相同的计算量常规数值积分结果的精确度要高于 Monte Carlo 数值积分的结果。那么，我们为何还需要用 Monte Carlo 来算数值积分呢？答案的关键在于，常规数值积分的精度直接取决于每个维
8、度上取点数量，维度增加了，但是每个维度上要取的点却不能减少。在多重积分中，随着被积函数维度增加，需要计算的函数值数量以指数速度递增。例如在一重积分中，只要沿着 x 轴取 N 个点；要达到相同大小的精确度，在 s 重积分中，仍然需要在每个维度上取 N 个点，s 个纬度的坐标相组合，共需要计算 Ns个坐标对应的 f（）函数值。取点越多，会占用计算机大量内存，也需要更长运算时间，最终导致这种计算方法不可行！Monte Carlo 方法却不同，不管是积分有多少重，取 N 个点计算的结果精确度都差不多。因此，即使在一重积分的情形下，Monte Carlo 方法的效率比不过常规数值积分，但随着积分维度增加
9、，常规数值积分的速度呈指数下降，Monte Carlo 方法的效率却基本不变。经验表明，当积分重数达到 4 重积分甚至更高时，Monte Carlo 方法将远远优于常规数值积分方法。现在回到金融产品定价，欧式期权理论定价公式只需要一重积分，此时 Monte Carlo 方法的效果不明显，但是如果我们考虑一个亚式期权：期限为 1 年期，期权价格基于此 1 年内每天某个时点时的价格，全年共 252 个交易日，这样此亚式期权理论定价公式是一个 252重积分。常规的数值积分方法，需要取 N252个点，这个数有多大，你自己去计算一下就知道了（注意：N 取值要远远大于 2），常规数值积分方法不可行，只能用
10、 Monte Carlo。综上，如果计算高维度多重积分，如路径依赖的 exotic options（奇异期权）等金融产品定价，我们一般用的方法都是 Monte Carlo。十十五、五、Monte CarloMonte Carlo 方法原理方法原理（选读选读）Monte Carlo 方法计算的结果收敛的理论依据来自于大数定律，且结果渐进地（Asymptotically）服从正态分布的理论依据是中心极限定理。以上两个属性都是渐进性质，要进行很多次抽样，此属性才会比较好地显示出来，如果Monte Carlo 计算结果的某些高阶距存在，即使抽样数量不太多，这些渐进属性也可以很快地达到。这些原理在理论上
11、意义重大，但由于我们一般遇上的 Monte Carlo 问题都是收敛的、结果也都是渐进正态分布，所以工作中使用时可以不加考虑。详细推导见相关书籍。第二章：随机数的生成第二章：随机数的生成讲课人：Xaero Chang|课程主页:http:/macro2.org/notes/intro2mc本章第一节会简要复习随机变量的一些概念，但学习本章最好要有一定的数学基础。第二节主要介绍如何生成一维概率分布的随机数，第三节介绍如何生成高维分布的随机数。最后略提及伪随机数问题的应对策略。由前文可知，Monte Carlo 积分解决的问题形如，f（x）值只需由 x 值决定，因此此处最重要的就是如何生成服从（x
12、）概率分布的随机数。可以说，正确生成随机数，Monte Carlo 方法就做完了一半。一、随机变量基本概念一、随机变量基本概念（一）随机变量（一）随机变量现实世界中有很多可以用数字来衡量的事物，站在当前时间点来看，它们在未来时刻的值是不确定的。例如，我们掷一骰子，在它停稳前，我们不可能知道掷出多少点（传说中的赌王除外，哈哈）；例如某只股票在明天的股价，没有人能准确知晓第二天股票的价格（不然他就发惨了！）。但是，我们却可以描述这些事物未来各种值的可能性。（二）离散型随机变量（二）离散型随机变量离散型随机变量最重要的是分布律，即每个取值的概率是多少。例如掷骰子，我们认为扔出任何一个点的概率都是 1
13、/6。那么掷骰子得到的点数的分布律如下表：骰子点数123456概率1/61/61/61/61/61/6（三）连续性随机变量（三）连续性随机变量连续型随机变量有两个重要的概念。概率密度函数（PDF）和累积概率分布函数（CDF），具体定义见数学书籍。PDF 函数本身不是概率，只有对 x 的某段区间中的 PDF 积分得到的数值才有概率的含义。CDF 是概率的意思，点 x 上 CDF 的值表示该随机变量可能取值小于 x 的概率的大小。如图是正态分布的 PDF 和 CDF（四）多元分布（四）多元分布在这个课程里面，我不将多维的随机变量拆分成多个一维的随机变量来表述，而是将各个维度组合成一个随机变量向量。
14、多元随机变量分布需要掌握联合分布、边缘分布和条件分布。联合分布和单变量 PDF 类似，如果将其对随机变量某个取值范围做积分就可以得到随机变量最终取值落在该区间内的概率。边缘分布是只考虑多维随机变量中的某一维，其他维度不考虑情况下的 PDF条件分布是固定随机变量在其他维度的取值，再考虑剩余那个维度上的 PDF。如图是一个二维正态分布（两个维度间相关系数为 0.3）的示意图。两个小图分别是第一和第二维度的边缘分布 PDF 图（都是标准正态分布 PDF）。右上角的大图是依据此二维正态联合分布生成的随机数。从随机数的疏密程度可以看出联合分布 PDF 函数在该区域的大小。研究这些分布不是我们的目的，只是
15、达到目的的手段。我们所需要做的事情是生成符合各种分布的随机数。由分布的不同类型连续型和离散型，常规型（Matlab 中有内置函数）和特殊型（Matlab 中无内置函数），一维分布和多维分布我们接下来就各种情况分别讲述如何在 Matlab 中生成各种各样的随机数。十六、一维随机数十六、一维随机数两个注意事项：1生成随机数有两种选择，可以每次只生成一个随机数，直接用此数计算 f（x）,然后循环重复此过程，最后求平均值；另一种方法是每次生成全部循环所需的随机数，利用 Matlab矩阵运算语法计算 f（x），不需要写循环，直接即可求平均值。前一种方法代码简单，但速度慢；后一种方法代码相对更难写。此处一
16、定要掌握如何生成随机数组成的向量和矩阵（单个随机数就是一个 1*1 的矩阵），在后面章节的例子里面一般会有两个版本的代码分别讲述这两种生成方法。2生成了一维的随机数后，可以用 hist（）函数查看这些数服从的大致分布情况。（一）（一）MatlabMatlab 内部函数内部函数a.基本随机数基本随机数Matlab 中有两个最基本生成随机数的函数。1rand（）rand（）生成（0,1）区间上均匀分布的随机变量。基本语法：rand（M,N,P.）生成排列成 M*N*P.多维向量的随机数。如果只写 M，则生成 M*M 矩阵；如果参数为M,N可以省略掉方括号。一些例子：rand（5,1）%生成 5 个
17、随机数排列的列向量，一般用这种格式rand（5）%生成 5 行 5 列的随机数矩阵rand（5,4）%生成一个 5 行 4 列的随机数矩阵生成的随机数大致的分布。x=rand（100000,1）;hist（x,30）;由此可以看到生成的随机数很符合均匀分布。（视频教程会略提及 hist（）函数的作用）2randn（）randn（）生成服从标准正态分布（均值为 0，方差为 1）的随机数。基本语法和 rand（）类似。randn（M,N,P.）生成排列成 M*N*P.多维向量的随机数。randn（5,1）%生成 5 个随机数排列的列向量，一般用这种格式randn（5）%生成 5 行 5 列的随机数
18、矩阵randn（5,4）%生成一个 5 行 4 列的随机数矩阵生成的随机数大致的分布。x=randn（100000,1）;hist（x,50）;由图可以看到生成的随机数很符合标准正态分布。b.连续型分布随机数连续型分布随机数如果你安装了统计工具箱（Statistic Toolbox），除了这两种基本分布外，还可以用 Matlab内部函数生成符合下面这些分布的随机数。3unifrnd（）unifrnd（）和 rand（）类似，这个函数生成某个区间内均匀分布的随机数。基本语法unifrnd（a,b,M,N,P,.）生成的随机数区间在（a,b）内，排列成 M*N*P.多维向量。如果只写 M，则生成
19、M*M矩阵；unifrnd（-2,3,5,1）%生成 5 个随机数排列的列向量，一般用这种格式unifrnd（-2,3,5）%生成 5 行 5 列的随机数矩阵unifrnd（-2,3,5,4）%生成一个 5 行 4 列的随机数矩阵%注：上述语句生成的随机数都在（-2,3）区间内.生成的随机数大致的分布。x=unifrnd（-2,3,100000,1）;由图可以看到生成的随机数很符合区间（-2,3）上面的均匀分布。4normrnd（）normrnd（）和 randn（）类似，此函数生成指定均值、标准差的正态分布的随机数。基本语法normrnd（mu,sigma,M,N,P,.）生成的随机数服从均
20、值为 mu，标准差为 sigma（注意标准差是正数）正态分布，这些随机数排列成 M*N*P.多维向量。normrnd（2,3,5,1）%生成 5 个随机数排列的列向量，一般用这种格式normrnd（2,3,5）%生成 5 行 5 列的随机数矩阵normrnd（2,3,5,4）%生成一个 5 行 4 列的随机数矩阵%注：上述语句生成的随机数所服从的正态分布都是均值为 2，标准差为 3.生成的随机数大致的分布。x=normrnd（2,3,100000,1）;如图，上半部分是由上一行语句生成的均值为 2，标准差为 3 的 10 万个随机数的大致分布，下半部分是用小节“randn（）”中最后那段语句生
21、成 10 万个标准正态分布随机数的大致分布。注意到上半个图像的对称轴向正方向偏移（准确说移动到 x=2 处），这是由于均值为 2的结果。而且，由于标准差是 3，比标准正态分布的标准差（1）要高，所以上半部分图形更胖（注意 x 轴刻度的不同）。5chi2rnd（）chi2rnd（）此函数生成服从卡方（Chi-square）分布的随机数。卡方分布只有一个参数：自由度 v。基本语法chi2rnd（v,M,N,P,.）生成的随机数服从自由度为 v 的卡方分布，这些随机数排列成 M*N*P.多维向量。chi2rnd（5,5,1）%生成 5 个随机数排列的列向量，一般用这种格式chi2rnd（5,5）%生
22、成 5 行 5 列的随机数矩阵chi2rnd（5,5,4）%生成一个 5 行 4 列的随机数矩阵%注：上述语句生成的随机数所服从的卡方分布的自由度都是 5生成的随机数大致的分布。x=chi2rnd（5,100000,1）;6frnd（）frnd（）此函数生成服从 F 分布的随机数。F 分布有 2 个参数：v1,v2。基本语法frnd（v1,v2,M,N,P,.）生成的随机数服从参数为（v1,v2）的卡方分布，这些随机数排列成 M*N*P.多维向量。frnd（3,5,5,1）%生成 5 个随机数排列的列向量，一般用这种格式frnd（3,5,5）%生成 5 行 5 列的随机数矩阵frnd（3,5,
23、5,4）%生成一个 5 行 4 列的随机数矩阵%注：上述语句生成的随机数所服从的参数为（v1=3,v2=5）的 F 分布生成的随机数大致的分布。x=frnd（3,5,100000,1）;从结果可以看出来，F 分布集中在 x 正半轴的左侧，但是它在极端值处也很可能有一些取值。7trnd（）trnd（）此函数生成服从t（Students t Distribution，这里Student不是学生的意思，而是Cosset.W.S.的笔名）分布的随机数。t 分布有 1 个参数：基本语法trnd（v,M,N,P,.）生成的随机数服从参数为 v 的 t 分布，这些随机数排列成 M*N*P.多维向量。trnd
24、（7,5,1）%生成 5 个随机数排列的列向量，一般用这种格式trnd（7,5）%生成 5 行 5 列的随机数矩阵trnd（7,5,4）%生成一个 5 行 4 列的随机数矩阵%注：上述语句生成的随机数所服从的参数为（v=7）的 t 分布生成的随机数大致的分布。x=trnd（7,100000,1）;可以发现 t 分布比标准正太分布要“瘦”，不过随着自由度 v 的增大，t 分布会逐渐变胖，当自由度为正无穷时，它就变成标准正态分布了。接下来的分布相对没有这么常用，同时这些函数的语法和前面函数语法相同，所以写得就简略一些在视频中也不会讲述，你只需按照前面那几个分布的语法套用即可，应该不会有任何困难时间
25、足够的话这是一个不错的练习机会。8betarnd（）betarnd（）此函数生成服从 Beta 分布的随机数。Beta 分布有两个参数分别是 A 和 B。下图是A=2,B=5 的 beta 分布的 PDF 图形。生成 beta 分布随机数的语法是：betarnd（A,B,M,N,P,.）9exprnd（）exprnd（）此函数生成服从指数分布的随机数。指数分布只有一个参数:mu,下图是 mu=3 时指数分布的 PDF 图形生成指数分布随机数的语法是：betarnd（mu,M,N,P,.）10gamrnd（）gamrnd（）生成服从 Gamma 分布的随机数。Gamma 分布有两个参数：A 和
26、B。下图是 A=2,B=5 Gamma 分布的 PDF 图形生成 Gamma 分布随机数的语法是：gamrnd（A,B,M,N,P,.）11lognrnd（）lognrnd（）生成服从对数正态分布的随机数。其有两个参数：mu 和 sigma，服从这个这样的随机数取对数后就服从均值为 mu，标准差为 sigma 的正态分布。下图是 mu=-1,sigma=1/1.2 的对数正态分布的 PDF 图形。生成对数正态分布随机数的语法是：lognrnd（mu,sigma,M,N,P,.）12raylrnd（）raylrnd（）生成服从瑞利（Rayleigh）分布的随机数。其分布有 1 个参数：B。下图是
27、 B=2 的瑞利分布的 PDF 图形。生成瑞利分布随机数的语法是：raylrnd（B,M,N,P,.）13wblrnd（）wblrnd（）生成服从威布尔（Weibull）分布的随机数。其分布有 2 个参数：scale 参数 A 和 shape 参数 B。下图是 A=3，B=2 的 Weibull 分布的 PDF 图形。生成 Weibull 分布随机数的语法是：wblrnd（A,B,M,N,P,.）还有非中心卡方分布（ncx2rnd），非中心 F 分布（ncfrnd），非中心 t 分布（nctrnd），括号中是生成服从这些分布的函数，具体用法用：help 函数名查找。c.离散型分布随机数离散型分
28、布随机数离散分布的随机数可能的取值是离散的，一般是整数。14unidrnd（）unidrnd（）此函数生成服从离散均匀分布的随机数。Unifrnd 是在某个区间内均匀选取实数（可为小数或整数），Unidrnd 是均匀选取整数随机数。离散均匀分布随机数有 1 个参数：n,表示从1,2,3,.N这 n 个整数中以相同的概率抽样。unidrnd（n,M,N,P,.）这些随机数排列成 M*N*P.多维向量。unidrnd（5,5,1）%生成 5 个随机数排列的列向量，一般用这种格式unidrnd（5,5）%生成 5 行 5 列的随机数矩阵unidrnd（5,5,4）%生成一个 5 行 4 列的随机数矩阵%注：上述语句生成的随机数所服从的参数为（10,0.3）的二项分布生成的随机数大致的分布。x=unidrnd（9,100000,1）;hi
