SAS和蒙物卡羅模擬技術(shù)

2024年2月13日發(fā)(作者：天真)

SAS和蒙特卡羅模擬_隨機數(shù)

(2020-11-10 16:24:47)

標(biāo)簽：

轉(zhuǎn)載自：

一、什么緣應(yīng)選擇SAS做蒙特卡羅模擬

什么緣故要用SAS做蒙卡第一，對我來講，我只會用SAS，而且打算用SAS完成我所有的工作。固然，其他一些通用的理由有(Fan, etc.,2002)：

1. 蒙卡是個計算密集的活，而SAS Ba、SAS Macro、SAS/IML壯大而靈活的編程能力能知足這一要求；

2. 做蒙卡時要用到大量的統(tǒng)計/數(shù)學(xué)技術(shù)，而SAS就內(nèi)置了大量的統(tǒng)計/數(shù)學(xué)函數(shù)(在

SAS/Stat和SAS/ETS)；用Fortran或C++固然也是超級好的主意，只是他們?nèi)鄙賰?nèi)置的統(tǒng)計函數(shù)，代碼就要冗長復(fù)雜很多。

二、什么是蒙卡一個啟發(fā)性例子

好，開始，什么是蒙卡了解它背景知識的最好方法固然是。蒙特卡羅是位于摩洛哥的一家賭場，二戰(zhàn)時，美國Los Alamos國家實驗室把它作為核裂變運算機模擬的代碼名稱。作為模擬方式，蒙卡以前就叫統(tǒng)計抽樣(statistical sampling)，咱們感愛好的結(jié)果因為輸入變量的不確信而不可知，但如果是能依概率產(chǎn)生輸入變量的樣本，咱們就能夠夠估量到結(jié)果變

量的散布。跟蒙卡對應(yīng)的，還有一種模擬技術(shù)叫系統(tǒng)模擬，包括排隊、庫存等模型，這些模型都跟從時刻推移而顯現(xiàn)的事件序列有關(guān)。下面舉個蒙卡的例子，來自Evans, etc.( 2001)的超級簡化版。

假設(shè)一家企業(yè)，利潤是其需求量的函數(shù)，需求是隨機變量。為了簡化討論，假定利潤確實是需求的兩倍。那個地址輸入變量確實是不可控的需求，結(jié)果變量確實是咱們感愛好的利潤。假設(shè)需求以相同的概率取10、20、30、40、50、60這六種情形。在如此的簡化下，咱們就能夠夠投一枚均勻的骰子來產(chǎn)生需求的樣本，若是點數(shù)為1，對應(yīng)得需求確實是10，點數(shù)為2，需求確實是20，以下類推。如此，咱們的模擬進程確實是：

1. 投骰子；

2. 依照骰子的點數(shù)確信需求量；

3. 依照需求量，求利潤。

擲10次骰子，假設(shè)咱們的模擬結(jié)果如下：

重復(fù)次數(shù)

1

2

3

4

5

6

7

骰子點數(shù)

5

3

3

6

1

3

4

需求量

50

30

30

60

10

30

40

利潤

100

60

60

120

20

60

40

8

9

10

5

2

5

50

20

50

100

20

50

如此通過模擬需求的樣本，咱們對結(jié)果利潤的散布也就有所知曉，比如平均利潤能夠算出確實是63。蒙卡一個重要的步驟確實是生成隨機數(shù)，那個地址咱們是用投骰子來完成。

三、又一個例子：利用蒙特卡羅模擬方式求圓周率∏(pi)

再舉個很出名的例子，確實是估量圓周率∏的值，來自Ross(2006)。那個實驗的思路正好能夠幫咱們溫習(xí)一下幾何概型的概念。咱們明白概率的古典概型，確實是把求概率的問題轉(zhuǎn)化為計數(shù)：樣本空間由n個樣本點組成，事件A由k個樣本點組成，那么事件A的概率確實是k/n。考慮到概率和面積在測度上具有某種共性，幾何概型的大體方式是把事件跟幾何區(qū)域相對應(yīng)，用面積來計算概率，其要點是：

1. 以為樣本空間Ω是平面上的某個區(qū)域，其面積記為υ(Ω)；

2. 向區(qū)域Ω隨機拋擲一點，該點落入Ω內(nèi)任何部份區(qū)域內(nèi)的可能性只與這部份區(qū)域的面積成比例，而與這部份區(qū)域的位置和形狀無關(guān)；

