金融工程与并行计算：第二章仿真法在财务工程的使用 Part 2

第三节随机数生成物件

在程序的设计中，当我们需要随机数时，通常我们可以由链接库中取得内建随机数。这类随机数可称之为假随机数(Pseudorandom Number)。这是因为我们是采用确定的方法(Deterministic Method)去模仿产生随机数，本质上他们并不是真正的随机数，只是刻意地使其看起来像是随机数而已。当然，他们具备了一些随机数该有的性质。

另外一类在程序设计中可充当随机数来源的是所谓的准随机数(Quasi-random Number)，这类随机数在本质上就不是以模仿随机数的方式去产生，他们是以数论的理论去产生均匀分布的数列(Low-discrepancy Sequence)，因此又称之为准随机数列(Quasi-random Sequence)。

事实上，如果要产生”真正”的随机数，那还是要求助于大自然，例如，原子核衰变过程中，放射出来的粒子的方位，便是一个天然的随机数源。网络上有网站提供由其产生的数列，这是一个有趣的网站。

一、随机数生成逻辑

计算机使用的假随机数通常是以迭代的方式去产生，我们可以表示如下，

，………………………………………………………...(2.3.1)

此假随机数都是以均等分配随机数的方式产生的。由于计算机的离散性质，上式中的数值x_i原则上都是整数，而起始值x₀由外生给定。在上一节中，我们便以1234为起始种子，初始化随机数对象。

以C#的随机数对象为例，当我们呼叫Next()方法，产生的随机数，是介于0与maxValue之间的整数。maxValue是32位整数的最大值，2,147,483,647，以十六进制表示为7FFFFFFF。将此产生的随机数除以maxValue便可产生范围介于[0,1]的均等随机数。成为我们产生其他随机数的来源。

实务上，线性同余法是最为常见的随机数生成器，其迭代的公式如下，

……………………………………………………...(2.3.2)

其中，a, c, m都是正常数，mod表取余数运算。在著名的IBM 360系统中，常数选择如下，

读者可以根据(2.3.2)式自行撰写自己的产生器。

二、Mersenne Twister随机数生成器

一个良好的随机数生成器应该具备下面三项性质，首先，产出的随机数要能均匀分布，其次，随机数数列的周期要够长，最后，产生的速度要够快。下面简单说明，

A. 均匀分布

理论上的均等分布随机数要能均匀的分布于范围之内，对于一维的随机数，这应该不是太大的问题。然而，如果我们要的是100维的随机数，则算法的质量就可能有差异了。一些算法在高维度时，会出现群聚(Clustering)的现象，造成高维空间分布不均匀的现象。通过一些统计检定，我们可以判断随机数均匀分布的性质。建议有兴趣的读者，可以参考Knuth Donald教授的巨著，The Art Of Computer Programming, Vol. 2:Semi-numerical Algorithms, Addison-Wesley, Reading, MA, 2^nd Ed, 1981。

B. 长周期

由(1.3)的跌代式可知，随机数的产生是一个有序的过程，最终将会回到源头，自我重复。然而，重点是多久后会自我重复，也就是随机数的周期有多长。在现代的商品模拟计算中，抽取上亿个随机数是相当常见的。如果随机数周期不能达到十亿以上的数量级，则仿真的质量自也堪虑。

C. 速度快

由于抽取上亿个随机数是相当常见的，如果算法的计算效率不高，则实用价值也就不大。通常为求加速，可以考虑尽可能使用位运算的方式，来完成算法。有时使用汇编语言来撰写程序代码，以求性能的要求。

实务上，一个在过去20年颇受欢迎的算法，已逐渐为各个主要数值链接库所采用。那是由日本学者Makoto Matsumoto (松本眞)与Takuji Nishimura(西村拓士)于1997年发表的Mersenne Twister (MT)随机数生成器。MT已实作于SPSS与SAS软件中，并被认为是较可信赖的产生器。Boost C++Library，GNU Scientific Library，NAG Numerical Library，Matlab，GAUSS，Python与其他软件都已实作MT。

