1、第十五章 蒙特卡罗模拟和自助法,主要内容,随机数的生成蒙特卡罗模拟重复抽样自助法,实验15-1： 随机数的生成,实验基本原理服从在(0,1)区间上均匀分布的随机变量的样本值，被称为“均匀分布随机数”，简称“随机数”（random number）。有了均匀分布的随机数后，就几乎可以产生任何分布的随机数了。然而，由电脑产生的所谓“随机数”只是“伪随机数”（pseudo random number），因为它仍然由递推公式产生。尽管如此，这些“伪随机数”都能通过独立性与均匀分布的统计检验，故可以作为“随机数”来使用。“随机数”初始值被称为“种子”（seed）。为了使随机样本具有可重复性，保证以后抽样或
2、别人抽样也能得到完全一样的样本，在生成随机数时通常需要设定“种子”。如果不给定“种子”，则Stata按照“电脑内置钟表”（computer clock）来自动选择“种子”，这样每次抽样的样本就会不同，模拟的结果也就不能完全复制。,实验内容及数据来源本实验中，我们会介绍如何生成均匀分布和正态分布的随机数，以及如何设定随机数的种子。讲解这些内容不需要使用具体的数据文件，一些简单的命令就可以实现。,实验操作指导1 随机数的生成最基本而常用的随机数序列是均匀分布随机数序列。生成均匀分布随机数序列的基本命令为：generate newvar=runiform()其中，generate为生成新变量的基本命
3、令，newvar为新变量的名称，runiform()是生成均匀分布于区间0, 1)随机数的函数。需要注意的是，runiform()中没有参数，但括号却必不可少。如果要生成位于其他区间的均匀分布，我们可以进行简单的变形。例如，要生成均匀分布于区间a, b)的随机数，相应的函数为：a+(b-a)* runiform()要生成均匀分布于区间a, b的随机数，相应的函数为：a+int(b-a+1)* runiform()其中，函数int()表示取整。生成标准正态分布的随机数的函数为：invnorm(uniform()生成均值为m、标准差为s的正态分布的随机数的函数则是：m+s*invnorm(runi
4、form(),例如，我们要生成一个均值为3、方差为5且服从正态分布的序列，可输入命令：set obs 80gen norm=3+5*invnorm(uniform()这里，第一步为设定观测值的个数为80；第二步我们生成一个均值为3、方差为5且服从正态分布的序列，并将新生成的变量命名为norm。需要说明的一点是，如果不设定观测值个数，则新变量的观测值个数会与原序列的观测值个数相同；而未打开任何数据文件时，原观测值个数显然为0。下面，我们看一下变量norm的描述统计量。输入命令：sum norm如果我们要作图看一下norm的分布，可输入命令：hist norm, normal这里，hist表示做直
5、方图，选项normal表示画出相应的正态分布。,2 种子的设定设定种子的基本命令为：set seed #|code其中，set seed是设定种子的基本命令。#是正整数，代表runiform()迭代的初始值， Stata的初始设置是123456789。code是字符串，代表随机数生成器runiform()的状态，我们可以用“c(seed)”来查看其当前值。例如，我们输入命令：disp c(seed)这里，disp表示“显示”，c(seed)表示种子的当前值。例如，我们先设定种子，然后生成随机数序列，可输入如下命令：set seed 1001set obs 60gen r=runiform()这
6、里，第一步设定种子为1001，第二步设定观测值个数为60，第三步生成位于0, 1)随机数序列，并将其命名为r。,实验15-2： 蒙特卡罗模拟,实验基本原理通过计算机模拟从已知分布的总体中抽取大量随机样本的计算方法被统称为“蒙特卡罗方法”（Monte Carlo Methods）。在计量经济学中，常使用蒙特卡罗法来确定统计量的小样本性质。我们知道，许多统计量的精确分布没有解析解。一种解决方法是使用大样本理论，用渐近分布来近似真实分布。然而，现实中的样本容量常常较小。,实验内容及数据来源本实验中，我们会介绍蒙特卡洛模拟的基本操作，以及如何编写程序并通过模拟获得真实的显著性水平。这里，我们不需要使用
7、具体的数据文件，通过命令就可以实现相应的模拟。,实验操作指导1 蒙特卡罗模拟的基本操作,2 用蒙特卡罗模拟获得真实显著性水平,为此，先使用命令program定义一个叫“reschi2”程序进行一次抽样与检验，然后用命令simulate来重复此程序1,000 次。程序的命令为：program reschi2, rclass （定义一个叫reschi2的程序。rclass表示结果以 r()形式储存）version 10 （设定所用程序版本为Stata 10）drop _all （删去内存中已有数据）set obs 50 （确定随机抽样的样本容量为50）gen double x = rchi2(1)
8、 （生成服从分布的解释变量x）gen y = 3+ 2*x + rchi2(1)-3 （生成被解释变量）reg y x （线性回归）return scalar t2 = (_bx-2)/_sex （计算t统计量的值，并将其命名为t2。这里，_bx代表x的系数，_sex代表系数的标准差，scalar表明t2是一个标量（数值）return scalar r2 = abs(return(t2)invttail(48,0.025) （是否拒绝原假设）end （程序结束）,其中，命令的倒数第二行的invttail(48, 0.025)表示自由度为48（样本容量50-参数个数2）、显著性水平为5%的双边检
9、验的t统计量的临界值。return(t2)代表我们前面模拟计算并保存的样本t统计量的值t2。如果t2的绝对值大于invttail(48, 0.025)，则应拒绝原假设；这时，r2的返回值为1。此外，由于程序的首行定义结果的形式为rclass，所以，如果我们要引用计算的r2的值，其形式应该是“r(r2)”。对于上面的程序语句，我们可以将其写入一个“Do-File Editor”，将程序保存为扩展名为“ado”的文件，并存放在Stata安装目录的ado文件夹下。这样，我们后面就可以直接使用该程序了。下面，我们进行1000次模拟，输入命令：simulate reject=r(r2), reps(10
10、00) nodots seed(101): reschi2这里，每次模拟的命令来自于程序“reschi2”，选项reps(1000)表明模拟次数为1000次，选项seed(101)表明设定种子为101，选项nodots表示不显示模拟过程的点。表达式“reject=r(r2)”表明我们将每次模拟的返回值r(r2)保存在reject中。因为我们在程序“reschi2”中，每次拒绝原假设时r2的返回值为1；这样，返回值为1的比例即为拒绝原假设的比例。也就是说，序列reject的均值反映了真实的显著性水平。为了获得reject的均值，我们输入命令：mean reject,事实上，我们在模拟之前，可以先
11、运行一次我们所编的程序，检查一下程序是否有问题。对于我们这个程序，可输入如下命令：set seed 2010reschi2其中，第一步设定种子为2010，第二步为运行程序。要知道这一次模拟是否拒绝原假设，我们可以输入命令：disp r(r2)另一种返回程序结果的命令语句为：return list,实验15-3： 重复抽样,实验基本原理,实验内容及数据来源本书附带光盘data文件夹下的“gender.dta”工作文件，是我们为了讲解重复抽样而编制的文件。该文件包括5810个观测值，只有一个变量gender，表示性别。其取值为1时表示女性（female），取值为0时表示男性（male）。该文件中，
12、女性观测值个数为3418，男性观测值个数为2392。此外，本实验中，我们还会利用本书附带光盘data文件夹下的数据文件“resample.dta”来介绍重复抽样部分选项的应用。该文件也是人为构造的，主要变量包括：group=分组变量（取值为A、B、C、D、E），strid=分层变量（取值为1或2，表明被调查者属于哪种类型），x=观测值。利用这些数据，我们来讲解重复抽样的基本命令和相关选项的应用。,实验操作指导1 重复抽样的基本命令,此外，默认情况下，命令bsample会将内存中的数据替换为抽样的观测值，但设定选项weight()会将抽取的样本频数存放在变量varname中，也就是说，这时只有v
13、arname的值改变，原数据不会改变。但选项weight()和选项idcluster()不能同时设定。另外，在bsample命令之后，选项weight(varname)中的varname可以用在Stata的其他命令中作为fweight（如果该命令允许设定fweight）。下面，我们结合“gender.dta”和“resample.dta”数据，对各选项做进一步的说明。,2 简单随机抽样对于“gender.dta”的数据，假设我们要采用简单随机抽样法抽取300个样本，可输入命令：bsample 300下面，我们来看一下数据文件现在的样本容量。输入命令：count如果要只对男性进行简单随机抽样，我
14、们可以利用条件语句。输入命令：bsample 200 if gender=0这里，我们对gender取值为0（表示male）的观测值抽取了200个随机样本。需要注意的是，条件语句中，gender后为两个等号。另外，需要说明的一点是，读者若要复制这个操作，则需要在前面简单随机抽样后重新打开一次数据；否则，这里的分层抽样就会在原来抽样结果的基础上进行。后面的操作也同理。下面，我们看一下抽取的样本的情况。输入命令：tab gender这里，命令tab表示显示变量gender的频数分布表。,3 分层抽样如果我们要令样本中包括100个女性和100个男性，可以采取分层抽样。输入命令：bsample 100
15、, strata(gender)这里，100表示每一层的样本容量都是100，选项strata(gender)表示按变量gender的不同取值来分层。下面，我们看一下分层抽样的样本情况。输入命令：tab gender,此外，注意到原来的数据中约3000个为女性，2000个为男性。这样，我们考虑按照原数据的男女比例进行抽样；也就是说，抽取300个女性，200个男性。要做到这一点，要先生成一个新变量，令gender为female时新变量的值为300，gender为male时新变量的值为200。输入命令：gen st = cond(gender,300,200)这里，我们将新变量命名为st。对于条件函
16、数cond(x, a, b)，其含义为：如果x为真（或取值不是0），则返回a的值；如果x为假（或取值为0），则返回b的值。对于本例，如果变量“gender”的值为female（1），则令变量st的值为300；如果变量“gender”的值为male（0），则令变量st的值为200。下面，我们就可以利用变量st作为“exp”来进行分层抽样。输入命令：bsample st, strata(gender)这句命令的含义为，按变量gender的值进行分层抽样，且变量gender各个取值对应的样本容量为st的值。也就是说，对于gender取值为female的观测值，对应的抽样样本容量为300；对于gend
17、er取值为male的观测值，对应的抽样样本容量为200。下面，我们看一下样本的情况。输入命令：tab gender,4 扩展样本容量首先要说明的一点是，如果不设定要抽取的样本容量，则样本容量会与原数据相同。而如果同时又设定了选项strata()，则各层中抽取的样本个数会与原数据相同。此外，如果要抽取的样本容量超过现有的观测值个数，我们可以先对样本容量进行扩展。输入命令：expand 2则我们将样本容量扩展为原来的2倍（当然，我们这里也可以选择将样本容量扩为原来的3倍等；但对于本实验，2倍就足够了）。这里，需要说明的是，该命令后也可以加条件语句，选择对分层变量的某个取值扩展样本容量。接着，我们重
18、新进行分层抽样，且不设定样本容量“exp”。输入命令：bsample, strata(gender)这里，选项strata(gender)表示按变量gender进行分层。下面，我们来看一下抽取的样本的频数情况。输入命令：tab gender,5 生成频率权重在前面的几个例子中，原样本都被新抽取的样本所替代。但我们也可以使用选项weight()，将抽取的样本频数放到一个新变量中去，并保留原来的样本。例如，我们还是要抽取200个男性样本，并将频数放到一个新变量中。为此，我们先生成一个变量，好用于存放频数。输入命令：gen fweight=.这里，我们生成一个新变量，将其命名为fweight，并令其
19、初始值缺失（即为“.”）。下面，我们进行抽样。输入命令：set seed 1111bsample 200 if gender=0, weight(fweight)这里，第一步命令为设定种子，这是为了复制结果的方便。第二步，我们对gender取值为0的观测值抽取了200个随机样本；选项weight(fweight)表明，我们将抽取的样本的频数存放在变量fweight中。下面，我们看一下样本的频数分布情况。输入命令：tab fweight gender这里，我们生成变量fweight和gender的二维交互表。,6 分组抽样对于分组数据，如果设定选项cluster()，则Stata会按组来进行抽样
20、。下面，我们利用“resample.dta”的数据，重点来讲一下要抽取的组多于原来的组时该如何处理。“resample.dta”的数据只有5组，而我们想要抽取7组的样本。显然，我们应该先对样本进行扩容。但需要注意的是，命令expand在这里是不恰当的。这是因为，expand命令不能为新生成的组设置其自身的标记；这样，Stata会认为数据还是只有5组，从而，抽取7组样本仍不可能。要解决这个问题，我们可以使用分组数据的复制命令。该命令的基本格式为：expandcl = exp if in, cluster(varlist) generate(newvar)其中，expandcl为复制分组数据的基本
21、命令（注意最后一个字母是小写字母el，而非数字1），表达式exp用于指定复制为原来的几份，if代表条件语句，in代表范围语句。选项cluster()用于指定分组变量，选项generate()用于生成新的变量来保存各个组的识别标志。,在对样本进行扩容之前，我们先看一下原来的数据情况，方便与后面进行对比。输入命令：tabstat x, stat(n mean) by(group)这里，tabstat表示用表格显示数据的描述统计量，x表示我们要得到变量x的统计量。选项stat(n mean)表示我们要获得x的观测值个数和均值，by(group)表示我们按组来获得相应的统计量。下面，我们对“resam
22、ple.dta”的数据进行扩容。输入命令：expandcl 2, generate(newg) cluster(group)这里，我们将数据复制为原来的两份。选项cluster(group)指定原来的分组变量为group，选项generate(newg)表示生成新变量newg用来保存复制后的各个组的识别标志。下面，我们生成一个新变量用于保存抽样频数，并进行分组抽样。输入命令：gen fweight= .set seed 1111bsample 7, cluster(newg) weight(fweight)这里，第一步生成用于保存频数的变量fweight，并将其值设定为缺失值。第二步设定种子为
23、1111，以方便结果的复制。第三步进行抽样，以newg为分组变量，抽取了7组，并将抽取的频数保存在变量fweight中。下面，我们看一下抽取的各组的频数情况。输入命令：tab fweight group,7 分层分组抽样对于“resample.dta”的数据，我们知道变量strid分为2层，且每层中有5组。但两层之间的分组变量不能被唯一地识别。要在每一层中获得7组样本，我们可以先使用expandcl命令进行扩展，并生成新变量唯一地识别每层的各个组。为了与后面的结果进行比较，我们先看一下各组各层的频数分布情况。输入命令：tab group strid下面，我们对数据进行扩展。输入命令：expan
24、dcl 2, generate(newg) cluster(group strid)注意，我们这里使用了选项cluster(group strid)，表明为每一个group和strid的组合设置一个组别的标识。下面，我们来进行分层分组抽样。输入命令：gen fweight= .set seed 1111bsample 7, cluster(newg) strata(strid) weight(fweight)其中，第三步的选项设定分组变量为newg，分层变量为strid。这样，我们会在每一层中抽取7个组。下面，我们看一个各层分别抽到了哪些组。输入命令：by strid: tabulate fw
25、eight group其中，“tabulate fweight group”表示绘制fweight和group的二维表，前缀“by strid”表示为各层分别绘制表格。,实验15-4： 自助法,实验基本原理,实验内容及数据来源本书附带光盘data文件夹下的“usaauto.dta”工作文件给出了美国的汽车业相关数据。利用该数据，我们来分析各因素对每加仑油所行驶的里程数的影响，并讲解自助估计的基本操作以及多种置信区间的获得等内容。,实验操作指导1 自助估计的基本操作进行自助法估计的基本命令为：bootstrap exp_list , options: command其中，bootstrap是自助
26、法估计的基本命令，exp_list表示要保存的结果的表达式，command指定一次抽样中执行的命令，options代表其他选项。需要说明的一点是，表达式exp_list的格式与蒙特卡洛模拟处相同；如果命令command每次改变的是估计的模型的系数，则exp_list为可选项且默认值为_b。此外，表15.3列出了主要的options选项。,下面，我们利用“usaauto.dta”的数据来讲解自助法估计的相关操作。首先，我们进行OLS回归，获取回归系数和标准差，以方便和后面的对比。输入命令：regress mpg weight gear foreign下面，我们采取自助法进行回归。输入命令：boo
27、tstrap, reps(100) seed(123): regress mpg weight gear foreign这里，reps(100)表示进行100次重复抽样， seed(123)表示设置种子为123，冒号后为回归的命令。,如果假定“usaauto.dta”的数据按照变量rep78（汽车在78年的修理次数）来聚类，我们可以使用cluster()选项来进行适当的修正。另外，我们这里想获得变量weight和gear_ratio的系数之差。这样，输入命令：keep if rep78 .bootstrap diff=(_bweight - _bgear), seed(123) cluster
28、(rep78): regress mpg weight gear foreign这里，第一步是删掉变量rep78的缺失值；这是因为命令bootstrap不允许选项cluster()中的变量存在缺失值。第二步进行自助法回归，我们将变量weight和gear_ratio的系数之差命名为diff，并设置种子为123，聚类变量为rep78。对于我们这个回归，下面的命令也可以得到相同的结果：bootstrap diff=(_bweight - _bgear), seed(123): regress mpg weight gear foreign, vce(cluster rep78)这里，我们在回归中设
29、置了选项vce(cluster rep78)，表明按变量rep78聚类；这样，bootstrap命令就会按照变量rep78的类别进行抽样。,2 多种置信区间的获得,对于前面的自助回归，我们要获得多种置信区间。这样，先输入如下命令：qui bootstrap, reps(1000) saving(bs) bca seed(123): regress mpg weight gear foreign这里，选项reps()表明抽样次数为1000；如我们前面所讲，要比较可靠地获得置信区间的估计值，抽样次数要大一些。saving(bs)表示将自助估计的结果保存到文件“bs.dta”中。qui表明不显示回归
30、结果。下面，我们来获取多种置信区间。输入命令：estat bootstrap, all这里，选项all表示显示所有可行的置信区间。前面，我们将自助估计的结果保存到文件“bs.dta”中。下面，我们可以利用该文件检验各系数的自助分布是否为渐进正态。输入命令：use bs, clear我们将打开数据文件“bs.dta”。下面，我们来看一下该文件的变量情况。输入命令：describe *这里，*表示各个变量。下面，我们可以对系数自助分布的正态性进行检验。例如，可输入命令：pnorm _b_weight这里，pnorm是通过绘制标准整体概率图来检验变量的正态性。,习题,1.模拟中心极限定理。提示：用命
31、令program定义一个程序，从均匀分布抽取一个容量为n（例如，60）的随机样本，并计算样本均值；然后命令simulate来重复此程序若干次（例如，10000次），得到若干个样本均值的取值；最后用命令histogram来画样本均值的直方图，并与正态分布比较。2.本书附带光盘data文件夹下的“lowbirth2.dta”工作文件给出了按母亲年龄进行匹配的56对婴儿体重的配对数据。主要变量包括：pairid=每组的编号，low=婴儿是否体重偏低（1表示体重偏低），age=母亲的年龄，lwt=母亲怀孕前体重，smoke=母亲怀孕期间是否抽烟（1表示抽烟），ptd=母亲之前是否有早产婴儿（1表示有）
32、，ht=母亲是否高血压（1表示是），ui=母亲是否有子宫炎（1表示有），race=母亲的种族（1表示白人，2表示黑人，3表示其他）。以pairid为分组变量，进行重复抽样。首先抽取30组样本；然后对数据进行扩容，抽取80组样本。要求设置频率权重。,3.仍然利用习题2的数据，利用自助法，分析lwt、smoke、ptd、ht、ui、race等各因素对婴儿体重是否偏低的影响，以pairid为分组变量。因为被解释变量二值变量，数据又是配对数据，考虑条件logit模型。这里需要注意的是，因为每组数据可能被多次抽样抽到，所以需要通过选项idcluster()为每组抽取的样本重新进行标识。（提示：可以采用命令 bootstrap, seed(1) cluster(pairid) idcluster(newpairid): clogit low lwt smoke ptd ht ui race, group(newpairid)）4.获取习题3中自助回归系数的多种置信区间。,
