proc (0)=hausman_ov(x,y,n,t);
/*
hausm_ov.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 Z_mu,P,Q;
local P_mu,P_lambda;
local Q_1,Q_2,Q_3,Q_4;
local k,z;
local sig_1;
local sig_v,sig_mu,sig_lam;
local lam_1,lam_2,lam_3,lam_4;
local theta_1,theta_2,theta_3;
local theta_SA;
local z_re,y_re,x_fe;
local x_ow,y_ow,z_ow;
local x_tw;
local res_u,urss;
local res_r,rrss;
local F_stat;

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

/* Transforming data using SWAR (one way model) */

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);

/* Estimating 1-way fixed individual effect */
x_fe=Q*x;

/* Calculate residual sum of squares for restricted and unrestricted model */
res_u=y_re-(z_re~x_fe)*(y_re/(z_re~x_fe));
urss=res_u'*res_u;

res_r=y_re-(z_re)*(y_re/z_re);
rrss=res_r'*res_r;

F_stat=((rrss-urss)/cols(x_fe))/(urss/(n*t-cols(z_re~x_fe)));


/*Calculate Hausman statistic */

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 "--------------------------------------OV version-----------------------------------------";
print "**************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " F_stat;
format /m1 /rd 8,4; print "p-value:           " cdffc(F_stat,cols(x_fe),n*t-cols(z_re~x_fe));
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);

/* Estimating 1-way fixed time effect */
x_fe=Q*x;


/* Calculate residual sum of squares for restricted and unrestricted model */
res_u=y_re-(z_re~x_fe)*(y_re/(z_re~x_fe));
urss=res_u'*res_u;

res_r=y_re-(z_re)*(y_re/z_re);
rrss=res_r'*res_r;

F_stat=((rrss-urss)/cols(x_fe))/(urss/(n*t-cols(z_re~x_fe)));


/*Calculate Hausman statistic */

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 time-unit FE ";
print "--------------------------------------OV version-----------------------------------------";
print "**************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " F_stat;
format /m1 /rd 8,4; print "p-value:           " cdffc(F_stat,cols(x_fe),n*t-cols(z_re~x_fe));
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);
k=cols(z_re);

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

/* Calculate residual sum of squares for restricted and unrestricted model */
res_u=y_re-(z_re~x_tw)*(y_re/(z_re~x_tw));
urss=res_u'*res_u;

res_r=y_re-(z_re)*(y_re/z_re);
rrss=res_r'*res_r;

F_stat=((rrss-urss)/cols(x_fe))/(urss/(n*t-cols(z_re~x_fe)));


/*Calculate F-statistic */

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 "----------------------------------------------------------------OV version----------------------------------------------------------------";
print "*************************************************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " F_stat;
format /m1 /rd 8,4; print "p-value:           " cdffc(F_stat,cols(x_tw),n*t-cols(z_re~x_tw));
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);

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

/* Calculate residual sum of squares for restricted and unrestricted model */
res_u=y_re-(z_re~x_tw)*(y_re/(z_re~x_tw));
urss=res_u'*res_u;

res_r=y_re-(z_re)*(y_re/z_re);
rrss=res_r'*res_r;

F_stat=((rrss-urss)/cols(x_fe))/(urss/(n*t-cols(z_re~x_fe)));

/*Calculate F-statistic */

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 an individual RE and a time FE under the null.   ";
print "----------------------------------------------------------------OV version----------------------------------------------------------------";
print "*************************************************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " F_stat;
format /m1 /rd 8,4; print "p-value:           " cdffc(F_stat,cols(x_tw),n*t-cols(z_re~x_tw));
print "*************************************************************************************************";


/*********** Two way test ***********/

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

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

y_re=y-theta_1*P_mu*y-theta_2*P_lambda*y+theta_3*P_mu*P_lambda*y;
z_re=z-theta_1*P_mu*z-theta_2*P_lambda*z+theta_3*P_mu*P_lambda*z;


/*************************** 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_fe=Q*x;

/* Calculate residual sum of squares for restricted and unrestricted model */
res_u=y_re-(z_re~x_fe)*(y_re/(z_re~x_fe));
urss=res_u'*res_u;

res_r=y_re-(z_re)*(y_re/z_re);
rrss=res_r'*res_r;

F_stat=((rrss-urss)/cols(x_fe))/(urss/(n*t-cols(z_re~x_fe)));


/*Calculate Hausman statistic */

print "*************************** Two way RE/FE test ***************************";
print "Hausman test for the null of fully RE against the alternative of fully FE";
print "--------------------------------------OV version-----------------------------------------";
print "*****************************************************************************";
format /m1 /rd 8,4; print "Test statistic:    " F_stat;
format /m1 /rd 8,4; print "p-value:           " cdffc(F_stat,cols(x_fe),n*t-cols(z_re~x_fe));
print "*****************************************************************************";

retp;
endp;