一般使用的MT随机数生成器产生32位的整数数列，它具有下面的优点。首先，它具有非常长的周期，2¹⁹⁹³⁷-1，足以满足一般的使用。这也是我们通常称其为MT19937的原因。一般早期的随机数生成器周期多为2³²。其次，MT可以达到623维的均匀随机数分布。第三，MT通过许多随机性的统计测试。

然而，MT也有一些缺点的，首先，算法中有很大的状态空间，其状态变量位数很长，造成CPU快取负荷很大。其次，MT的计算速率不高，除非使用平行运算版本SFMT。第三，MT并没有通过最严格的随机性统计测试，TestU01。

MT已有多个语言的实作版本，下面这个C#版本是由Mitil Ooyama所改写，网址为http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937ar.cs。程序代码如下。

#001 public class MersenneTwister

#002 {

#003 #region "PrivateParameter"

#004 /* Periodparameters */

#005 private const Int16 N= 624;

#006 private const Int16 M= 397;

#007 /* constant vector a */

#008 private const UInt32 MATRIX_A= (UInt32)0x9908b0df;

#009 /* most significant w-r bits */

#010 private const UInt32 UPPER_MASK = (UInt32)0x80000000;

#011 /* least significant r bits */

#012 private const UInt32 LOWER_MASK = (UInt32)0x7fffffff;

#013 private UInt32[] mt; /* thearray for the state vector */

#014 private UInt16 mti; /* mti==N+1means mt[N] is not initialized */

#015 private UInt32[] mag01;

#016 #endregion

#017

#018 #region"Constructor"

#019 public MersenneTwister(UInt32 s)

#020 {

#021 MT();

#022 init_genrand(s);

#023 }

#024

#025 // coded by Mitil.2006/01/04

#026 public MersenneTwister()

#027 {

#028 MT();

#029 // autogenerate seed for .NET

#030 UInt32[] seed_key = new UInt32[6];

#031 Byte[] rnseed = new Byte[8];

#032 seed_key[0]=(UInt32)System.DateTime.Now.Millisecond;

#033 seed_key[1]=(UInt32)System.DateTime.Now.Second;

#034 seed_key[2]=(UInt32)System.DateTime.Now.DayOfYear;

#035 seed_key[3]=(UInt32)System.DateTime.Now.Year;

#036 System.Security.Cryptography.RandomNumberGenerator rn

#037 = new System.Security.Cryptography.RNGCryptoServiceProvider();

#038 rn.GetNonZeroBytes(rnseed);

#039 seed_key[4]=((UInt32)rnseed[0]<<24)|((UInt32)rnseed[1]<<16)

#040 |((UInt32)rnseed[2]<<8)|((UInt32)rnseed[3]);

#041 seed_key[5]=((UInt32)rnseed[4]<<24)|((UInt32)rnseed[5]<<16)

#042 |((UInt32)rnseed[6]<<8)|((UInt32)rnseed[7]);

#043 init_by_array(seed_key);

#044 rn=null;

#045 seed_key=null;

#046 rnseed=null;

#047 }

#048

#049 public MersenneTwister(UInt32[] init_key)

#050 {

#051 MT();

#052 init_by_array(init_key);

#053 }

#054

#055 private void MT()

#056 {

#057 mt = new UInt32[N];

#058 mag01 = new UInt32[] {0,MATRIX_A};

#059 /* mag01[x]= x * MATRIX_A for x=0,1 */

#060 mti = N+1;

#061 }

#062 #endregion

#063

#064 #region "Destructor"

#065 ~MersenneTwister()

#066 {

#067 mt = null;

#068 mag01 = null;

#069 }

#070 #endregion

#071

#072 #region "seed init"

#073 /* initializesmt[N] with a seed */

#074 private void init_genrand(UInt32s)

#075 {

#076 mt[0]= s;

#077 for (mti=1; mti<N;mti++)

#078 {

#079 mt[mti] = ((UInt32)1812433253*(mt[mti-1]^(mt[mti-1]>>30))+mti);

#080 }

#081 }

#082

#083 /* initialize byan array with array-length */

#084 /* init_key is thearray for initializing keys */

#085 /* key_length isits length */

#086 /* slight changefor C++, 2004/2/26 */

#087 private void init_by_array(UInt32[]init_key)

#088 {

#089 UInt32 i, j;

#090 Int32 k;

#091 Int32key_length=init_key.Length;

