金融工程与并行计算：第二章 仿真法在财务工程的使用 Part 2
第三节 随机数生成物件
在程序的设计中，当我们需要随机数时，通常我们可以由链接库中取得内建随机数。这类随机数可称之为假随机数(Pseudorandom Number)。这是因为我们是采用确定的方法(Deterministic Method)去模仿产生随机数，本质上他们并不是真正的随机数，只是刻意地使其看起来像是随机数而已。当然，他们具备了一些随机数该有的性质。
另外一类在程序设计中可充当随机数来源的是所谓的准随机数(Quasi-random Number)，这类随机数在本质上就不是以模仿随机数的方式去产生，他们是以数论的理论去产生均匀分布的数列(Low-discrepancy Sequence)，因此又称之为准随机数列(Quasi-random Sequence)。
事实上，如果要产生”真正”的随机数，那还是要求助于大自然，例如，原子核衰变过程中，放射出来的粒子的方位，便是一个天然的随机数源。网络上有网站提供由其产生的数列，这是一个有趣的网站。
一、随机数生成逻辑
计算机使用的假随机数通常是以迭代的方式去产生，我们可以表示如下，
，
………………………………………………………...(2.3.1)
此假随机数都是以均等分配随机数的方式产生的。由于计算机的离散性质，上式中的数值xi原则上都是整数，而起始值x0由外生给定。在上一节中，我们便以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, 2nd 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位的整数数列，它具有下面的优点。首先，它具有非常长的周期，219937-1，足以满足一般的使用。这也是我们通常称其为MT19937的原因。一般早期的随机数生成器周期多为232。其次，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，
其中，p1, p2,…, pd为前d个质数。
(二)Halton范例
以2,3,5,7为质数之四维HaltonSequence，第37个点由下产生，
Base n=3710
2 1001012 0.1010013 =0.640625
3 11013 0.10113 =0.382716
5 1225 0.2215 =0.488000
7 527 0.257 =0.387755
(三)Sobol数列
首先，针对每一维度需要一组方向数(direction number) vi与一个质数多项式。令vi为一二元分数，可表示为
选择正确的起始vi是非常重要的。对mi的正式要求只有mi为奇数，且
在Sobol产生器中MAXBIT= 30。
每一维度需要一个质数多项式如下，
对每一多项式，需要q个起始方向数vi。一但起始方向数vi选好后，vi可由下是递归产生。
其中，⊕为位XOR运算，上式亦可表示为
。
在求得所有的vi后，可由下式Sobol数列，
c为n以二元表示下，最右方0位的位置。
(四)Sobol范例
以下面质数多项式产生第四维之Sobol数列，
。前三个方向数如下表
i | 1 | 2 | 3 |
mi | 1 | 1 | 5 |
vi(Binary) | 0.1 | 0.01 | 0.101 |
可得
可得
i | 4 | 5 | 6 |
mi | 3 | 15 | 51 |
vi(Binary) | 0.0011 | 0.01111 | 0.110011 |
可得
x0 = 0, n = 02, c = 1
x1 = x0⊕v1 = 0.0⊕0.1 = 0.12 =0.5,
n= 012, c = 2
x2 = x1⊕v2 = 0.1⊕0.01 = 0.112=0.75,
n= 102, c = 1
x3 = x2⊕v1 = 0.11⊕0.1 = 0.012=0.25,
n= 112, c = 3
x4 = x3⊕v3 = 0.01⊕0.011 = 0.1112=0.875,
n= 1002, c = 1
x5 = 0.375
x6 = 0.125
x7 = 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数列的平均数与变异数更加接近均等分布的理论值。