3. 設(shè)事件A是Ω的某個區(qū)域，面積為υ(A)，那么向區(qū)域Ω上隨機拋擲一點，該點落在區(qū)域A的可能性稱為事件A的概率，P(A)=υ(A)/υ(Ω)。

扯遠(yuǎn)了，回到用蒙卡估量圓周率∏的實驗思路。假設(shè)樣本空間是一個邊長為2面積為4的正方形，咱們感愛好的

事件是正方形內(nèi)的一個半徑為1面積為∏的圓，因此向正方形內(nèi)隨機拋擲一點，落在圓里面的概率為∏/4。實驗的思路如下：

1. 1）生成隨機數(shù)——生成n個均勻落在正方形內(nèi)的點；

2. 2）對落在正方形內(nèi)的n個點，數(shù)一數(shù)正好落在圓里面的點的個數(shù)，假設(shè)為k（另外n-k個點就落在圓外面的正方形區(qū)域內(nèi)）。

3. 3）k/n就能夠夠大致以為是圓的面積與正方形的面積之比，另其等于∏/4，就能夠夠求出圓周率∏的估量值。

又，提供了利用電子表格Excel求解以上問題的詳細(xì)進程，有愛好能夠看看。另外，也有一個詳細(xì)的說明，用Mathematica實現(xiàn)。

四、隨機數(shù)很重要（和下期預(yù)報）

在第一個例子中，我強調(diào)了骰子若是均勻的。那里骰子確實是咱們的隨機數(shù)生成器，骰子的均勻程度確實是咱們隨機數(shù)的“隨機”成份。在上面的10次實驗中，投出了3次點數(shù)“3”和3次點數(shù)“5”，然后點數(shù)“1”、“2”、“4”、“6”各顯現(xiàn)一次——因為只是重復(fù)了10次，如此咱們對骰子的均勻程度不行評估，但如果是重復(fù)無數(shù)次——假設(shè)是十萬次，若是還顯現(xiàn)類似比例的結(jié)果，那么就有理由疑心這粒骰子的均勻程度了。若是骰子灌了鉛，比如投出點數(shù)“6”的可能性從六分之一上升到6分之三，那么依照這粒骰子來做模擬，咱們就有高估需求和利潤的危險。

下次就講講隨機數(shù)的生成原理，固然結(jié)合SAS來講，或也會提一下Excel。

五、其他軟件包

在商業(yè)領(lǐng)域，基于Excel的插件比較興盛，比如應(yīng)用就較為普遍，有試用版能夠下載。其他競爭產(chǎn)品有、等。此刻聽說也開始興盛起來，大有后發(fā)先至之勢。他的開發(fā)者Johnathan

Mun還在Wiley Finance系列有一本書，。

S語言（R或S-Plus）由于其面對對象的特性，加上豐碩的內(nèi)置函數(shù)和諸多用戶提供的庫，使得R或S-Plus也是蒙卡研究的得力工具——或比SAS更有前景。那個我不熟，略之。

Random Numbers

最簡單、最大體、最重要的隨機變量是在[0,1]上均勻散布的隨機變量。一樣地，咱們把[0,1]上均勻散布隨機變量的抽樣值稱為隨機數(shù)，其他散布隨機變量的抽樣都是借助于隨機數(shù)來實現(xiàn)的。以下談的都是所謂“偽隨機數(shù)”(Pudo Random Numbers)。產(chǎn)生隨機數(shù)，能夠通過物理方式取得（好久好久以前，蘭德公司就曾以隨機脈沖源做信息源，利用電子旋轉(zhuǎn)輪來產(chǎn)生隨機數(shù)表），但現(xiàn)今最為普遍的乃是在運算機上利用數(shù)學(xué)方式產(chǎn)生隨機數(shù)。這種隨機數(shù)依照特定的迭代公式計算出來，初值確信后，序列就能夠夠預(yù)測出來，因此不能算是真正的隨機數(shù)（就成為“偽隨機數(shù)”）。只是，在應(yīng)用中，只要產(chǎn)生的偽隨機數(shù)列能通過一系列統(tǒng)計查驗，就能夠夠把它們當(dāng)做“真”隨機數(shù)來用。

此刻大多軟件包內(nèi)置的隨機數(shù)產(chǎn)生程序，都是利用同余法(Congruential Random

Numbers Generators)。“同余”是數(shù)論中的概念。

0.預(yù)備知識：同余

撿回小學(xué)一年級的東西："4/2=2"讀作“4除以2等于2”，或，“2除4等于2”。還有求模的符號mod(number,divisor)，其中，number是被除數(shù)（在上式中，為4），divisor