#092

#093 init_genrand(19650218);

#094 i=1; j=0;

#095 k = (N>key_length ?N : key_length);

#096

#097 for (; k>0; k--)

#098 {

#099 mt[i] =(mt[i]^((mt[i-1]^(mt[i-1]>>30))*(UInt32)1664525))

#100 + init_key[j] + (UInt32)j; /* non linear*/

#101 i++; j++;

#102 if (i>=N) { mt[0] =mt[N-1]; i=1; }

#103 if (j>=key_length)j=0;

#104 }

#105 for (k=N-1;k>0; k--)

#106 {

#107 mt[i] =(mt[i]^((mt[i-1]^(mt[i-1]>>30))*(UInt32)1566083941))

#108 - (UInt32)i; /* non linear*/

#109 i++;

#110 if (i>=N) { mt[0] =mt[N-1]; i=1; }

#111 }

#112

#113 /* MSB is1; assuring non-zero initial array */

#114 mt[0] = 0x80000000;

#115 }

#116 #endregion

#117

#118 #region "Get Unsigned Int32bit number"

#119 /* generates arandom number on [0,0xffffffff]-Interval */

#120 public UInt32 genrand_Int32()

#121 {

#122 UInt32 y;

#123 if (mti >= N)

#124 { /* generate N words at one time */

#125 Int16 kk;

#126 if (mti == N+1) /* if init_genrand() has not been called, */

#127 init_genrand(5489); /* a default initial seed is used */

#128 for (kk=0;kk<N-M;kk++)

#129 {

#130 y =((mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK)) >> 1 ;

#131 mt[kk] = mt[kk+M] ^mag01[mt[kk+1] & 1] ^ y;

#132 }

#133 for (;kk<N-1;kk++)

#134 {

#135 y =((mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK)) >> 1 ;

#136 mt[kk] = mt[kk+(M-N)] ^mag01[mt[kk+1] & 1] ^ y;

#137 }

#138 y = ((mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK))>> 1 ;

#139 mt[N-1] = mt[M-1] ^mag01[mt[0] & 1] ^ y;

#140

#141 mti = 0;

#142 }

#143 y = mt[mti++];

#144

#145 /*Tempering */

#146 y ^= (y >> 11);

#147 y ^= (y << 7)& 0x9d2c5680;

#148 y ^= (y << 15)& 0xefc60000;

#149 y ^= (y >> 18);

#150 return y;

#151 }

#152 #endregion

#153

#154 #region "Get Int31number"

#155 /* generates arandom number on [0,0x7fffffff]-Interval */

#156 public UInt32 genrand_Int31()

#157 {

#158 return (genrand_Int32()>>1);

#159 }

#160 #endregion

#161

#162 #region "Get type'double'number"

#163 /* generates arandom number on [0,1]-real-Interval */

#164 public double genrand_real1()

#165 {

#166 return genrand_Int32()*((double)1.0/4294967295.0);

#167 /* dividedby 2^32-1 */

#168 }

#169

#170 /* generates arandom number on [0,1)-real-Interval */

#171 public double genrand_real2()

#172 {

#173 return genrand_Int32()*((double)1.0/4294967296.0);

#174 /* dividedby 2^32 */

#175 }

#176

#177 /* generates arandom number on (0,1)-real-Interval */

#178 public double genrand_real3()

#179 {

#180 return (((double)genrand_Int32())+0.5)*((double)1.0/4294967296.0);

#181 /* dividedby 2^32 */

#182 }

#183

#184 /* generates arandom number on [0,1) with 53-bit resolution*/

#185 public double genrand_res53()

#186 {

#187 UInt32a=genrand_Int32()>>5, b=genrand_Int32()>>6;

#188 return((double)a*67108864.0+b)*((double)1.0/9007199254740992.0);

#189 }

#190 /* These realversions are due to Isaku Wada, 2002/01/09 added */

#191 #endregion

#192 }

程序行表2.2

我们可以将程序行表2.1中的#004行，更改调用MT的对象如下，

#004 MersenneTwister Rnd = new MersenneTwister(1234);

#009行呼叫产生随机数的方法改为MT的方法，

#009 X[i] = Rnd.genrand_real1();

