蒙特卡罗算法与matlab（精品教程）

第一章：Monte Carlo方法概述

讲课人：Xaero Chang | 课程主页: http://macro2.org/notes/intro2mc

本章主要概述Monte Carlo的一些基础知识，另外包括一个最简单的用Monte Carlo方法计算数值积分的例子。

一、Monte Carlo历史渊源

Monte Carlo方法的实质是通过大量随机试验，利用概率论解决问题的一种数值方法，基本思想是基于概率和体积间的相似性。它和Simulation有细微区别。单独的Simulation只是模拟一些随机的运动，其结果是不确定的；Monte Carlo在计算的中间过程中出现的数是随机的，但是它要解决的问题的结果却是确定的。

历史上有记载的Monte Carlo试验始于十八世纪末期（约1777年），当时布丰（Buffon）为了计算圆周率，设计了一个“投针试验”。（后文会给出一个更加简单的计算圆周率的例子）。虽然方法已经存在了200多年，此方法命名为Monte Carlo则是在二十世纪四十年，美国原子弹计划的一个子项目需要使用Monte Carlo方法模拟中子对某种特殊材料的穿透作用。出于保密缘故，每个项目都要一个代号，传闻命名代号时，项目负责人之一von Neumann灵犀一点选择摩洛哥著名赌城蒙特卡洛作为该项目名称，自此这种方法也就被命名为Monte Carlo方法广为流传。

十一、Monte Carlo方法适用用途 （一）数值积分

计算一个定积分，如，如果我们能够得到f(x)的原函数F(x)，那么直接由表达式: F(x1)-F(x0)可以得到该定积分的值。但是，很多情况下，由于f(x)太复杂，我们无法计算得到原函数F(x)的显示解，这时我们就只能用数值积分的办法。如下是一个简单的数值积分的例子。

数值积分简单示例

如图，数值积分的基本原理是在自变量x的区间上取多个离散的点，用单个点的值来代替该小段上函数f(x)值。

常规的数值积分方法是在分段之后，将所有的柱子（粉红色方块）的面积全部加起来，用这个面积来近似函数f(x)（蓝色曲线）与x轴围成的面积。这样做当然是不精确的，但是随着分段数量增加，误差将减小，近似面积将逐渐逼近真实的面积。

1

Monte Carlo数值积分方法和上述类似。差别在于，Monte Carlo方法中，我们不需要将所有方柱的面积相加，而只需要随机地抽取一些函数值，将他们的面积累加后计算平均值就够了。通过相关数学知识可以证明，随着抽取点增加，近似面积也将逼近真实面积。

在金融产品定价中，我们接触到的大多数求基于某个随机变量的函数的期望值。考虑一个欧式期权，假定我们已经知道在期权行权日的股票服从某种分布（理论模型中一般是正态分布），那么用期权收益在这种分布上做积分求期望即可。

（五）随机最优化

Monte Carlo在随机最优化中的应用包括：模拟退火(Simulated Annealing)、进化策略(Evolution strategy)等等。一个最简单的例子是，已知某函数，我们要求此函数的最大值，那么我们可以不断地在该函数定义域上随机取点，然后用得到的最大的点作为此函数的最大值。这个例子实质也是随机数值积分，它等价于求此函数的无穷阶范数（-Norm）在定义域上的积分。

由于在金融产品定价中，这部分内容用的相对较不常见，所以此课程就不介绍随机最优化方法了。

十二、Monte Carlo形式与一般步骤 （一）积分形式

做Monte Carlo时，求解积分的一般形式是：

