proc (0)=hausman_mat(x,y,n,t);
/*
hausm_m.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,P;
local Z_mu;
local P_mu,P_lambda,Q_lambda;
local Q_1,Q_2,Q_3,Q_4;
local lam_1,lam_2,lam_3,lam_4;
local theta_1,theta_2,theta_3;
local theta_SA;
local sig_v,sig_1,sig_lam,sig_mu;
local omega_hat;
local x_tw,y_tw,z,z_re,y_re;
local x_ow,z_ow,y_ow;
local var_within,beta_within,beta_fgls,var_fgls;
local h_stat;
local q1,var_q1;
local k,nvar;

nvar=cols(x);

/* ********** One-way individual effect test ********** */

Z_mu=eye(n).*.ones(T,1);
P=Z_mu*inv(Z_mu'*Z_mu)*Z_mu'; /* Define P and Q matrices */
Q=eye(n*t)-P;

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

/* Estimating FGLS using SWAR */

sig_v=sqrt((y'*Q*y-y'*Q*x*inv(x'*Q*x)*x'*Q*y)/(n*(t-1)-k));
sig_1=sqrt((y'*P*y-y'*P*z*inv(z'*P*z)*z'*P*y)/(n-k-1));

theta_SA=1-sig_v/sig_1;

y_re=y-theta_SA*(P*y);
z_re=z-theta_SA*(P*z);
beta_fgls=y_re/z_re;
k=cols(z_re);
beta_fgls=beta_fgls[2:k];
omega_hat=sig_1^2*P+sig_v^2*Q;
var_fgls=inv(z'*inv(omega_hat)*z);
var_fgls=var_fgls[2:k,2:k];

/* Estimating 1-way fixed individual effect */
x_tw=Q*x;
y_tw=Q*y;
beta_within=y_tw/x_tw;
var_within=sig_v^2*inv(x'*Q*x);

/*Calculate Hausman statistic */
q1=beta_fgls-beta_within;
var_q1=var_within-var_fgls;
h_stat=q1'*inv(var_q1)*q1;

print "********************* One way RE/FE test *********************";
print "Hausman test for the null of RE against the alternative of FE";
print " The one way fixed effects model contains an individual FE ";
print "-----------------------------------Matrix version-------------------------------------";
print "**************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " h_stat;
format /m1 /rd 8,4; print "p-value:           " cdfchic(h_stat,nvar);
print "**************************************************************";


/* ********** One-way time effect test ********** */

P=ones(n,n)/n.*.eye(t); 
Q=eye(n*t)-P; /* Sweeps averages over different individuals for a specific time period. */

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

/* Estimating FGLS using SWAR */

sig_v=sqrt((y'*Q*y-y'*Q*x*inv(x'*Q*x)*x'*Q*y)/(n*(t-1)-k));
sig_1=sqrt((y'*P*y-y'*P*z*inv(z'*P*z)*z'*P*y)/(n-k-1));

theta_SA=1-sig_v/sig_1;

y_re=y-theta_SA*(P*y);
z_re=z-theta_SA*(P*z);
beta_fgls=y_re/z_re;
k=cols(z_re);
beta_fgls=beta_fgls[2:k];
omega_hat=sig_1^2*P+sig_v^2*Q;
var_fgls=inv(z'*inv(omega_hat)*z);
var_fgls=var_fgls[2:k,2:k];

/* Estimating 1-way fixed time effect */
x_tw=Q*x;
y_tw=Q*y;
beta_within=y_tw/x_tw;
var_within=sig_v^2*inv(x'*Q*x);

/*Calculate Hausman statistic */
q1=beta_fgls-beta_within;
var_q1=var_within-var_fgls;
h_stat=q1'*inv(var_q1)*q1;

print "********************* One way RE/FE test *********************";
print "Hausman test for the null of RE against the alternative of FE";
print " The one way fixed effects model contains a time-unit FE ";
print "-----------------------------------Matrix version-------------------------------------";
print "**************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " h_stat;
format /m1 /rd 8,4; print "p-value:           " cdfchic(h_stat,nvar);
print "**************************************************************";



/* ********** Two-way individual FE and time RE test ********** */
/*********** Two way test H0={mu_FE,lam_RE} H1={mu_FE,lam_FE}***********/

Z_mu=eye(n).*.ones(T,1);
P=Z_mu*inv(Z_mu'*Z_mu)*Z_mu'; /* Define P and Q matrices */
Q=eye(n*t)-P;

x_ow=Q*x;
y_ow=Q*y;

k=cols(x);
z_ow=ones(rows(x),1)~x_ow;

/* Estimating FGLS using SWAR */

P=ones(n,n)/n.*.eye(t); 
Q=eye(n*t)-P; /* Sweeps averages over different individuals for a specific time period. */

sig_v=sqrt((y_ow'*Q*y_ow-y_ow'*Q*x_ow*inv(x_ow'*Q*x_ow)*x_ow'*Q*y_ow)/(t*(n-1)-k));
sig_1=sqrt((y_ow'*P*y_ow-y_ow'*P*z_ow*inv(z_ow'*P*z_ow)*z_ow'*P*y_ow)/(t-k-1));

theta_SA=1-sig_v/sig_1;

y_re=y_ow-theta_SA*(P*y_ow);
z_re=z_ow-theta_SA*(P*z_ow);
beta_fgls=y_re/z_re;
k=cols(z_re);
beta_fgls=beta_fgls[2:k];
omega_hat=sig_1^2*P+sig_v^2*Q;
var_fgls=inv(z'*inv(omega_hat)*z);
var_fgls=var_fgls[2:k,2:k];

/* Estimating 2-way fixed effect */
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_tw=Q*x;
y_tw=Q*y;
beta_within=y_tw/x_tw;
var_within=sig_v^2*inv(x'*Q*x);

/*Calculate Hausman statistic */
q1=beta_fgls-beta_within;
var_q1=var_within-var_fgls;
h_stat=q1'*inv(var_q1)*q1;

print "************************************** Two way RE/FE test *************************************";
print "Hausman test for the null of time RE and indiv. FE against the alternative of fully FE model";
print " The two way fixed effects model contains a time RE and an individual FE under the null.   ";
print "--------------------------------------------------------------Matrix version-------------------------------------------------------------";
print "**************************************************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " h_stat;
format /m1 /rd 8,4; print "p-value:           " cdfchic(h_stat,nvar);
print "**************************************************************************************************";






/* ********** Two-way individual RE and time FE test ********** */
/*********** Two way test H0={mu_RE,lam_FE} H1={mu_FE,lam_FE}***********/

P=ones(n,n)/n.*.eye(t); 
Q=eye(n*t)-P; /* Sweeps averages over different individuals for a specific time period. */

x_ow=Q*x;
y_ow=Q*y;

k=cols(x);
z_ow=ones(rows(x),1)~x_ow;

/* Estimating FGLS using SWAR */

Z_mu=eye(n).*.ones(T,1);
P=Z_mu*inv(Z_mu'*Z_mu)*Z_mu'; /* Define P and Q matrices */
Q=eye(n*t)-P;

sig_v=sqrt((y_ow'*Q*y_ow-y_ow'*Q*x_ow*inv(x_ow'*Q*x_ow)*x_ow'*Q*y_ow)/(t*(n-1)-k));
sig_1=sqrt((y_ow'*P*y_ow-y_ow'*P*z_ow*inv(z_ow'*P*z_ow)*z_ow'*P*y_ow)/(t-k-1));

theta_SA=1-sig_v/sig_1;

y_re=y_ow-theta_SA*(P*y_ow);
z_re=z_ow-theta_SA*(P*z_ow);
beta_fgls=y_re/z_re;
k=cols(z_re);
beta_fgls=beta_fgls[2:k];
omega_hat=sig_1^2*P+sig_v^2*Q;
var_fgls=inv(z'*inv(omega_hat)*z);
var_fgls=var_fgls[2:k,2:k];

/* Estimating 2-way fixed effect */
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_tw=Q*x;
y_tw=Q*y;
beta_within=y_tw/x_tw;
var_within=sig_v^2*inv(x'*Q*x);

/*Calculate Hausman statistic */
q1=beta_fgls-beta_within;
var_q1=var_within-var_fgls;
h_stat=q1'*inv(var_q1)*q1;

print "************************************** Two way RE/FE test *************************************";
print "Hausman test for the null of indiv RE and time FE against the alternative of fully FE model";
print " The two way fixed effects model contains a individual RE and an time FE under the null.   ";
print "--------------------------------------------------------------Matrix version-------------------------------------------------------------";
print "**************************************************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " h_stat;
format /m1 /rd 8,4; print "p-value:           " cdfchic(h_stat,nvar);
print "**************************************************************************************************";



/*********** Two way test H0={mu_RE,lam_RE} H1={mu_FE,lam_FE}***********/

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 */ 

Q_1=(eye(n)-ones(n,n)/n).*.(eye(t)-ones(t,t)/t);
Q_2=(eye(n)-ones(n,n)/n).*.(ones(t,t)/t);
Q_3=(ones(n,n)/n).*.(eye(t)-ones(t,t)/t);
Q_4=(ones(n,n)/n).*.(ones(t,t)/t);

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

/* Estimating FGLS using SWAR */
/* Swamy and Arora */

sig_v=sqrt((y'*Q_1*y-y'*Q_1*x*inv(x'*Q_1*x)*x'*Q_1*y)/((n-1)*(t-1)-k+1));
lam_1=sig_v^2;

lam_2=(y'*Q_2*y-y'*Q_2*x*inv(x'*Q_2*x)*x'*Q_2*y)/((n-1)-k+1);
sig_mu=sqrt((lam_2-sig_v^2)/t);

lam_3=(y'*Q_3*y-y'*Q_3*x*inv(x'*Q_3*x)*x'*Q_3*y)/((t-1)-k+1);
sig_lam=sqrt((lam_3-sig_v^2)/n);

lam_2=T*maxc(0|(lam_2-lam_1)/T)+lam_1;
lam_3=N*maxc(0|(lam_3-lam_1)/N)+lam_1;
lam_4=lam_2+lam_3-lam_1;

lam_2; lam_3;

theta_1=1-sig_v/sqrt(lam_2); 
theta_2=1-sig_v/sqrt(lam_3); 
theta_3=theta_1+theta_2+sig_v/sqrt(lam_4)-1;



/*
sig_v=sqrt((y'*Q_1*y-y'*Q_1*x*inv(x'*Q_1*x)*x'*Q_1*y)/((n-1)*(t-1)-k+1));
lam_1=sig_v^2;

lam_2=(y'*Q_2*y-y'*Q_2*x*inv(x'*Q_2*x)*x'*Q_2*y)/((n-1)-k+1);
sig_mu=sqrt((lam_2-sig_v^2)/t);

lam_3=(y'*Q_3*y-y'*Q_3*x*inv(x'*Q_3*x)*x'*Q_3*y)/((t-1)-k+1);
sig_lam=sqrt((lam_3-sig_v^2)/n);

lam_4=lam_2+lam_3-lam_1;

theta_1=1-sig_v/sqrt(lam_2); if theta_1<0; theta_1=0; endif;
theta_2=1-sig_v/sqrt(lam_3); if theta_2<0; theta_2=0; endif;
theta_3=theta_1+theta_2+sig_v/sqrt(lam_4)-1; if theta_3<0; theta_3=0; endif;
*/

omega_hat=lam_1*Q_1+lam_2*Q_2+lam_3*Q_3+lam_4*Q_4;

y_re=y-theta_1*P_mu*y-theta_2*P_lambda*y+theta_3*ones(rows(y),1)*meanc(y);
z_re=z-theta_1*P_mu*z-theta_2*P_lambda*z+theta_3*ones(rows(y),cols(z)).*meanc(z)';

beta_fgls=y_re/z_re;
/*var_fgls=inv(z'*inv(omega_hat)*z);*/
var_fgls=(y_re-z_re*beta_fgls)'*(y_re-z_re*beta_fgls)/(n*t-cols(z))*inv(z_re'*z_re);
beta_fgls=beta_fgls[2:k];
var_fgls=var_fgls[2:k,2:k];

/* Estimating 2-way fixed effect */
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_tw=Q*x;
y_tw=Q*y;
beta_within=y_tw/x_tw;
var_within=sig_v^2*inv(x_tw'*x_tw);

/*Calculate Hausman statistic */
q1=beta_fgls-beta_within;
var_q1=var_within-var_fgls;
h_stat=q1'*inv(var_q1)*q1;

print;print;
print "*************************** Two way RE/FE test ***************************";
print "Hausman test for the null of fully RE against the alternative of fully FE";
print "-----------------------------------Matrix version-------------------------------------";
print "*****************************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " h_stat;
format /m1 /rd 8,4; print "p-value:           " cdfchic(h_stat,nvar);
print "*****************************************************************************";



retp;
endp;