第五节 R语言的仿真实作

下面为以R撰写的Black-Scholes程序下，阳春型选择权的仿真程序。

NPath = 1024;

MStep = 365;

dt = 1.0/365.0;

S0 = 100.0;     # Spot price

K = 100.0;      # Strike price

TTM = 1.0;      # Maturity

sigma = 0.3;    # Heston parameter : volatility of variance

r = 0.08;       # Risk free rate

q = 0.04;       # Dividend yield

drift = ((r-q)-(sigma*sigma)/2.0)*dt;

volterm = sigma*sqrt(dt);

VST = numeric(NPath);

VCT = numeric(NPath);

VPT = numeric(NPath);

set.seed(1234);

for(i in 1:NPath)

{

    St= S0;

   for(j in 1:MStep)

    {

       n = rnorm(1);

       St = St*exp(drift+volterm*n);

    }

   VST = St;

   VCT = max(St-K, 0.0);

   VPT = max(K-St, 0.0);

}

CallValue = mean(VCT)*exp(-r*TTM);

PutValue = mean(VPT)*exp(-r*TTM);

print(CallValue);

print(PutValue);

程序4.5.1(BSSimulation.R)

输出结果为，

> print(CallValue);

[1] 13.77969

> print(PutValue);

[1] 9.128698

圖4.5.1

第六节 R与C语言的沟通

使用R、Python或是Matlab这类解译式的语言，一个很大的问题是执行效率不高，通常，C是首选的语言，一方面C的执行效能高，另一方面C在计算器的领域已是市场的标准。在执行模拟法中，有效产生大量高质量的随机数，是一个很大的计算负荷，因此我们打算以C语言撰写MT的链接库，再由R来呼叫，产生需要的大量随机数，以加快模拟的效能。此范例使用VisualStudio 2010做为开发工具。

一、使用Visual Studio 产生动态链接库(DLL)

圖4.6.1

首先，以新增项目的方式，产生RUseDLL方案以及MTDLL项目。选取Win32控制台应用程序。按确定，出现Win32应用程序精灵。按下一步。

圖4.6.2

圖4.6.3

在Win32应用程序精灵中，应用程序类型选取DLL，其他选项选取空项目。

在项目中加入原始程序文件，MT19937.cpp。

圖4.6.4

查阅项目属性，组态类型为动态链接库(.dll)。

圖4.6.5

在C/C++的一般页面中，其他Include目标要含入R的程序开发所需的文件头，此目录是安装R时的设定。

圖4.6.6

二、MT程序代码

#001 //*******************************************************************

#002 #include <R.h>

#003

#004 /* Periodparameters */  

#005 #define N 624

#006 #define M 397

#007 #define MATRIX_A0x9908b0dfUL   /*constant vector a */

#008 #define UPPER_MASK0x80000000UL /* most significant w-r bits */

#009 #define LOWER_MASK0x7fffffffUL /* least significant r bits */

#010

#011 static unsigned long mt[N]; /* the array for the state vector  */

#012 static int mti=N+1; /* mti==N+1means mt[N] is not initialized */

#013

#014//*******************************************************************

#015 extern "C" __declspec(dllexport)

#016 /*initializes mt[N] with a seed */

#017     voidinit_genrand(unsigned longs)

#018 {

#019     mt[0]= s & 0xffffffffUL;

#020     for(mti=1; mti<N; mti++) {

#021         mt[mti] =

#022         (1812433253UL * (mt[mti-1] ^ (mt[mti-1]>> 30)) + mti);

#023         /* See KnuthTAOCP Vol2. 3rd Ed. P.106 for multiplier. */

#024         /* In theprevious versions, MSBs of the seed affect  */

#025         /* only MSBs of the array mt[].                        */

#026         /*2002/01/09 modified by Makoto Matsumoto             */

#027         mt[mti] &= 0xffffffffUL;

#028         /* for>32 bit machines */

#029     }

#030 }

#031

#032 /* initializeby an array with array-length */

#033 /* init_keyis the array for initializing keys */

#034 /* key_lengthis its length */

#035 /* slightchange for C++, 2004/2/26 */

#036 void init_by_array(unsigned longinit_key[], int key_length)

#037 {

#038     int i, j,k;

#039     init_genrand(19650218UL);

#040     i=1; j=0;

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

#042     for (; k;k--) {

