/*
** Lesson 9.3: Breusch-Pagan and White Tests
** for Heteroscedasticity
*/
use gpe2;
output file = gpe\output9.3 reset;
load greene[52,3]= gpe\greene.txt;
data = greene[2:52,.];
               
data = miss(data,"NA"); @ NA to missing value @
data = packr(data);   @ deletes row w/missing value @
spending = data[.,2];
income = data[.,3]/10000;

/* Breusch-Pagan and White Tests */
call reset;
_names = {"spending","income","income^2"};
_bjtest = 1;
_bptest = 1;

call estimate(spending,income~income^2);

/* do bptest manually */
e2=__e^2;
v=meanc(e2);

_bjtest=0;
_bptest=0;
call estimate(e2,income~income^2);

e2hat=e2-__e;
ess=sumc((e2hat-meanc(e2hat))^2);
bp=ess/(2*v^2);
kb=ess/meanc((e2-v)^2);
print "Breusch-Pagan"  bp;
print "Koenkar-Basset" kb;

end;