X为自变量，它应该是随机的，定义域为(x0, x1)，f(x)为被积函数，ψ(x)是x的概率密度。在计算欧式期权例子中，x为期权到期日股票价格，由于我们计算期权价格的时候该期权还没有到期，所以此时x是不确定的（是一随机变量），我们按照相应的理论，假设x的概率密度为ψ(x)、最高可能股价为x1(可以是正无穷）、最低可能股价为x0（可以是0），另外，期权收益是到期日股票价格x和期权行权价格的函数，我们用f(x)来表示期权收益。

（二）一般步骤

我将Monte Carlo分为三加一个步骤：

1．依据概率分布ψ(x)不断生成随机数x, 并计算f(x)

由于随机数性质，每次生成的x的值都是不确定的，为区分起见，我们可以给生成的x赋予下标。如xi表示生成的第i个x。生成了多少个x，就可以计算出多少个f(x)的值

2．将这些f(x)的值累加，并求平均值

例如我们共生成了N个x，这个步骤用数学式子表达就是

3．到达停止条件后退出

常用的停止条件有两种，一种是设定最多生成N个x，数量达到后即退出，另一种是检测计算结果与真实结果之间的误差，当这一误差小到某个范围之内时退出。

有趣的类比：积分表达式中的积分符合类比为上式中累加符号，dx类比为1/N（数学知识告诉我们积分实质是极限意义下的累加；f(x)还是它自己，积分中的ψ(x)可类比为依据ψ(x)生成随机数 4．误差分析

2

Monte Carlo方法得到的结果是随机变量，因此，在给出点估计后，还需要给出此估计值的波动程度及区间估计。严格的误差分析首先要从证明收敛性出发，再计算理论方差，最后用样本方差来替代理论方差。在本课程中我们假定此方法收敛，同时得到的结果服从正态分布，因此可以直接用样本方差作区间估计。详细过程在例子中解释。

这个步骤的理论意义很重要，但在实际应用中，它的重要性有所淡化，倘若你的老板不太懂这些知识，你报告计算结果时可以只告诉他点估计即可。

注意，前两大步骤还可以继续细分，例如某些教科书上的五大步骤就是将此处的前两步细分成四步。

十三、最简单的例子

举个例子:

计算从函数从0到2的定积分值。

=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在后面小节详细讲）。看Matlab代码：

N=100;x=unifrnd(0,2,N,1);mean(2*exp(x))

同样的，通过增大N，这种方法得到的结果也将越来越接近真实值。

解释

这个例子要求的积分形式是：， 还不完全是形式，我们先做变

换，，这里是f(x)；1/2是ψ(x)，它表示，在取值范围(0,2)区间内，x服从均匀分布。

前一例子共三条语句，逐句解释如下：

N=100;

设定停止条件，共做N次Monte Carlo模拟。 x=unifrnd(0,2,N,1);

按照(0,2)区间均匀分布概率密度对x随机抽样，共抽取N个xi。此句相当于第一个步骤中的前半部分。

mean(2*exp(x))

2*exp(x)作用是对每个xi计算f(xi)的值,共可得到N个值，这个相当于第一个步骤后半部分；Mean()函数的作用是将所有的f(xi)加起来取平均值，相当于第二个步骤。

3

这段代码中的停止条件隐含于N值设定中，它一次性生成N个x值，完成此次计算后整个程序就结束了。

十四、Monte Carlo方法的优点

对比前面常规数值积分和Monte Carlo数值积分代码，同样数量的N值——也就意味这几乎相同的计算量——常规数值积分结果的精确度要高于Monte Carlo数值积分的结果。那么，我们为何还需要用Monte Carlo来算数值积分呢？

答案的关键在于，常规数值积分的精度直接取决于每个维度上取点数量，维度增加了，但是每个维度上要取的点却不能减少。在多重积分中，随着被积函数维度增加，需要计算的

函数值数量以指数速度递增。例如在一重积分中，只要沿着x轴取N个点；

要达到相同大小的精确度，在s重积分

s

中，仍然需要在每个维度上取N个点，s个纬度的坐标相组合，共需要计算N个坐标对应的f()函数值。取点越多，会占用计算机大量内存，也需要更长运算时间，最终导致这种计算方法不可行！

Monte Carlo方法却不同，不管是积分有多少重，取N个点计算的结果精确度都差不多。因此，即使在一重积分的情形下，Monte Carlo方法的效率比不过常规数值积分，但随着积分维度增加，常规数值积分的速度呈指数下降，Monte Carlo方法的效率却基本不变。经验表明，当积分重数达到4重积分甚至更高时，Monte Carlo方法将远远优于常规数值积分方法。

现在回到金融产品定价，欧式期权理论定价公式只需要一重积分，此时Monte Carlo方法的效果不明显，但是如果我们考虑一个亚式期权：期限为1年期，期权价格基于此1年内每天某个时点时的价格，全年共252个交易日，这样此亚式期权理论定价公式是一个252

252

重积分。常规的数值积分方法，需要取N个点，这个数有多大，你自己去计算一下就知道了（注意：N取值要远远大于2），常规数值积分方法不可行，只能用Monte Carlo。

综上，如果计算高维度多重积分，如路径依赖的exotic options（奇异期权）等金融产品定价，我们一般用的方法都是Monte Carlo。

十五、Monte Carlo方法原理(选读)

Monte Carlo方法计算的结果收敛的理论依据来自于大数定律，且结果渐进地（Asymptotically）服从正态分布的理论依据是中心极限定理。

以上两个属性都是渐进性质，要进行很多次抽样，此属性才会比较好地显示出来，如果Monte Carlo计算结果的某些高阶距存在，即使抽样数量不太多，这些渐进属性也可以很快地达到。

这些原理在理论上意义重大，但由于我们一般遇上的Monte Carlo问题都是收敛的、结果也都是渐进正态分布，所以工作中使用时可以不加考虑。

详细推导见相关书籍。

4

第二章：随机数的生成

讲课人：Xaero Chang | 课程主页: http://macro2.org/notes/intro2mc

本章第一节会简要复习随机变量的一些概念，但学习本章最好要有一定的数学基础。第二节主要介绍如何生成一维概率分布的随机数，第三节介绍如何生成高维分布的随机数。最后略提及伪随机数问题的应对策略。

由前文可知，Monte Carlo积分解决的问题形如，f(x)值只需由x值决定，因此此处最重要的就是如何生成服从ψ(x)概率分布的随机数。可以说，正确生成随机数，Monte Carlo方法就做完了一半。

一、随机变量基本概念 （一）随机变量

现实世界中有很多可以用数字来衡量的事物，站在当前时间点来看，它们在未来时刻的值是不确定的。例如，我们掷一骰子，在它停稳前，我们不可能知道掷出多少点（传说中的赌王除外，哈哈）；例如某只股票在明天的股价，没有人能准确知晓第二天股票的价格（不然他就发惨了！）。但是，我们却可以描述这些事物未来各种值的可能性。

（二）离散型随机变量

离散型随机变量最重要的是分布律，即每个取值的概率是多少。例如掷骰子，我们认为扔出任何一个点的概率都是1/6。那么掷骰子得到的点数的分布律如下表： 骰子点数 1 2 3 4 5 6 概率 1/6 1/6 1/6 1/6 1/6 1/6 （三）连续性随机变量 连续型随机变量有两个重要的概念。概率密度函数(PDF)和累积概率分布函数(CDF)，具体定义见数学书籍。

PDF函数本身不是概率，只有对x的某段区间中的PDF积分得到的数值才有概率的含义。CDF是概率的意思，点x上CDF的值表示该随机变量可能取值小于x的概率的大小。

如图是正态分布的PDF和CDF

5

生成Gamma分布随机数的语法是： gamrnd(A,B,[M,N,P,...])

11．lognrnd()

生成服从对数正态分布的随机数。其有两个参数：mu和sigma，服从这个这样的随机数取对数后就服从均值为mu，标准差为sigma的正态分布。下图是mu=-1, sigma=1/1.2的对数正态分布的PDF图形。

生成对数正态分布随机数的语法是： lognrnd(mu,sigma,[M,N,P,...])

11

12．raylrnd()

生成服从瑞利（Rayleigh）分布的随机数。其分布有1个参数：B。下图是B=2的瑞利分布的PDF图形。

生成瑞利分布随机数的语法是： raylrnd(B,[M,N,P,...])

13．wblrnd()

生成服从威布尔（Weibull）分布的随机数。其分布有2个参数：scale 参数 A和shape 参数 B。下图是A=3，B=2的Weibull分布的PDF图形。

12

生成Weibull分布随机数的语法是： wblrnd(A,B,[M,N,P,...])

还有非中心卡方分布(ncx2rnd)，非中心F分布(ncfrnd)，非中心t分布（nctrnd)， 括号中是生成服从这些分布的函数，具体用法用：

