function [result,epshat,se_epshat]=heterotest(X,GG,r,maxcvH,skips,DEMEAN);

 if DEMEAN ==2;
 x=standard(X);
 G=standard(GG);
 end;

if DEMEAN ==1;
 x=demean(X);
 G=demean(GG);
 end;
if DEMEAN==0;
  x=X;
  G=GG;
  end;




[T,N]=size(x);
nz=cols(G);
skip=ceil(skips*T);
up=T-skip;
dn=skip+1;


Ntest=zeros(T,nz);
Mtest=zeros(T,nz);
Ttest=zeros(1,nz);
Ftest=zeros(1,nz);
R2=zeros(1,nz);
se=zeros(T,nz);
CIplus=zeros(T,nz);
CIminus=zeros(T,nz);
ci=zeros(T,nz);
maxci=zeros(T,nz);
cover=zeros(1,nz);
maxcover=zeros(1,nz);
maxtest=zeros(1,nz);
abstest=zeros(1,nz);
epshat=zeros(T,nz);
NIR=zeros(1,nz);
sig2eps=zeros(1,nz);




xx=x*x';                 
[A1,A2,A3]=svd(xx);

fhat=A1(:,1:r)*sqrt(T);

lamhat=x'*fhat/T;
VNT=A2(1:r,1:r)/(N*T);
chat=fhat*lamhat';
ehat=x-chat;
ehat2=ehat.*ehat;

for i=1:r;
if fhat(1,i)<0; fhat(:,i)=-fhat(:,i);end;
  end;


reg= fhat;


% for each observed factor G(:,j), regress them on the estimated factors.
% check if the fitted value is close to the actual G(:,j)


for j=1:nz;
delta=reg\G(:,j);                                % delta is r by 1
Ghat=reg*delta;
epshat(:,j)=G(:,j)-Ghat;
R2(j)=Ghat'*Ghat/(G(:,j)'*G(:,j));
Ftest(j)=((T-r)/r)*R2(j)/(1-R2(j));
%pval(j)=1-cdf('F',Ftest(j),r,T-r);
sig2eps(j)=epshat(:,j)'*epshat(:,j)/T;



for t=1:T;
dum=0;
for i=1:N;
  dum=dum+lamhat(i,:)'*lamhat(i,:)*ehat2(t,i);
end;

VF=dum/N;                                        % VF is r by r
VCV= (delta'*inv(VNT)*VF*inv(VNT)*delta)*(N*T)/(N*T-r*(N+T));  
se(t,j)=sqrt(VCV/N);
Ntest(t,j)=(Ghat(t)-G(t,j))/se(t,j);
CIplus(t,j)=Ghat(t)+1.96*se(t,j);
CIminus(t,j)=Ghat(t)-1.96*se(t,j);
ci(t,j)= (G(t,j) > CIminus(t,j) & G(t,j)< CIplus(t,j));
s2t=(fhat(t,:)*fhat(t,:)')*sig2eps(j)+(T/N)*VCV;;
se_epshat(t,j)=1.96*sqrt(s2t)/sqrt(T);

end;  % end T loop

Ttest(j)= mean(abs(Ntest(:,j)) > 1.96);
cover(j)=mean(ci(:,j));
NIR(j)=(epshat(:,j)'*epshat(:,j))/(Ghat'*Ghat);


[dum1,Index]=sort(Ntest(:,j),1);
maxtest(j)=dum1(up)  ;


[dum1,Index]=sort(abs(Ntest(:,j)),1);



abstest(j)= dum1(up) ;
d=Index(up);
maxCIplus=abs(Ghat(d))+maxcvH*se(j);
maxcover(j)=(abs(G(d,j)) < maxCIplus) & (-abs(G(d,j)) > -maxCIplus);



end;  % end j loop


result.t=Ttest;
result.max=maxtest;
result.abs=abstest;
result.r2=R2;
result.nir=NIR;


[T,k]=size(fhat);
[T,m]=size(G);

for j=1:nz;
results =canonical(fhat,G(:,j));
result.r2(j)=results.r2;
result.r2up(j)= results.r2+results.se;
result.r2dn(j)= results.r2-results.se;
end;



results=canonical(fhat,G);
cr2=results.r2;
cr2se=results.se;