其他不变，编译执行，输出如下。

三、Quasi-Random Number

准随机数又称之为低差异数列(Low-discrepancy Sequence, LDS)，一数列之所以称为低差异是指，此数列中的点落于特定区域B的比率，与该区域B的空间测度成正比。空间测度可以有不同的定义，但在一维空间中，常用测度为距离；在二维空间中，常用测度为面积；在三维空间中，常用测度为体积。因此，此一性质也可见于均匀分布数列之中。产生LDS的方法有很多种，下面举两种方式为例说明之。

(一)Halton数列

对任意k≧1，可以用唯一的方式如下表示之，

其中，ai为整数，。r满足下面条件，。

定义radicalinverse function in base如下，

上式为一个介于[0,1)的有理数，藉由将k反射，以m为基底之分数而产生之。

我们可以如下定义一个d维之HaltonPoints，

其中，p₁, p₂,…, p_d为前d个质数。

(二)Halton范例

以2,3,5,7为质数之四维HaltonSequence，第37个点由下产生，

Base n=37₁₀

2 100101₂ 0.101001₃ =0.640625

3 1101₃ 0.1011₃ =0.382716

5 122₅ 0.221₅ =0.488000

7 52₇ 0.25₇ =0.387755

(三)Sobol数列

首先，针对每一维度需要一组方向数(direction number) vi与一个质数多项式。令v_i为一二元分数，可表示为

选择正确的起始vi是非常重要的。对mi的正式要求只有mi为奇数，且

在Sobol产生器中MAXBIT= 30。

每一维度需要一个质数多项式如下，

对每一多项式，需要q个起始方向数vi。一但起始方向数vi选好后，vi可由下是递归产生。

其中，⊕为位XOR运算，上式亦可表示为

。

在求得所有的vi后，可由下式Sobol数列，

c为n以二元表示下，最右方0位的位置。

(四)Sobol范例

以下面质数多项式产生第四维之Sobol数列，

。前三个方向数如下表

i	1	2	3
m_i	1	1	5
v_i(Binary)	0.1	0.01	0.101

可得

i	4	5	6
m_i	3	15	51
v_i(Binary)	0.0011	0.01111	0.110011

可得

x⁰ = 0, n = 0₂, c = 1

x¹ = x⁰⊕v₁ = 0.0⊕0.1 = 0.1₂ =0.5,

n= 01₂, c = 2

x² = x¹⊕v₂ = 0.1⊕0.01 = 0.11₂=0.75,

n= 10₂, c = 1

x³ = x²⊕v₁ = 0.11⊕0.1 = 0.01₂=0.25,

n= 11₂, c = 3

x⁴ = x³⊕v₃ = 0.01⊕0.011 = 0.111₂=0.875,

n= 100₂, c = 1

x⁵ = 0.375

x⁶ = 0.125

x⁷ = 0.625

实务上，LDS数列的主要缺点是，在高维度时会发生群聚的现象，使的随机数分布产生不均匀的结果。相对Halton数列，Sobol数列的品质较佳。因此，我们以一个Sobol数列的实作，做为应用范例。

#001 public class Sobseq

#002 {

#003 public Sobseq()

#004 {

#005 ix = new ulong[MAXDIM];

#006 iu = new ulong[MAXBIT,MAXDIM];

#007 iv = new ulong[MAXBIT* MAXDIM];

#008 iv[0] = 1; iv[1] =1; iv[2] = 1; iv[3] = 1;

#009 iv[4] = 1; iv[5] =1; iv[6] = 3; iv[7] = 1; iv[8] = 3;

#010 iv[9] = 3; iv[10] = 1;iv[11] = 1; iv[12] = 5; iv[13] = 7;

#011 iv[14] = 7; iv[15] = 3;iv[16] = 3; iv[17] = 5; iv[18] = 15;

#012 iv[19] = 11; iv[20] = 5; iv[21] = 15; iv[22] = 13; iv[23] = 9;

#013 }

#014 private int MAXBIT =30;

#015 private int MAXDIM =6;

#016 private double fac;

#017 private ulong in0;

#018 private ulong[] ix;

#019 private ulong[,] iu;

#020 private ulong[] mdeg= { 1, 2, 3, 3, 4, 4 };