是除數(shù)（上式中的2）。如此的約定對SAS和Excel都通用，如mod(4,13)=4，mod(13,4)=1。

此刻咱們能夠開講“同余”了。設(shè)m是正整數(shù)，用m去除整數(shù)a、b，取得的余數(shù)相同，那么稱“a與b關(guān)于模m同余”。上面的概念能夠讀寫成，對整數(shù)a、b和正整數(shù)m，假設(shè)mod(a,m)=mod(b,m)，那么稱“a與b關(guān)于模m有相同的余數(shù)（同余）”，記做a≡b(mod m)（這確實是同余式）。舉個例子，mod(13,4)=1，mod(1,4)=1，那么讀成13和1關(guān)于模4同余，記做13≡1(mod 4)。固然，同余具有對稱性，上式還能夠?qū)懗?=13(mod 4)。

a≡b(mod m)的一個充要條件是a=b+mt，t是整數(shù)，比如13=1+3*4。a=b+mt能夠?qū)懗?a-b)/m=t，即m能整除(a-b)。

這些就夠了。更多基礎(chǔ)性的介紹，能夠參考《》。

1.乘同余法

同余法是一大類方式的統(tǒng)稱，包括加同余法、乘同余法等。因為這些方式中的迭代公式都能夠?qū)懗缮厦嬖蹅円娺^的同余式形式，故統(tǒng)稱同余法。經(jīng)常使用的確實是下面的乘同余法(Multiplicative Congruential Generator.)。符號不行敲，做些約定，如R(i)確實是R加一個下標(biāo)i。

乘同余法隨機數(shù)生成器的同余形式如下：R(i+1)=a*R(i) (mod m)。那個迭代式能夠?qū)懗筛庇^的形式，R(i+1)=mod[a*R(i),m]，其中初值R(0)稱為隨機數(shù)種子。因為mod(x,m)老是等于0到m-1的一個整數(shù)，因此最后把R(i+1)那個隨機序列都除以m，就能夠夠取

得在[0,1]上均勻散布的隨機數(shù)。下面用電子表格演示一遍，假設(shè)隨機數(shù)種子R(0)=1，a=4，m=13：

1

2

3

4

A

i

0

1

2

B

R(i)

1

=D2

=D3

C

a*R(i)

D E

mod[a*R(i),m] mod[a*R(i),m]/m

=B2*4 =mod(C2,13) =D2/13

=B3*4 =mod(C3,13) =D3/13

=B4*4 =mod(C4,13) =D4/13

上面顯示的是公式（Excel中，公式與計算結(jié)果的轉(zhuǎn)換，用快捷鍵ctrl+~實現(xiàn)），結(jié)果看起來是如此的，E列確實是咱們想要的在[0,1]上均勻散布的隨機數(shù)：

1

2

3

4

A

i

0

1

2

B

R(i)

1

4

3

C D E

a*R(i) mod[a*R(i),m] mod[a*R(i),m]/m

4

16

12

4

3

12

上表中，a*R(1)=16，R(2)=3，m=13，一個同余式確實是R(2)=a*R(1) (mod m)，3=16(mod 13)，或說，16-3能夠被13整除——同余法確實是這意思了。說“乘同余法”，要點在于a*R(i)是乘法形式。在SAS系統(tǒng)中，a=397,204,094，m = 2^31-1=47（一個素數(shù)），隨機種子R(0)能夠取1到m-1之間的任何整數(shù)。

2.偽隨機數(shù)的查驗

此刻內(nèi)置于各大軟件包的隨機數(shù)生成器都通過了完全的查驗，咱們固然能夠安心地利用那個黑箱。或那么咱們也能夠多了解些。一個好的“偽”隨機數(shù)列，應(yīng)該看起來就跟從[0,1]上均勻散布的隨機數(shù)列中隨機抽掏出來的一樣。兩個比較直觀的查驗有：

1. 均勻性查驗——所有的數(shù)是不是均勻地散布在[0,1]區(qū)間；

2. 獨立性（或不相關(guān)）查驗——序列中的數(shù)不存在序列相關(guān)，每一個數(shù)都跟其前后顯現(xiàn)的數(shù)獨立無關(guān)。一個例子是如、、、、，那個地址每一個接踵的數(shù)都正比如它前面的一個數(shù)大，如此那個數(shù)列就似乎有了某種“格局”。

具體的查驗方式就略之了，更多請參見徐鐘濟（1985）。

3.生成其他概率散布的隨機數(shù)

