对付界截面数据的异方差的stata命令是什么啊

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

<P>如题 </P>
<P>谢谢你</P>

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：stata命令 Stata 截面数据 tata 异方差

white estimator

reg y x, robust

for others, see

http://jinhe.xjtu.edu.cn/bbs/dispbbs.asp?BoardID=5&id=422&replyID=1654&star=2&skin=

or

http://arlion.8j.cn
这里是我以前做得一份答案，我做了几处修改，所以可能会出现错误。

*--------------------------------- version 1 -----------------------------

//Greene chp12 example12.7 pp515
* Given by arlion, jinhe center, xi'an Jiaotong University
* 10.8 2005

clear
cap log close
cap cmdlog close
log using Greene_exp12_7,text replace
cmdlog using Greene_exp12_7_cmd, replace
use "D:\stata8\ado\Examples\others\Greene4th\TBL5-1.DTA", clear

//=======================================================
//=======================Example12.7=====================
//=======================================================

label var x3 "Avg. monthly credit card expenditure"
label var x1 "Age in years + 12ths of a year"
label var x2 "Income, divided by 10,000"
label var x4 "OwnRent, individual owns (1) or rents (0) home"
label var x5 "Self employed (1=yes, 0=no)"
label var y1 "Number of derogatory reports"
label var y2 "Credit card application accepted (1=yes)"
rename x1 age
rename x2 income
rename x3 expenditure
rename x4 ownrent
rename x5 self_emp
rename y1 badreport
rename y2 card_apply
gen income2 = income^2
gen cons = 1
keep if exp>0 /* Now,we have only 72 obs */

mkmat exp,mat(y)
mkmat cons age ownrent income income2,mat(x)

// Row1 : sigma2_i = sigma2
// Method: OLS
mat beta_hat = inv(x'*x)*(x'*y)
mat beta_hat = beta_hat'
mat list beta_hat

// Row2 : sigma2_i = sigma2*income
// Method: WLS
gen weight2 = income
mkmat weight2,mat(weight)
matrix omiga = diag(weight)
mat beta_hat = inv(x'*inv(omiga)*x)*(x'*inv(omiga)*y)
mat beta_hat = beta_hat'
mat list beta_hat

// Row2-1 : sigma2_i = sigma2*Income
// e2 on Zi = (Income),noconstant
// Purpose: using FGLS to test Row2
// Method : FGLS
mat b1 = inv(x'*x)*(x'*y)
mat b0 = b1*10 /*initial*/
mkmat income,mat(z)
local dif = mreldif(b1,b0)
local i = 1