help 函数名

查找。

c. 离散型分布随机数

离散分布的随机数可能的取值是离散的，一般是整数。 14．unidrnd()

此函数生成服从离散均匀分布的随机数。Unifrnd是在某个区间内均匀选取实数（可为小数或整数），Unidrnd是均匀选取整数随机数。离散均匀分布随机数有1个参数：n, 表示从{1, 2, 3, ... N}这n个整数中以相同的概率抽样。基本语法：

unidrnd(n,[M,N,P,...])

这些随机数排列成M*N*P... 多维向量。如果只写M，则生成M*M矩阵；如果参数为[M,N]可以省略掉方括号。一些例子：

unidrnd(5,5,1) %生成5个随机数排列的列向量，一般用这种格式 unidrnd(5,5) %生成5行5列的随机数矩阵

unidrnd(5,[5,4]) %生成一个5行4列的随机数矩阵

%注：上述语句生成的随机数所服从的参数为(10,0.3)的二项分布 生成的随机数大致的分布。 x=unidrnd(9,100000,1); hist(x,9);

可见，每个整数的取值可能性基本相同。 15．binornd()

此函数生成服从二项分布的随机数。二项分布有2个参数：n,p。考虑一个打靶的例子，每枪命中率为p，共射击N枪，那么一共击中的次数就服从参数为（N,p）的二项分布。注意p要小于等于1且非负，N要为整数。基本语法：

binornd(n,p,[M,N,P,...])

生成的随机数服从参数为(N,p)的二项分布，这些随机数排列成M*N*P... 多维向量。如果只写M，则生成M*M矩阵；如果参数为[M,N]可以省略掉方括号。一些例子：

binornd(10,0.3,5,1) %生成5个随机数排列的列向量，一般用这种格式 binornd(10,0.3,5) %生成5行5列的随机数矩阵

binornd(10,0.3,[5,4]) %生成一个5行4列的随机数矩阵

%注：上述语句生成的随机数所服从的参数为(10,0.3)的二项分布 生成的随机数大致的分布。 x=binornd(10,0.45,100000,1); hist(x,11);

我们可以将此直方图解释为，假设每枪射击命中率为0.45，每论射击10次，共进行10万轮，这个图就表示这10万轮每轮命中成绩可能的一种情况。

16．geornd()

此函数生成服从几何分布的随机数。几何分布的参数只有一个：p。几何分布的现实意义可以解释为，打靶命中率为p，不断地打靶，直到第一次命中目标时没有击中次数之和。注意p是概率，所以要小于等于1且非负。基本语法：

geornd(p,[M,N,P,...])

13

这些随机数排列成M*N*P... 多维向量。如果只写M，则生成M*M矩阵；如果参数为[M,N]可以省略掉方括号。一些例子：

geornd(0.4,5,1) %生成5个随机数排列的列向量，一般用这种格式 geornd(0.4,5) %生成5行5列的随机数矩阵

geornd(0.4,[5,4]) %生成一个5行4列的随机数矩阵

%注：上述语句生成的随机数所服从的参数为(0.4)的二项分布 生成的随机数大致的分布。 x=geornd(0.4,100000,1); hist(x,50);

17．poissrnd()

此函数生成服从泊松(Poisson)分布的随机数。泊松分布的参数只有一个：lambda。此参数要大于零。基本语法：

geornd(p,[M,N,P,...])

这些随机数排列成M*N*P... 多维向量。如果只写M，则生成M*M矩阵；如果参数为[M,N]可以省略掉方括号。一些例子：

poissrnd(2,5,1) %生成5个随机数排列的列向量，一般用这种格式 poissrnd(2,5) %生成5行5列的随机数矩阵