前面提到過，咱們把[0,1]上均勻散布隨機變量的抽樣值稱為隨機數(shù)，其他散布隨機變量的抽樣都是借助于隨機數(shù)來實現(xiàn)的。對其他持續(xù)型散布，（積存）散布函數(shù)為F(x)，r是在[0,1]上均勻散布的隨機變量，另r=F(x)，解出的x確實是該持續(xù)型散布的隨機抽樣。由于x能夠?qū)懗蓃的反函數(shù)，該變換就稱作“逆變換”(Inver Transformation method)。對離散型散布，思路類似，只是繁瑣些，具體能夠見徐鐘濟（1985）、Ross（2006）和Evans等（2001）。大多數(shù)相關(guān)軟件包都會提供各類散布的隨機數(shù)生成函數(shù)，明白有這回事就能夠夠。最重要的，是咱們陳述一類問題時，明白需要用哪一種概率散布來描述Evans等（2001）：

經(jīng)常使用的持續(xù)散布

1. 均勻散布(Uniform Distribution)描述在某最小值和最大值之間所有結(jié)果等可能的隨機變量的特性。

2. 正態(tài)散布(Normal Distribution)是對稱的，具有中位數(shù)等于均值的性質(zhì)，確實是咱們熟悉的鐘形形狀。各類類型的誤差常常是正態(tài)散布的。最后，作為中心極限定理的推論，大量具有任意散布的隨機變量的平均數(shù)也呈正態(tài)散布。

3. 三角形散布(Triangular Distribution)由三個參數(shù)來概念，最大值、最小值和最可能值。臨近最可能值的結(jié)果比那些位于端點的結(jié)果有更大的顯現(xiàn)機遇。

4. 指數(shù)散布(Exponential Distribution)經(jīng)常使用來構(gòu)建在時刻上隨機重現(xiàn)的事件的模型，如顧客抵達(dá)效勞系統(tǒng)的時刻距離、機械元件失效前的工作時刻等。它的要緊性質(zhì)是“無經(jīng)歷性”(Memoryless) ，即當(dāng)前時刻對以后結(jié)果沒有阻礙。例如，不論機械原先已經(jīng)運轉(zhuǎn)了多長時刻，它繼續(xù)運轉(zhuǎn)直至顯現(xiàn)故障的時刻距離總有一樣的散布。

5. 對數(shù)正態(tài)散布(Lognormal Distribution)：假設(shè)隨機變量x的自然對數(shù)是正態(tài)的，那么x就服從對數(shù)正態(tài)散布。對數(shù)正態(tài)散布的常見例子是股票價錢。

經(jīng)常使用的離散散布

1. 貝努利散布(Bernoulli Distribution)描述的是只有兩個以常數(shù)概率顯現(xiàn)的可能結(jié)果的隨機變量的特性；

2. 二項散布(Binomial Distribution)給出每次實驗成功概率為p的n次獨立重復(fù)貝努利實驗的模型；

3. 泊松散布(Poisson Distribution)用于成立某種氣宇單位內(nèi)發(fā)生次數(shù)的模型的一種離散散布，如某時段內(nèi)某事件發(fā)生的次數(shù)。

其他有效的散布如伽瑪散布(Gamma Distribution)、威布爾散布(Weibull Distribution)、貝塔散布(Beta Distribution)略之。

4.下期預(yù)報

上次落了些東西。說到蒙卡與C++，在數(shù)量金融領(lǐng)域，不能不提的一本書確實是Mark Joshi的C++ Design Patterns and Derivatives Pricing (Cambridge Uni. Press, 2004)，面試中必然要說自己讀過的，如這人家以為你至少會用C++做一個蒙特卡羅模擬。那個系列只講SAS，下次先講如何用SAS生成隨機數(shù)，然后確實是具體的項目。

SAS隨機數(shù)函數(shù)及CALL子程序

**************************************************************************************************************************

《》——簡介，通過例子成立起蒙卡的直觀概念；參考軟件包及書目

《》——隨機數(shù)入門；參考書目

**************************************************************************************************************************************************

一、SAS隨機數(shù)函數(shù)和CALL子程序

SAS系統(tǒng)產(chǎn)生隨機數(shù)，兩種方式，利用SAS函數(shù)(Functions)或CALL子程序(CALL

Routines)，它們的語法格式是（具體的區(qū)別容后討論）：

方式

函數(shù)

代碼 說明

var=name(ed,) var為存儲隨機數(shù)列的變量，name為特定的分布函數(shù)形式，ed為隨機數(shù)種子，為特定分布要求的參數(shù)（可選）

CALL子程序 call 同上,記得ed=0, ±1,±2, , ±

