I decided to check out the code and it has many problems. First, tempvars
are generated and then replaced to equal a constant and then a scalar with
the same name (declared thru tempname) is created to hold the constant
(confusing at best). In many places, the code takes several steps to do
something that could be done much more simply and the potential presence of
missing values in the grouping variables is handled inconsistently.
But the question of negative chi sq values is most likely due to the fact
that the formula appears to be implemented incorrectly in that the number of
groups is used instead of the number of observations when creating the first
term in the LR statistic (I'm not sure about how the code should perhaps be
different for commands other than -reg-).
If you substitute the code below for everything inside the original quietly
block and change the tempvar and tempname declarations to
tempvar e xbgls nindg sigg arg perg
tempname sig sumg LR
then it may work. You can also drop the restore at the end of the program,
since Stata does that automatically. This code may still need further
tweaks and could probably be further simplified, but this may be a good
start to fixing it. When I tried this with OLS it appeared to provide
feasible results, but I haven't checked it any further than that.
* I'll assume we want to exclude missing values of grouping var
keep if e(sample) & `1'<.
bys `1': gen `nindg'=_N
local ntot=_N
bys `1': gen byte `perg'=_n==1
sum `perg', meanonly
local nind=r(sum)
replace `e'=`e'^2
sum `e', meanonly
scalar `sig'=r(sum)/`ntot'
egen `sigg'=sum(`e'), by (`1')
replace `sigg'=`sigg'/`nindg'
bysort `1': keep if _n==_N
gen `arg'=log(`sigg')*`nindg'
sum `arg', meanonly
scalar `sumg'=r(sum)
scalar `LR'=`ntot'*log(`sig')-`sumg'
return scalar lr= `LR'
return scalar pval= chiprob(`nind'-1,`LR')
Michael Blasnik
[email protected]
----- Original Message -----
From: "Clive Nicholas" <[email protected]>
To: <[email protected]>
Sent: Saturday, March 13, 2004 7:43 PM
Subject: st: Problems with -hetgrot-
> All,
>
> After plenty of 'behind-the-curtains' discussions with Richard Williams,
> we agreed that it would be a good idea to bring this to the attention of
> the list.
>
> Essentially, there appear to be problems in getting sensible chi-squared
> estimates out of -hetgrot-, a user-written procedure testing for groupwise
> heteroscedascity after -reg-, -areg-, -xtreg- and -xtgls- (I've been
> calling it after the last of these) which is based on a proposal by Greene
> (1993: 395).
>
> Statalisters will be aware that -hetgrot- as originally formulated was
> deleting cases after execution. With help from Kit Baum, Richard was able
> to solve this problem by carefully positioning -preserve- and -restore- in
> the right places. The revised code runs thus:
>
> --------------------------------------------------------------------------
--
> prog def hetgrot, rclass
>
> ***by Luca Nunziata, Nuffield College, Oxford
> ***Version 1.0, 1/4/00
> version 6
>
> preserve
>
> tempvar nind nindg e xbgls temp1 temp2 sig sigg a sumg arg
> tempname sig sumg LR
>
> if "`1'"=="" {
> display in red "ERROR, Name of group variable required as 1st argument"
> exit
> }
>
> if e(cmd)=="" {
> display in red "ERROR, need reg, xtreg or areg as command"
> exit 498
> }
>
> if e(cmd)== "regress" | e(cmd)== "areg" {
> predict `e', res
> }
>
> if e(cmd)=="xtreg" {
> predict `e', e
> }
>
> if e(cmd)=="xtgls" {
> predict `xbgls'
> gen `e'=$S_E_depv-`xbgls'
> }
>
> if e(cmd)~="regress" & e(cmd)~="areg" & e(cmd)~="xtreg" & e(cmd)~="xtgls"{
> display in red "ERROR, need reg, xtreg, xtgls or areg as command"
> exit 301
> }
>
> qui {
> gen `temp1'=1 if e(sample)==1
>
> egen `nindg'= count(`temp1'), by (`1')
>
> drop if e(sample)==0
>
>
> *determination of number of groups
> egen `nind' = group(`1')
>
> replace `nind'=`nind'[_N]
> scalar `nind'=`nind'[_N]
>
>
> *calculation of sigma
> replace `e'=`e'^2
> gen `sig'=sum(`e')
> replace `sig'=`sig'[_N]/`nind'
>
> scalar `sig'=`sig'[_N]/`nind'
>
> *calculation of sigman
> egen `sigg'=sum(`e'), by (`1')
> replace `sigg'=`sigg'/`nindg'
>
>
> ***calculation of LR test statistic
> *sum of groups var
> sort `1'
> by `1': gen `a'=1 if _n==_N
> replace `sigg'=log(`sigg')
> qui gen `arg'=`sigg'*`a'*`nindg'
> egen `sumg'=sum(`arg')
> scalar `sumg'=`sumg'
>
> scalar `LR'=`nind'*log(`sig')-`sumg'
> return scalar lr= `LR'
> return scalar pval= chiprob(`nind'-1,`LR')
>
> }
>
> *display results
>
> display""
> display in green "Testing for Groupwise heteroscedasticity"
> display""
> display "Ho: homoscedasticity"
> display "H1: groupwise heteroscedasticity by " "`1'"
> display""
> di _col(10) in gr "chi2(" in ye `nind'-1 in gr ")" /*
> */ _col(23) "=" in ye /*
> */ _col(25) %9.2f `LR'
> di _col(10) in gr "Prob > chi2 =" in ye _col(30) /*
> */ %6.4f chiprob(`nind'-1,`LR')
>
> restore
>
> end
> --------------------------------------------------------------------------
--
>
> The problem now is that -hetgrot- keeps producing _negative_ chi-squared
> values after -xtgls-. I've used it about 30 times (before and since the
> recoding by Kit and Richard) and it has always done this.
>
> I've not tried it against -areg- (as I little about the procedure), but
> calling
>
> . hetgrot [panelvar]
>
> consistently produces negatives after -xtgls- and -xtreg- (all forms).
> Richard mentioned to me that some papers from -hetgrot-'s originator, Luca
> Nunziata, does report using the procedure and comes out with positive
> numbers. I cannot replicate this with my data!
>
> I should say that I have tried to contact Luca, but no response has been
> forthcoming. Any help on this problem from anybody on the list would be
> most appreciated. Thanks.
>
> CLIVE NICHOLAS |t: 0(44)191 222 5969
> Politics Building |e: [email protected]
> School of Geography, |f: 0(44)870 126 2421
> Politics & Sociology |
> University of |
> Newcastle-upon-Tyne |
> Newcastle-upon-Tyne |
> NE1 7RU |
> United Kingdom |http://www.ncl.ac.uk/geps
>
> Greene WH (2nd ed, 1993) ECONOMETRIC ANALYSIS, New York: Prentice Hall.
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
