function [Qtests,dfQtests,pQtests,Ttests,pTtests]=test_condate(Y,X,T,r);
% Function that tests Zero Conditional ATE, Constant Conditional ATE, 
% and Zero ATE hypotheses as in Crump, R. K.,V. J. Hotz, G. W. Imbens 
% and O. A. Mitnik (2008) "Nonaparametric Tests for Treatment Effect 
% Heterogeneity", Review of Economics and Statistics, 
% August, 90(3): 389-405.
%
%
% Inputs:
%    Y      = Outcome vector over which ATE will be tested
%    X      = Matrix of conditioning covariates 
%    T      = Treatment indicator vector
%    r      = Indicator variable to indicate which var-cov matrix 
%             used to estimate the tests. If:  
%             =1 heteroskedastic-robust var-cov matrix used 
%             =0 homoskedastic var-cov matrix used 
%                     
% Outputs:
% Qtests    = Chi-Squared test statistics of the zero conditional ATE,
%             constant conditional ATE and zero ATE null hypotheses
%             (1x3 row vector) in this order
% dfQtests  = Degrees of freedom for statistics in Qtests
% pQtests   = p-value associated with statistics in Qtests
% Ttests    = Normal test statistic of the zero conditional ATE, 
%             constant conditional ATE and zero ATE null hypotheses
%             (1x3 row vector) in this order
% pTtests   = p-value associated with statistics in Ttests (one-sided 
%             except the one for zero ATE)


[N,K]=size(X);
N1=sum(T);
N0=sum(1-T);

