% Estimators
function [beta, VC, rho, rho_se, sgmu2]=sht_het(Y,XU,XC,ZU,ZC,W,T); 
% Spatial Hausman-Taylor estimator under heteroskedastic innovations
% OUTPUT
% beta:   parameter estimates on time-varying and time-invariant variables
% VC:     Variance-Covariance matrix of beta
% rho:	  spatial autocorrelation parameter
% rho_se: standard error of the spatial autorcorrelation parameter
% sgmu2:  sigma_mu^2
% INPUT
% T:  total number of time units
% Y:  dependent variable (NTx1) where N denotes the total number of
%       individual units and NT the total number of observations
% XU: matrix of time-varying exogenous covariates(NTxku) where ku denotes
%       the number of exogenous time-varying covariates 
% XC: matrix of time-varying endogenous covariates(NTxkc) where kc denotes
%       the number of endogenous time-varying covariates 
% ZU: matrix of time-invariant exogenous covariates(NTxru) where ru denotes
%       the number of exogenous time-invariant covariates 
% ZC: matrix of time-invariant endogenous covariates(NTxrc) where rc denotes
%       the number of endogenous time-invariant covariates 
% W:  normalized spatial weights matrix (NTxNT)


% Preliminaries
NT=length(Y); N=NT/T; JT=ones(T); Q1=kron(JT/T, eye(N)); Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));
Xbar=Q1*[XU XC]; Ybar=Q1*Y; Xtil=Q0*[XU XC]; Ytil=Q0*Y;

% Step 1: Estimate bet0 with fixed effects 
XtilXtil=Xtil'*Xtil; XtilYtil=Xtil'*Ytil;
bet0=XtilXtil\XtilYtil;
dhat=Ybar-Xbar*bet0;

