% Hausman and Taylor Estimator
% Unknown Author
% This function does not allow time-invariant variable Z
% Updated by Wonho Song, May 2014
% E-mail: whsong@cau.ac.kr, whsong73@hotmail.com

function [bht,seht,HT,EHT,R2h]=ht(y,x,p1,n)

[nt,p]=size(x);
t=nt/n;
v=p-p1;


iota=ones(nt,1);
Piota=iota*inv(iota'*iota)*iota';
%Qiota=eye(nt)-Piota;                      % eye(nt): out of memory problem encountered 
y=y-Piota*y;
x=x-Piota*x; 


xbar = reshape(mean(reshape(x,t,n*p)),n,p);
ybar = mean(reshape(y,t,n))';

xtilde = x - kron(xbar,ones(t,1));
ytilde = y - kron(ybar,ones(t,1));

ttilde = inv(xtilde'*xtilde)*xtilde'*ytilde;

utilde = ytilde - xtilde*ttilde;
stilde2= utilde'*utilde/(n*(t-1));
sbar=ybar - xbar*ttilde;

sbarm = mean(sbar);
ssbar = sbar - sbarm;
sb2 = ssbar'*ssbar/n;
theta = (stilde2/(t*sb2))^(.5);

z = x(:,p-v+1:p);
zbar = xbar(:,p-v+1:p);

nx = x(:,1:p1);
nxbar= xbar(:,1:p1);

nxtilde = xtilde(:,1:p1);
ztilde = z - kron(ones(t,1),zbar);

% Hausman - Taylor
xqv = [kron(nxbar,ones(t,1)) xtilde];

% x2hat = xqv*inv(xqv'*xqv)*xqv'*z;  
                    %  we do not use this line because it does not work if the matrix is too large. 
                    % Instead, we use the following lines.

tmp1=inv(xqv'*xqv);
tmp2=xqv'*z;
x2hat = xqv*tmp1*tmp2;

xx = [nx x2hat ones(n*t,1)];
[temp,xp]=size(xx);


xxbar   =   reshape(mean(reshape(xx,t,n*xp)),n,xp);

xstar = xx - kron(xxbar,ones(t,1))+kron((theta*xxbar),ones(t,1));
ystar = y - kron(ybar,ones(t,1)) + kron((theta*ybar),ones(t,1));

bht = inv(xstar'*xstar)*xstar'*ystar;
tglscov = stilde2*inv(xstar'*xstar);

bht = bht(1:p);
seht = sqrt(diag(tglscov));
seht=seht(1:p);

residual=y-x*bht;
residual=residual-mean(residual);

eh=residual;

% Calculation of R2 
ey=y-mean(y);
R2h=1-(eh'*eh)/(ey'*ey);
R2h=[R2h;1-(1-R2h)*(nt-1)/(nt-p-n)];



HT=mean(reshape(residual,t,n))';
HT=HT-mean(HT);
EHT=exp(HT-max(HT));