% Estimate regression coefficients for controls and treated
% Makes sure that the same variables that are used in estimating the regression function
% for controls, are used for treated (i.e. re-estimates until the variables are the same)
w=ones(N,1);
xcols=[1:K];
fl=0;
while fl==0;
	[b0,se_hom0,se_rob0,sigma20,var_hom0,var_rob0,Kxx0,inx0]=wols_rob(Y(T==0),[ones(N0,1),X(T==0,xcols)],w(T==0));
	[b1,se_hom1,se_rob1,sigma21,var_hom1,var_rob1,Kxx1,inx1]=wols_rob(Y(T==1),[ones(N1,1),X(T==1,xcols)],w(T==1));
	if Kxx0==Kxx1,
		b0=b0(inx0==1);
		b1=b1(inx1==1);
		K=Kxx0-1;
		fl=1;
	else,
		inx00=zeros(K,1);
		inx00(xcols')=inx0(2:end);
		inx11=zeros(K,1);
		inx11(xcols')=inx1(2:end);
		xcols=find(inx00'==1 & inx11'==1); % Variables to be used in both regression functions
	end,
end;

% Calculate Q test for zero conditional ATE hypothesis
Vhat_rob=var_rob0+var_rob1;
Vhat_hom=var_hom0+var_hom1;

if r==1;
	Vhat=Vhat_rob;
elseif r==0;
	Vhat=Vhat_hom;
else;
	disp('Error: you need to specify whether to use homoskedastic or robust var-cov matrix')
end;

Qzcte=(b1-b0)'*inv(Vhat)*(b1-b0); % statistic 
dfQzcte=K+1; % dof
pQzcte=1-chi2cdf(Qzcte,dfQzcte); % p-value

% Calculate T test for zero conditional ATE hypothesis
Tzcte=(Qzcte-(K+1))/sqrt(2*(K+1)); % statistic
pTzcte=1-normcdf(Tzcte); % p-value

% Calculate Q test for constant conditional ATE hypothesis
if r==1;
	Vhat=var_rob0(2:K+1,2:K+1)+var_rob1(2:K+1,2:K+1);
elseif r==0;
	Vhat=var_hom0(2:K+1,2:K+1)+var_hom1(2:K+1,2:K+1);
end;
Qccte=(b1(2:K+1)-b0(2:K+1))'*inv(Vhat)*(b1(2:K+1)-b0(2:K+1)); % statistic 
dfQccte=K; % dof
pQccte=1-chi2cdf(Qccte,dfQccte); % p-value

% Calculate T test for constant conditional ATE hypothesis
Tccte=(Qccte-K)/sqrt(2*K); % statistic
pTccte=1-normcdf(abs(Tccte)); % p-value (1-sided)

% Calculate T test for zero ATE (difference of means)
Tzate=(mean(Y(T==1))-mean(Y(T==0)))/sqrt((var(Y(T==1))/N1)+(var(Y(T==0))/N0));
pTzate=(1-normcdf(abs(Tzate)))*2; % p-value
Qzate=Tzate^2;
pQzate=1-chi2cdf(Qzate,1); % p-value

% Collect results (all are 1x3 row vectors with 1st, 2nd and 3rd column representing the 
% corresponding result for the zero conditional ATE, constant conditional ATE and
% zero ATE null hypotheses in this order
Qtests=[Qzcte,Qccte,Qzate];
dfQtests=[dfQzcte,dfQccte,1];
pQtests=[pQzcte,pQccte,pQzate];
Ttests=[Tzcte,Tccte,Tzate];
pTtests=[pTzcte,pTccte,pTzate];


function [b,se_hom,se_rob,sigma2,var_hom,var_rob,Kxx,inx]=wols_rob(y,x,w);

% Inputs:
%    y          = Dependent variable vector
%    X          = Matrix of covariates 
%    w          = Vector of weights for WOLS (if vector of ones, estimates OLS)
                     
% Outputs (coeff and se vectors filled with zeros for variables dropped due to colinearity):
%    b          = OLS or WOLS coefficients (depending on W vector)
%    se_hom     = Homoskedastic standard errors
%    se_rob     = Heteroskedastic robust standard errors (Eicker-White)
%    sigma2     = Sigma squared from homoskedastic SE formula
%    var_hom    = Homoskedastic var-cov matrix
%    var_rob    = Heteroskedastic robust var-cov matrix
% NOTE: both var-cov matrices will have the size of the non-dropped vector of coefficients due to colinearity
%    Kxx        = Number of coefficients estimated after dropping colinear variables (including constant 
%                 if there's one)
%    inx        = Vector of ones (for variables in the regression) and zeros (for variables dropped due
%                 to colinearity)

cc=0.001;   % cutoff point for eigenvalues of covariance matrix of regressors
[N,K]=size(x);
xx=x;
inx=ones(K,1);
if min(eig(xx'*xx))<cc,
   inx=[1;zeros(K-1,1)];
   xx=x(:,1);
   for k=2:K,
       min(eig([xx,x(:,k)]'*[xx,x(:,k)]));
       if min(eig([xx,x(:,k)]'*[xx,x(:,k)]))>cc,
          xx=[xx,x(:,k)];
          inx(k,1)=1;
          end,
       end,
    end,   
[N,Kxx]=size(xx);

yw=y.*sqrt(w);
xxw=xx.*(sqrt(w)*ones(1,Kxx));

% Estimated coefficients (zeros for coefficients not estimated because of colinearity)
beta=inv(xxw'*xxw)*(xxw'*yw);
b=zeros(K,1);
b(inx==1,1)=beta;

% Homoskedastic standard errors and var-cov matrix
eps=y-xx*beta;
sigma2=eps'*eps/(N-Kxx);
var_hom=inv(xx'*xx)*sigma2;
se1=sqrt(diag(var_hom));
se_hom=zeros(K,1);
se_hom(inx==1,1)=se1;

% Heteroskedastic robust standard errors and var-cov matrix
epsx=xx.*(eps*ones(1,Kxx));
var_rob=inv(xx'*xx)*(epsx'*epsx)*inv(xx'*xx)*N/(N-sum(inx));  % Multiply by N/N-K (true K) for finite-sample correction
%var_rob=inv(xx'*xx)*(epsx'*epsx)*inv(xx'*xx);  % Robust var-cov matrix without finite-sample correction (asymptotic formula)
se2=sqrt(diag(var_rob));
se_rob=zeros(K,1);
se_rob(inx==1,1)=se2;