計量分析與stata應(yīng)用高級章節(jié)b9 mc

-on2011-3-5123456789-on2011-3-5123456789*********計量分析與STATA應(yīng)用*郵arlionn@163.com頁:::高級部分********計量分析與Stata應(yīng)第九講setschemes2color*JK的基本思想JK應(yīng)用實例JK與BSBSMonteCarloSimulation(蒙特卡羅模擬)MC的基本思想：圓周率的估算檢驗統(tǒng)計量的Power和********************** 簡*偽隨機數(shù)的產(chǎn)生問題的來源：無論是BS還是MC，都需要隨機抽樣，而抽樣的基礎(chǔ)是隨機setobsgenx1=uniform()histogramx1genz=invnormal(uniform())histogramz,normalgenx2=uniform()listx1x2in[0,1seedx3=uniform()seed123x4=listx3x4inQ1:Printedon2011-3-5-Printedon2011-3-5******(j=1,2,X_j=(a*X_{j-1}+c)modX_j=(9*X_{j-1}+3)mod2^9(j=1,2,local=種子值*/上限*/local = =setobsgenid=_ngenz=`seed'forvaluesi=dis_cin}z=mod((9*z[`i'-1]+3),`bound')ifid==`i'dropin/*第一個觀察值是種子tabzlocalbound=genzz=z/`bound'sumzzhistogram/*[0,1)均勻分布的隨機數(shù)*genxinvnormal(zz/**/sumxttestx=0histogramx,normal*N=100； bound=2^4tssetidgenz_lag=scatterz*N=10000；bound=2^9tssetidgenz_lag=L.zscatterz************************a]at,GeorgeMarsaglia(1994)32-bitpseudorandomnumberKISS(KeepItSimpleStupid)====(69069X_j-1+1234567)mod65184(z_j-1mod2^16)+int(z_j-1/2^16)63663(w_j-1mod2^16)+int(w_j-1/2^16)R_j=(x_jy_j+z_j+x_0=y_0=z_0=w_0=/*這就是所謂的“種子(seed)”,setseed改變的就是該值(1)每次使用seed命令改變的是x_0為了保證隨機數(shù)的質(zhì)量，Stata會丟棄前100個隨機數(shù)R_j/2^32便可得到服從均勻分布[0,1)Printedon2011-3-5-Printedon2011-3-5**dis/**setobsgeny1=uniform()disc(seed)geny2=uniform()disc(seed)dis%10.0f2^(31)-1若將種子值設(shè)定為負值，則stata會自動將其轉(zhuǎn)換為正setobssetseed-123genx=uniform()setseed123geny=uniform()listxy為何要設(shè)定種子？seeGentle(2003保證MC或BS*******gen=genxa(b-a)*uniform()[a,b區(qū)間上的隨機正整數(shù)，公式為：genx=a+int((b-setobs/*U[0,1)genx1=genx2=5+(12-/*U[5,12)genx3=5+int((12-5)*uniform())/*intU[5,12)*/listin1/10histogramx1histogramx2histogram正態(tài)分布: 標準正態(tài)分布的隨機數(shù)，即，X--N(0,1)，公式為 genx=mu+setobsgenx1=genx2=10+5*invnormal(uniform())sumx1x2,detailhistogramx1,normalhistogramx2,setobsgenu1=uniform()genu2=genx1=sqrt(-2*ln(u1))*cos(2*_pi*u2)genx2=sqrt(-*******x1--N(0,1);x2--N(0,1);pwcorrx1x2,Printedon2011-3-5-Printedon2011-3-5sumx1x2,detailhistogramx1,normalhistogramx2,其他分布：t(df)分布F(k1,k2)分布，Chi2(df)分布，gamma分布weibull*AppliedRobustStatistics,Chapter*常用分布的數(shù)值生成方法參見*Microeoconometric:MethodsBrndtobsdfrndchiobsdfrndfobsdf_ndf_drndlgnobsmeanvarrndpoiobsmeanrndpoix[mu]rndbinobsprobnumbrndbinx[prob]rndgamobsshapescalerndgamx[mu],s(#)rndivgobsmeansigma***************[mu],rndexpobsBetarndweiobsshapescalerndbbobsdenomprob*chi2*-計算公式：Sum^{k}Sum_1^{k}*-即，krndchi10005histogramxc,normal/*Chi2(5分布rndchi10000/*Chi2(20)twoway(histogramxc)(kdensityxc)sumxc,detail****-***E[x]=0;Var[x]=k/(k-2)(forrndt100005sumxthistogramxt,rndt1000030histogramt(5)分布t(30分布t(120)分布rndt10000histogramxt,****時，t(k分布具有“尖峰厚尾”k—>無窮大時，t(k**-計算公式：-rndexp100003sumxetwoway(histogramrndexp10000120twoway/*Exp(3分布/*E[xe]=Simga[xe]=lamda*/xe)(kdensityxe)/*Exp(120分布xe)(kdensity***xlamdaExp(lamda)，則E[xlamdaVar[xlamda^2*on2011-3-5-on2011-3-5*如，在隨機邊界(SFA)模型中,常用于描述非效率成分的分布特征*F(k1k2**計算公式：rndf3/*F(3,20分布histogramrndf100020*邏輯分布:*r*計算公式：ab*ln其中，r--**E[x]=a1-rsetobsgenr=genx=3+5*ln(r/(1-sumhistogramx*-泊松分布**rndpoi10040histogramxpforx=*-*Cameron(2005)Microeoconometric:Methodsand*AppendixB,截斷型半正態(tài)分布（truncatedstandardnormalhelp這在模擬分析隨機邊界模型(StochasticFrontierModel)時非常有用Olive(2005)AppliedRobustStatisticsChapter4***setobs1000gentrunx,left(0)sumxhistogramgentruny,right(-0.5)sumyhistogram多元正態(tài)分布(bivariatenormal*變量x1和x2構(gòu)成向量[x1x2N(Mu或(X1,X2N(m1m2s1^2s2^2,**matrixm=(2,3)matrixsd=(.5,2)drawnormxy,n(2000)means(m)sds(sd)corryxscatteryx相關(guān)系數(shù)為0.8**matrixmmatrixCsumxycorrx==x(0.3,(1,0.8\0.8,y,n(1000)means(m)Printedon2011-3-5-Printedon2011-3-5scatteryx,三元正態(tài)分布序列:設(shè)定方差-*X(x1x2x3)N(M,S)(M和S為矩陣則Var[XE[X*X'E(x1x2x3)'*(x1x2********[Var(x1)Var[X]=E[[Cov(x1,x3)Cov(x2,x3) matM=5,-6,matV=matlistmatlistVdrawnormx1x2\\x3,n(1000)cov(V)correlatex1x2x3,(1)分解方差-協(xié)方差矩陣：V=A'A(cholesky分解)(2產(chǎn)生K個N(0,1)序列，XN(0,1)，X=(x1,x2,...,xk)'(3YMA'X，則YN(MK****setobsmkbilognx1x2,r(0.8)m1(1)s1(2)m2(3)s2(4)genln_x1=ln(x1)genln_x2=ln(x2)sumln_x1ln_x2corrln_x1ln_x2histogramln_x1,normalxline(1,lw(thick))scatterln_x1ln_x2,msymbol(x)SeeDekking(2005,DividsonandMickinnon(2004,p.159)** 簡*=Boostrap********************x123456個觀察值x_m=4,請問x_m的標準誤，即se(x_m)我們可以采用重復(fù)抽樣的方法得到一系列的經(jīng)驗樣本(empirical如,j=10011233234451224434359====100xi_mx_mBootstrap的基本思想x_m替換為其他統(tǒng)計量，如OLS回歸中的系數(shù)估計值b，Printedon2011-3-5-Printedon2011-3-5**********從“觀測樣本”中進行多次可重復(fù)抽樣可以近似復(fù)制出“母體因此，我們可以通過分析這些“經(jīng)驗樣本”的特征來考察“母體”的分布特征23459saveB9_intro_x,useB9_intro_x,clear50次抽樣的useB9_intro_x,localr=matB=J(`r',7,.)matL=J(6,1,1)forvaluesj1/`r'{preserve/*抽樣次數(shù)mkmatx,mat(s)matB[`j',1]=s'matmean=matB[`j',7]=}matlistsvmatB,names(m)dropm1-m6renamem7meansummeanhistogrammean可見，BS得到的mean值與4非常接近se(x_m)~=***-練習(xí)：將抽樣次數(shù)更改為500*-問題：6*答案：462statauseB9_intro_x,bootstrapx_m=r(mean),reps(50):sumbsx_m=r(mean),reps(500)nodotsnowarnsum********bootstrapbs；reps(saving()選項可將抽樣數(shù)據(jù)保存至指定文件中nodot選項可不在屏幕上打點；nowarn選項可不在屏幕上提示；*Printedon2011-3-5-Printedon2011-3-5usebs_intro1.dta,clearusebs_intro1.dtaclearsample8countBootstrap抽樣：可重usebs_intro1.dta,clearbsample8usebs_intro1.dtaclear/*抽樣具有隨機性*/bsample9************觀測樣本：13451213132233445244359x1_m=x2_m=x3_m=xj_m=4******BSBS若采用BS獲得置信區(qū)間，500-1000次即可若采用BS獲得多個描述樣本分布的統(tǒng)計量，1000次以上；多數(shù)情況下，1000次可重復(fù)抽樣即可獲得非常穩(wěn)定的結(jié)果。*****usebs_panel.dta,cleartssetidyearusebs_panel.dta,clearbsample5,cluster(id)sortidyearusebs_panel.dtaclear/*抽樣具有隨機性*/bsample5,cluster(id)sortidyeartab*usebs_panel.dta,bsample5,cluster(id)idcluster(id_bs)tssetid tssetid_bsyear/*可能產(chǎn)生錯誤*--組合樣本的個數(shù)：為什么采用BS可以近似模擬母體分布2011-3-5-2011-3-5***(((m) (其中，m=2n- ncapprogramdropbs_Nprogramdefinebs_Nargsnlocalm=2*`n'-localN=comb(`m',`n')/*helpmathfunctions*/dising"n="iny`n'ing"N="iny%-20.0fbs_N/*運行程序n=2為例：Xx1,x2)則BS樣本將有3個：(x1,x1)(x2,x2)(x1,x2)***以n=3為例，XBS10(x1,x1,x1);(x2,x2,x2);(x1,x2,x2);(x2,x3,x3);(x2,x2,x3)*******表樣本數(shù)與BS樣本組合個***********nN23453*BS*Thestandarderrorofastatisticisstandardofitsbootstrapdistribution.Itmeasureshowmuchthestatisticvariesrandom*OLSbinv(X'X)*(X'y分別對這些實證樣本進行OLS估計，可以得到一系列OLS估計值(b1,b2,...,b1000)；使用這1000個估計值便可近似得到b的分布特征。******參數(shù)BS與非參數(shù)參數(shù)；否則，稱為“非參數(shù)*************模型:y_i=xi*b+(i=若假設(shè)e_iN(0,s^2)N(0,s^2分布中抽取隨機數(shù)，作為模型(1)的干擾項；因為要生成服從N(0,s^2的隨機數(shù)，必須事先確定參數(shù)若不對e_i的分布情況做任何假設(shè)，則我們可以直接對原始樣本進行自抽樣；Printedon2011-3-5-Printedon2011-3-5非參數(shù)Bootstrap(Non-pramatric*seeEfron(1993,*模型:y_i=xi*b+(i=*Y=Xb+*-方法1:Bootstrapping*BS(Y_bs1,step1step2K次，比如，K=300，將得到b1,b2,...,b300;****計算s.d.(b1,b2,...,b300)blocalreps=300 matB=J(`reps',1,.)forvaluesi=1(1)`reps'{qui/*自抽樣quiregresspricematB[`i',1]=}svmatB,sumlocalse_bsr(sd)/*自抽樣標準誤*/dis`se_bs'與OLS標準誤對比sysuseauto,clearregpriceweight*tz*localb_ols_b[weight]/*OLS估計系數(shù)*/localt_bs=`b_ols`se_bs'dis"t_bs="計算p**localp_z=1-normprob(`t_bs')dis"p-value="%6.4f`p_z'小樣本，tlocalp_t=1-tprob(74-2,*dis"p-value="%6.4f*reps501000*-方法2:Bootstrapping**stataregpriceweight,vce(bootstrap,*step1:regYonX,得到估計系數(shù)b_ols，和殘差序列/*采用真實數(shù)據(jù)*step2BS(E_bs1),構(gòu)造BS樣本，Y_bs1X*b_ols/*X的觀察值保持不變計算s.d.(b1,b2,...,b300)便可得到系數(shù)b的標準誤 sysuseauto,clearregpriceweipredictres,residualkeepresPrintedon2011-3-5-Printedon2011-3-5saveres_data,localreps=/*抽樣次數(shù)matB=forvaluesi=*useres_data,cleargenid=_nsortid/*便于與主數(shù)據(jù)合并tempfilesave"`e_data'",*構(gòu)造BS樣sysuseauto,clearregresspriceweightpredicty_hatkeepy_hatweightgenid=_nsortmergeidusing"`e_data'"geny_bs=y_hat+res對BS樣本進行OLS估計，獲得BS系regy_bsweightmatb=e(b)matB[`i',1]=*}}svmatB,names(b)sumb1localse_bs=/*自抽樣標準誤*****至于t_bs和p說明：(1)tempfile命令的具體使用方法，參見“A5-編程初步*sysuseauto,bootstrap,reps(300)nodotsnoheader:regprice*regpriceweight,vce(bootstrap,*兩種方法比較BSResidualsBSPairs********模型yxbeE(e|X)=0即模型的設(shè)定是正確的。本質(zhì)上講，BSP意味著X_i在抽樣過程中是可變的，BSR則意味著X_iB>1000兩種方法得到的標準誤會非常接近。文獻中使用較多的是第二種方法：bootstrapresiduals 參數(shù)Bootstrap(Pramatric*Efron(1993,sysuseauto,regpriceweilenmpgforeignpredicteresidual/*OLS殘差*/sume,detailhistograme,normalsktesteswilk我們可以仍然假設(shè)e服從正態(tài)分布，并采用BSPrintedon2011-3-5-Printedon2011-3-5**********regyonxs^2；生成服從N(0,s^2)的隨機數(shù)e_n；對e_n進行自抽樣，得到殘差序列e_bs1；構(gòu)造BS樣本，即，Y_bs1X*be_bs1；采用OLS估計BS樣本，得到b_bs1；(b_bs1b_bs2b_bs300)(b_bs1b_bs2b_bs300的"標準差"(s.d.)，即為b的BS"標準誤"(s.e.)sysuseauto,quiregpriceweilenmpgeretlistlocals2=e(rss)/(74-locals=gene_n=`s'*invnorm(uniform())/*生成服從正態(tài)分布的殘差序列res即可*下面的步驟與上例相似，這是把這里的e_n*keepsaveres_data,localreps=1000matB=J(`reps',1,.)/*抽樣次數(shù)forvaluesi=*useres_data,cleargenid=_nsortid/*便于與主數(shù)據(jù)合并tempfilesave"`e_data'",*構(gòu)造BS樣sysuseauto,regresspriceweilenmpgforeignpredicty_hatkeepy_hatweilenmpgforeigngenid=_nsortmergeidusing"`e_data'"geny_bs=y_hat+e_n對BS樣本進行OLS估計，獲得BS系regy_bsweilenmpgforeignmatb=e(b)matB[`i',1]=*}}svmatB,names(b)sumb1localse_bs=/*自抽樣標準誤*OLSsysuseauto,regpriceweilenmpgforeign,小結(jié)：(1)采用這種方法得到的結(jié)果與根據(jù)OLSVar(b)=s^2*(X'X)^{- *****StataBS*Printedon2011-3-5-Printedon2011-3-5sysuseauto,*regresspriceweirep78eststorem_bs*BSregresspriceweivce(bootstrap,**regresspriceweiregresspriceweirep78mpgforeign,regresspriceweirep78mpgforeign,**sysuseauto,cleardropifrep78==bootstrap,cluster(rep78):regpriceweirep78mpgforeigneststorem_bs_clusteresttabm_bsm_bs_cluster,mtitle(m_bs*bootstrap:regpriceweirep78mpgforeign,工具變量法(IVregresshelp******hsngvalusehsng2.dta,ivregrentpcturban(hsngvaleststoreivivregrentpcturban(hsngvaleststoreiv_bs*=famincreg1-=famincreg1-reg4),esttabiviv_bs,*在這種情況下，也可以采用vce()選項得到BS標準誤分位數(shù)回歸(helpsysuseauto,regpriceweightlengthforeigneststoreolsregpriceweightlengthforeign,vce(bs,reps(200))eststoreols_bsqregpriceweightlengthforeigneststorep50bsqregpriceweightlengthforeign,reps(200)eststorep50_bsqregpriceweightlengthforeign,quantile(.75)eststorep75bsqregpriceweightlengthforeign,quantile(.75)reps(200)eststorep75_bslocalmodelolsols_bsp50p50_bsp75p75_bsesttab`model',mtitle(`model')compress*-評論**(3采用BS*Printedon2011-3-5-Printedon2011-3-5useinvest2.dta,xtregmarketinveststock,feeststorefextregmarketinveststock,fevce(bs,nodotsreps(200)seed(135))eststorefe_bslocalmfefe_bsesttab`m',mtitle(`m')useinvest2.dta,cleartsset,clear**bootstrap,reps(200)seed(135)cluster(id)idcluster(new_id):xtregmarketinveststock,fei(new_id)eststore**localmfe_bsfe_bs0esttab`m',mtitle(`m')*********cluster(id選項，idcluster(new_id)選項必須設(shè)定，在BS過程中，statatssetnew_idt命令；cluster(id),idcluster(new_id)和i(new_id) *BS檢驗其它統(tǒng)計量*usextcs.dtaclear/*中國上市公司資本結(jié)構(gòu)數(shù)據(jù)*/sortsizegengroupSize=group(3)tabgroupSize,gen(dsize)dropdsize2genxsize1dsize1*size/*小規(guī)模公司交乘項*/genxsize3dsize3*size/*大規(guī)模公司交乘項xtregtldsize1dsize3xsize1xsize3ndtstangnprtobin,*方法1：傳統(tǒng)的t檢驗testxsize1xsize3tsset,clear*reps(200)seed(123)cluster(code)xtregtldsize1dsize3xsize1xsize3ndtstangnprtobin,fe*xtregtldsize1dsize3xsize1xsize3ndtstangnprtobin,fei(code)///vce(bootstrap(_b[xsize1]-_b[xsize3]),reps(200)seed(123))********t檢驗的結(jié)果和BS的結(jié)果可以發(fā)現(xiàn)，但t檢驗的結(jié)果更加顯著(可能存在過度拒絕lf reps(1000)，看看結(jié)果有何差異？（BS結(jié)論非常穩(wěn)健*(2比較不同Tobin說明：Statavce(以便確定是否支持vce(bs)選項。常用命令包括：probitlogitcloglitglmivregxtregreg*****Printedon2011-3-5-Printedon2011-3-5*BS*seeDekking(2005,p.350);Efron(1993,Chp12-****************xN(muse^2)x--x-其Prob{z(0.025)<=<=z(1-0.025)}=1-0.05x-Prob{-1.96<=<=1.96}=1-0.05=do************************Prob{mu_in_[x-z(0.975)*se,x-z(0.05)*se]}=因此，mu95[x-z(0.975)*se,x-z(0.025)*sez(0.975z(1-0.975-x+(-)[x-1.96*se,x+1.96*se含義：在95%的情況下，上述區(qū)間都會包含x的真實值muCI_a=[x-z(1-a)*se,x-z(a)*se其中，xa表示置信水平，如z(a表示標準正態(tài)分布在第100*a百分位上的數(shù)值，如z(0.025)=-tt_k(a)，k表示tse計算過程解析：cisysuseauto,clearciprice,level(90)/*ci:*cilocalk=r(N)-localt_a=invttail(`k',0.05)t(k分布的第5*/t(73)分布的5%臨界值*/dis`t_a'*pricesumlocalse_price=sqrt(r(Var)/r(N))dis"se(price)="`se_price'*計算置信區(qū)locallb=localub=disbound-+`t_a'*`se_price'"`lb'"disbound*ciciprice,Printedon2011-3-5-Printedon2011-3-5*基于Bootstrap*see[U]bootstrapp.217-*Bootstrapx采用se_bs代替(1)式中的se，即可得到置信區(qū)間，即CI_a_bsx-z(1-a)*se_bsx-z(a)*se_bs(2)*******Normal-based[90ConfInterval]sysuseauto,clearbootstrapmean_bs=r(mean),saving(bs_ci01,replace)quisumlocalmean=r(mean)quiusebs_ci01,clearquisummean_bs/*均值采用price的均值localse_bs=/*標準誤由BS平均值的標準差計算得到localz_95invnorm(0.95)/*正態(tài)分布第95百分位數(shù)值*/localz_5=invnorm(0.05)/*正態(tài)分布第5百分位數(shù)值locallb`mean`z_95'*`se_bs'/*下限*/localub`mean`z_5*`se_bs'/*上限*/dis"90CI[`lb`ub']"*基于BS*see[U]bootstrapp.217-setobssetseedgenz=invnorm(uniform())histogramz,normal*******z90=[p5,(p5表示第5百分位上的數(shù)值本例 CI90=[-1.337876,z,就是x由小到大排序后的第5個觀察值(或第5個與第6個觀察值的平均值就是x由小到大排序后的第95個觀察值(或第95個與第96個觀察值的平均值sortdisdisp5*/p95**quisumz,detaillocalp5=r(p5)localp95=r(p95)histogramz,text(0.44`=`p5'+0.6'"")text(0.440.2"90%text(0.44`=`p95'-0.6' sysuseauto,bsmeanr(meanreps(1000nodotssaving(bs_mp,replacelevel(90)sumpriceestatbootstrap,all/*列出所有三種置信區(qū)間*/*注意：(P)**usebs_mp.dtaclear/*bs_mp.dta中存儲了1000個BSmean的觀察值*/quisummean,detailPrintedon2011-3-5-Printedon2011-3-5locallb_bs=r(p5)localub_bs=r(p95)dising"90%percentileCI=["iny%8.3f`lb_bs'","http:///iny%8.3f`ub_bs'ing"]"summean,*sysuseauto,clearsumpricecipricelevel(90)/*傳統(tǒng)方法*/locallb=r(lb)localub=usebs_mp.dta,twowayhistogrammean,xline(`lb'xline(`lb_bs'`ub_bs',lcolor(blue))caption("red:cibound""blue:Bootstrap*說明：紅線表示采用傳統(tǒng)的t分布得到的10%Bootstrap10置信區(qū)間；****例2exp(meansetobsgenx=invnormal(uniform())setseed123456bstheta=exp(r(mean)),reps(1000)saving(bs_n01,replace)level(90):estatbootstrap,allxusebs_n01,clearsumtheta,detailhistogramtheta由于theta的分布具有向右拖尾的特征采用傳統(tǒng)的ci方法得到的CI可能過citheta,即基于BS得到的百分位數(shù)值確定置信區(qū)間***********model：y=B0*(1-e^(-adopath+D:\stata9\ado\personal\Net_course\B9_BS_MCviewsourcenlexpgr.adouseB9_production.dta,cleargenx=capital*NLSnlexpgr采用BS獲得標準誤（僅需估計模型系數(shù)標準誤時推薦采用該方法nlexpgrlnout,vce(bs,reps(200)setseedbootstrapb_x=_b[B1],reps(200):nlexpgr**Printedon2011-3-5-Printedon2011-3-5*B1bootstraplog_b_x=ln(_b[B1]),reps(200):nlexpgr*--*********y_i=x_i*b1+e_i(group1:i=1,2,...,N1)y_i=x_i*b2+e_i(group2:b2Ho:Bootstrap在Ho將N1和N2個樣本混合后得到的估計系數(shù)b與*****************步驟：參見N1+N2的“混合樣本”；將前N1group1,Diff(bs_j)=b1-b2,計算“實證P值P_bs=#{Diff(bs_j)>=Diff(0)}/2模型：price=a0+a1*weight+a2*length+a3*mpg+問題：在國產(chǎn)車和進口車兩個子樣本中，汽車重量(weight)對價格影響相同嗎***sysuseauto,regpriceweightmpgifforeign==1/*進口車localb1=regpriceweightlengthmpgifforeign==0/*國產(chǎn)車*/localb2=_b[weight]scalardiff0=`b1'-dis"diff0_weight="*BootstrapHob1localreps=matDJ(`reps'3/*存儲結(jié)果的矩陣*/forvaluesj=1/`reps'{quisysuseauto,clearquiregpriceweilenmpginlocalb1=_b[weight]quiregpriceweilenmpginlocalb2=_b[weight]localdiff=`b1'- /*前22個觀察值，視為進口車23/74/*后52個觀察值，視為國產(chǎn)車matD[`j',1]=(`b1',`b2',}svmatD,names(d)renamed1b1renamed2b2renamed3diff*計算“實證P值countifdiff>diff0|diff==diff0localp=`r(N)'/`reps'dis"實證P值Printedon2011-3-5-Printedon2011-3-5*localdiff0=histogramdiff,xline(`diff0',lw(thick))localdiff0=scalar(diff0)kdensitydiff,lw(thick)Boostrapb1b2sumb1*sysuseauto,regpriceweightlengthmpgifforeign==1,noheader/*進口車*/**-擴展I：如何計算模型中所有系數(shù)的“實證P值 sysuseauto,clearregpriceweilenmpgeretlistmatlistsysuseautoclearquiregpriceweightmatb1=quiregpriceweightmatb2=e(b)*/*進口車/*國產(chǎn)車matD0=b1-*采用Bootstrap得到實證差異和實證Plocalreps=matDJ(`reps'4/*存儲結(jié)果的矩陣*/forvaluesj=1/`reps'{quisysuseauto,clearquiregpriceweilenmpgin /*前22個觀察值，視為進口車matrixb1=quiregpriceweilenmpgmatrixb2=e(b)matrixdiff=b1-b2matD[`j',1]=diffin23/74/*后52個觀察值，視為國產(chǎn)車}matJ(4,2,./*記錄系數(shù)真實差異和實證P值的矩陣forvaluesj=localdiff0_`j'=quicountif(diff`j'>`diff0_`j''|diff`j'==`diff0_`j'')localp=`r(N)'/1000matP[`j',1]=(`diff0_`j'',}matcolnamesP=系數(shù)真實差 實證PmatrownamesP=weightlengthmatlistPmpg**圖示：mpg變量的BS系數(shù)差異的分histogramdiff3,xline(-161.1735,*擴展II：如何得到面板數(shù)據(jù)模型的“實證P值解決思路：按個體進行抽樣(每個面板中的截面都是一個個體see[XT]p.215TechnicalPrintedon2011-3-5-Printedon2011-3-5例：成長機會(Tobin)對資本結(jié)構(gòu)的影**Ho:b1=(b_large=usextcs.dtaclear/*中國上市公司資本結(jié)構(gòu)數(shù)據(jù)*/tssetcodeyearegensize_meanmean(sizeby(code/*保證每家公司都分在同一組listcodeyearsizesize_meaninsortgengroup=group(3)tabgroupdropifgroup==2xtregtlsizeeststoresmallxtregtlsizendts/*去掉中間組，保證最大和最小兩個組的差異足夠顯著tangnprtobinifgroup==1,fe/*小規(guī)模公司tangnprtobinifgroup==3,fe/*大規(guī)模公司eststoreesttabsmalllarge,mtitle(small可見，在不同規(guī)模的公司中，Tobin對TL的影響是不同的下面采用BS進行檢驗*tangnprtobinifgroup==1,fe/*小規(guī)模公司matb1=xtregtlsizendtstangnprtobinifgroup==3,fe/*大規(guī)模公司mat=matD0=b1-matD0,title(系數(shù)估計值的真實差異*采用Bootstrap得到實證差異和實證P*usextcs.dta,clearsetseed12345bsample,cluster(code)idcluster(code_bs)tssetcode_bsyearlistcode_bscodeyearin1/35if*可見：公司000002000012都被抽中了2*code_bslocalreps=localvarlist"tlsizendtstangnpr/*存儲結(jié)果的矩陣matD=J(`reps',6,.)forvaluesj=1/`reps'{usextcs.dta,cleargengg=quitssetcode_bs/*將樣本等分為三組quixtreg`varlist'matrixb1=e(b)quixtreg`varlist'gg==1,fe/*第一組，視為小規(guī)模公司gg==3,fe/*第三組，視為大規(guī)模公司matrixb3=matrixdiff=b1-b3matD[`j',1]=diff}svmatD,*diff1-diff6sizendtstangnprtobinmatPJ(6,2,.)/*記錄系數(shù)真實差異和實證P值的矩陣*/forvaluesj=1/6{localdiff0_`j'=quicountif(diff`j'>`diff0_`j''|diff`j'==`diff0_`j'')&diff`j'!=localp=`r(N)'/`reps'matP[`j',1]=(`diff0_`j'',}matcolnamesP系數(shù)真實差異實證PmatrownamesP=sizendtstangnprtobin.Printedon2011-3-5-Printedon2011-3-5matlist*結(jié)論：除SIZE外，其它變量的組間系數(shù)估計值并不存在顯著差異*Reps5000次的結(jié)果(大約需要20分鐘 實證P* * *PermutationTests（組合檢驗*Efron(1993,*rndallo.adoZ(z1z2z_n),分布為Y(y1y2 分布為****************Ho：F=(例如：mean(ZZYn+mH；從n+m個樣本中隨機抽取n個觀察值(不重復(fù))，計算均值計算d1(m1將上述三個步驟重復(fù)B=1000次，得到d_i>d_0的次數(shù)，假設(shè)為Q次，則permutetest的“p值”為：p-value=Q/B其中，d_0例1：檢驗兩組數(shù)字的均值是否相*做了兩組實驗，每組三個樣本，考察treatment是否有助于提高老鼠的存活時inputytreatment 11101214saveB9_permu.dta,sumbysorttreatment:sumsetseedpermuteymean=r(meanrightnodropnowarnsumyiftreatment/*默認100次*/permuteymean=r(mean),rightreps(20)nodropnowarn:sumyiftreatment雖然精確的概率是0.1，但當Permute傳統(tǒng)tttesty,PermutestuseB9_permu.dta,clearbysorttreatment:sum*d_0=-3=9Printedon2011-3-5-Printedon2011-3-5useB9_permu.dta,genid_n/*便于后續(xù)合并數(shù)據(jù)*/sortidsavedata_all,*-Permunationsetmatsize8000localB=6000matD=J(`B',1,.)forvaluesi=useB9_permu.dta,genrdm=uniform()sortrdmsumyin1/3localm1=r(mean)sumyin4/6localm2=r(mean)locald=`m1'-`m2'matD[`i',1]=/*隨機組合/*前三個觀察值/*后三個觀察值}}D,if(diff<-3|diff==-3)Q=r(N)p_value=`Q'/dis"p-value="結(jié)論：二者的差異在10%(5%圖示diff1的分布情histogramdiff1,xline(-3,lw(thick))****(((N)n (其中，m=2n-**其中，N為總樣本數(shù)，n為第一組的樣本數(shù)，N-n為第二組的樣本matA=localj=foreachNofnumlist46810203050100localnint(`N'/2)/*23,4看看結(jié)果如何？*/localpermute=comb(`N',`n')matA[`j',1]=(`N',`n',`permute')localj=`j'+1}matcolnamesA=NnPermutematlistA,format(%20.0f)*結(jié)論：對于N>10的樣本，采用permutetest*而對于N<10*******拒絕標準(Efron1993可以拒絕（borderlingASL<ASL<ASL<ASL（reasonablystrong（veryevidence）拒絕原假設(shè)*Permutation（Efron,p.211,Tab****Printedon2011-3-5-Printedon2011-3-5*1000PermutationTest與BootstrapTests“同”，二者的基本思想非常相似在原假設(shè)Ho:m1=m2“異”，二者的抽樣方式不同Bootstrap為可重復(fù)抽樣(samplingwithreplacement)Permutationtest為不可重復(fù)抽樣(samplingwithoutreplacement)******例2：進口車和國產(chǎn)車的“平均價格”相等嗎capprogramdropmy_diffprogramdefinemy_diff,rclassargsvgvquisum`v'if`gv'localmean1=r(mean)quisum`v'if!`gv'localmean2=r(mean)retscalardiff=`mean1'-sysusepermutepricediff=r(diff):permutepricediff=my_diffpricereps(2000)nodotsnowarn:my_diffpricepermutepricemean=*tttestprice,reps(1000)right*例3：進口車和國產(chǎn)車“價格的波動”(方差)相等嗎capprogramdropmy_diff2programdefinemy_diff2,rclassargsvgvquisum`v'if`gv'localvar1=r(Var)quisum`v'if!`gv'localvar2=retscalardiff_var=ln(`var1'/`var2')*sysusepermutepriced_var=r(diff_var),reps(5000)nodotsnowarn:my_diff2pricerperm,permtest1,**面板數(shù)據(jù)的Permutation*listcodegenrr=bysortcode:egenrm=mean(rr)tssetrmyear*xtnewid/*連玉君自行編寫的命令，產(chǎn)生連續(xù)編號的公司代碼egencode_new=Printedon2011-3-5-Printedon2011-3-5listcodeyear多次運行上述代碼，理解“隨機排序”的含接下來，我們就可以根據(jù)code_new**采用rpermusextcs_simple.dta,clearlistcodeyear*finditrpermcode,gen(code_new)cluster(code)listcodeyearcode_new/*重新組合排序tssetcode_new*xtnewidegenid=group(code_new)listcodeyearcode_newid*至此，我們可以根據(jù)id*例：成長機會(Tobin)對資本結(jié)構(gòu)的影****Ho:b1=(b_large=注意：與此前Bootstrap檢驗的差異是什么usextcs.dta,bysortcode:egenmsize=mean(size)sortmsizegengsize=group(3)dropifgsize==2*sizendtstangnprtobinifgsize==1,fe/*小規(guī)模公司=sizendtstangnprtobinifgsize==3,fe/*大規(guī)模公司=e(b)b1-D0,title(系數(shù)估計值的真實差異matxtregmatb3matD0=matlocalreps=/*模擬次數(shù)localvarlist"tlsizendtstangnprmatDJ(`reps'6/*存儲結(jié)果的矩陣*/forvaluesj=1/`reps'{usextcs.dta,cleargenrr=uniform()bysortcode:egenrm=mean(rr)quitssetrmyearegencode_new=/*隨機排序quitssetcode_new*dislocalg3=*disquixtreg`varlist'matrixb1=e(b)quixtreg`varlist'/*第三組的分界點，第292家公司以后的公司code_new<`g1',fe/*第一組，視為小規(guī)模公司code_new>`g3',fe/*第三組，視為大規(guī)模公司matrixb3=matrixdiff=b1-b3matD[`j',1]=diff}svmatD,*diff1-diff6sizendtstangnprtobinmatPJ(6,2,.)/*記錄系數(shù)真實差異和實證P值的矩陣*/forvaluesj=1/6{localdiff0_`j'=quicountif(diff`j'>`diff0_`j''|diff`j'==`diff0_`j'')&diff`j'!=.localp=`r(N)'/`reps'matP[`j',1]=(`diff0_`j'',Printedon2011-3-5-Printedon2011-3-5}matcolnamesP系數(shù)真實差異實證PmatrownamesP=sizendtstangnprtobinconsmatlistP*模擬1000* 實證P* .11700291-.04144458-.08150638-.01324722*-.62137931*Jackknife法(刀切法*************簡在四十年代末，五十年代初提出，見多用于勘查“離群值是因為采用Bootstrap計算標準誤和置信區(qū)間更為有效和穩(wěn)健*************給定樣本(x1,x2,...,x100)，記為x_i只計算剩下99個觀察值的均值，得到100個均值，100個均值，便可計算其標準誤和置信區(qū)間了。m1=( m2= ...m100=(x1+x2+x3+...+x99)***偽數(shù)值*XM_2=* (去除第二個觀察值后的樣本均值Printed2011-3-5則-Printed2011-3-5則************************和XM_2(N-1)XM_2+x2XM_a=N((N-1)XM_j+x_jXM_a=N由[2x_j=N*XM_a-(N-(x_j表示第該值在Jackknife檢驗中稱為“偽數(shù)值設(shè)Q為任意統(tǒng)計量(如標準差、方差等Q_j為刪除第j則第j個觀察值對應(yīng)的“偽數(shù)值”為q_j=N*Q-(N-**********Q的一階無偏估計量；的標準誤是統(tǒng)計量Q的標準誤；={{{1Nq_seSUM(q_j-q_mean)^2N(N-1)}例1：解析Jackknife的計算過程setobsgenid_n/*樣本序號*/genx=_n+10genpx /*x的偽值forvaluesj=quisumQlocalm_0=sumxif(id!=`j')Q_jlocalm_j=r(mean)replacepx=10*`m_0'-in/*計算偽值}sumxStatajackknifem_jk=r(mean),rclass計算偽sumreplacepx=forvaluesj=1/10{quisumxlocalm_0=sumxif(id!=`j')localm_j=r(Var)replacepx=10*`m_0'-(10-}sumx**statain/*計算偽值replacem_jk=jackknifem_jk=r(Var),rclasskeep:sumPrintedon2011-3-5-Printedon2011-3-5例2：JK得到的標準誤和置信區(qū)間與ci命令的結(jié)果相同(傳統(tǒng)方法sysuseauto,jackknifer(mean),doublesumpriceciprice/*參見Bootstrap小節(jié)*/*這里，double用于設(shè)定數(shù)值的精度，即“雙精度例3：Jackknife計算標準jackknifesd=r(sd),rclasskeep:sum*x的“標準差”為*x的“標準差”的95置信區(qū)間為[-x的“標準差”的JK標準誤為我們發(fā)現(xiàn)，最后一個觀察值對應(yīng)的“偽值”很大，明顯不同于其他偽值至于11個觀察值是否為離群值，則視具體情況而定本例中，這11個觀察值來自于指數(shù)分布，因此它并非離群值histogram例4：離群值的查sysuseauto,jackknifesd=r(sd),rclasskeep:sumsumsd,listpricemakesdif例5：同時計算多個統(tǒng)計量的標準sysuseauto,jackknifesd=r(sd)skew=r(skewness)var=r(Var),rclass:例6：計算估計系數(shù)的標準*-截面資sysuseauto,jackknife:regpriceweilenmpg*regpriceweileneststorer_jkregpriceweileneststorer_bsregpriceweiforeign,foreign,vce(bootstrap,eststorelocalmm"r_olsr_bsr_jk"esttab`mm',mtitle(`mm')*結(jié)論：相對于傳統(tǒng)的t檢驗，BS和JK*同時，BS和JK得到的標準誤非常接近*-面板數(shù)usextcs.dta,clearkeepifcode<1000Printedon2011-3-5-Printedon2011-3-5jackknife,eclasscluster(code)idcluster(code_new):///xtregtlsizendtstangtobin,festorextregtlsizendtstangtobin,storelocalmm"fe_jkesttab`mm',mtitle(`mm')/*標準誤esttab`mm',mtitle(`mm')/*t值*結(jié)論：相對于傳統(tǒng)的t檢驗，JK*說明：類似于BS，在采用JK*cluster(idcluster()例7：Jackknifeinstrumentalvariables“刀切法”IV對該方法的介紹和相關(guān)的MonteCarlo過程，參見usehsng2.dta,jiverentpcturban(hsngval=famincreg2-ivregrentpcturban(hsngval=famincreg2-Journal2006(3):304-*********JackknifeBootstrapseeEfron(1993,Chp11,雖然從表面上看，Jackknife似乎只利用了非常有限的樣本信息換言之，Jack的準確程度決定于統(tǒng)計量與其線性展開的接近程度*例：比較BS和JK*參見setobslocalB=/*模擬次數(shù)matR /*存儲模擬結(jié)果的矩陣localstatis/*統(tǒng)計量*localstatis"r(mean)^2*統(tǒng)計量2：[mean(x)]^2*/forvaluesi=1/`B'{tempvargen`x'=1+invnorm(uniform())/*x-tempfilequibsm_bs=r(mean),reps(20)nowarnsaving("`bsdata'"):quiuse"`bsdata'",clearquisummatR[`i',1]=`statis'capdropm_jkquijackknifem_jk=r(mean),rclasskeep:sum`x'quisumm_jkmatR[`i',2]=sum}svmatR,names(r)renamer1r_bsrenamer2r_jksumr*graphboxr_bssetobsPrintedon2011-3-5-Printedon2011-3-5genx=.genmx=*localjk"in2/100"forvaluesi=1/100{quireplacex=1+invnorm(uniform())quisumx`jk'quireplacemx=r(mean)in}summx****Jackknife不適用的情seeEfron(1993,原始樣本具有“平滑性即，數(shù)據(jù)的微小變動(如，刪除一個觀察值)，僅會導(dǎo)致統(tǒng)計量的微小變動在滿足這一假設(shè)的前提下，JK結(jié)果會非常接近于BS反之，JK的結(jié)果將是有偏的(biased)*-例1 JK中位 *==問題：變量xuseB9_jk_fail.dta,clear*s.e.usingjackknifemed=r(p50),rclasskeep:sums.e.usingbootstrapmed=r(p50),reps(100)nodotsnowarn:sumx,**-我們發(fā)現(xiàn)，采用JK得到的標準誤明顯偏小JK**采用JK *在一個邊長為r的正方形內(nèi)均勻投點********area= ==>pi=|||4Printedon2011-3-5-Printedon2011-3-5**doviewsource/***實例：rcapprogramdroppiprogramdefinepi,rclassversion9.2scalarx=uniform()scalary=scalard=sqrt((x-0.5)^2+(y-0.5)^2)locals=d<=0.5retscalars=`s'simulates=r(s),reps(100000):piscalarpi=dising"pi="iny*例2*genx=1+(3-=sumlocalI=1/(3-1)*r(mean)dis"I=`I'"setobs10000genR=uniform()genX=-ln(R)genF=X^(1.9-1)sum*setobsgenR1=uniform()genR2=uniform()genX=-ln(R1)*R2genF=X^(0.1)sumFMC和BS換言之，MCBS則無需事先給定分布函數(shù)*****“對數(shù)正態(tài)分布”的期望和方capprogramdroplnsimprogramdefinelnsim,rclassversion9.2drop_allsetobs100genz=exp(invnormal(uniform()))/*ln(z)--N(0,1)*/sumzreturnscalarmean=r(mean)returnscalarVar=r(Var)Printedon2011-3-5-Printedon2011-3-5simulatemean=r(mean)Var=r(Var),reps(1000):lnsim僅執(zhí)行一次的效果MC則是將該過程重復(fù)多次retsetseed1234ret**一個相對完整的設(shè)定：改變“對數(shù)-正態(tài)分布”的參capprogramdroplnsim2programdefinelnsim2,rclassversionsyntax[,obs(integer100)mu(real0)sigma(real1)]drop_allsetobs`obs'tempvarzgen`z'=exp(`mu'summarize`z'returnscalarmeanreturnscalarVar+==simulatemean=r(mean)var=r(Var),reps(1000):lnsim,obs(100)simulatemean=r(mean)var=r(Var),reps(1000):lnsim,obs(50)mu(-3)*例2：內(nèi)生性問題和異方差對OLS**Insomecasesitispossibletocalculatethesampling*fromtheeconometricmodel.Butsometimes,especiallyforfinite*samples,thisiseithernotpossibleorverydifficult.Inthesecases*Carloexperimentsareanintuitivewaytoobtaininformationabout*samplingdistributionandhenceaboutthe“quality”ofthe*模型：y_i12*x_i*u_i=z_i+*顯然：Var[u_i*Corr(u_i,x_i)*-數(shù)據(jù)生成drop_allgenx=invnormal(uniform())geny_true=1+2*xsavetrue_data,*-capprogramdrophetero1programdefinehetero1version9.2argsusetrue_data,genz=invnormal(uniform())geny=y_true+z+`c'*xregyxsimulate_b_se,reps(10000):hetero13*事實上，_b_x=5=*capprogramdrophetero2programdefinehetero2Printedon2011-3-5-Printedon2011-3-5versionsyntaxvarname[,c(real3)]tempvarygen`y'=y_true+(invnormal(uniform())+`c'*`varlist')reg`y'xusetrue_data,simulate_b_se,reps(1000)nodots:usetrue_data,simulate_b_se,reps(1000)nodots:*usetrue_data,genz=invnormal(uniform())simulate_b_se,reps(1000)nodots:*結(jié)論：異方差并不會影響x系數(shù)OLS(1(2*例3：遺漏變量和增加無關(guān)變量對OLS估計的影對于模型ya0a1*x_1a2*x_2u--x_1a_2OLS估計將是有偏的；x_3a_1a_2y=0.5+x_1+2*x_2setobsgenx1=invnormal(uniform())genx2=invnormal(uniform())genx3=invnormal(uniform())geny=0.5+1*x1+2*x2savemyomit_data,replace*******capprogramdropmyomitprogramdefinemyomit,eclassversion9.2tempvarugen`u'=gen`yy`u'/*讓干擾項變動是模擬的基礎(chǔ)*/reg`y'`varlist'*simulate_b_se,*myomitx1*usemyomit_data,simulate_b_se,usemyomit_data,simulate_b_se,*x1x2::結(jié)論***Printedon2011-3-5-Printedon2011-3-5*例4：采用MC比較兩個正態(tài)分布序列中位數(shù)的比******y--N(mu,{N(mu1,sigma_1^2)y--{{N(mu2,n1n2*capprogramdropmyroatioprogramdefinemyratio,rclassversion9.2syntax[,n1(int5)mu1(int3)sigma1(int1)n2(int10)mu2(int3)sigma2(intdroplocalN=`n1'+setobsgen`y'=`y'=`mu1'+`y'*`sigma1'if_n<=`n1'/*y1--`y'=`mu2'+`y'*`sigma2'if_n> /*y2--*replace`y'=cond(_n<=`n1',`mu1'+`y'*`sigma1',`mu2'+`