#043         mt = (mt ^ ((mt[i-1] ^ (mt[i-1]>> 30)) * 1664525UL))

#044           + init_key[j] + j; /* non linear */

#045         mt &= 0xffffffffUL; /* for WORDSIZE > 32 machines */

#046         i++; j++;

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

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

#049     }

#050     for(k=N-1; k; k--) {

#051         mt = (mt ^ ((mt[i-1] ^ (mt[i-1]>> 30)) * 1566083941UL))

#052           - i; /*non linear */

#053         mt &= 0xffffffffUL; /* for WORDSIZE > 32 machines */

#054         i++;

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

#056     }

#057

#058     mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */

#059 }

#060

#061 /* generatesa random number on [0,0xffffffff]-interval */

#062 unsigned long genrand_int32(void)

#063 {

#064     unsigned long y;

#065     static unsigned longmag01[2]={0x0UL, MATRIX_A};

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

#067

#090     /* Tempering */

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

#092     y ^= (y << 7) & 0x9d2c5680UL;

#093     y ^= (y << 15) & 0xefc60000UL;

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

#095

#096     return y;

#097 }

.............................

#134

#135//*******************************************************************

#136 extern "C" __declspec(dllexport)

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

#138     voidgenerate(int* dim, double*result)

#139 {

#140     unsigned long a, b;

#141     for(int i=0; i<*dim; i++)

#142     {

#143          a=genrand_int32()>>5, b=genrand_int32()>>6;

#144          *(result+i) = (a*67108864.0+b)*(1.0/9007199254740992.0);

#145     }

#146 }

程序行表4.6.1

在MT19937.cpp中输入列表4.1程序代码。#002先把R的挡头含入。#015的宣告extern "C" __declspec(dllexport)，是针对#017 voidinit_genrand(unsigned longs)。由于CPP的编译程序会对函数做复杂的修饰，为了要使R认的此函数，我们要告知编译程序要以C语言的格式产出此函数的链接库。#136的宣告 extern "C" __declspec(dllexport)也是针对#138 void generate(int* dim, double*result)同样的作用。如此，我们便可在R中，引用这两个函数。读者如果需要由R中呼叫其他函数，也可类似处理。

建置此项目，输出结果如下。其中MTDLL.dll是R要引用的链接库。将其拷贝的R的工作目录，D:\RProjects下，以备R的呼叫使用。

圖4.6.7

R的工作目录，D:\RProjects下的内容。

圖4.6.8

三、R使用DLL

#005 > setwd("D:\\RProjects")

#006 > dyn.load("MTDLL.dll")

#007 > seed = 1234

#008 > dim = 1024

#009 > .C("init_genrand", as.integer(seed))

#010 [[1]]

#011 [1] 1234

#012

#013 > Res <- .C("generate", as.integer(dim), result =double(dim))

#014 > Res

#015 [[1]]

#016 [1] 1024

#017

#018 $result

#019     [1] 0.9981249604 0.42447415090.5900033008 0.4061420923 0.3762355761 0.4631217649 0.8263664039

#020     [8] 0.1999571176 0.31256047810.6980994841 0.5451827569 0.7955019851 0.8777663181 0.5317340976

#021      …

#022      …

#023      …

#026 [1016] 0.0050764715 0.77529226810.2824727317 0.0860086472 0.2160509267 0.1308943564 0.2769393093

#027 [1023] 0.3684407598 0.5566065340

#028

#029 > v1 = Res[[1]]

#030 > v2 = Res[[2]]

#031 > mean(v2)

#032 [1] 0.498358

#033 > var(v2)

#034 [1] 0.08669075

#035 > dyn.unload("MTDLL.dll")

程序行表4.6.2

在R的交谈环境中输入程序行表4.2，#005设定工作目录到D:\RProjects。注意#005setwd(“D:\\RProjects”)中的双斜线，这是控制程序中跳脱字符的方式。#006动态加载MTDLL.dll链接库。#009呼叫init_genrand(intseed)函数。#013呼叫generate(int dims, double* result)函数，并将结果传回给 Res。#014显示Res的内容，一共有1024个均等分布[0, 1)随机数。Res为一个列表的数据结构，内有两个成分。第一个成分为dims，第二个成分为result。#029将dims传给v1，#030将result传给v2。#031与#033分别计算v2的平均数与变异数。