poissrnd(2,[5,4]) %生成一个5行4列的随机数矩阵

%注：上述语句生成的随机数所服从的参数为(2)的泊松分布 生成的随机数大致的分布。 x=poissrnd(2,100000,1); hist(x,50);

其他离散分布还有超几何分布(Hyper-geometric, 函数是hygernd)等，详细见Matlab帮助文档。

（六）特殊连续分布

这里我将Matlab中没有对应函数的分布称为特殊分布。有多种方法可以用于生产服从这些分布的随机数。这里主要介绍两种最常见的。

1．逆CDF函数法

如果我们已知某特定一维分布的CDF函数，经过如下几个步骤即可生成符合该分布的随机数。（其中数学推导等在此处略去，详见相关数学书籍）

1. 计算CDF函数的反函数：

2. 生成服从(0,1)区间上均匀分布的初始随机数a

3. 令x=，则x即服从我们需要的特定分布的随机数。

为了更形象解说这种方法，这里选取柯西（Cauchy）分布作为例子。有时也称其为洛仑兹分布或者 Breit-Wigner 分布。柯西分布有一大特点就是，它是肥尾（Fat-tail，又译作胖尾）分布。在金融市场中，肥尾分布越来越受到重视，因为在传统的正态分布基本不考虑像当前次贷危机等极端情况，而肥尾分布则能很好地将很极端的情形考虑进去。

14

上图是Cauchy 分布和标准正态分布PDF图对比，看看是不是Cauchy分布的尾巴（x轴两端）更“胖”一点？

柯西分布的PDF函数是：

简化起见我们只考虑x0=0, γ=1情形。此时PDF函数是：

PDF函数对x作积分，就得到CDF函数（推导过程略）：

现在我们套用这三个步骤来生成服从Cauchy分布的随机数：

1. 计算得到Cauchy分布CDF函数的反函数为：

2. 使用rand()函数生成（0，1）区间上均匀分布的初始随机数。我习惯一次生成一堆这种随机数。

original_x=rand(1,100000);

3. 将初始随机数代入CDF反函数即可得到我们需要的Cauchy随机数。。 cauchy_x=tan((original_x-1/2)*pi);

上面这两句代码结合起来就生成了10万个服从参数为(x0=0, γ=1)Cauchy分布的随机数。

这种方法生成随机数与Cauchy分布有多大相似之处呢？这里有一个图可以说明：

15

蓝色的图形就是用hist画出的随机数的样本分布情况，红色线条是Cauchy分布理论的PDF图形。由此可看出生成的随机数挺符合Cauchy分布。

注意：上图中，我略去了x轴小于-12.5和大于12.5部分的图形——因为Cauchy是胖尾分布，会生成出的不少取值很大的随机数，而那些很大的值使得我们不方便用hist函数来画随机数分布图。

注意，这种方法本身虽然很简单，效率也很高，但有如下受限之处：

1．它有个可能会出错的地方，有的CDF函数的反函数在0或者1处的值是正/负无穷，例如此处的Cauchy分布就是这样，倘若用(0,1)均匀分布产生的初始随机数中包含0或者1，那么这个程序会出错。幸运的是，迄今为止，我用Matlab的rand()函数生成的随机数中还没有出现过0或者1。但不同版本的Matlab的这种情况也许会改变。此处提醒你，如果程序出错，不要忘记检查是否是这个错误。

2．CDF函数必须严格单调递增，这也就意味着，PDF函数在x定义域内必须处处严格大于零，否则CDF的反函数不存在。

3．即使CDF函数存在，如果它太复杂，可能导致计算速度太慢，甚至无法计算的后果。

2．接受/拒绝法 Acceptance-Rejection Method

Accelptence-Rejection方法的精髓在于“形似”，可以形象地将其比喻为制作冰雕——二者相同之处在于都要首先堆砌出雏形，然后再用将多出的部分削去。用此法生成服从f(x)分布的随机数，分为如下几大步骤：

1．首先，选用某个分布，如pdf为g(x)的分布，此时要计算一个常数c, 使得,对x定义域内任意的x都成立——这相当于使cg(x)图形完全“覆盖”住f(x)图形，容易理解，做冰雕时，最初堆出来的那堆冰块要比最终得到的雕塑大。

2．生成服从pdf为g(x)分布的随机数，假设生成的随机数为x0。 3．再生成一个服从（0，1）间的均匀分布的随机数y 4．如果的随机数。

，丢弃生成的x0；反之，生成的x0就是我们需要的、服从f(x)分布

16

下面用一个例子结合图形解释这种方法，假设我们要生成的分布是:

，此pdf图形如下图的蓝色曲线。

1．我们选用（0，2）之间的均匀分布作为原始分布，即g(x)=0.5，此分布的pdf图见

下图中的绿色线。由条件：无论哪个x，都要成立，我们计算得到c要大于等于10.8。这种情况下，我们一般选择c=1.875。因为c选得越大，意味着我们堆砌的原始雏形越大，需要削去的部分越多，效率越低，所以我们要使得c尽量地小。

2．生成服从（0，2）之间的均匀分布的随机数，设它为x0

X0=unifrnd(0,2);

3．然后再生成一个服从（0，1）间的均匀分布的随机数y Y=rand;

4．如果，丢弃生成的x0，重新生成；反之，生成的x0就是我们需要的、服从f(x)分布的随机数，用于做后续计算。