while `dif'>e^(-6){
mat e = y - x*b`i'
mat e2 = hadamard(e,e)
mat beta_z= inv(z'*z)*(z'*e2)
mat omiga = diag(z*beta_z)
local j= `i' + 1
mat b`j' = inv(x'*inv(omiga)*x)*(x'*inv(omiga)*y)
local dif = abs(mreldif(b`j', b`i'))

mat c`j' = b`j''
dis _n in g " =====itering times: " in y `i' in g "====== "
mat list c`j'
local i = `i'+1
}

// Row3 : sigma2_i = sigma2*income2
// Method: WLS
gen weight3 = income2
mkmat weight3,mat(weight)
matrix omiga = diag(weight)
mat beta_hat = inv(x'*inv(omiga)*x)*(x'*inv(omiga)*y)
mat beta_hat = beta_hat'
mat list beta_hat

// Row4 : sigma2_i = sigma2*(1,logIncome)
// e2 on Zi = (1,logIncome)
// Method: FGLS
gen logincome = ln(income)
mat b1 = inv(x'*x)*(x'*y)
mkmat cons logincome,mat(z)

forvalues i=1(1)20{
mat e = y - x*b`i'
mat e2 = hadamard(e,e)
mat beta_z= inv(z'*z)*(z'*e2)
mat omiga = diag(z*beta_z)
local j = `i' + 1
mat b`j' = inv(x'*inv(omiga)*x)*(x'*inv(omiga)*y)
mat c`j' = b`j''
dis _n in g " =====itering times: " in y `i' in g "====== "
mat list c`j'
}

// Row5 : sigma2_i = sigma2*(1,logIncome)^2
// |e| on Zi = (1,logIncome)
// Method: FGLS
cap gen logincome = ln(income)
mat b1 = inv(x'*x)*(x'*y)
mat b0 = b1*10
mkmat cons logincome,mat(z)
local dif = mreldif(b1,b0)
local i = 1
while `dif'>e^(-4){
mat e = y - x*b`i'
svmat e , names(res)
replace res = abs(res)
mkmat res,mat(e)
mat beta_z= inv(z'*z)*(z'*e)
mat omiga = diag(z*beta_z)
local j= `i' + 1
mat b`j' = inv(x'*inv(omiga)*x)*(x'*inv(omiga)*y)
local dif = mreldif(b`j', b`i')
mat c`j' = b`j''
dis _n in g " =====itering times: " in y `i' in g "====== "
mat list c`j'
drop res
local i = `i'+1
}

// Row6 : sigma2_i = sigma2*exp(z'*alpha)
// log_e2 on Zi = (1,logIncome)
// Method : FGLS

mat b1 = inv(x'*x)*(x'*y)
mat b0 = b1*10 /* initializtion */
local dif = mreldif(b1,b0)

cap gen logincome = ln(income)
mkmat cons logincome,mat(z)

local i = 1

while `dif'>e^(-2.3){
mat e = y - x*b`i'
mat e2 = hadamard(e,e)
svmat e2,names(res2)
replace res2 = ln(res2)
mkmat res2,mat(ln_e2)

mat beta_z= inv(z'*z)*(z'*ln_e2)
mat log_e2_hat = z*beta_z
svmat log_e2_hat,names(log_e2_hat)
replace log_e2_hat = exp(log_e2_hat)
mkmat log_e2_hat,mat(e2_hat)
mat omiga = diag(e2_hat)
local j= `i' + 1
mat b`j' = inv(x'*inv(omiga)*x)*(x'*inv(omiga)*y)
local dif = abs(mreldif(b`j', b`i'))
mat c`j' = b`j''
dis _n in g " =====itering times: " in y `i' in g "====== "
mat list c`j'
local i = `i'+1
drop res2 log_e2_hat
}
dis _n _s(6) in g "alpha = " in y beta_z[2,1]
/* we can get alpha value in matirx beta_z*/

// Row8 : sigma2_i = sigma2*(x'b)
// Method: WLS using (x'b)^2
mat b = inv(x'*x)*(x'*y)
mat y_hat = x*b
qui reg exp age own income income2
qui predict yhat
cap gen yhat2 = yhat^2
mkmat yhat2,mat(weight)
matrix omiga = diag(weight)
mat beta_hat = inv(x'*inv(omiga)*x)*(x'*inv(omiga)*y)
mat beta_hat = beta_hat'
mat list beta_hat

* --------------------------------------------------------------------------------

* -----------------------------------version 2 采用一个下载的命令---------------------------------------------

* Lecture 4-1: 2 Greene chp12 exercises
* GLS estimation for Heteroscedasticity
* Given by Arlion, Jinhe Center, Xi'an Jiaotong University.

* see Chapter 8 of Chung's Lecture; Chapter 12 of Greene(2000)
** The main reference is Chapter 4 of Arlion's notes,
** which can be download at http://arlion.8j.cn
* Suggestions are welcome, E-mail: arlionn@163.com

* Greene chp12 exercises
* Including eg-12.1, eg-12.4, eg-12.5, eg-12.7

cap log close
cap cmdlog close
log using Greene12,text replace
cmdlog using Greene12_cmd, replace
use "D:\stata8\ado\Examples\others\Greene4th\TBL5-1.DTA", clear

//=======================================================
//==================------Example12.1------pp.500========
//=======================================================

label var x3 "Avg. monthly credit card expenditure"
label var x1 "Age in years + 12ths of a year"
label var x2 "Income, divided by 10,000"
label var x4 "OwnRent, individual owns (1) or rents (0) home"
label var x5 "Self employed (1=yes, 0=no)"
label var y1 "Number of derogatory reports"
label var y2 "Credit card application accepted (1=yes)"
rename x1 age
rename x2 income
rename x3 expenditure
rename x4 ownrent
rename x5 self_emp
rename y1 badreport
rename y2 card_apply
gen income2 = income^2

reg expenditure age ownrent income income2 if exp>0
* Testing for HSK
* Method1: plotting
predict e,resid
gen e2 = e^2

rvpplot income , ytitle("Residuals") xtitle("Income") saving(chp4_income.gph,replace)
rvpplot age , ytitle("Residuals") xtitle("Age") saving(chp4_age.gph,replace)
graph combine chp4_income.gph chp4_age.gph

/*
twoway scatter e income, ///
title("Fig 12.1_1: for Greene chp12 exp12.1") ///
note("Note: we can see the residual is heteroscedastic with income")
twoway scatter e2 income if e2<1000000, ///
title("Fig 12.1_2: plot between e2 to income") ///
note("Note: we can see here the heteroscedastic is rather obvious")
twoway scatter e age, ///
titl("Fig12.1_3: plot between e to age") ///
note("Note: we can see that no HTS between age to e")
*/

*=======================================================
*==================------Example12.5------ PP.500-511 ==
*=======================================================

* ---------------------------------------------
* Method1: =========Goldfeld-Quandt test==========
* ---------------------------------------------
sort income
gen obs = _n
qui reg expend age ownrent income income2 if obs<50&exp>0
scalar sse1 = e(rss)
local N1 = e(N)
qui reg expend age ownrent income income2 if obs>50&exp>0
scalar sse2 = e(rss)
local N2 = e(N)
scalar k = e(df_m)
local n1 = `N1'-k
local n2 = `N2'-k
local F = sse2/sse1 /*the lager one is on the upper side*/
local p = Ftail(`n1',`n2',`F')

* Display the results of test
di _n in g "Goldfeld-Quandt test for heteroskedasticity : " ///
_n _s(10) in g "H0: Constant variance H1: otherwise" ///
_n _n _s(10) in g "F("in y %2.0f `n1' in g "," in y `n2' in g ")= " ///
in y %6.4f `F' ///
_n _s(10) in g "P-value = "in y %7.4f `p'

* 采用Arlion编写的命令——gqhet.ado
gqhet expend age ownrent income income2, sortvar(income)

* ---------------------------------------------
* Method2: =========Breush-Pagan Test==========
* ---------------------------------------------

* Using the commmand given by stata [see] help for hettest
qui reg exp age ownrent income income2 if exp>0
hettest income
hettest income income2
/*
rvfplot , title("Fig12.2 Breush-Pagan test") ///
note(" Note: variables using to test is : income income2")
*/
hettest age
hettest ownrent

* ===written by arlion===version1 Refer to Greene(2000, pp.509) Eq.(12-13)
quietly reg expend age ownrent income income2 if exp>0
tempvar ei touse
gen byte `touse' = e(sample)
quietly predict `ei' if `touse' ,resid
quietly replace `ei' = `ei'^2/(e(rss)/e(N))
qui reg `ei' income income2 if `touse'
local testvar "income income2"
scalar chi2 = 0.5*e(mss)
scalar df = e(df_m)
scalar p = chi2tail(e(df_m),e(mss)/2)
dis in g "The Breusch-Pagan test for heteoscedasticity: " _n ///
in g "The variables used in testing is : " in y "`testvar'" _n _n ///
in g _s(6) "Chi2(" in y (e(df_m)) in g ") = " in y %8.4f chi2 _n ///
in g _s(6) "P-value = " in y %8.4f p

* ===written by arlion===version2: A robust test of B-Pagan
* Refer to Greene pp510 Koenker(1981) suggest LM test
* Note: this test is robust for not assumping normality
qui reg exp age ownrent income income2 if exp>0
tempvar token
gen byte `token' = e(sample)
keep if `token'
local n = e(N)
cap drop e e2
qui predict e if `token', resid
gen e2=e^2
mkmat e,mat(e)
mkmat e2 , mat(e2)
matrix u = e2
mat sum_e2 = e'*e
scalar u_bar = sum_e2[1,1]/`n'
gen u_bar = u_bar
gen V_i = (e2 - u_bar)^2
qui sum V_i
scalar V = r(mean)
matrix u_bar = J(`n',1,u_bar) /*get matrix u_bar*i*/
gen cons = 1
mkmat cons income income2 if `token' , mat(z)
matrix LM = (1/V)*(u - u_bar)'*z*inv(z'*z)*z'*(u - u_bar)
scalar LM = LM[1,1]
local testvar "income income2"
local freedom=2
scalar p_value = chi2tail(`freedom',LM)
dis _n _n ///
in g "The Robust version of Breusch-Pagan test for heteoscedasticity: " _n ///
in g "Suggested by Koenker(1981) and Bassett(1982) Refer to Greene PP510" _n ///
in g _s(3) "Given by arlion" _n ///
in g "The variables used in testing is : " in y "`testvar'" _n _n ///
in g _s(6) "LM = Chi2(" in y `freedom' in g ") = " in y %8.4f LM _n ///
in g _s(6) " P-value = " in y %8.4f p

* Method3: =========White's Test==========

local i=1
local expvars " age ownrent income income2 "
foreach pp of local expvars{
foreach qq of local expvars{
gen `pp'_`qq' = `pp'*`qq'
}
}
drop income2_age income2_ownrent income2_income /// /* 3 */
income_age income_ownrent income_income /// /* 3 */
ownrent_age ownrent_ownrent /* 2 */
qui reg e2 age ownrent income income2 ///
age_* ownrent_* income_* income2_* if exp>0 /*5+4+2+1+1 */
scalar white = e(N)*e(r2)
scalar p = chi2tail(e(df_m), white)
dis _n _n ///
in g "The White's test for heteoscedasticity: " _n ///
in g "Refer to Greene PP508-510" _n ///
in g _s(3) "Given by arlion" _n _n ///
in g _s(6) " Chi2(" in y e(df_m) in g ") = " in y %8.4f white _n ///
in g _s(6) " P>chi2 = " in y %8.4f p

