proc (0)=hsced_test(x,y,n,t,j);
/*
hsced.src is written by Kristian Jönsson (November 12, 2002), Department of Economics, Lund University.
Contact info:  kristian.jonsson@nek.lu.se

See separate documentation and 'Econometric Analysis of Panel Data' 2nd edition, Baltagi (2001).

Note:
This code can be used freely as long as proper reference is given.
No performance guarantee is made.
Bugreports are welcome.
*/
local Q;
local x_fe,y_fe,y_hat,y_fit;
local u,u_sq;
local res_r,res_u;
local rrss,urss;
local F_stat,pv;
local i,k;
local res_mat;
local var_i,var_tot;
local Bartlett,pv_Bartlett;

k=cols(x);

Q=eye(n).*.eye(t)-eye(n).*.(ones(t,t)/t)-(ones(n,n)/n).*.eye(t)+(ones(n,n)/n).*.(ones(t,t)/t);

x_fe=Q*x;
y_fe=Q*y;

u=y_fe-x_fe*(y_fe/x_fe);

y_hat=y-u; /* Or: xb */
y_fit=y_hat;

for i (2,j,1);
y_fit=y_fit~y_fit^i;
endfor;

y_fit=ones(rows(y_hat),1)~y_fit;
u_sq=u^2;
res_r=u_sq-ones(rows(y_hat),1)*(u_sq/ones(rows(y_hat),1));
rrss=res_r'*res_r;

res_u=u_sq-y_fit*(u_sq/y_fit);
urss=res_u'*res_u;

F_stat=((rrss-urss)/(cols(y_fit)-1))/(urss/(rows(y)-cols(y_fit)));
pv=cdffc(F_stat,cols(y_fit)-1,rows(y)-cols(y_fit));

print "**********************************************************************";
print "Bickel's version of the B-P test for the null of no heteroskedasticity";
print "**********************************************************************";
print "F-statistic:  " F_stat;
print "    p-value:  " pv;
print "**********************************************************************";


/* Bartlett's test */
/* Uses the same residuals as above */

res_mat=reshape(u,n,t)';
var_i=sumc(res_mat.*res_mat)/(t-1);
var_tot=u'*u/(rows(u)-n-k);

Bartlett=((t-1)*(n*log(var_tot)-sumc(log(var_i))))/(1+(n+1)/(3*t-3));
pv_Bartlett=cdfchic(Bartlett,n-1);

print "*****************************************************";
print "Bartlett's test for the null of no heteroskedasticity";
print "*****************************************************";
print "Test-statistic:  " Bartlett;
print "    p-value:  " pv_Bartlett;
print "*****************************************************";
print;
print "****************************************************************";
print "Note: The Hsced tests are conducted in a two-way FE model";
print "****************************************************************";

retp;
endp;