以上步骤每次只能处理一个随机数，效率较低，下面这段代码可以一次性生成一堆随机数。

N=400000;c=1.875;gx=0.5 x0=unifrnd(0,2,1,N); y=rand(1,N);

fx0=(x0-0.5).*(x0-0.5)/2.4; final_x=x0(y<=fx0./c/gx);

在视频教程中我会逐句解释每句含义，如果没听懂，一般是因为你对Matlab向量运算不熟悉，请参照Matlab基础教程学习此部分的内容，后面章节会有很多地方用得上。

这段语句生成的变量final_x即为服从f(x)分布的随机数组成的一个行向量。我们可以用hist查看这些随机数大致的分布。

hist(final_x,50);title('f(x)=(x-0.5)^2/2.4');

17

如图所示，生成的随机数挺符合f(x)分布。

这种方法很简单，也不需要计算CDF函数的反函数，但它也有如下受限之处： 1．由于我们用随机数y来控制是否削去某个随机数x0，所以我们无法准确预知最终得到的随机数数量多少。

2．选择合适的g(x)分布是此方法最关键的技巧所在。g(x)的选择原则是在完全覆盖f(x)的前提下尽可能与f(x)形似，二者形状越相似，需要削去的部分就越少，这种方法的效率就越高。需要记住的是：很多时候，人们不选用这种方法的原因几乎都在于它的效率过低。

（七）特殊离散分布

离散分布关键在获得它的分布律，有了分布律我们计算骰子掷出点数小于等于某个数字的累积概率分布。一个简单的例子，假设我们有一个不均匀的骰子，获得六个点数的概率分别是：

点数 1 2 3 4 5 6 概率 .1 累积点数 1 累积概率 .1 0.3 2 0.4 3 0.6 4 0.8 5 06 1 0.2 0.1 0.2 0.2 0.2 0生成符合该分布随机数的步骤是：

1．生成一个（0，1）间均匀分布的随机数x0。

2．依据x0介于累积概率哪个区间来决定掷出骰子的点数x。如0

代码是

x0=rand; if x0<0.1 x=1;

18

elseif x0<0.3 x=2;

elseif x0<0.4 x=3

elseif x0<0.6 x=4

elseif x0<0.8 x=5 else x=6 end

这段语句能生成一个服从上表中离散分布的随机数x，如果生成多个x，可以用循环语句，也可以考虑将上述代码向量化。下图是我用上述代码生成10万个随机数所画出的分布直方图，可见这些随机数很符合上表中的分布律。

十七、生成多维联合分布随机数

一维随机变量是标量（也就是指单独的一个数字），而多维随机变量是一个向量。一个n维随机变量x是有n个分量的向量，(X0,X1,...,Xn)，用f(X0,X1,...,Xn)表示联合分布，用fk(Xk)表示第k维的边缘分布，用fk(Xk|X1=x1,X2=x2,....Xk-1=xk-1, Xk+1=xk+1,...,Xn=xn)表示当分量X1=x1,X2=x2,....Xk-1=xk-1, Xk+1=xk+1,...,Xn=xn时第k个分量xk的分布。这里大写X表示随机变量某个维度上的分量，小写x表示具体的数值。关于边缘分布、条件分布、联合分布一定要明白，这些都是基础数学知识，非本课程内容。如果手头没有书，通过google搜索或上维基百科临阵磨枪也是可以的。

各种生成多维分布随机数的方法一般步骤都是，逐个维度生成随机数分量，最后将这些分量依次组合起来——如先生成x0,再x1,...,最后xn,，最终写成(x0,x1,...,xn)。

在详细讲如何生成这些分量前，我们讲讲如何储存生成的随机数。 如果一次生成一个n维的随机数向量，可以用n变量来储存这个随机数的n个分量，也可以将这n个分量按照次序（次序不能乱）存于一个1*n的行向量中。如果一次生成随机数的数量很多，例如N个随机数，前面两种办法都可以用，即可用n个变量来储存这些随机数的每个分量，此时每个变量是N*1的列向量；也可以只用一个N*n矩阵储存随机数所有

19

分量，这个矩阵每一行是一个服从规定的联合分布的随机数，共有N行即表示共储存N个这样的随机数，矩阵的每一列表示这N个随机数的一个维度上的分量，共有n个维度。

（一）最简单的——各维度独立

各维度独立的联合分布随机数的生成最为方便，由于联合分布函数就是每个维度边缘分布函数的直接乘积，所以只要分别生成每个维度的随机数分量，然后组合成随机数向量即可得到服从该联合分布的随机数。

例子1，生成一个在0≤x≤2,0≤y≤2，正方形区域上的二维均匀分布。二维均匀分布在每个维度上都是均匀分布（即两个维度的边缘分布都是(0,2)上的均匀分布），且两个维度互相独立。

用第一种存储方法，

x=unifrnd(0,2); y=unifrnd(0,2);

则每个维度上分别生成一个服从(0,2)均匀分布并分别储存在x,y这两个变量中。如果一次生成多个随机数，如N个，可用

N=400;

x=unifrnd(0,2,N,1); y=unifrnd(0,2,N,1);

这里x，y都是N*1大小的列向量，分布存储着这N个随机数的第一维和第二维两个分量。我们看看这些随机数是否很好得符合二维均匀分布特性。

scatter(x,y);

接着，我们看看用上述第二种存储方法， X=[x,y];