* 采用stata命令
qui reg expend age ownrent income income2
whitetst /*using STATA command : findit whitetst to download this command online.*/
imtest

*=======================================================
*============------Example12.4------===PP.506 ==========
*=======================================================
* ols
reg expend age ownrent income income2
est store ols

* D&M(1) (sigmai=n.ei^2/(n-K))
reg expend age ownrent income income2, robust
est store dm1

* white
matrix d = vecdiag(e(V))
matrix v = cholesky(diag(d))
matrix s = sqrt((72-5)/72)*vecdiag(v)
matrix V = diag(s)
matrix V = V*V'
matrix b = e(b)
ereturn post b V
ereturn display

* D&M(2) (sigmai=ei^2/m_ii)
reg expend age ownrent income income2, hc2
est store dm2

* D&M(3) (sigmai=ei^2/m_ii^2)
reg expend age ownrent income income2, hc3
est store dm3

*=======================================================
*============------Example12.7------===PP.515===========
*=======================================================

* --- Method 1 ---
* using option WLS：[aweight] 选项

reg expend age ownrent income income2 [aweight=1/income] /* 12.3a */
reg expend age ownrent income income2 [aweight=1/income2] /* 12.3b */

* FGLS: twostage regression
qui reg expend age ownrent income income2
cap drop e e2
predict e , res
gen ee = e^2
qui reg ee income income2, noconstant
qui predict p1
reg expend age ownrent income income2 [aweight=1/p1] /* 12.3c */

