/*
** Lesson 9.4: Multiplicative Heteroscedasticity
** Greene [1997], Chap. 12.5
*/
use gpe2;
output file=gpe\output9.4 reset;
load data[52,3]=gpe\greene.txt;
data=data[2:52,2]~(data[2:52,3]/10000); @ scale data @
data=packr(miss(data,"NA")); @ take care of missing obs @
b=data[.,1]/(ones(rows(data),1)~data[.,2]~(data[.,2]^2));

call reset;
_method=4;
_iter=100;
_restart=10;
_b=b|2.0;
_nlopt=1;

call estimate(&llf,data);

end;

proc llf(data,b);
    local n,y,x,e,h,ll;
    y=data[.,1]; @ public school spending @
    x=data[.,2]; @ income @
    n=rows(y);
    h=x^b[4];    @ multiplicative hetero @
/*
    h=exp(b[4]*x);  
*/
    e=(y-b[1]-b[2]*x-b[3]*(x^2))./sqrt(h);
    ll=-0.5*n*(1+ln(2*pi)+ln(e'e/n))-0.5*sumc(ln(h));
    retp(ll);
endp;