#021 private ulong[] ip ={ 0, 1, 1, 2, 1, 4 };

#022 private ulong[] iv;

#023

#024 public void sobseq(int n, double[] x)

#025 {

#026 int j, k, l;

#027 ulong i, im, ipp;

#028

#029 if (n < 0)

#030 {

#031 for (k = 0; k < MAXDIM; k++) ix[k] =0;

#032 in0 = 0;

#033 if (iv[0] != 1) return;

#034 fac = 1.0 / (1L << MAXBIT);

#035 for (j = 0; j < MAXBIT; j++)

#036 for (k = 0; k < MAXDIM; k++)

#037 iu[j, k] = iv[MAXDIM * j + k];

#038 for (k = 0; k < MAXDIM; k++)

#039 {

#040 for (j = 0; j < (int)mdeg[k];j++)

#041 iu[j, k] <<= (MAXBIT - j - 1);

#042 for (j = (int)mdeg[k];j < MAXBIT; j++)

#043 {

#044 ipp = ip[k];

#045 i = iu[j - (int)mdeg[k], k];

#046 i ^= (i >> (int)mdeg[k]);

#047 for (l = (int)mdeg[k]- 1; l >= 1; l--)

#048 {

#049 if ((ipp & 1) != 0) i ^= iu[j - l,k];

#050 ipp >>= 1;

#051 }

#052 iu[j, k] = i;

#053 }

#054 }

#055 }

#056 else

#057 {

#058 im = in0++;

#059 for (j = 0; j < MAXBIT; j++)

#060 {

#061 if ((im & 1) == 0) break;

#062 im >>= 1;

#063 }

#064 if (j > MAXBIT - 1)

#065 {

#066 try

#067 {

#068 throw newException();

#069 }

#070 catch (Exception)

#071 {

#072 throw newApplicationException(

#073 "MAXBIT too small in sobseq invalidmethod");

#074 }

#075 }

#076 im = (ulong)(j * MAXDIM);

#077 for (j = 0; j < MAXBIT; j++)

#078 {

#079 for (k = 0; k < MAXDIM; k++)

#080 iv[MAXDIM * j + k] = iu[j, k];

#081 }

#082 for (k = 0; k < Math.Min(n, MAXDIM); k++)

#083 {

#084 ix[k] ^= iv[(int)im + k];

#085 x[k] = ix[k] * fac;

#086 }

#087 }

#088 }

#089 }

程序行表2.3

下面是使用Sobol数列对象的测试范例，#008产生对象，#010为此对象的初始化。初始化方法的第一个参数若小于0，则程序会自动产生MAXBIT位的方向向量。因此，此处我们输入-1便可。由于我们产生的数列是二维的，#013呼叫的产生方法，分别传入维数2以及储存的2维数组D2。我们以#011~#015循环产生十万个随机数，并取出第一维的随机数，存入X数组之中，并计算平均数与变异数。

#001 class Program

#002 {

#003 static void Main(string[] args)

#004 {

#005 const intOBS = 100000;

#006 double[] X = newdouble[OBS];

#007 double[] D2 = newdouble[2];

#008 Sobseq obj = newSobseq();

#009

#010 obj.sobseq(-1, D2);

#011 for (inti = 0; i < OBS; i++)

#012 {

#013 obj.sobseq(2, D2);

#014 X[i] = D2[0];

#015 }

#016

#017 double Sum = 0.0;

#018 for (inti = 0; i < OBS; i++)

#019 {

#020 Sum = Sum + X[i];

#021 }

#022 double Mean = Sum / OBS;

#023

#024 Sum = 0.0;

#025 for (inti = 0; i < OBS; i++)

#026 {

#027 Sum = Sum + (X[i] - Mean) * (X[i] - Mean);

#028 }

#029 double Var = Sum / (OBS - 1);

#030

#031 Console.WriteLine("Mean of Uniform(0, 1) = " +Mean.ToString());

#032 Console.WriteLine("Variance of Uniform(0, 1) = " +Var.ToString());

#033 Console.ReadKey();

#034 }

#035 }

程序行表2.4

输出结果可以发现，相对于之前的Pseudo随机数，Sobol数列的平均数与变异数更加接近均等分布的理论值。