* sigma_i^2 = sigma^2exp(zi'a) zi = (1 ln_income)
gen e2 = ln(e^2)
qui reg e2 income
predict p2
qui replace p2 = exp(p2)
reg expend age ownrent income income2 [aweight=1/p2] /* 12.3d */

* --- Method 2 ---
* using wls0.ado command

* Table 12.3 in Greene(2000, pp.515)
* Note: the wls0 command should be download by typing: findit wls0

*===============WLS0================
wls0 exp age ownrent income income2 , wvar(income) type(abse) noconst /* 12.3a */
est store Tab12_3a
wls0 exp age ownrent income income2 , wvar(income2) type(abse) noconst /* 12.3b */
est store Tab12_3b
wls0 exp age ownrent income income2 , wvar(income income2) type(e2) noconst /* 12.3c */
est store Tab12_3c
wls0 exp age ownrent income income2 , wvar(income income2) type(abse) noconst /* 12.3d */
est store Tab12_3d
wls0 exp age ownrent income income2 , wvar(income income2) type(loge2) /* 12.3e */
est store Tab12_3e
wls0 exp age ownrent income income2 , wvar(income income2) type(xb2) /* 12.3h */
est store Tab12_3f

est table Tab*, stats(r2) b(%6.3f) se(%6.3f)

*========================OVER===========================

Wow; didn't expect it to be so complicated; thought would be just a command or two. Need time to digest. But thanks a looooooot!

BTW, can I just use WLS ---------- reg Y X (aweight *** )?

Maybe it is not accurate enough?

这是Greene的书上的例子，介绍如何和处理检验异方差。

你说的 reg Y X (aweight *** ) 也有例子在里面。

Greene的书估计大家都有，对着书看例子就可以了。

http://www.ats.ucla.edu/stat/stata/examples/greene/greene12.htm

112326.rar (21.85 KB) 本附件包括：

• Stata Econometric Analysis by William Greene, Chapter 12.mht

2007-4-30 13:22:00 上传

对付界截面数据的异方差的stata命令是什么啊

Whoa, I'm so reading these stuff now! Thx!

能不能提供Greene的第四版呢

我找了半天 都是第5版的

想看一下原题 以对照你的程序

 i want to know too