【正文】 蒙特卡羅模擬與 編程 張曉峒 南開大學(xué)數(shù)量經(jīng)濟(jì) 研究所所長、 博士生導(dǎo)師 東北財(cái)經(jīng)大學(xué) 兼職教授 中國數(shù)量經(jīng)濟(jì)學(xué)會(huì)常務(wù)理事 、 天津市數(shù)量經(jīng)濟(jì)學(xué)會(huì)理事長 n k e vi e w s@ya h om .c n 蒙特卡羅模擬與 編程 1 ．蒙特卡羅（ M on t e C a r l o ）模擬和自舉（ B oo t s t r a p ）的命名過程 2 ．蒙特卡羅模擬和自舉 原理 3 ．設(shè)計(jì)計(jì)算過程 4 ． 生成 服從某種分布的隨機(jī)數(shù) ，生成各種類型的隨機(jī) 序列 5 ． 模擬模型回歸系數(shù)有限樣本估計(jì)量和檢驗(yàn)統(tǒng)計(jì)量的分布特征 6 ．估計(jì)響應(yīng)面函數(shù) 7 ． 極大似然估計(jì) 編程 8 ． 模擬中應(yīng)注意的問題 蒙特卡羅模擬與 編程 大家知道，只有當(dāng)回歸模型滿足 O L S 法所有的假定條件時(shí)，參數(shù)的估計(jì)量才具有最佳線性無偏特性，同時(shí)也具有一致性。 當(dāng)假定條件不成立時(shí)（比如存在異方差、自相關(guān)等），所采用的廣義最小二乘法， 以及對聯(lián)立方程模型的估計(jì)，動(dòng)態(tài)分布滯后模型的估計(jì)，向量自回歸模型的估計(jì)所得參數(shù)的估計(jì)量 只具 有漸近特性（漸近無偏性、一致性）。這意味著，只有當(dāng)樣本容量相當(dāng)大時(shí)，漸近特性才起作用。而當(dāng)樣本容量不是很大，甚至很小時(shí)，仍然不知道估計(jì)量的有限樣本分布特征。 另外通過對非平穩(wěn)過程的研究知道，單位根檢驗(yàn)式和用非平穩(wěn)變量建立的回歸函數(shù)的參數(shù)和 t 統(tǒng)計(jì)量都不服從正態(tài)分布。他們都是漸近地服從 W i e n e r過程的泛函 。參數(shù)估計(jì)量和統(tǒng)計(jì)量的有限樣本特性不能用解析的方法求解。 蒙特卡羅模擬與 編程 對于上述兩種情形，若要研究這些估計(jì)量和統(tǒng)計(jì)量的有限樣本分布特 征，通常采用兩種方法。 一種為 數(shù)值計(jì)算法 。也稱為 有限樣本近似法 （ f i n i t e sa m p l e ap p r ox i m at i o n ）。這種方法要用到許多數(shù)學(xué)知識(shí)，專業(yè)性很強(qiáng)，使沒有受過專門訓(xùn)練的人員運(yùn)用此方法受到限制。 另一種為 蒙特卡羅模擬方法 。又稱 隨機(jī)模擬法 。本章主要介紹蒙特卡羅模擬 與計(jì)算機(jī)編程 。 1 ．蒙特卡羅（ M on t e C ar l o ）模擬和自舉（ B oot s t r a p ） 的命名過程 蒙特卡羅（ M on t e C a r l o ）模擬 是一種通過設(shè)定隨機(jī)過程（數(shù)據(jù)生成系統(tǒng)），反復(fù)生成 隨機(jī) 序列，并計(jì)算參數(shù)估計(jì)量和統(tǒng)計(jì)量，進(jìn)而研究其分布特 征的方法。 “蒙特卡羅模擬” 這個(gè)術(shù)語是 美國 物理學(xué)家 M e t r op o l i s 在 第 2 次世界大戰(zhàn)時(shí)期執(zhí)行曼哈頓計(jì)劃 ( M a n h a t t a n Pr o j e c t ) 過程中 提出的。 作為地名， 蒙特卡羅在歐洲的摩那哥 （ M on ac o ） ，以著名賭城而得名。 若再晚些時(shí)候，蒙特卡羅模擬也許就稱作 L as V e gas （在美國的 N e v ad a 州，著名賭城）模擬方法了。 作為地名，蒙特卡羅在歐洲的摩納哥（ Monaco），以著名賭城而得名。 蒙特卡羅 自舉與蒙特卡羅模擬既有聯(lián)系，又不相同。自舉（ B o ot s t r ap ）這個(gè)名詞是E f r on 在 1979 年提出的。 “ 自舉 ” 一詞來源于 童 話 故事。 指一個(gè)人落水時(shí)，試圖用自提鞋扣兒的方法自 救。 自舉（ B oo t s t r ap ）有人翻譯成 “靴襻” 不恰當(dāng)。 自舉，即采用從 樣本 中反復(fù)抽取 子 樣本的方法計(jì)算參數(shù)估計(jì)量的值，置信區(qū)間或相應(yīng)統(tǒng)計(jì)量的值并估計(jì)這些量的分布。 這里介紹的遠(yuǎn)不是 蒙特卡羅模擬 的全 部 ，而是參數(shù)估計(jì)方面的應(yīng)用。 因?yàn)檫@些方法的實(shí)現(xiàn)是以高容量和高速度的計(jì)算機(jī)為前提條件，所以只是在近年才得到 廣泛 推廣。 20 世紀(jì) 80 ， 90 年代發(fā)展很快 ，現(xiàn)在已很普及 。 2 ．蒙特卡羅模擬和自舉 原理 進(jìn)行蒙特卡羅模擬和 自舉 首先要設(shè)定數(shù)據(jù)生成系統(tǒng)。而設(shè)定數(shù)據(jù)生成系統(tǒng)的關(guān)鍵是要產(chǎn)生大量的隨機(jī)數(shù)。 例如模擬樣本 容量 為 100 的一元線性回歸模型中參數(shù)的分布，若試驗(yàn) 1 萬次，則需要生成 200 萬個(gè)隨機(jī)數(shù)。 計(jì)算機(jī)所生成的隨機(jī)數(shù)并不是“真隨機(jī)數(shù)”，而是具有某種相同統(tǒng)計(jì)性質(zhì)的隨機(jī)數(shù)。計(jì)量經(jīng)濟(jì)學(xué)中蒙特卡羅模擬和 自舉 所用到的隨機(jī)數(shù)一般是服從N ( 0,1 ) 分布或均勻分布的隨機(jī)數(shù)。計(jì)算機(jī)生成的隨機(jī)數(shù)稱作“偽隨機(jī)數(shù)”（ p s e u d o r an d om n u m b e r ）（以下簡稱 隨機(jī)數(shù) ） 。生成的隨機(jī)數(shù)的程序稱作“偽隨機(jī)數(shù)生成系統(tǒng)”。實(shí)際上計(jì)算機(jī)不可能生成真隨機(jī)數(shù)。 3 ．設(shè)計(jì)計(jì)算過程。 蒙特卡羅模擬和 自舉 的實(shí)現(xiàn)要通過計(jì)算機(jī)編程來實(shí)現(xiàn)。常用的高級編程語言軟件有 M at h e m at i c a ， G au s s ， Ox ， E V i e w s ， S P l u s ， M at l ap ， r a t s ， s t a t a ， R 等。 這里主要介紹利用 E V i e w s 做蒙特卡羅模擬與編程。編程也是一個(gè)很寬的概念，這里只介紹與蒙特卡羅模擬有關(guān)的編程。 編程像一門藝術(shù)，需要經(jīng)驗(yàn)和技巧。只有多 思，熟記各種命令，才能編出最優(yōu)的程序來。 否則會(huì)浪費(fèi)大量時(shí)間。 首先需要把自己需要計(jì)算的問題設(shè)計(jì)出合理的計(jì)算流程框圖，然后轉(zhuǎn)化為正確的計(jì)算機(jī)語言。 提出研究的問題 設(shè)計(jì)計(jì)算流程 框圖 轉(zhuǎn)換為計(jì)算機(jī)語言 4 ． 生成 服從某種分布的隨機(jī)數(shù) ，生成各種類型的隨機(jī) 序列 生成 服從某種分布的隨機(jī)數(shù) 、 累計(jì)分布函數(shù)、概率密度函數(shù)、 p 分位數(shù)函數(shù)的表達(dá)語言見表 1 。 表 1 ： 生成 18 種 分布的隨機(jī)數(shù) 、 累計(jì)分布函數(shù)、概率密度函數(shù)、 p 分位數(shù)函數(shù)的 用 語 分布 類型 生成隨機(jī)數(shù) 累計(jì)分布函數(shù) 概率密度函數(shù) p 分位數(shù)函數(shù) B e t a 分布 @ r b e t a( a, b ) @c be t a ( x , a , b ) @dbe t a ( x , a , b ) @qbe t a ( p , a , b ) 二項(xiàng)式分布 @ r b i n o m ( n , p ) @c bi no m ( x, n , p) @db i no m ( x, n , p) @qb i no m ( s , n , p) 卡方分布 @ r c h i s q ( v ) @c c hi s q( x , v ) @dc hi s q( x , v ) @qc hi s q( p , v ) 指數(shù)分布 @ r e xp ( m ) @c e xp( x , m ) @de x p( x , m ) @qe x p( p , m ) 極值分布 I @ r e xt r e m e @c e xt r e m e ( x) @de x t r e m e ( x) @qe x t r e m e ( p) F 分布 @ r f d i s t ( v 1, v 2) @c f di s t ( x , v 1, v 2) @df di s t ( x, v 1, v 2) @qf di s t ( p, v 1, v 2) G am m a 分布 @ r gam m a( b , r ) @c g a m m a ( x, b, r ) @dg a m m a ( x , b, r ) @qg a m m a ( p , b, r ) 廣義誤差分布 @ r ge d ( r ) @c g e d( x, r ) @dg e d( x , r ) @qg e d( p , r ) 拉普拉斯分布 @ r l ap l a c e @c l a p l a c e ( x ) @dl a p l a c e ( x ) @ql a p l a c e ( p ) l o gi s t i c 分布 @ r l o gi s t i c @c l o g i s t i c ( x ) @dl o g i s t i c ( x ) @ql o g i s t i c ( p ) 對數(shù)正態(tài)分布 @ r l o gn o r m ( m , s ) @c l o g n o r m ( x, m , s ) @dl o g no r m ( x, m , s ) @ql o g no r m ( p, m , s ) 負(fù)二項(xiàng)分布 @ r n e gb i n ( n , p) @c ne g bi n( x, n , p) @dne g bi n( x , n , p ) @qne g bi n( s , n , p) 標(biāo)準(zhǔn)正態(tài)分布 @ r n o r m , n r n d @c no r m ( x) @dno r m ( x ) @qno r m ( p ) 泊松分布 @ r p o i s s o n ( m ) @c po i s s o n( x , m ) @dpo i s s o n( x , m ) @qpo i s s o n( p , m ) 帕雷托分布 @ r p ar e t o ( a, k ) @c pa r e t o ( x , a , k ) @dpa r e t o ( x, a , k) @qpa r e t o ( a , k ) t 分布 @ r t d i s t ( v ) @c t di s t ( x , v ) @dt di s t ( x , v ) @qt di s t ( p , v ) 均勻分布 @ r u n i f ( a, b ) @c uni f ( x , a , b ) @dun i f ( x , a , b ) @qun i f ( p , a , b ) 威布爾分布 @ r w e i b ( m , a ) @c w e i b( x, m , a ) @dw e i b( x, m , a ) @qw e i i b ( p , m , a ) （ 1 ） 生成 服從某種分布的隨機(jī)數(shù) 序列 【例】 生成 標(biāo)準(zhǔn)正態(tài)分布 、 指數(shù)分布、 p oi s s on 分布、 t 分布 的隨機(jī)數(shù)序列 E V i e w s 程序如下： 39。 生成 服從某種分布的隨機(jī)數(shù) 序列 （ T= 10 00 ） ，并存入 r an d om 1 文件。 w or k f i l e r an d o m 1 u 1 10 00 s e r i e s Y 1= @r n or m s e r i e s Y 2 = n r n d s e r i e s Y 3 = @r e x p ( 3 ) s e r i e s Y 4 = @r p oi s s on ( 3 ) s e r i e s Y 5 = @r t d i s t ( 3 ) 6420242 5 0 5 0 0 7 5 0 1 0 0 0Y 1標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù) 01020302 5 0 5 0 0 7 5 0 1 0 0 0Y 3指數(shù)分布隨機(jī)數(shù) 02468102 5 0 5 0 0 7 5 0 1 0 0 0Y 43個(gè)自由度的 poisson分布隨機(jī)數(shù) 1 5 1 0505102 5 0 5 0 0 7 5 0 1 0 0 0Y 5df=3的 t分布隨機(jī)數(shù) （ 2 ） 生成 任 意（非標(biāo)準(zhǔn)） 正態(tài)分布 隨機(jī)數(shù) 序列 用 n r n d 和 @r n or m 可以 生成 標(biāo)準(zhǔn)正態(tài)分布 隨機(jī)數(shù)，那么，任何參數(shù)的正態(tài)分布 隨機(jī)數(shù)都可以生成。 用 Z 表示標(biāo)準(zhǔn)正態(tài)分布的 隨機(jī)數(shù)，用 X 表示任 一 參數(shù)的正態(tài)分布的偽隨機(jī)數(shù)，那么，把 X 標(biāo)準(zhǔn)化的公式是 ZsXXX??)(? N ( 0,1 ) 反過來，若用 Z 表示 X ，則公式為 X = Z ? s( X ) + X 按上公式， E V i e w s 可以生成 任 意 正態(tài)分布隨機(jī)數(shù) 。 【例】 生成 均值為 50 ，標(biāo)準(zhǔn)差為 的 正態(tài)分布隨機(jī)數(shù) X ? N( 5 0, 2) 。 se r i e s Z =n r n d se r i e s X= Z * 0. 5+ 50 4 8 . 44 8 . 84 9 . 24 9 . 65 0 . 05 0 . 45 0 . 85 1 . 25 1 . 65 2 . 050 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 4 0 0 4 5 0 5 0 0 5 5 0 6 0 0 6 5 0 7 0 0 7 5 0 8 0 0 8 5 0 9 0 0 9 5 0 1 0 0 0X（ 3 ）在 生成 服從某種分布的隨機(jī)數(shù)序列的基礎(chǔ)上求