紧接着第一种存储方法中的语句，我们将生成的两个分量组合起来储存到一个变量中。当然这里还有一种取巧的办法，由于两个维度的边缘分布都相同且独立，我们只需用unifrnd函数一次性生成一个N*n大小的矩阵就可以了。

X=unifrnd(0,2,N,2);

例子2，我们要生成的随机数服从一个三维联合分布，其第一维边缘分布服从标准正太，第二维边缘分布是自由度为4的t分布，第三维边缘分布是自由度为(7,8)的F分布，各个维度边际的边缘分布之间相互独立。我们只要用如下代码：

x1=rand x2=trnd(4) x3=frnd(7,8) x=[x1,x2,x3]

x1,x2,x3分布储存三个维度的分量（第一种存储方法），将这些分量组合起来存入x中（第二种存储方法）。

如果要一次就能生成一堆这样的随机数。可以用如下的代码：

N=1000;

x=[rand(N,1),trnd(4,[N,1]),frnd(7,8,[N,1])];

这段代码略过了中间过程，直接生成第二种存储方法所说的矩阵，这个矩阵大小为N*3， 我们可以大致观察该联合分布在每个区域内的概率密度的大小。 scatter3(x(:,1),x(:,2),x(:,3),'marker','.','sizedata',1);

注：x(:,1)表示将x矩阵的第一列（也即随机数第一维上的分量）提取出来。

20

点越密集的地方，该联合分布概率密度函数的值越大。你可以在Matlab中旋转图形更细致观察该分布。特别的，分别从x,y,z三个轴的角度看次图形的横截面图可以分别看到此三维度边缘分布的大致分布图形。

（二）最常用的——协方差阵生成多元正态分布

一个n维的随机变量，其协方差矩阵为一个n*n大小的矩阵，该矩阵对角线上的元素是随机变量各个分量的方差，矩阵其他位置的元素是各维分量两两之间的协方差；此随机变量的相关系数矩阵也是一个n*n大小的矩阵，该矩阵对角线上的元素都是1，其他位置的元素是各维分量两两之间的相关系数。这两个矩阵关系非常密切。更详细的内容请参考数学书籍。

协方差矩阵及相关系数矩阵揭示了不同纬度之间的线性相关关系，它只是高维随机变量的部分性质，一般而言，仅仅依据协方差矩阵或相关系数矩阵再加上各维度的边缘分布信息，我们还不能确定此随机变量的联合分布。如下例子：

下表是一个两维的离散型随机变量。两个维度的边缘分布都是（1/4,1/2,1/4)，两维间相关系数是0. 维度1\\维度2 -1 0 1 此维边缘分布 -1 0 1 此维边缘分布 0 1/4 0 1/4 1/4 0 1/4 1/2 0 1/4 0 1/4 1/4 1/2 1/4 下表也是一个两维的离散型随机变量。两个维度的边缘分布都是（1/4,1/2,1/4)，两维间相关系数是0。但是这两个表的联合分布是不同的，即二者是两个不同的随机变量。这个例子说明边缘分布加相关系数并不能完全代表多维分布的所有信息。 维度1\\维度2 -1 0 1 此维边缘分布 -1 0 1 此维边缘分布 1/8 0 1/8 1/4 0 1/2 0 1/2 1/8 0 1/8 1/4 1/4 1/2 1/4 但是有一类特殊的分布：多元正态分布，它的全部信息可以浓缩成边缘分布加相关系数。多元正态分布的边缘分布都是正太分布，只要我们知道每个维度上的边际正太分布的均值和标准差，再加上相关系数矩阵，我们就可以得到整个联合分布。

此节所述方法的步骤：

1．依照给定的边缘分布的均值和标准差，分别独立地生成各个维度上的符合正态分布的随机数。并依次序组合成一个向量。

2．将相关系数矩阵作Cholesky分解

3．用分解得到的矩阵乘以第一步中生成的向量即可得到我们需要的随机数。

例子1：假定我们要生成一个三维的多元正态分布。各个维度均值标准差如下表：

维度 均值 标准差 1 2 3 相关系数矩阵如下表： 1 0.3 0.4 详细做法：

21

2 -1 0 0.3 1 0.2 3 2 1 0.4 0.2 1 1．生成各维度上的独立的正态分布随机数。注：此处代码一次性生成10万个三维正态分布随机数，这些数组成了一个100000*3大小的矩阵。

N=100000;

x0=[normrnd(2,3,N,1),normrnd(-1,2,N,1),normrnd(0,1,N,1)]; 2．将系数矩阵R做Cholesky分解得到矩阵L。 R=[1,0.3,0.4;0.3,1,0.2;0.4,0.2,1]; L=chol(R);

3．计算x0*L，即可得到10万个符合上述二表中条件要求的多元正态分布随机数，这些随机数被存储在一个100000*3大小矩阵中。

x=x0*L;

注意，最后一个语句是矩阵乘法，L和x0的次序不能颠倒，否则会出错。

例子2：生成本章第一节第三小节“多元分布”中所述相关系数为0.3的二维正太分布随机数。

1．这个二维分布在每个维度上的边缘分布都是标准正太。所以我们先生成由标准正太随机数组成的N*2矩阵

N=1000;

x0=randn(N,2);

2．由两个维度间相关系数为0.3，可以知道其相关系数矩阵，并作Cholesky分解 R=[1,0.3;0.3,1]; L=chol(R);