name(ed,,var) (2**31-2)

SAS可用的隨機數(shù)函數(shù)列表如下（能夠參見SAS Help and Documentation-SAS

Products-Ba SAS-SAS Language Dictionary-Functions and CALL

Routines-Functions and CALL Routines by Category）：

分布 SAS函數(shù) 說明

n為獨立實驗的次數(shù)，p為成功概率

二項分布(Binomial) ranBin(ed,n,p)

指數(shù)分布(Exponential)

ranExp(ed)

正態(tài)分布(Normal) ranNor(ed),or

normal(ed)

ranNor和normal是同質(zhì)的，但normal沒有相對應(yīng)的CALL子程序

泊松分布(Poisson) ranPoi(ed,m) m為均值

均勻分布(Uniform) ranUni(ed),or

uniform(ed)

ranUni和uniform是同質(zhì)的，但uniform沒有相對應(yīng)的CALL子程序

柯西分布(Cauchy) ranCau(ed)

伽瑪分布(Gamma) ranGam(ed,a)

a>0為形狀參數(shù)

由分布律表格決定的ranTbl(ed,p1,p2,...pn) ∑p=1

離散分布(tabled

probability

distribution)

三角分布(Triangular)

ranTri(ed,h) h為斜邊（最可能值）

上面的隨機函數(shù)，除normal和uniform，都能夠由直接相應(yīng)的CALL子程序挪用。

二、SAS隨機數(shù)函數(shù)和CALL子程序：比較

用SAS隨機數(shù)函數(shù)同時創(chuàng)建的多個隨機數(shù)變量其實都屬于同一個隨機數(shù)列。這話費解，一個例子先，創(chuàng)建兩個隨機數(shù)變量，各包括3個記錄，其中x1的種子為123，x2的種子為456：

data ranuni(drop=i);

retain ed1 123 ed2 456;

do i=1 to 3;

x1=ranuni(ed1);

x2=ranuni(ed2);

output;

end;

run;

proc print data=ranuni;run;

結(jié)果為：

Obs ed1 ed2 x1 x2

1 123 456

2 123 456

3 123 456

仿佛沒什么異樣。咱們把上面的x1增加為6個記錄：

data ranuni2(drop=i);

retain ed1 123;

do i=1 to 6;

x1=ranuni(ed1);

output;

end;

run;

proc print data=ranuni2;run;

結(jié)果如下，把上下用紅色標(biāo)出的數(shù)字對照看一看：

Obs ed1 x1

1 123

2 123

3 123

4 123

5 123

6 123

什么意思在第一段代碼中，x2的種子全然不起作用，把x2的記錄安插到x1里，看起來確實是用種子123產(chǎn)生的隨機數(shù)列加長了罷了。x2的第一個值并非是由種子456產(chǎn)生的，而是產(chǎn)生第一個x1后的下一個值，x1、x2其實屬于同一個隨機數(shù)列，盡管x2的種子被指定為456，而x1的被指定為123。此刻就能夠夠重復(fù)上面的一句話：用SAS隨機數(shù)函數(shù)同時創(chuàng)建的多個隨機數(shù)變量其實都屬于同一個隨機數(shù)列。

用CALL子程序就能夠夠同時產(chǎn)生獨立的隨機數(shù)列。

data ranuni3(drop=i);

retain ed3 123 ed4 456;

do i=1 to 3;

call ranuni(ed3,x3);

call ranuni(ed4,x4);

output;

end;

run;

proc print data=ranuni3;run;

結(jié)果如下：

Obs ed3 ed4 x3 x4

1 28 6

2 6 4

3 4 0

以上x3確實是x1。x1和x3的初始種子都是123，但x3那個結(jié)果顯示的種子是當(dāng)前種子值。要在SAS隨機數(shù)函數(shù)語句中顯示隨機種子的當(dāng)前值，能夠由以下代碼給出：

data ranuni4(drop=i);

retain ed1 123;

do i=1 to 3;

x1=ranuni(ed1);

ed=x1*(2**31-1);

output;

end;

run;

proc print data=ranuni4;

var ed1 ed x1;

run;

結(jié)果如下，能夠跟上面由CALL子程序得出的結(jié)果對照：

Obs ed1 ed x1

1 123 28

2 123 6

3 123 4

---------參考資料---------

1. Xitao Fan, etc..Monte Carlo Studies: A Guide for Quantitative

Rearchers. SAS Institute Inc.,2002

2. 朱世武《SAS編程技術(shù)與金融數(shù)據(jù)處置》，北京：清華大學(xué)出版社，2003