SAS和蒙特卡罗模拟（3 ） – 模型编程

关键词：蒙特卡罗模拟、蒙特卡罗模型、蒙特卡罗模拟法、蒙特卡罗模拟 matlab

一、SAS随机数函数和CALL子程序
SAS系统产生随机数，两种方式，利用SAS函数(Functions)或者CALL子程序(CALL Routines)，它们的语法格式是（具体的区别容后讨论）：

方式	代码	说明
函数	var=name(seed,)	var为存储随机数列的变量，name为特定的分布函数形式，seed为随机数种子，为特定分布要求的参数（可选）
CALL子程序	call name(seed,,var)	同上,记得seed=0, ±1,±2, , ± (2**31-2)

SAS可用的随机数函数列表如下（可以参见SAS Help and Documentation-SAS Products-Base SAS-SAS Language Dictionary-Functions and CALL Routines-Functions and CALL Routines by Category）：

分布	SAS函数	说明
二项分布(Binomial)	ranBin(seed,n,p)	n为独立实验的次数，p为成功概率
指数分布(Exponential)	ranExp(seed)
正态分布(Normal)	ranNor(seed),or normal(seed)	ranNor和normal是同质的，但normal没有相对应的CALL子程序
泊松分布(Poisson)	ranPoi(seed,m)	m为均值
均匀分布(Uniform)	ranUni(seed),or uniform(seed)	ranUni和uniform是同质的，但uniform没有相对应的CALL子程序
柯西分布(Cauchy)	ranCau(seed)
伽玛分布(Gamma)	ranGam(seed,a)	a>0为形状参数
由分布律表格决定的离散分布(tabled probability distribution)	ranTbl(seed,p1,p2,…pn)	∑p=1
三角分布(Triangular)	ranTri(seed,h)	h为斜边（最可能值）

上面的随机函数，除了normal和uniform，都可以由直接相应的CALL子程序调用。
二、SAS随机数函数和CALL子程序：比较
用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。这话费解，一个例子先，创建两个随机数变量，各包含3个记录，其中x1的种子为123，x2的种子为456：

data ranuni(drop=i);
retain seed1 123 seed2 456;
do i=1 to 3;
x1=ranuni(seed1);
x2=ranuni(seed2);
output;
end;
run;
proc print data=ranuni;run;

结果为：

Obs    seed1    seed2       x1         x2
1      123      456     0.75040    0.32091
2      123      456     0.17839    0.90603
3      123      456     0.35712    0.22111

好像没什么异样。我们把上面的x1增加为6个记录：

data ranuni2(drop=i);
retain seed1 123;
do i=1 to 6;
x1=ranuni(seed1);
output;
end;
run;
proc print data=ranuni2;run;

结果如下，把上下用红色标出的数字对照看一看：

Obs    seed1       x1
1      123     0.75040
2      123     0.32091
3      123     0.17839
4      123     0.90603
5      123     0.35712
6      123    0.22111

什么意思？在第一段代码中，x2的种子根本不起作用，把x2的记录安插到x1里，看起来就是用种子123产生的随机数列加长了而已。x2的第一个值并不是由种子456产生的，而是产生第一个x1后的下一个值，x1、x2其实属于同一个随机数列，尽管x2的种子被指定为456，而x1的被指定为123。现在就可以重复上面的一句话：用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。
用CALL子程序就能够同时产生独立的随机数列。

data ranuni3(drop=i);
retain seed3 123 seed4 456;
do i=1 to 3;
call ranuni(seed3,x3);
call ranuni(seed4,x4);
output;
end;
run;
proc print data=ranuni3;run;

结果如下：

Obs       seed3        seed4         x3         x4
1     1611463328    736440516   0.75040    0.34293
2      689153326    774069794    0.32091    0.36045
3      383088854    686944750    0.17839    0.31988

以上x3就是x1。x1和x3的初始种子都是123，但x3那个结果显示的种子是当前种子值。要在SAS随机数函数语句中显示随机种子的当前值，可以由以下代码给出：

data ranuni4(drop=i);
retain seed1 123;
do i=1 to 3;
x1=ranuni(seed1);
seed=x1*(2**31-1);
output;
end;
run;
proc print data=ranuni4;
var seed1 seed x1;
run;

结果如下，可以跟上面由CALL子程序得出的结果对照：

Obs    seed1       seed          x1
1      123     1611463328   0.75040
2      123      689153326    0.32091
3      123      383088854    0.17839

———参考资料———

• Xitao Fan, etc..Monte Carlo Studies: A Guide for Quantitative Researchers. SAS Institute Inc.,2002
• 朱世武《SAS编程技术与金融数据处理》，北京：清华大学出版社，2003

转载请注明：数据分析 » SAS和蒙特卡罗模拟3 -模型编程_蒙特卡罗模拟