3．最后计算x0*L，即可得到10万个符合上述二表中条件要求的多元正态分布随机数，这些随机数被存储在一个100000*2大小矩阵中。

x=x0*L;

我们可以将这10万个随机数画在二维平面上。用如下语句： scatter(x(:,1),x(:,2),'marker','.','sizedata',1)

（三）最一般的——由联合分布生成多维分布随机数（选读）

我们知道，联合分布函数包含多维分布随机数所有信息，所以我们直接从联合分布函数出发，通过相关的技巧生成随机数应该可以解决任意形式分布的问题。的确，只要给出联合分布函数pdf，无论此联合分布如何诡异，我们用此节所述方法都有可能将服从此分布的随机数向量生成出来。但是这种方法涉及到很多计算，相当麻烦，故使用地很少。

这种方法的总体原则是：对联合分布PDF函数积分，计算出某个维度的边缘分布，用其生成随机数；再将已生成的这个维度的随机数代回联合分布函数，得到这个维度分量数值给定条件下的新的联合分布函数，不断重复上述过程，直到所有维度的值都确定为止。最后将各个维度的值组合起来即可得到我们所需的随机数。 假定联合分布的PDF函数为：步骤

，生成符合此分布的随机数需要如下的

1．选择一个维度，计算该维度上的边缘分布，例如可以选x1，计算得边缘分布为2．生成一个服从上述边缘分布的随机数3．令

，计算得到剩余维度的概率密度函数

4．再选择一个维度，依据计算得边缘分布为

计算该维度上的边缘分布，例如可以选x2，

。

22

5．生成一个服从上述边缘分布的随机数

6．重复进行3，4，5三个步骤，直到所有维度的随机数都生成。

7．最后将各维上的随机数组合成一个n维向量，该向量就是服从分布的随机数向量。

生成服从某边缘分布随机数的方法要用到前一节讲的生成一维分布的方法。如果是常规分布，可以直接用Matlab内置函数，如果是特殊的分布，可以用逆CDF法或者Acceptance-Rejection法。

在视频教程中将以一个相对简单的例子进一步说明此方法。此例子是一个二维的随机数向量，联合密度函数是：

PDF：，定义域：

（八）最流行的——Copula

现在Copula很热门，书籍很多，当然它也确实很有用。如果要详细讲Copula，估计要一本挺厚的书才行。这里我们只简要介绍一点Copula的基础知识。需要指出，本课程的简单例子并没有用到Copula，所以这个部分你完全可以跳过。

我们知道，多维随机变量联合分布一般都很复杂，为了将这个复杂的问题简化，我们常常考虑，各个维度上的边缘分布。现在的问题在于，联合分布与边缘分布之间的关系也很复杂。

最简单的情形是各个维度独立，则联合分布函数就是各个维度上边缘分布函数的直接乘积

稍微复杂一点是多元正态分布，我们知道各维度正态分布的均值方差（也即相当于知道了边缘分布），再加上两两维度之间相关系数的信息，我们就可以得到整个多元正态分布的联合分布函数。

但是，除此之外，我们就很难找到联合分布与边缘分布间比较简单的关系了。这种困难孕育了Copula。Copula是一种函数，这种函数揭示了联合分布与各维度边缘分布的关系。下面式子是一个二维分布的Copula函数：

左侧是联合累积概率分布(多维随机变量的CDF函数)，右侧的C表示Copula函数，C()内部分布是x维度上的边缘分布和y维度上的边缘分布。这个式子的意义在于，我们只要知道边缘分布和Copula函数，整个联合分布就可以计算出来。

Copula的名称也随着C()函数形式的不同而改变。常见的有高斯Copula，t-Copula，Archimedean Copula等等。

Copula的一个作用在于：我们从现实中能直接观察到的信息主要是边缘分布，利用这些边缘分布，然后选择合适的Copula函数类型，用这些信息去拟合现实数据，最后可以确定Copula函数中待定参数。至此，我们就能得到联合分布函数。

高版本Matlab附带的统计工具箱中有几个常用的Copula命令，使用很方便。这里我也顺带提一下如何自写代码生成n维的Gauss-Copula。

1．Gauss-Copula的参数是一个n*n相关系数矩阵，将此矩阵做Cholesky分解得到L 2．生成一个n维、各维为独立标准正太的多元正态分布随机数向量x0 3．计算x0*L，并将结果的各个分量求其标准正态分布CDF函数的值

4．用各个分量分别代换到各维的CDF函数的反函数中，计算得到各维上的数值

23

5．最后将这些数值按照顺序组合成一个向量，此向量即为我们要的随机数向量。 代码及Matlab自带Copula命令见视频教程的PPT。

（九）最熟悉的——Acceptance-Rejection Method

方法原理和前面一维的情况下完全一样。只要将那一小节中的所有步骤中的x改为x向量，就变成了多维情形下的Acceptance-Rejection Method。此方法请自行研究。

多维情形下这种方法效率过低的弱点更加明显，所以选择一个好的g(x)尤为重要。但问题在于，维数越高，联合分布越复杂，选择g(x)就越困难。

十八、伪随机数的诅咒

很多书介绍了计算机生成的随机数都不是真正意义上的随机数这个问题，而且还详细讲明伪随机数会带来很多可能的危害。这个问题的确不可忽视，但我们也不要杞人忧天。以我们现在的技术水平，要想生成比Matlab更“真”的随机数基本不可能，所以，不要花太多时间在此方面。我建议注意以下两点：