% Step 2: Estimate thet0 with 2SLS of dhat on Z using instrument set HI=[XU ZU];
HI=[XU ZU]; 
PHI=HI*inv(HI'*HI)*HI';
ZPHIZ=[ZU ZC]'*PHI*[ZU ZC]; ZPHIdhat=[ZU ZC]'*PHI*dhat;
thet0=ZPHIZ\ZPHIdhat;
[bethet0]=[bet0;thet0];
u=Y-[XU XC ZU ZC]*bethet0;

% Step 3: Estimate rho and sigma_mu^2
[rho, rho_se,sgmu2]= gmhetht(u,N,T,W,Q0,Q1,eTIN,XU,XC,ZU,ZC);

% Step 4: 
% Spatial Cochrane Orcutt transformation
WY=W*Y; WXU=W*XU; WXC=W*XC; WZU=W*ZU; WZC=W*ZC;
Ys=Y-rho*WY; XUs=XU-rho*WXU; XCs=XC-rho*WXC; ZUs=ZU-rho*WZU; 
ZCs=ZC-rho*WZC;

% Estimate transformed model
As=[XUs XCs ZUs ZCs]; Xstil=Q0*[XUs XCs]; XUsbar=Q1*XUs;
HHTS=[Xstil XUsbar ZUs];
HHTS2=HHTS'*HHTS; temp=(HHTS/HHTS2); PAHT=temp*HHTS';
APA=As'*PAHT*As; APY=As'*PAHT*Ys;
beta=APA\APY;

u=Y-[XU XC ZU ZC]*beta;
[SgHR]= sigmaHR(u,N,T,W,Q0,eTIN,rho);
Fmu=(kron(ones(1,T),eye(N)))*HHTS; Fnu=HHTS;
Omg_mu=1/NT*sgmu2*(Fmu'*Fmu);
Omg_nu=1/NT*(Fnu'*SgHR*Fnu);
Omg=Omg_mu+Omg_nu;
HA=HHTS'*As; AH=As'*HHTS;
P1=HHTS2\HA; temp=AH/HHTS2; P2=temp*(1/NT*HA); D=P1/P2;
VC=1/NT*(D'*Omg*D);


function [rho,rho_se,sgmu2]= gmhetht(u,N,T,W,Q0,Q1,eTIN,XU,XC,ZU,ZC);
% GM procedure for the spatial Hausman-Taylor model under heteroskedastic
% innovations
% OUTPUT
% rho:	  spatial autocorrelation parameter
% rho_se: standard error of the spatial autorcorrelation parameter
% sgmu2:  sigma_mu^2
% INPUT
% N:      total number of individual units 
% T:      total number of time units
% u: 	  estimated residual (NTx1)
% W:      normalized spatial weights matrix (NTxNT)
% Q0:     within-transformation matrix (NTxNT)
% Q1:     between-transformation matrix (NTxNT)
% eTIN:   eTIN=kron(ones(1,T),eye(N))
% XU:     matrix of time-varying exogenous covariates(NTxku) where ku 
%           denotes the number of exogenous time-varying covariates 
% XC:     matrix of time-varying endogenous covariates(NTxkc) where kc 
%           denotes the number of endogenous time-varying covariates 
% ZU:     matrix of time-invariant exogenous covariates(NTxru) where ru 
%           denotes the number of exogenous time-invariant covariates 
% ZC:     matrix of time-invariant endogenous covariates(NTxrc) where rc 
%           denotes the number of endogenous time-invariant covariates 

% Preliminaries
NT=N*T; WpW=W(1:N,1:N)'*W(1:N,1:N); trWpW=trace(WpW); ITWW=kron(eye(T),WpW);
Q1ITWW=Q1*ITWW;

% GM Elements
ubar=W*u;
ubarbar=W*ubar;

ga=zeros(4,1); GA=zeros(4,3);
ga(2,1)=1/(N*(T-1))*(ubar'*Q0*u);
ga(3,1)=1/N*(ubar'*Q1*ubar-trace(bsxfun(@times,Q1ITWW,u.^2)));

GA(2,1)=1/(N*(T-1))*(ubarbar'*Q0*u+ubar'*Q0*ubar);
GA(2,2)=-1/(N*(T-1))*(ubarbar'*Q0*ubar);
GA(2,3)=0;
GA(3,1)=2/N*(ubar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,(ubar.*u))));
GA(3,2)=-1/N*(ubarbar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,ubar.^2)));
GA(3,3)=(T-1)/N*trWpW;

gmfun0=@(aa)(ga(2,1)-GA(2,1:2)*[aa(1);aa(1)^2]);  
a0=0; 

% Estimate for rho
[ai] = lsqnonlin(gmfun0,a0);
rho=ai(1,1); 

% Estimate for sigma_mu^2
aa13=ga(3,1)-GA(3,1)*rho-GA(3,2)*rho^2;  
sgmu2=aa13*N/((T-1)*trWpW);

% Standard error of rho 
J=GA(2,1:2)*[1; 2*rho];
JI=inv(J);
XZ=[XU XC ZU ZC];
Xtil=Q0*[XU XC]; XUbar=Q1*XU;
HHT=[Xtil XUbar ZU];
HHT2=HHT'*HHT; 
HXZ=HHT'*XZ; XZH=XZ'*HHT;
P1=HHT2\HXZ; temp=XZH/HHT2; P2=temp*(1/NT*HXZ);
P=P1/P2;
A2=Q0*W;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho*W(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-XZ'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*HHT;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
[SgHR]= sigmaHR(u,N,T,W,Q0,eTIN,rho);
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=R*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
cc=T/((T-1)^2*N);
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho_se=sqrt(Omrho);


function [SgHR]= sigmaHR(u,N,T,W,Q0,eTIN,rho);
% INPUT
% N:      total number of individual units 
% T:      total number of time units
% u:      estimated residual (NTx1) 
% W:      normalized spatial weights matrix (NTxNT)
% Q0:     within-transformation matrix (NTxNT)
% eTIN:   eTIN=kron(ones(1,T),eye(N))
% rho:    spatial autocorrelation parameter

rho=(eye(N*T)-rho*W);
ee=rho*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=eTIN*nutil2;
tempB=kron(ones(T,1),temp);
tempF=1/((T-1)*(T-2))*tempB;
nutil2=T/(T-2)*nutil2;
nu2HR=nutil2-tempF;
[SgHR]=diag(nu2HR);

function [beta, VC, rho, rho_se]=sfe_het(Y,X,W,T); 
% Spatial fixed effects estimator under heteroskedastic innovations
% OUTPUT
% beta:   parameter estimates
% VC:     Variance-Covariance matrix of beta
% rho:	  spatial autocorrelation parameter
% rho_se: standard error of the spatial autorcorrelation parameter
% INPUT
% T:      total number of time units
% Y:      dependent variable (NTx1) where N is the total number of
%               individual units and NT the total number of observations
% X:      matrix of time-varying covariates (NTxk) where k is the number of
%               time-varying covariates
% W:      normalized spatial weights matrix (NTxNT)


% Preliminaries
NT=length(Y); N=NT/T;
JT=ones(T); Q1=kron(JT/T, eye(N)); Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));

Xtil=Q0*X; Ytil=Q0*Y;
XtilXtil=Xtil'*Xtil; XtilYtil=Xtil'*Ytil;

% Initial estimator
bet0=XtilXtil\XtilYtil;
u=Y-X*bet0;

% Estimate and standard error of rho
[rho, rho_se] = gmhetfe(u,N,T,W,Q0,eTIN,X);

% Spatial Cochrane Orcutt and within transformation
WY=W*Y; WX=W*X;
Ys=Y-rho*WY; Xs=X-rho*WX; 
Xstil=Q0*Xs; Ystil=Q0*Ys;

% Estimate transformed model
XstilXstil=Xstil'*Xstil; XstilYstil=Xstil'*Ystil;
beta=XstilXstil\XstilYstil;

us=Y-X*beta;
[SgHR] = sigmaHR(us,N,T,W,Q0,eTIN,rho);
Fnu=Q0*Xs;
Omg=1/NT*(Fnu'*SgHR*Fnu);
P_fes=(1/NT*(Fnu'*Fnu))\eye(size(X,2));
VC=1/NT*(P_fes'*Omg*P_fes);


function [rho, rho_se]= gmhetfe(u,N,T,W,Q0,eTIN,X);
% GM procedure for the spatial fixed effects model under heteroskedastic
% innovations
% OUTPUT
% rho:	  spatial autocorrelation parameter
% rho_se: standard error of the spatial autorcorrelation parameter
% INPUT
% N:      total number of individual units 
% T:      total number of time units
% u: 	  estimated residual (NTx1)
% W:      normalized spatial weights matrix (NTxNT)
% Q0:     within-transformation matrix (NTxNT)
% eTIN:   eTIN=kron(ones(1,T),eye(N))
% X:      matrix of time-varying covariates (NTxk) where k is the number 
%           of time-varying covariates

% Preliminaries
NT=N*T;

% GM Elements
ubar=W*u;
ubarbar=W*ubar;

GA=zeros(1,2);
ga=1/(N*(T-1))*(ubar'*Q0*u);
GA(1,1)=1/(N*(T-1))*(ubarbar'*Q0*u+ubar'*Q0*ubar);
GA(1,2)=-1/(N*(T-1))*(ubarbar'*Q0*ubar);

gmfun0=@(aa)(ga-GA*[aa(1);aa(1)^2]);  
a0=0; 

% Estimate for rho
[ai] = lsqnonlin(gmfun0,a0);
rho=ai(1,1);

% Standard error of rho 
J=GA*[1; 2*rho];
JI=inv(J);
Xfe=Q0*X;
P=(1/NT*(Xfe')*Xfe)\eye(size(Xfe,2)); 
A2=Q0*W;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*(A2A2t);
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho*W(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-Xfe'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*Xfe;
anu=Fnu*P*alpha;
[SgHR]= sigmaHR(u,N,T,W,Q0,eTIN,rho);
Enu=1/NT*(anu'*SgHR*anu);
ee=R*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
cc=T/((T-1)^2*N);
Enuc=cc*temp;
E=Enu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho_se=sqrt(Omrho);

function [beta, VC, rho, rho_se, sgmu2]=sre_het(Y,X,W,T); 
% Spatial random effects estimator under heteroskedastic innovations
% OUTPUT
% beta:   parameter estimates on time-varying and time-invariant variables
% VC:     Variance-Covariance matrix of beta
% rho:	  spatial autocorrelation parameter
% rho_se: standard error of the spatial autorcorrelation parameter
% sgmu2:  sigma_mu^2
% INPUT
% T:      total number of time units
% Y:      dependent variable (NTx1) where N is the total number of
%           individual units and NT the total number of observations
% X:      matrix of time-varying and time-invariant covariates (NTx(k+r))
%           where k denotes the number of time-varying and r the number of
%           time-invariant covariates, respectively.
% W:      normalized spatial weights matrix (NTxNT)


% Preliminaries
NT=length(Y); N=NT/T;
JT=ones(T); Q1=kron(JT/T, eye(N)); Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));

% Initial estimator
XX=X'*X; XY=X'*Y;
bols=XX\XY;
u=Y-X*bols;

% Estimate and standard error of rho and sigma_mu^2
[rho, rho_se,sgmu2]= gmhetre(u,N,T,W,Q0,Q1,eTIN,X);

% Spatial Cochrane Orcutt transformation
WY=W*Y; WX=W*X; Ys=Y-rho*WY; Xs=X-rho*WX; 

% Estimate transformed model
XsXs=Xs'*Xs; XsYs=Xs'*Ys;
beta=XsXs\XsYs;

u=Y-X*beta;
[SgHR]= sigmaHR(u,N,T,W,Q0,eTIN,rho);
Fmu=eTIN*Xs; Fnu=Xs;
Omg_mu=1/NT*sgmu2*(Fmu'*Fmu);
Omg_nu=1/NT*(Fnu'*SgHR*Fnu);
Omg=Omg_mu+Omg_nu;
P_res=(1/NT*XsXs)\eye(size(X,2));
VC=1/NT*(P_res'*Omg*P_res);


function [rho,rho_se,sgmu2]= gmhetre(u,N,T,W,Q0,Q1,eTIN,X);
% GM procedure for the spatial random effects model under heteroskedastic
% innovations
% OUTPUT
% rho:	  spatial autocorrelation parameter
% rho_se: standard error of the spatial autorcorrelation parameter
% sgmu2:  sigma_mu^2
% INPUT
% N:      total number of individual units 
% T:      total number of time units
% u: 	  estimated residual (NTx1)
% W:      normalized spatial weights matrix (NTxNT)
% Q0:     within-transformation matrix (NTxNT)
% Q1:     between-transformation matrix (NTxNT)
% eTIN:   eTIN=kron(ones(1,T),eye(N))
% X:      matrix of time-varying and time-invariant covariates (NTx(k+r))
%           where k denotes the number of time-varying and r the number of
%           time-invariant covariates, respectively.
% 

% Preliminaries
NT=N*T;
WpW=W(1:N,1:N)'*W(1:N,1:N); trWpW=trace(WpW);
ITWW=kron(eye(T),WpW);
Q1ITWW=Q1*ITWW;

% GM Elements
ubar=W*u;
ubarbar=W*ubar;

ga=zeros(4,1); GA=zeros(4,3);
ga(2,1)=1/(N*(T-1))*(ubar'*Q0*u);
ga(3,1)=1/N*(ubar'*Q1*ubar-trace(bsxfun(@times,Q1ITWW,u.^2)));

GA(2,1)=1/(N*(T-1))*(ubarbar'*Q0*u+ubar'*Q0*ubar);
GA(2,2)=-1/(N*(T-1))*(ubarbar'*Q0*ubar);
GA(2,3)=0;
GA(3,1)=2/N*(ubar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,(ubar.*u))));
GA(3,2)=-1/N*(ubarbar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,ubar.^2)));
GA(3,3)=(T-1)/N*trWpW;

gmfun0=@(aa)(ga(2,1)-GA(2,1:2)*[aa(1);aa(1)^2]);  
a0=0;

% Estimate for rho
[ai] = lsqnonlin(gmfun0,a0);
rho=ai(1,1); 

% Estimate for sigma_mu^2
aa13=ga(3,1)-GA(3,1)*rho-GA(3,2)*rho^2;  
sgmu2=aa13*N/((T-1)*trWpW);

% Standard error of rho
J=GA(2,1:2)*[1; 2*rho];
JI=inv(J);
P=(1/NT*(X')*X)\eye(size(X,2)); 
A2=Q0*W;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho*W(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-X'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*X;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
[SgHR]= sigmaHR(u,N,T,W,Q0,eTIN,rho);
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=R*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
cc=T/((T-1)^2*N);
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho_se=sqrt(Omrho);


function [beta, VC, rho1, rho1_se, rho2, rho2_se, rho3, rho3_se, sgmu2]=sht_het_multr(Y,XU,XC,ZU,ZC,W1,W2,W3,W,T); 
% Spatial Hausman-Taylor estimator under heteroskedastic innovations with
%   three spatial weights matrices
% OUTPUT
% beta:   parameter estimates on time-varying and time-invariant variables
% VC:     Variance-Covariance matrix of beta
% rho1, rho2, rho3:	spatial autocorrelation parameters
% rho1_se, rho2_se, rho3_se: standard errors of the spatial autorcorrelation parameters
% sgmu2:  sigma_mu^2
% INPUT
% T:  total number of time units
% Y:  dependent variable (NTx1) where N denotes the total number of
%       individual units and NT the total number of observations
% XU: matrix of time-varying exogenous covariates(NTxku) where ku denotes
%       the number of exogenous time-varying covariates 
% XC: matrix of time-varying endogenous covariates(NTxkc) where kc denotes
%       the number of endogenous time-varying covariates 
% ZU: matrix of time-invariant exogenous covariates(NTxru) where ru denotes
%       the number of exogenous time-invariant covariates 
% ZC: matrix of time-invariant endogenous covariates(NTxrc) where rc denotes
%       the number of endogenous time-invariant covariates 
% W1, W2, W3, W: normalized spatial weights matrices (NTxNT)


% Preliminaries
NT=length(Y); N=NT/T;
JT=ones(T); Q1=kron(JT/T, eye(N)); Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));

Xbar=Q1*[XU XC]; Ybar=Q1*Y; Xtil=Q0*[XU XC]; Ytil=Q0*Y;

% Step 1: Estimate bet0 with fixed effects 
XtilXtil=Xtil'*Xtil; XtilYtil=Xtil'*Ytil;
bet0=XtilXtil\XtilYtil;
dhat=Ybar-Xbar*bet0;

% Step 2: Estimate thet0 with 2SLS of dhat on Z using instrument set HI=[XU ZU];
HI=[XU ZU]; 
PHI=HI*inv(HI'*HI)*HI';
ZPHIZ=[ZU ZC]'*PHI*[ZU ZC]; ZPHIdhat=[ZU ZC]'*PHI*dhat;
thet0=ZPHIZ\ZPHIdhat;
[bethet0]=[bet0;thet0];
u=Y-[XU XC ZU ZC]*bethet0;

% Step 3: Estimate rho1, rho2, rho3 and sigma_mu^2
[rho1, rho1_se, rho2, rho2_se, rho3, rho3_se, sgmu2]= gmhetht_multr(u,N,T,W1,W2,W3,W,Q0,Q1,eTIN,XU,XC,ZU,ZC);

% Step 4: 
% Spatial Cochrane Orcutt transformation
R=eye(NT)-rho1*W1-rho2*W2-rho3*W3;
Ys=R*Y; XUs=R*XU; XCs=R*XC; ZUs=R*ZU; ZCs=R*ZC;

% Estimate transformed model
As=[XUs XCs ZUs ZCs]; Xstil=Q0*[XUs XCs]; XUsbar=Q1*XUs;
Hs=[Xstil XUsbar ZUs]; Hs2=Hs'*Hs; temp=(Hs/Hs2); PAHT=temp*Hs';
APA=As'*PAHT*As; APY=As'*PAHT*Ys;
beta=APA\APY;

u=Y-[XU XC ZU ZC]*beta;
[SgHR]= sigmaHR_multrho(u,T,Q0,eTIN,R);
Fnu=Hs; Fmu=eTIN*Hs; 
Omg_mu=1/NT*sgmu2*(Fmu'*Fmu);
Omg_nu=1/NT*(Fnu'*SgHR*Fnu);
Omg=Omg_mu+Omg_nu;
HAs=Hs'*As; AHs=As'*Hs;
P1=Hs2\HAs; temp=AHs/Hs2; P2=temp*(1/NT*HAs);
D=P1/P2;
VC=1/NT*(D'*Omg*D);

function [rho1, rho1_se, rho2, rho2_se, rho3, rho3_se, sgmu2]= gmhetht_multr(u,N,T,W1,W2,W3,W,Q0,Q1,eTIN,XU,XC,ZU,ZC);
% GM procedure for the spatial Hausman-Taylor model under heteroskedastic
% innovations with three spatial weights matrices
% OUTPUT
% rho1, rho2, rho3:	  spatial autocorrelation parameters
% rho1_se, rho2_se, rho3_se: standard error of the spatial autorcorrelation
% parameters
% INPUT
% N:      total number of individual units 
% T:      total number of time units
% u: 	  estimated residual (NTx1)
% Q0:     within-transformation matrix (NTxNT)
% Q1:     between-transformation matrix (NTxNT)
% eTIN:   eTIN=kron(ones(1,T),eye(N))
% XU:     matrix of time-varying exogenous covariates(NTxku) where 
%            ku denotes the number of exogenous time-varying covariates 
% XC:     matrix of time-varying endogenous covariates(NTxkc) where kc 
%           denotes the number of endogenous time-varying covariates 
% ZU:     matrix of time-invariant exogenous covariates(NTxru) where ru 
%           denotes the number of exogenous time-invariant covariates 
% ZC:     matrix of time-invariant endogenous covariates(NTxrc) where rc 
%           denotes the number of endogenous time-invariant covariates 
% W1, W2, W3, W: normalized spatial weights matrices (NTxNT)

% Preliminaries
NT=N*T; WpW=W(1:N,1:N)'*W(1:N,1:N); trWpW=trace(WpW); ITWW=kron(eye(T),WpW);
Q1ITWW=Q1*ITWW;

% GM Elements
u1bar=W1*u;
u1barbar=W1*u1bar;
u2bar=W2*u;
u2barbar=W2*u2bar;
u3bar=W3*u;
u3barbar=W3*u3bar;

ga1(1,1)=1/(N*(T-1))*(u1bar'*Q0*u);
GA1(1,1)=1/(N*(T-1))*(u1barbar'*Q0*u+u1bar'*Q0*u1bar);
GA1(1,2)=-1/(N*(T-1))*(u1barbar'*Q0*u1bar);

ga2(1,1)=1/(N*(T-1))*(u2bar'*Q0*u);
GA2(1,1)=1/(N*(T-1))*(u2barbar'*Q0*u+u2bar'*Q0*u2bar);
GA2(1,2)=-1/(N*(T-1))*(u2barbar'*Q0*u2bar);

ga3(1,1)=1/(N*(T-1))*(u3bar'*Q0*u);
GA3(1,1)=1/(N*(T-1))*(u3barbar'*Q0*u+u3bar'*Q0*u3bar);
GA3(1,2)=-1/(N*(T-1))*(u3barbar'*Q0*u3bar);

gmfun1=@(aa)(ga1-GA1*[aa(1);aa(1)^2]);                   
gmfun2=@(aa)(ga2-GA2*[aa(1);aa(1)^2]);  
gmfun3=@(aa)(ga3-GA3*[aa(1);aa(1)^2]);  
                    
ubar=W*u;
ubarbar=W*ubar;

ga=zeros(4,1); GA=zeros(4,3);
ga(2,1)=1/(N*(T-1))*(ubar'*Q0*u);
ga(3,1)=1/N*(ubar'*Q1*ubar-trace(bsxfun(@times,Q1ITWW,u.^2)));

GA(2,1)=1/(N*(T-1))*(ubarbar'*Q0*u+ubar'*Q0*ubar);
GA(2,2)=-1/(N*(T-1))*(ubarbar'*Q0*ubar);
GA(2,3)=0;
GA(3,1)=2/N*(ubar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,(ubar.*u))));
GA(3,2)=-1/N*(ubarbar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,ubar.^2)));
GA(3,3)=(T-1)/N*trWpW;

gmfun0=@(aa)(ga(2,1)-GA(2,1:2)*[aa(1);aa(1)^2]);  
a0=0;

% Estimates for rho1, rho2, rho3
rho1 = lsqnonlin(gmfun1,a0);
rho2 = lsqnonlin(gmfun2,a0);
rho3 = lsqnonlin(gmfun3,a0);

% Estimate for sigma_mu^2
[ai] = lsqnonlin(gmfun0,a0);
aa13=ga(3,1)-GA(3,1)*ai-GA(3,2)*ai^2;  
sgmu2=aa13*N/((T-1)*trWpW);

% Standard errors of rho1, rho2, rho3
Xtil=Q0*[XU XC]; XUbar=Q1*XU;
H=[Xtil XUbar ZU];
A=[XU XC ZU ZC];
HH=H'*H; HA=H'*A; AH=A'*H;
P1=HH\HA; temp=AH/HH; P2=temp*(1/NT*HA); P=P1/P2;
RR=eye(NT)-rho1*W1-rho2*W2-rho3*W3;
[SgHR]= sigmaHR_multrho(u,T,Q0,eTIN,RR);
cc=T/((T-1)^2*N);

% rho1
J=GA1*[1; 2*rho1]; JI=inv(J);
A2=Q0*W1;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho1*W1(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-A'*C2*u);
R0=R0\eye(N);                       
RI=kron(eye(T),R0);
Fnu=RI*H;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=RR*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho1_se=sqrt(Omrho);

% rho2
J=GA2*[1; 2*rho2]; JI=inv(J);
A2=Q0*W2;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho2*W2(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-A'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*H;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=RR*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho2_se=sqrt(Omrho);

% rho3
J=GA3*[1; 2*rho3]; JI=inv(J);
A2=Q0*W3;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho2*W2(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-A'*C2*u);
R0=R0\eye(N);                       
RI=kron(eye(T),R0);
Fnu=RI*H;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=RR*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho3_se=sqrt(Omrho);


function [SgHR]= sigmaHR_multrho(u,T,Q0,eTIN,R);
% INPUT
% T:      total number of time units
% u:      estimated residual (NTx1) where N is the total number of
%           individual units and NT the total number of observations
% Q0:     within-transformation matrix (NTxNT)
% eTIN:   eTIN=kron(ones(1,T),eye(N))
% R:      R=eye(NT)-rho1*W1-rho2*W2-rho3*W3 where rho1, rho2, and rho3 
%           denote the spatial autocorrelation parameter estimates and W1, W2,
%           and W3 the respective NTxNT normalized spatial weights matrices

ee=R*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=eTIN*nutil2;
tempB=kron(ones(T,1),temp);
tempF=1/((T-1)*(T-2))*tempB;
nutil2=T/(T-2)*nutil2;
nu2HR=nutil2-tempF;
[SgHR]=diag(nu2HR);



function [beta, VC, rho1, rho1_se, rho2, rho2_se, rho3, rho3_se]=sfe_het_multr(Y,X,W1,W2,W3,T); 
% Spatial fixed effects estimator under heteroskedastic innovations with
%   three spatial weights matrices
% OUTPUT
% beta:                      parameter estimates
% VC:                        Variance-Covariance matrix of beta
% rho1, rho2, rho3:	         spatial autocorrelation parameters
% rho1_se, rho2_se, rho3_se: standard errors of the spatial autorcorrelation parameters
% INPUT
% T:            total number of time units
% Y:            dependent variable (NTx1) where N is the total number of
%               individual units and NT the total number of observations
% X:            matrix of time-varying covariates (NTxk) where k is the 
%                 number of time-varying covariates
% W1, W2, W3,:  normalized spatial weights matrices (NTxNT)


% Preliminaries
NT=length(Y); N=NT/T;
JT=ones(T); Q1=kron(JT/T, eye(N)); Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));

Xtil=Q0*X; Ytil=Q0*Y;
XtilXtil=Xtil'*Xtil; XtilYtil=Xtil'*Ytil;

%Initial estimator
bfe=XtilXtil\XtilYtil;
u=Y-X*bfe;

%Estimate and standard error of rho1, rho2, rho3
[rho1, rho1_se, rho2, rho2_se, rho3, rho3_se]= gmhetfe_multr(u,N,T,W1,W2,W3,Q0,eTIN,X);

%Spatial Cochrane Orcutt and within transformation
R=eye(NT)-rho1*W1-rho2*W2-rho3*W3;
Ys=R*Y; Xs=R*X; 
Xstil=Q0*Xs; Ystil=Q0*Ys;

% Estimate transformed model
XstilXstil=Xstil'*Xstil; XstilYstil=Xstil'*Ystil;
beta=XstilXstil\XstilYstil;

u=Y-X*beta;
[SgHR]= sigmaHR_multrho(u,T,Q0,eTIN,R);
Fnu=Q0*Xs;
Omg=1/NT*(Fnu'*SgHR*Fnu);
P=(1/NT*(Fnu'*Fnu))\eye(size(X,2));
VC=1/NT*(P'*Omg*P);

function [rho1, rho1_se, rho2, rho2_se, rho3, rho3_se]= gmhetfe_multr(u,N,T,W1,W2,W3,Q0,eTIN,X);
% GM procedure for the spatial fixed effects model under heteroskedastic
% innovations with three spatial weights matrices
% OUTPUT
% rho1, rho2, rho3:	  spatial autocorrelation parameters
% rho1_se, rho2_se, rho3_se: standard errors of the spatial autorcorrelation
% parameters
% INPUT
% N:          total number of individual units 
% T:          total number of time units
% u: 	      estimated residual (NTx1)
% W1, W2, W3: normalized spatial weights matrices (NTxNT)
% Q0:         within-transformation matrix (NTxNT)
% eTIN:       eTIN=kron(ones(1,T),eye(N))
% X:          matrix of time-varying covariates (NTxk) where k is the 
%               number of time-varying covariates

% Preliminaries
NT=N*T;

% GM Elements 
u1bar=W1*u;
u1barbar=W1*u1bar;
u2bar=W2*u;
u2barbar=W2*u2bar;
u3bar=W3*u;
u3barbar=W3*u3bar;

ga1(1,1)=1/(N*(T-1))*(u1bar'*Q0*u);
GA1(1,1)=1/(N*(T-1))*(u1barbar'*Q0*u+u1bar'*Q0*u1bar);
GA1(1,2)=-1/(N*(T-1))*(u1barbar'*Q0*u1bar);

ga2(1,1)=1/(N*(T-1))*(u2bar'*Q0*u);
GA2(1,1)=1/(N*(T-1))*(u2barbar'*Q0*u+u2bar'*Q0*u2bar);
GA2(1,2)=-1/(N*(T-1))*(u2barbar'*Q0*u2bar);

ga3(1,1)=1/(N*(T-1))*(u3bar'*Q0*u);
GA3(1,1)=1/(N*(T-1))*(u3barbar'*Q0*u+u3bar'*Q0*u3bar);
GA3(1,2)=-1/(N*(T-1))*(u3barbar'*Q0*u3bar);

gmfun1=@(aa)(ga1-GA1*[aa(1);aa(1)^2]);                   
gmfun2=@(aa)(ga2-GA2*[aa(1);aa(1)^2]);  
gmfun3=@(aa)(ga3-GA3*[aa(1);aa(1)^2]);  
                    
%Estimates for rho1, rho2, rho3
a0=0;
rho1 = lsqnonlin(gmfun1,a0);
rho2 = lsqnonlin(gmfun2,a0);
rho3 = lsqnonlin(gmfun3,a0);

% Standard errors of rho1, rho2, rho3
Xfe=Q0*X;
P=(1/(N*T)*(Xfe')*Xfe)\eye(size(Xfe,2)); 
RR=eye(NT)-rho1*W1-rho2*W2-rho3*W3;
[SgHR]= sigmaHR_multrho(u,T,Q0,eTIN,RR);
cc=T/((T-1)^2*N);

% rho1
J=GA1*[1; 2*rho1]; JI=inv(J);
A2=Q0*W1;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*(A2A2t);
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho1*W1(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-Xfe'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*Xfe;
anu=Fnu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
ee=RR*u; nutil=Q0*ee; nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho1_se=sqrt(Omrho);

% rho2
J=GA2*[1; 2*rho2]; JI=inv(J);
A2=Q0*W2;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*(A2A2t);
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho2*W2(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-Xfe'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*Xfe;
anu=Fnu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
ee=RR*u; nutil=Q0*ee; nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho2_se=sqrt(Omrho);

% rho3
J=GA3*[1; 2*rho3]; JI=inv(J);
A2=Q0*W3;
A2A2t=A2+A2';
A2nu=1/(2*(T-1))*(A2A2t);
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho3*W3(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-Xfe'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*Xfe;
anu=Fnu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
ee=RR*u; nutil=Q0*ee; nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho3_se=sqrt(Omrho);


function [beta, VC, rho1, rho1_se, rho2, rho2_se, rho3, rho3_se, sgmu2]=sre_het_multr(Y,X,W1,W2,W3,W,T); 
% Spatial random effects estimator under heteroskedastic innovations with
%       three spatial weights matrices
% OUTPUT
% beta:   parameter estimates on time-varying and time-invariant variables
% VC:     Variance-Covariance matrix of beta
% rho1, rho2, rho3:	         spatial autocorrelation parameters
% rho1_se, rho2_se, rho3_se: standard errors of the spatial autorcorrelation parameters
% sgmu2:  sigma_mu^2
% INPUT
% T:     total number of time units
% Y:     dependent variable (NTx1) where N is the total number of
%          individual units and NT the total number of observations
% X:     matrix of time-varying and time-invariant covariates (NTx(k+r))
%          where k denotes the number of time-varying and r the number
%          of time-invariant covariates, respectively.
% W1,W2,W3,W:  normalized spatial weights matrices (NTxNT)


% Preliminaries
[NT]=length(Y); N=NT/T;
JT=ones(T); Q1=kron(JT/T, eye(N)); Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));

% Initial estimator
XX=X'*X; XY=X'*Y;
bols=XX\XY;
u=Y-X*bols;

% Estimate and standard error of rho and sigma_mu^2
[rho1, rho1_se, rho2, rho2_se, rho3, rho3_se, sgmu2]= gmhetre_multr(u,N,T,W1,W2,W3,W,Q0,Q1,eTIN,X);

% Spatial Cochrane Orcutt transformation
R=eye(NT)-rho1*W1-rho2*W2-rho3*W3;
Ys=R*Y; Xs=R*X; 

% Estimate transformed model
XsXs=Xs'*Xs; XsYs=Xs'*Ys;
beta=XsXs\XsYs;
u=Y-X*beta;
[SgHR]= sigmaHR_multrho(u,T,Q0,eTIN,R);
Fmu=eTIN*Xs; Fnu=Xs;
Omg_mu=1/NT*sgmu2*(Fmu'*Fmu);
Omg_nu=1/NT*(Fnu'*SgHR*Fnu);
Omg=Omg_mu+Omg_nu;
P=(1/NT*XsXs)\eye(size(X,2));
VC=1/NT*(P'*Omg*P);


function [rho1, rho1_se, rho2, rho2_se, rho3, rho3_se, sgmu2]= gmhetre_multr(u,N,T,W1,W2,W3,W,Q0,Q1,eTIN,X);
% GM procedure for the spatial random effects model under heteroskedastic
% innovations with three spatial weights matrices
% OUTPUT
% rho1, rho2, rho3:	  spatial autocorrelation parameters
% rho1_se, rho2_se, rho3_se: standard errors of the spatial autorcorrelation parameters
% sgmu2:  sigma_mu^2
% INPUT
% N:      total number of individual units 
% T:      total number of time units
% u: 	  estimated residual (NTx1)
% W1, W2, W3, W:   normalized spatial weights matrices (NTxNT)
% Q0:     within-transformation matrix (NTxNT)
% Q1:     between-transformation matrix (NTxNT)
% eTIN:   eTIN=kron(ones(1,T),eye(N))
% X:      matrix of time-varying and time-invariant covariates (NTx(k+r))
%           where k denotes the number of time-varying and r the number of
%           time-invariant covariates, respectively.

% Preliminaries
NT=N*T; WpW=W(1:N,1:N)'*W(1:N,1:N); trWpW=trace(WpW);
ITWW=kron(eye(T),WpW); Q1ITWW=Q1*ITWW;

% GM Elements
u1bar=W1*u;
u1barbar=W1*u1bar;
u2bar=W2*u;
u2barbar=W2*u2bar;
u3bar=W3*u;
u3barbar=W3*u3bar;

ga1(1,1)=1/(N*(T-1))*(u1bar'*Q0*u);
GA1(1,1)=1/(N*(T-1))*(u1barbar'*Q0*u+u1bar'*Q0*u1bar);
GA1(1,2)=-1/(N*(T-1))*(u1barbar'*Q0*u1bar);

ga2(1,1)=1/(N*(T-1))*(u2bar'*Q0*u);
GA2(1,1)=1/(N*(T-1))*(u2barbar'*Q0*u+u2bar'*Q0*u2bar);
GA2(1,2)=-1/(N*(T-1))*(u2barbar'*Q0*u2bar);

ga3(1,1)=1/(N*(T-1))*(u3bar'*Q0*u);
GA3(1,1)=1/(N*(T-1))*(u3barbar'*Q0*u+u3bar'*Q0*u3bar);
GA3(1,2)=-1/(N*(T-1))*(u3barbar'*Q0*u3bar);

gmfun1=@(aa)(ga1-GA1*[aa(1);aa(1)^2]);                   
gmfun2=@(aa)(ga2-GA2*[aa(1);aa(1)^2]);  
gmfun3=@(aa)(ga3-GA3*[aa(1);aa(1)^2]);  
                    
ubar=W*u;
ubarbar=W*ubar;

ga=zeros(4,1); GA=zeros(4,3);
ga(2,1)=1/(N*(T-1))*(ubar'*Q0*u);
ga(3,1)=1/N*(ubar'*Q1*ubar-trace(bsxfun(@times,Q1ITWW,u.^2)));

GA(2,1)=1/(N*(T-1))*(ubarbar'*Q0*u+ubar'*Q0*ubar);
GA(2,2)=-1/(N*(T-1))*(ubarbar'*Q0*ubar);
GA(2,3)=0;
GA(3,1)=2/N*(ubar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,(ubar.*u))));
GA(3,2)=-1/N*(ubarbar'*Q1*ubarbar-trace(bsxfun(@times,Q1ITWW,ubar.^2)));
GA(3,3)=(T-1)/N*trWpW;

gmfun0=@(aa)(ga(2,1)-GA(2,1:2)*[aa(1);aa(1)^2]);  
a0=0;

% Estimate for rho
a0=0;
rho1 = lsqnonlin(gmfun1,a0);
rho2 = lsqnonlin(gmfun2,a0);
rho3 = lsqnonlin(gmfun3,a0);

% Estimate for sigma_mu^2
[ai] = lsqnonlin(gmfun0,a0);
aa13=ga(3,1)-GA(3,1)*ai-GA(3,2)*ai^2; 
sgmu2=aa13*N/((T-1)*trWpW);

% Standard error of rho1, rho2, rho3
P=(1/NT*X'*X)\eye(size(X,2)); 
RR=eye(NT)-rho1*W1-rho2*W2-rho3*W3;
[SgHR]= sigmaHR_multrho(u,T,Q0,eTIN,RR);
cc=T/((T-1)^2*N);

% Rho1
J=GA1*[1; 2*rho1]; JI=inv(J);
A2=Q0*W1; A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho1*W1(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-X'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*X;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=RR*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho1_se=sqrt(Omrho);

% Rho2
J=GA2*[1; 2*rho2]; JI=inv(J);
A2=Q0*W2; A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho2*W2(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-X'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*X;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=RR*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho2_se=sqrt(Omrho);

% Rho3
J=GA3*[1; 2*rho3]; JI=inv(J);
A2=Q0*W3; A2A2t=A2+A2';
A2nu=1/(2*(T-1))*A2A2t;
A2nu2=A2nu.*A2nu;
R0=eye(N)-rho3*W3(1:N,1:N);
R=kron(eye(T),R0);
C2=1/(2*(T-1))*(R'*A2A2t*R);
alpha=2/N*(-X'*C2*u);
R0=R0\eye(N);                        
RI=kron(eye(T),R0);
Fnu=RI*X;
Fmu=eTIN*Fnu;
anu=Fnu*P*alpha; 
amu=Fmu*P*alpha;
Enu=1/NT*(anu'*SgHR*anu);
Emu=1/N*sgmu2*(amu'*amu);
ee=RR*u;
nutil=Q0*ee;
nutil2=nutil.^2;
temp=nutil2'*A2nu2*nutil2;
Enuc=cc*temp;
E=Enu+Emu+Enuc;
Omrho=JI*E*JI;
Omrho=1/N*Omrho;
rho3_se=sqrt(Omrho);

% Tests
function [HTS, df, pv]= hausmantest_het(bres,VCres,rhore,bfes,VCfes,rhofe,Y,X,W,N,T);
% Spatial Hausman Test under heteroskedastic innovations
% OUTPUT
% HTS: Hausman test statistic
% df:  degrees of freedom
% pv:  p-value
% INPUT
% bres:  beta of the spatial random effects estimator
% VCres: variance-covariance matrix of the spatial random effects estimator
% bfes:  beta of the spatial fixed effects estimator
% VCfes: variance-covariance matrix of the spatial fixed effects estimator
% rhore: spatial autocorrelation parameter of the spatial random effects
%        estimator
% rhofe: spatial autocorrelation parameter of the spatial fixed effects
%        estimator
% N:     total number of individual units  
% T:     total number of time units
% Y:     dependent variable (NTx1) 
% X:     matrix of time-varying and time-invariant covariates (NTxp) where
%           p denotes the number of covariates
% W:     normalized spatial weights matrix (NTxNT)

% Preliminaries
NT=N*T; k=size(bfes,1);
JT=ones(T); Q1=kron(JT/T, eye(N)); 
Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N)); WX=W*X; 

% VC matrix
Omega=zeros(2*k,2*k);
Omega(1:k,1:k)=VCres(1:k,1:k);
Omega(k+1:2*k,k+1:2*k)=VCfes(1:k,1:k);
Fnu_res=X-rhore*WX;
Fnu_res=Fnu_res(1:NT,1:k);
P_res=1/NT*(Fnu_res'*Fnu_res)\eye(k);
Xs=X-rhofe*WX;
Fnu_fes=Q0*Xs;
P_fes=1/NT*(Fnu_fes'*Fnu_fes)\eye(k);
ufes=Y-X*bfes;
[SgHRfes]= sigmaHR(ufes,N,T,W,Q0,eTIN,rhofe);
VCresfes=1/NT*P_res'*(1/NT*Fnu_res'*SgHRfes*Fnu_fes)*P_fes;
Omega(1:k,k+1:2*k)=VCresfes;
Omega(k+1:2*k,1:k)=VCresfes';

% Test statistic
bres=bres(1:k,1);
q=bres-bfes;
B=[eye(k) -eye(k)];
HTS=q'*pinv(B*(Omega)*B')*q;
if HTS<0;
temp2=VCfes(1:k,1:k)-VCres(1:k,1:k);
HTS=q'*pinv(temp2)*q;
end
df=k;
pv=1-chi2cdf(HTS,df);


function [HTTS, df, pv]= hausmantaylortest_het(bhts,VChts,rhoht,bfes,VCfes,rhofe,Y,XU,XC,ZU,ZC,W,N,T);
% Spatial Hausman-Taylor Test under heteroskedastic innovations
% OUTPUT
% HTTS: Hausman-Taylor test statistic
% df:   degrees of freedom
% pv:   p-value
% INPUT
% bhts:  beta of the spatial Hausman-Taylor estimator
% VChts: variance-covariance matrix of the spatial Hausman-Taylor estimator
% bfes:  beta of the spatial fixed effects estimator
% VCfes: variance-covariance matrix of the spatial fixed effects estimator
% rhoht: spatial autocorrelation parameter of the spatial Hausman-Taylor
%        estimator
% rhofe: spatial autocorrelation parameter of the spatial fixed effects
%        estimator
% N:     total number of individual units  
% T:     total number of time units
% Y:     dependent variable (NTx1) 
% XU:    matrix of time-varying exogenous covariates(NTxku)
% XC:    matrix of time-varying endogenous covariates(NTxkc)
% ZU:    matrix of time-invariant exogenous covariates(NTxru)
% ZC:    matrix of time-invariant endogenous covariates(NTxrc)
% W:     normalized spatial weights matrix (NTxNT)

% Preliminaries
NT=N*T; k=size(bfes,1);
JT=ones(T); Q1=kron(JT/T, eye(N)); 
Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));
ku=size(XU,2);
rc=size(ZC,2);

% VC matrix
Omega=zeros(2*k,2*k);
Omega(1:k,1:k)=VChts(1:k,1:k);
Omega(k+1:2*k,k+1:2*k)=VCfes(1:k,1:k);
WXU=W*XU; WXC=W*XC; WZU=W*ZU; WZC=W*ZC;
XUs=XU-rhoht*WXU; XCs=XC-rhoht*WXC;
ZUs=ZU-rhoht*WZU; ZCs=ZC-rhoht*WZC; 
As=[XUs XCs ZUs ZCs];
Xstil=Q0*[XUs XCs]; XUsbar=Q1*XUs;
HHTs=[Xstil XUsbar ZUs]; 
HHTs2=HHTs'*HHTs; HA=HHTs'*As; AH=As'*HHTs;
P1=HHTs2\HA; temp=AH/HHTs2; P2=temp*(1/NT*HA);
P_hts=P1/P2;
P_hts=P_hts(1:k,1:k);
Fnu_hts=HHTs(1:NT,1:k);
X=[XU XC]; WX=[WXU WXC];
Xs=X-rhofe*WX;
Fnu_fes=Q0*Xs;
P_fes=1/NT*(Fnu_fes'*Fnu_fes)\eye(k);
ufes=Y-[XU XC]*bfes;
[SgHRfes]= sigmaHR(ufes,N,T,W,Q0,eTIN,rhofe);
VChtsfes=1/NT*P_hts'*(1/NT*Fnu_hts'*SgHRfes*Fnu_fes)*P_fes;
Omega(1:k,k+1:2*k)=VChtsfes;
Omega(k+1:2*k,1:k)=VChtsfes';

% Test statistic
bhts=bhts(1:k,1);
q=bhts-bfes;
B=[eye(k) -eye(k)];
HTTS=q'*pinv(B*(Omega)*B')*q;
if HTTS<0;
temp2=VCfes(1:k,1:k)-VChts(1:k,1:k);
HTTS=q'*pinv(temp2)*q;
end
df=ku-rc;
pv=1-chi2cdf(HTTS,df);

function [HTS, df, pv]= hausmantest_het_multr(bres,VCres,Rre,bfes,VCfes,Rfe,Y,X,N,T);
% Spatial Hausman Test under heteroskedastic innovations with three spatial weights
% matrices
% OUTPUT
% HTS: Hausman test statistic
% df:  degrees of freedom
% pv:  p-value
% INPUT
% bres:  beta of the spatial random effects estimator
% VCres: variance-covariance matrix of the spatial random effects estimator
% bfes:  beta of the spatial fixed effects estimator
% VCfes: variance-covariance matrix of the spatial fixed effects estimator
% N:     total number of individual units  
% T:     total number of time units
% Y:     dependent variable (NTx1) 
% X:     matrix of time-varying and time-invariant covariates (NTxp) where
%           p denotes the number of covariates
% Rfe:   Rfe=eye(NT)-rhofe1*W1-rhofe2*W2-rhofe3*W3 where rhofe1, rhofe2, 
%           and rhofe3 denote spatial autocorrelation parameter estimates 
%           and W1, W2, and W3 the respective NTxNT normalized spatial weights matrices


% Preliminaries
NT=N*T; k=size(bfes,1);
JT=ones(T); Q1=kron(JT/T, eye(N)); 
Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N)); 

% VC matrix
Omega=zeros(2*k,2*k);
Omega(1:k,1:k)=VCres(1:k,1:k);
Omega(k+1:2*k,k+1:2*k)=VCfes(1:k,1:k);
Fnu_res=Rre*X;
Fnu_res=Fnu_res(1:NT,1:k);
P_res=1/NT*(Fnu_res'*Fnu_res)\eye(k);
Xs=Rfe*X;
Fnu_fes=Q0*Xs;
P_fes=1/NT*(Fnu_fes'*Fnu_fes)\eye(k);
ufes=Y-X*bfes;
[SgHRfes]= sigmaHR_multrho(ufes,T,Q0,eTIN,Rfe);
VCresfes=1/NT*P_res'*(1/NT*Fnu_res'*SgHRfes*Fnu_fes)*P_fes;
Omega(1:k,k+1:2*k)=VCresfes;
Omega(k+1:2*k,1:k)=VCresfes';

% Test statistic
bres=bres(1:k,1);
q=bres-bfes;
B=[eye(k) -eye(k)];
HTS=q'*pinv(B*(Omega)*B')*q;
if HTS<0;
temp2=VCfes(1:k,1:k)-VCres(1:k,1:k);
HTS=q'*pinv(temp2)*q;
end
[df]=k;
[pv]=1-chi2cdf(HTS,df);


function [HTTS, df, pv]= hausmantaylortest_het_multr(bhts,VChts,Rht,bfes,VCfes,Rfe,Y,XU,XC,ZU,ZC,N,T);
% Spatial Hausman-Taylor Test under heteroskedastic innovations three spatial weights
% matrices
% OUTPUT
% HTTS: Hausman-Taylor test statistic
% df:   degrees of freedom
% pv:   p-value
% INPUT
% bhts:  beta of the spatial Hausman-Taylor estimator
% VChts: variance-covariance matrix of the spatial Hausman-Taylor estimator
% bfes:  beta of the spatial fixed effects estimator
% VCfes: variance-covariance matrix of the spatial fixed effects estimator
% N:     total number of individual units  
% T:     total number of time units
% Y:     dependent variable (NTx1) 
% XU:    matrix of time-varying exogenous covariates(NTxku)
% XC:    matrix of time-varying endogenous covariates(NTxkc)
% ZU:    matrix of time-invariant exogenous covariates(NTxru)
% ZC:    matrix of time-invariant endogenous covariates(NTxrc)
% Rfe:   Rfe=eye(NT)-rhofe1*W1-rhofe2*W2-rhofe3*W3 where rhofe1, rhofe2, 
%           and rhofe3 denote spatial autocorrelation parameter estimates 
%           and W1, W2, and W3 the respective NTxNT normalized spatial weights matrices

% Preliminaries
NT=N*T; k=size(bfes,1);
JT=ones(T); Q1=kron(JT/T, eye(N)); 
Q0=eye(NT)-Q1;
eTIN=kron(ones(1,T),eye(N));
ku=size(XU,2);
rc=size(ZC,2);

% VC matrix
Omega=zeros(2*k,2*k);
Omega(1:k,1:k)=VChts(1:k,1:k);
Omega(k+1:2*k,k+1:2*k)=VCfes(1:k,1:k);
XUs=Rht*XU; XCs=Rht*XC; ZUs=Rht*ZU; ZCs=Rht*ZC; 
As=[XUs XCs ZUs ZCs];
Xstil=Q0*[XUs XCs]; XUsbar=Q1*XUs;
HHTs=[Xstil XUsbar ZUs]; 
HHTs2=HHTs'*HHTs; HA=HHTs'*As; AH=As'*HHTs;
P1=HHTs2\HA; temp=AH/HHTs2; P2=temp*(1/NT*HA);
P_hts=P1/P2;
P_hts=P_hts(1:k,1:k);
Fnu_hts=HHTs(1:NT,1:k);
Xs=Rfe*[XU XC];
Fnu_fes=Q0*Xs;
P_fes=1/NT*(Fnu_fes'*Fnu_fes)\eye(k);
ufes=Y-[XU XC]*bfes;
[SgHRfes]= sigmaHR_multrho(ufes,T,Q0,eTIN,Rfe);
VChtsfes=1/NT*P_hts'*(1/NT*Fnu_hts'*SgHRfes*Fnu_fes)*P_fes;
Omega(1:k,k+1:2*k)=VChtsfes;
Omega(k+1:2*k,1:k)=VChtsfes';

% Test statistic
bhts=bhts(1:k,1);
q=bhts-bfes;
B=[eye(k) -eye(k)];
HTTS=q'*pinv(B*(Omega)*B')*q;
if HTTS<0;
temp2=VCfes(1:k,1:k)-VChts(1:k,1:k);
HTTS=q'*pinv(temp2)*q;
end
[df]=ku-rc;
[pv]=1-chi2cdf(HTTS,df);


