proc (0)=HAC_2wRe(x,y,n,t,theta,names);
/*
hac_2wre.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.
*/

/* Estimates Arellano's HAC VC-matrix */
/* See notes p 51 */
/* Est.Asy.Var(b)=inv(x'x)*sum_i(x_i'e_i*e_i'*x_i)*inv(x'x) */
/* x_i is a kx1 vector containing the i:th observation of the k variables */
/* A degrees of freedom correction is performed to obtain results identical to E-views (Note however vc is derived asy) */

local beta,res;
local S0,temp_sum,var_w;
local nobs,k;
local i,j;
local z_out,fmat,mask;
local name;
local Q, z_re,y_re,x_tr,z;
local P_mu,P_lambda;
local th_method;

th_method= "Wallace & Hussain (1969)" | "Amemiya (1971)" | "Swamy    and Arora (1972)" | "Nerlove (1971)";

P_mu=eye(n).*.ones(t,t)/t;     /* Time average for different individuals */
P_lambda=ones(n,n)/n.*.eye(t);  /* Individual average for different time-periods */ 

z=ones(rows(x),1)~x;


for j (1,4,1);

y_re=y-theta[j,1]*P_mu*y-theta[j,2]*P_lambda*y+theta[j,3]*P_mu*P_lambda*y;
z_re=z-theta[j,1]*P_mu*z-theta[j,2]*P_lambda*z+theta[j,3]*P_mu*P_lambda*z;

beta=y_re/z_re;
k=cols(z_re);
res=y_re-z_re*beta;
nobs=rows(x);

x_tr=z_re;

temp_sum=zeros(k,k);
for i (1,n,1);
temp_sum=temp_sum+x_tr[(i-1)*t+1:(i-1)*t+t,.]'*(res[(i-1)*t+1:(i-1)*t+t,.]*res[(i-1)*t+1:(i-1)*t+t,.]')*x_tr[(i-1)*t+1:(i-1)*t+t,.];
endfor;
S0=temp_sum;

var_w=inv(x_tr'*x_tr)*S0*inv(x_tr'*x_tr);


mask = ones(1,2); mask[1]=0; /*zero for character output and 1 for numeric output */
fmat=zeros(2,3); fmat[1,.]="*.*s"~8~8; fmat[2,.]="*.*lf"~12~3;

print;print;print;

print "*******************************************************************";
print "th_method: " j "[W&H(1), Amemiya(2), S&A(3), Nerlove(4)]";
print "*******************************************************************";
print "************ Arellano's HAC SEs in 2-way RE model ************";
z_out=printfm(names~sqrt(diag(var_w)),mask,fmat);
print "*******************************************************************";

endfor;

retp;
endp;