1．由于同一台机器生成的随机数的模式固定，且有一定的规律，所以可以在一台机器上长时间运行同一个程序，以此来测试此程序的稳健性。因为程序中可能有些很微小的错误需要运行很长一段时间、落入某个特殊“陷阱”时才能显现。

2．不同机器随机数生成的机制会有差别，所以可以将同一套程序在多台配置不同的机器上运行，如果不同机器上得到的结果差不多，说明此程序质量很好，不受特定的随机数生成机制影响；否则，如果结果悬殊，则证明你的程序很可能有问题。

第三章：随机过程模拟

讲课人：Xaero Chang | 课程主页: http://macro2.org/notes/intro2mc

我们不可能确切知道未来经济变量——如，股票价格、收益率、债券利率、汇率等等——确切走势，一般我们采用与布朗运动有关的随机过程来描述这些变量的未来趋势。这里我们遇上的问题就是如何用计算机代码来描述这些随机的运动。

注意：这里只讲基于布朗运动的一维随机过程，它是学习其他深层次知识的基础。 基于布朗运动的随机过程是连续的，我们如何使用它们呢？依据计算任务的不同，一般分为如下这些情况：如果金融产品价格只基于到期日时标的资产价格，那么我们只需要由随机过程推导出标的资产在到期日价格的分布，然后用Monte Carlo方法，例如欧式期权就是此类；如果金融产品价格基于在到期日前标的资产在离散时间点上的市场价格，我们就用随机过程推导在每个盯市时间点上标的资产的价格的分布，然后用Monte Carlo方法计算，例如离散时间盯市的亚式期权；如果金融产品价格基于在到期日前标的资产在连续时间内的市场价格，这种情况比较麻烦，详见下一段。

由于我们不可能将时间无限细分，所以我们无法用计算机直接模拟连续情况下的随机过程，一般的做法都是将连续时间近似看做离散时间的版本。例如要用Monte Carlo方法计算连续时间盯市的亚式期权，我们只能将其转化为离散时间盯市的亚式期权。这一步转化必然要带来误差，这一误差不可避免，但随着离散化时将盯市时间间隔区段越短，所得结果越精确，当然计算时间也越长。具体将区段划分到何种程度，以及转化的误差有多大，可以用复杂的数学方法进行分析，这里推荐的一个简单的办法是，用不同长度的时间间隔区段分别做

24

几次，将所得结果做一比较，即可大略知道区段划分会带来多大的误差影响，只要这个误差在可以接受范围之内，则一般可以将区段划分到该程度即可。

下面介绍三个简单的一维随机过程的模拟，更复杂的内容不太适宜在这个入门级别的课程中详细讲述，如果你想深入学习，可以参考相关书籍。

讲义和PPT的内容部分重叠，但表述方法不同，侧重点各异，故先简要讲解讲义，再详细讲解PPT内容。

一、标准布朗运动

标准布朗运动，又称维纳过程(Weiner Process)，用W(t)表示。它的性质比较独特： 1．W(0)=0， 即我们定义初始时刻的点为0点。

2．W(t)~N(0,t)，在t时刻的位置服从均值为0，方差为t的正态分布

3．W(s)-W(t)~N(0,s-t)，从时刻t走到时刻s（s要大于t），位置的变化服从均值为0，方差为s-t的正态分布，且该分布与W(t)独立，但不与W(s)独立，因为W(s)=W(t)+(W(s)-W(t))。

有了这些性质，我们可以模拟标准布朗运动（准确说是一维的标准布朗运动）。分两种情况：

如果我们只要看终点时刻T时的位置，则有性质2，我们可知，其在时刻T服从N(0,T)分布，直接生成服从该分布的随机数即可。

如果我们需要模拟此运动从时间0到时点T整个过程中的运动路径。首先需要明确一点，我们只能模拟离散时间点的位置。此处我们用t1,t2,...tn表示即将用于模拟的时间点。由上面的性质2和3，我们有:

??W(t1)~N(0,t1)

??W(t2)-W(t1)~N(0,t2-t1) ??W(t3)-W(t2)~N(0,t3-t2) ??.....

??W(T)-W(tn)~N(0,T-tn)

所以我们要做的是分别生成服从N(0,t1)，N(0,t2-t1)，N(t3-t2)，??，N(0,T-tn)这些随机数，在第i个时间点ti时此运动的位置X(ti)就是将服从N(0,ti-t(i-1))以及在这个分布之前的所有分布的随机数累加起来。

详细例子见视频教程中展示的PPT

十九、带漂移的布朗运动

标准布朗运动的表述是：

, 其微分表述是：

一般形式的布朗运动的表述是：

。这里μ是漂移系数，σ表示波动率。我们通过如下转化：

就得到了标准布朗运动。

这个转化告诉我们，如果要生成一般形式的布朗运动，我们只要先生成一个标准布朗运动，然后由

具体例子见ppt

即可得到所需的一般形式。

二十、几何布朗运动

几何布朗运动的作用是用来模拟股价的变动。它的好处在于，一般形式布朗运动中取值可能为负数，而几何布朗运动取值永远不小于0，这一点符合股价永远不为负的特征。

几何布朗运动微分形式的表述。或者称SDE（随机微分方程）形式：

25

本文来源：https://www.bwwdw.com/article/5lm2.html